运用广义线性混合模型分析随机区组重复测量的试验资料
张久权, 闫慧峰, 褚继登, 李彩斌
作物学报
2021, 47 ( 2):
294-304.
DOI: 10.3724/SP.J.1006.2021.04085
重复测量试验对同一受试对象进行多次测量, 各时间点数据间存在自相关性, 进行方差分析和均值比较时需要进行特殊处理。虽然此方法在农业等研究领域运用十分广泛, 但目前有效地相关统计方法鲜见。为了建立操作简单、实用性强、结果可靠的统计分析方法, 本研究采用SAS的广义线性混合模型(Generalized Linear Mixed Models, GLIMMIX), 以随机区组重复测量试验资料为例, 说明了协方差结构筛选、方差分析和均值比较的具体方法。结果表明, 用传统的裂区设计、多变量统计等方法会造成资料信息浪费, 统计功效降低, 缺区无法处理等问题, 甚至会导致错误的结论。GLIMMIX能很好地处理自相关问题, 功能强大, 结果可靠, 使用简单, 允许缺区, 是进行重复测量试验资料方差分析和均值比较的理想方法。目前在国内将其运用到农学类试验数据的统计分析的相关报道鲜见, 该文在本领域具有很强的实用性和创新性。
| Rain | N | Block | Core | Times | Y | | 1 | 1 | 1 | 1 | 6 | 20.20 | | 1 | 1 | 1 | 1 | 7 | 7.37 | | 1 | 1 | 1 | 1 | 8 | 5.75 | | 1 | 1 | 1 | 1 | 9 | 4.25 | | 1 | 1 | 1 | 1 | 10 | 2.42 | | 1 | 1 | 1 | 1 | 11 | 4.69 | | 1 | 1 | 1 | 1 | 12 | 5.13 | | 1 | 1 | 2 | 2 | 6 | 20.45 | | 1 | 1 | 2 | 2 | 7 | 12.33 | | … | … | … | | … | … | | 3 | 3 | 3 | 27 | 12 | 6.16 |
View table in article
表1
数据Microsoft Excel输入表示例
正文中引用本图/表的段落
1.1.1 土柱试验 该试验为室内模拟试验二因子(降水强度、N肥水平)三重复全因子设计, 小区排列方式为随机区组。土壤为中国农业科学院西南试验基地(四川西昌市)收集的紫色土。降水强度分别为0.44、0.68、1.20 mm min-1, 降水量均为50 mm 次-1。设N1 (常规施肥, CK)、N2 (减量30%)、N3 (减量30%+生物炭)共3个氮肥水平。对土柱进行了12次模拟降雨, 第6次开始收集到下渗淋洗液, 一直到第12次降雨完成, 每个土柱收集到7次淋洗液, 时间间隔相等。采用碱性过硫酸钾消解紫外分光光度法测定淋洗液中的全氮[14]。7次淋洗被当作重复测量因素。采用SAS 9.4的GLIMMIX进行统计分析, 数据输入示例见表1。
程序说明: 输入代码时注意分号为英文字符。数据存放在E:\数据\N.xls sheet1里, 格式见表1。程序(1)读取Excel文件的数据。程序(2)~(8)内容基本相同, 仅在“type =”后的协方差结构选项不同。对某些协方差结构, 包括一阶自回归[autoregressive(1), AR(1)]、循环相关(toeplitz, TOEP)、一阶前依赖[antedependence, ANTE(1)]等协方差结构, 需要考虑对象间误差, 因此需要增加random block core(Rain N)语句[15]。程序(2)~(8)调用GLIMMIX过程, 分别采用方差分量结构(variance components, VC)、复合对称结构(compound symmetry, CS)、不规则结构(unstructured, UN)、空间幂相关结构[space power, SP(POW)]、一阶自回归[AR(1)]、循环相关结构(TOEP)、一阶前依赖结构[ANTE(1)]模型进行方差分析。Class语句列出所有的分类变量。Model语句根据上文(1)式编写, Model语句后仅需要列出固定效应[4], 注意此与mixed的不同。Rain|N|times表示3个因素的主效应、2级和3级交互作用效应的线性组合。采用KR法对标准误和自由度进行修正[16]。random_residual_语句相当于mixed模型中的repeated语句[5]。选项type指定协方差结构类型。选项sub指定数据集中的受试对象(subject), 本例中为土柱(core), 即小区, 若为单因素试验, 直接指定对象名称; 若为多因素试验, 则在对象名称后加因素名称, 并加“()”。ods output语句输出模型拟合信息fitstatistics。程序(9)建立宏rename, 对数据集中已有变量名value更名, 否则原有数据列value的数据会被覆盖[17]。程序(10)对不同协方差结构拟合参数进行合并, 以便程序(11)输出拟合结果, 为挑选最佳协方差结构提供依据。
这里, rate为生物炭用量, block为区组, plot为小区, 即试验对象, time为取样时间。由于本田间试验取样时间间隔不相等, 对于SP(POW)协方差结构, 需要计算时间间隔变量, 本列以天数days表示, 并增加语句 “random _residual_/type=SP(POW)(days) sub=plot;”。对于需要增加random语句的协方差结构, 包括SP(POW)、ANTE(1)、AR(1)+RE、TOEP, random语句写法为“random plot block block*time;”。
从表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)。可以继续进行两处理水平在各特定时间点的差异比较, 找出最优处理组合。
本文的其它图/表
|
