关于算法:R语言因子实验设计nlme拟合非线性混合模型分析有机农业施氮水平

61次阅读

共计 3343 个字符,预计需要花费 9 分钟才能阅读完成。

原文链接:http://tecdat.cn/?p=24134 

测试非线性回归中的交互作用

因子试验在农业中十分广泛,它们通常用于测试试验因素之间相互作用的重要性。例如,能够在两种不同的施氮程度(例如高和低)下进行基因型评估,以理解基因型的排名是否取决于营养的可用性。对于那些不太理解农业的人,我只会说这样的评估是相干的,因为咱们须要晓得咱们是否能够举荐雷同的基因型,例如,在传统农业(高氮可用性)和有机农业中农业氮的可用性。

让咱们思考一个试验,在该试验中,咱们在残缺的区组因子设计中以两种氮含量(“高”和“低”)测试了三种基因型(为了简便起见,咱们称它们为 A、B 和 C),并进行四次反复。在八个不同的工夫(收获后天数:DAS)从 24 个地块中的每一个中取出生物量子样本,以评估生物量随工夫的增长。

加载数据并将“Block”变量转换为一个因子。

head(dataset)

数据集由以下变量组成:

  1. ‘Id’:察看的数字代码
  2. ‘DAS’:即收获后的天数。这是采集样本的时刻
  3. ‘Block’, ‘Plot’, ‘GEN’ 和 ‘N’ 别离代表每个察看的块、图、基因型和氮程度
  4. “产量”代表播种的生物量。

查看察看到的增长数据,如下图所示。

咱们看到增长是对称的(大略是逻辑的)并且察看的方差随着工夫的推移而减少,即方差与冀望因变量成正比。

问题是:咱们如何剖析这些数据?

模型

咱们能够凭教训假如生物量和工夫之间的关系是逻辑的:

其中 Y 是第 i 个基因型、第 j 个氮程度、第 k 个区块和第 l 个小区在 X 工夫察看到的生物量产量,d 是工夫进入无穷大时的最大渐进生物量程度,b 是拐点处的斜率,而 e 是生物量产量等于 d / 2 时的工夫。咱们次要对参数 d 和 e 感兴趣:第一个参数形容基因型的产量后劲,而第二个参数给出成长速度的测量。
 

每个小区都有反复测量,因而,模型参数可能显示出一些变动,取决于基因型、氮程度、区块和小区。特地是,假如 b 是相当恒定的,并且独立于上述因素,而 d 和 e 可能依据以下公式发生变化,这是能够承受的。

其中,对于每个参数,μμ 是截距,gg 是第 i 个基因型的固定效应,NN 是第 j 个氮程度的固定效应,gNgN 是固定交互效应,θ 是区块的随机效应,而 ζζ 是区块边疆块的随机效应。这两个方程齐全等同于通常用于线性混合模型的方程,在双因素因子区块设计的状况下,其中 ζζ 是残差误差项。事实上,原则上,咱们也能够思考两步法的拟合程序,即咱们。

  1. 将逻辑模型拟合到每个图的数据并取得 d 和 e 的估计值
  2. 应用这些预计来拟合线性混合模型

咱们不会在这里谋求这种两步法,咱们将专一于一步拟合。

谬误的办法

如果察看是独立的(即没有块和没有反复测量),这个模型能够通过应用传统的非线性回归来拟合。

编码报告如下。产量 “ 是(∼)DAS 的函数,通过一个三参数的 Logistic 函数。对于基因型和氮程度的不同组合必须拟合不同的曲线(id = N:GEN),只管这些曲线应该局部地基于独特的参数值(’models = …)。model” 参数须要一些补充阐明。它必须是一个矢量,其元素数与模型中的参数数一样多(在本例中是三个:B、D 和 E)。每个元素代表一个变量的线性函数,并按字母程序指向参数,即第一个元素指 b,第二个指 d,第三个指 e。参数 b 不依赖于任何变量(’~1’),因而在不同的曲线上拟合出一个常数;d 和 e 依赖于基因型和氮程度的齐全因子组合(~N*GEN = ~N + GEN + N:GEN)。最初,咱们应用参数 ’bcVal = 0.5’ 来指定咱们打算应用转换两边办法,即对方程的两边进行对数转换。这对于思考异方差是必要的,但它不影响参数估计。

rm(Yield ~ DAS, dta =daas,
           id = GEN:N,
            model = c(~ 1,  ~ N\*GEN,  ~ N\*GEN))

