运用广义线性混合模型分析随机区组重复测量的试验资料
张久权, 闫慧峰, 褚继登, 李彩斌
作物学报
2021, 47 ( 2):
294-304.
DOI: 10.3724/SP.J.1006.2021.04085
重复测量试验对同一受试对象进行多次测量, 各时间点数据间存在自相关性, 进行方差分析和均值比较时需要进行特殊处理。虽然此方法在农业等研究领域运用十分广泛, 但目前有效地相关统计方法鲜见。为了建立操作简单、实用性强、结果可靠的统计分析方法, 本研究采用SAS的广义线性混合模型(Generalized Linear Mixed Models, GLIMMIX), 以随机区组重复测量试验资料为例, 说明了协方差结构筛选、方差分析和均值比较的具体方法。结果表明, 用传统的裂区设计、多变量统计等方法会造成资料信息浪费, 统计功效降低, 缺区无法处理等问题, 甚至会导致错误的结论。GLIMMIX能很好地处理自相关问题, 功能强大, 结果可靠, 使用简单, 允许缺区, 是进行重复测量试验资料方差分析和均值比较的理想方法。目前在国内将其运用到农学类试验数据的统计分析的相关报道鲜见, 该文在本领域具有很强的实用性和创新性。
效应
Effect | 方差分量
VC | 复合对称
CS | 不规则
UN | 空间幂相关
SP | 一阶自回归AR(1) | 循环相关
TOEP | 一阶前依赖ANTE(1) | | Rain | <0.0001 | <0.0001 | <0.0001 | <0.0001 | <0.0001 | <0.0001 | <0.0001 | | N | <0.0001 | <0.0001 | <0.0001 | <0.0001 | <0.0001 | 0.0001 | <0.0001 | | Rain*N | 0.2375 | 0.2609 | 0.2809 | 0.3023 | 0.3023 | 0.3161 | 0.2822 | | Times | <0.0001 | <0.0001 | <0.0001 | <0.0001 | <0.0001 | <0.0001 | <0.0001 | | Rain*Times | 0.0033 | 0.0038 | 0.6167 | 0.0059 | 0.0059 | 0.0172 | 0.5533 | | N*Times | <0.0001 | <0.0001 | 0.0017 | <0.0001 | <0.0001 | <0.0001 | 0.0003 | | Rain*N*Times | 0.8526 | 0.8537 | 0.9457 | 0.8584 | 0.8584 | 0.9225 | 0.9727 |
View table in article
表2
采用不同协方差结构模型时F检验P值(III型检验)
正文中引用本图/表的段落
SAS程序I输出了F检验结果, 表2是对土柱试验数据7种协方差结构模型F检验结果的汇总。可以看出, 协方差结构模型对F检验的结果影响较大。降雨强度与施肥方式之间的交互作用(Rain*N)虽然P值都大于0.05, 但差异较大, 最低的为0.2375, 最高的为0.3161, 相差33.09%。降雨强度与淋洗次数之间的交互作用效应(Rain*Times) UN和ANTE(1)不显著, 而其余4种模型显示极显著。因此, 选用恰当的协方差结构模型进行F检验很关键, 否则有可能得出错误的结论。
根据统计学理论, 选用协方差结构模型时, 可参考赤池信息准则(akaike information criterion, AIC)、为小样本修正的赤池信息准则(akaike information corrected criterion, AICC)、修正的赤池信息准则(corrected akaike information criterion, CAIC)、贝叶斯信息准则(bayesian information criterion, BIC)、汉南-奎因信息准则(hannan-quinn information criterion, HQIC)、-2残差对数似然值准则(-2 res log likelihood, -2logL)[1,4-5,15,17]等, 值越小表示拟合性越好, 如果相近, 可通过χ2检验[17]并结合试验本身的特点进行判断。另外, 协方差结构模型需要估算的参数个数越少越好。从表3可以看出, 土柱试验UN和ANTE(1)模型各准则值明显低于其余5种协方差模型的值, 应优先考虑。UN模型协方差矩阵中需要估算的参数个数在所有模型中最多, 为n(n-1)/2, n为重复测量次数, 计算时可能无法收敛, 估算无法完成[5]。本例中n=7, 参数个数为21个; ANTE(1)模型需要估算的参数为2n-1, 本例中为15。综合考虑, 选用ANTE(1)模型进行F检验和均值比较。
从表2可以看出, 第6次淋洗时, N2比N1的总氮淋失量低8.05 g (图1), 接近5%差异显著水准(P = 0.078), 此为特定时间(第6次淋洗)两处理水平(N1 vs. N2)间的比较; 语句(5)输出的结果表明, N2处理总N淋失量第6次比第8次淋洗高10.94 g, 差异达1%极显著水平(P = 0.0025), 我们还可以检验第6次与第8、9、10、11、12次淋洗的差异显著性, 也可以对其他时间进行差异显著性比较, 确定总N淋失量随时间的变化规律。语句(6)输出的结果表明, N1第7次淋洗比N2处理第9次淋洗全氮淋失量多6.38 g氮(表5和图1), 差异极显著(P = 0.0008)。可以继续进行两处理水平在各特定时间点的差异比较, 找出最优处理组合。
语句(7)输出结果表明, 3种氮水平平均全氮淋失量第8次与第10次淋洗相差0.94 g, 未达到5%差异显著水平。由于N*Times交互作用效果显著, 这种比较并不合适, 此处仅仅用来说明方法。语句(8)输出的结果表明, 所有淋洗平均全氮淋失量N2比N1处理低1.78 g, 差异达5%显著水平。同样原因, 这种比较在此例中不合适, 仅用于方法说明。田间试验均值间比较思路类似, 不再赘述。
本文的其它图/表
|
