如何逐步建立线性混合效应模型

Question

Im模拟了一种方法，该方法可以使用R中的nmle((根据固定的顺序将项组添加到线性混合效应模型中。作者说，她把变量分为四类，每一类都有一些主要效应，它们的二次项和线性项，以及它们之间的一阶相互作用。像"A,B,C,D, A^2,B^2,C^2,D^2" A+B,B+C,C+D A*B,B*C,C*D

https://besjournals.onlinelibrary.wiley.com/action/downloadSupplement?doi=10.1111%2F1365-2664.12478；file=jpe12478-sup-0001-SuppInfo.docx
您可以在表S2中查看详细信息。

在表S2中有一些交互术语，我不知道我是否应该使用：或*来表示这种关系，也不知道+在该表中的含义是什么。这是我厌倦的代码

mod1<-lme(C50_0.45~1,random = ~  1|sitecode,data=practice,method="ML")
mod2<-lme(C50_0.45~MAT,random = ~  1|sitecode,data=practice,method="ML")
mod3<-lme(C50_0.45~MAP,random = ~  1|sitecode,data=practice,method="ML")
mod4<-lme(C50_0.45~MGSL,random = ~  1|sitecode,data=practice,method="ML")
mod5<-lme(C50_0.45~MGDD,random = ~  1|sitecode,data=practice,method="ML")
mod6<-lme(C50_0.45~MAT+I(MAT^2),random = ~  1|sitecode,data=practice,method="ML")
mod7<-lme(C50_0.45~MAP+I(MAP^2),random = ~  1|sitecode,data=practice,method="ML")
mod8<-lme(C50_0.45~MGSL+I(MGDD^2),random = ~  1|sitecode,data=practice,method="ML")
mod9<-lme(C50_0.45~MGDD+I(MGSL^2),random = ~  1|sitecode,data=practice,method="ML")
mod10<-lme(C50_0.45~MAT+MAP,random = ~  1|sitecode,data=practice,method="ML")
mod11<-lme(C50_0.45~MGDD+MAP,random = ~  1|sitecode,data=practice,method="ML")
mod12<-lme(C50_0.45~MGSL+MAP,random = ~  1|sitecode,data=practice,method="ML")
mod13<-lme(C50_0.45~MAT*MAP,random = ~  1|sitecode,data=practice,method="ML")
mod14<-lme(C50_0.45~MGDD*MAP,random = ~  1|sitecode,data=practice,method="ML")
mod15<-lme(C50_0.45~MGSL*MAP,random = ~  1|sitecode,data=practice,method="ML")
AIC(mod1)-AIC(mod2)
AIC(mod1)-AIC(mod3)
AIC(mod1)-AIC(mod4)
...
anova(mod1,mod2)
anova(mod1,mod3)
anova(mod1,mod4)
...

这就是我试图在AIC((中进行测试的方式，直到找到在卡方似然比删除测试(LRT(中降低AIC和显著的项，然后将这些项组合成一个新的模型。在这些测试之后；然后添加代表对功能的不同控制的另一组变量，并重复该过程">

但我无法回放作者在表S3-6中所做的结果

这是数据https://datadryad.org/stash/dataset/doi:10.5061/dryad.s7867(Manning等人，2015(。

FIRST例如

mod3<-lme(C50_0.45~MAP,random = ~  1|sitecode,data=practice,method="ML")
mod7<-lme(C50_0.45~MAP+I(MAP^2),random = ~  1|sitecode,data=practice,method="ML")
mod10<-lme(C50_0.45~MAT+MAP,random = ~  1|sitecode,data=practice,method="ML")
mod11<-lme(C50_0.45~MGDD+MAP,random = ~  1|sitecode,data=practice,method="ML")
mod12<-lme(C50_0.45~MGSL+MAP,random = ~  1|sitecode,data=practice,method="ML")
mod13<-lme(C50_0.45~MAT*MAP,random = ~  1|sitecode,data=practice,method="ML")
mod14<-lme(C50_0.45~MGDD*MAP,random = ~  1|sitecode,data=practice,method="ML")
mod15<-lme(C50_0.45~MGSL*MAP,random = ~  1|sitecode,data=practice,method="ML")
AIC(mod1)-AIC(mod3)
[1] 6.281076
AIC(mod1)-AIC(mod7)
[1] 4.440664
AIC(mod1)-AIC(mod10)
[1] 4.319208
AIC(mod1)-AIC(mod11)
[1] 4.390488
AIC(mod1)-AIC(mod12)
[1] 4.332252
AIC(mod1)-AIC(mod13)
[1] 3.261881
AIC(mod1)-AIC(mod14)
[1] 3.205841
AIC(mod1)-AIC(mod15)
[1] 4.132244

所有这些模型在LRT中都很重要，但如果我在模型中保留这些项，mod3中的MAP与mod7、10、11、12、13、14、15中的MAP相同，并且是像这个这样的一类最终模型

mods1<-lme(C50_0.45 ~MAP+MAP+I(MAP^2)+MAT+MAP+MGDD+MAP+MGDD+MAP+MGSL+MAP+MAT:MAP+MGDD*MAP+MGDD*MAP+MGSL*MAP,random = ~  1|sitecode,data=practice,method="ML")

我不知道启动LRT是否正确，因为例如，删除(MGSL+MAP(或(MAT+MAP(就像只删除MGSL或MAT 一样

第二个正如你在数据中看到的，我不知道如何按照作者在表S2中描述的方式添加类型，这可能是我无法复制正确的最终模型的原因

关于我的模型，由于这些问题，我已经尝试过了，但你知道，有时会删除必要的术语，或者只是无法删除一些不必要的术语。在C50_250或总共C中，我可以复制与表中相同的最终模型，但在其他模型中不起作用。

Answer 1

我能够使用下面的模型重现表S3中的值。也许这将帮助你找出如何复制他们的其余方法。data是指从您提供的链接中读入的原始表。

library(dplyr)
library(lme4)
practice <- dplyr::filter(data, datatype=="Initial")
model <- lme4::lmer(TotalC ~ MAT + I(MAT^2) + moist + I(moist^2) + PH + I(PH^2) + (1|site), data=practice)

我假设对于表S3至S6中提到的每个模型的响应是不同的。S3似乎使用TotalC，但推测C250、C50-250等是其他模型的响应。

在回答您关于指定交互项的问题时，此链接对于理解如何在R中指定公式非常有用：https://conjugateprior.org/2013/01/formulae-in-r-anova/

简而言之，A*B相当于A+B+A:B