这个模型对于其余状况(无区块和无反复测量)可能是有用的,但在咱们的例子中是谬误的。事实上,观测值在区块和地块内是聚在一起的;如果疏忽这一点,咱们就违反了模型残差独立的假如。残差与拟合值的图显示,不存在异方差的问题。

 

思考到上述情况,咱们必须在这里应用不同的模型,只管我将证实这种拟合可能会很有用。

非线性混合模型拟合

为了解释察看的类,咱们切换到非线性混合效应模型(NLME)。一个不错的抉择是 ’nlme()’ 函数(Pinheiro 和 Bates,2000),只管有时语法可能很麻烦。咱们须要指定以下内容:

  1. 模型参数的线性函数。nlme’ 函数中的 ’fixed’ 参数与下面函数中的 ’models’ 参数十分类似,即它须要一个列表,其中每个元素都是变量的线性函数。惟一的区别是,参数名称须要在函数的左侧指定。
  2. 模型参数的随机效应。这些是通过应用 “ 随机 “ 参数来指定的。在这种状况下,参数 d 和 e 无望在一个区块内的不同区块和不同地块之间显示随机变动。为了简略起见,因为参数 b 不受基因型和氮程度的影响,咱们也心愿它在区块和地块之间不显示任何随机变动。
  3. 模型参数的起始值。咱们须要指定模型参数的初始值。在这种状况下,我决定应用下面非线性回归的输入。

方程两边的变换。

nlme(sqtYld ~ srtLS.L3(DAS, b, d, e)
tTable

从上图中,咱们看到整体拟合良好。随机效应的固定效应和方差重量按如下形式取得:

summary(mdnle1)

当初,让咱们回到咱们最后的目标:测试 “ 基因型 x 氮 “ 交互作用的显著性。事实上,咱们有两个可用的测试:一个是参数 d,一个是参数 e。首先,咱们对两个 “ 简化 “ 模型进行编码。为此,咱们把固定效应从 ’~ N*GEN’ 改为 ’~ N + GEN’。同样在这种状况下,咱们应用非线性回归拟合来取得模型参数的起始值,用于上面的 NLME 模型拟合。

 

mdNave2 <- Yied ~ DAS
            modes = c( ~ 1,  ~ N + GEN,  ~ N * GEN

mdnle2 <-sqrt(Yeld) ~ sqrt(NS.3(DAS, b, d, e
mdaie3 <- Yied ~ DAS
            crid = N:GEN
            modls = c( ~ 1,  ~ N*GEN,  ~ N + GEN

mdne3 <- sqrt(Yild) ~ sqrt(NS.L3(DAS, b, d, e
                    ranom = d + e~ 1|Blck/Pot
                    fixd = lit(b ~ 1 d ~ N*GN, e ~ N + GN)

让咱们思考第一个放大的模型 ’modnlme2’。在这个模型中,” 基因型 x 氮 “ 的交互作用曾经被移除,用于 d 参数。咱们能够通过应用似然比测验来比拟这个模型和残缺的模型 “modnlme1″。

aova(mode1, mdne2)

该测验是显著的,但两个模型的 AIC 值十分靠近。思考到混合模型中的 LRT 通常比拟宽松,应该能够得出结论,” 基因型 x 氮素 “ 的交互作用不显著,因而,用 d 参数掂量的基因型在产量后劲方面的排名应该与氮素程度无关。

当初让咱们思考第二个简化模型“modnlme3”。在第二个模型中,“e”参数的“基因型 x 氮”交互作用已被删除。咱们还能够应用似然比测验将此简化模型与残缺模型“modnlme1”进行比拟:

anva(mole1, mdle3)

在这第二次测验中,” 基因型 x 氮素 “ 交互作用的不显著性仿佛比第一次测验更可信。


 

最受欢迎的见解

1. 基于 R 语言的 lmer 混合线性回归模型

2.R 语言用 Rshiny 摸索 lme4 狭义线性混合模型(GLMM)和线性混合模型(LMM)

3.R 语言线性混合效应模型实战案例

4.R 语言线性混合效应模型实战案例 2

5.R 语言线性混合效应模型实战案例

6. 线性混合效应模型 Linear Mixed-Effects Models 的局部折叠 Gibbs 采样

7.R 语言 LME4 混合效应模型钻研老师的受欢迎水平

8.R 语言中基于混合数据抽样 (MIDAS) 回归的 HAR-RV 模型预测 GDP 增长回归的 HAR-RV 模型预测 GDP 增长 ”)

9. 应用 SAS,Stata,HLM,R,SPSS 和 Mplus 的分层线性模型 HLM

正文完
 0