拓端tecdat|R语言nlme、nlmer、lme4用（非）线性混合模型non-linear mixed model分析藻类数据实例

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

原文出处：拓端数据部落公众号

序言

混合线性模型，又名多层线性模型(Hierarchical linear model)。它比较适合处理嵌套设计(nested)的实验和调查研究数据。此外，它还特别适合处理带有被试内变量的实验和调查数据，因为该模型不需要假设样本之间测量独立，且通过设置斜率和截距为随机变量，可以分离自变量在不同情境中(被试内设计中常为不同被试)对因变量的作用。

简单的说，混合模型中把研究者感兴趣的自变量对因变量的影响称为固定效应，把其他控制的情景变量称为随机效应。由于模型中包括固定和随机效应，故称为混合线性模型。无论是用方差分析进行差异比较，还是回归分析研究自变量对因变量的影响趋势，混合线性模型比起传统的线性模型都有更灵活的表现。
 

非线性混合模型就是通过一个连接函数将线性模型进行拓展，并且同时再考虑随机效应的模型。

非线性混合模型常常在生物制药领域的分析中会用到，因为很多剂量反应并不是线性的，如果这个时候数据再有嵌套结构，那么就需要考虑非线性混合模型了。

 本文中我们用（非）线性混合模型分析藻类数据。这个问题的参数是：已知截距（0日值）在各组和样本之间是相同的。

数据

用lattice和ggplot2绘制数据。

xyplot(jitter(X)~Day, groups=Group)

ggplot版本有两个小优势。1. 按个体和群体*均数添加线条[用stat_summary应该和用xyplot的type="a "一样容易])；2.调整点的大小，使重叠的点可视化。(这两点当然可以用自定义的 panel.xyplot 来实现 ...)

1.  
   ## 必须用手进行汇总
2.  
   ggplot(d,aes(x=Day,y=X,colour=Group))

从这些图片中得出的主要结论是：（1）我们可能应该使用非线性模型，而不是线性模型；（2）可能存在一些异方差（在较低的*均值上有较大的方差，好像在 X=0.7的数据有一个 "天花板"）；看起来可能存在个体间的变化（特别是基于t2的数据，其中个体曲线*乎*行）。然而，我们也将尝试线性拟合来说明问题。

使用nlme

用lme的线性拟合失败。

LME <- lme(X ~ 1, random = ~Day|Individual, data=d)

如果我们用control=lmeControl(msVerbose=TRUE))运行这个程序，就会得到输出，最后是。 

可以看到考虑到组*日效应的模型也失败了。 

LME1 <- lme(X ~ Group*Day, random = ~Day|Individual, data=d)

我试着用SSfpl拟合一个非线性模型，一个自启动的四参数Logistic模型（参数为左渐*线、右渐*线、中点、尺度参数）。这对于nls拟合来说效果不错，给出了合理的结果。 

1.  
   nlsfit1 <- nls(X ~ SSfp)
2.  
   coef(nlsfit1)

可以用gnls来拟合组间差异（我需要指定起始值

我的第一次尝试不太成功。 

1.  
   gnls(
2.  
   X ~ SSfpl)

 但如果我只允许asymp.R在各组之间变化，就能运行成功。 

params=symp.R~Group

绘制预测值。

1.  
    
2.  
   g1 + geom_line()

这些看起来很不错（如果能得到置信区间就更好了--需要使用delta法或bootstrapping）。 

1.  
   dp <- data.frame(d,res=resid(gnlsfit2),fitted=fitted(gnlsfit2))
2.  
   (diagplot1 <- ggplot(dp,aes(x=factor(Individual),
3.  
   y=res,colour=Group))+
4.  
   geom_boxplot(outlier.colour=NULL)+
5.  
   scale_colour_brewer(palette="Dark2"))
6.  
    

除了7号样本外，没有很多证据表明个体间的变异......如果我们想忽略个体间的变异，可以用 

anova(lm(res~Individual))

大的（p）值可以接受个体间不存在变异的无效假设...

更一般的诊断图--残差与拟合，同一个体的点用线连接。可以发现，随着*均数的增加，方差会逐渐减小。

plot(dp,(x=fitted,y=res,colour=Group))

我不能用nlme来处理三个参数因组而异模型，但如果我只允许asymp变化，就可以运行。 

nlme(model=list(fixed=with(c(asymp.R,xmid,scale,asymp.L),...)

右侧渐*线中的方差估计值是非零的。

加入随机效应后，参数根本就没有什么变化。 

最大的比例差异是3.1%（在比例参数中）。 

1.  
   nlmefit2 <- update(list(asyR+xmd+scal+asp ~1),
2.  
   start )

我们可以通过AIC或似然比检验来比较模型

AICtab(nlmefit1,nlmefit2,weights=TRUE)

anova(nlmefit1,nlmefit2)

可以做一个F测试而不是 LRT（即考虑到有限大小的修正）。
 

pchisq(iff,df=2,lower.tail=FALSE)

 

1.  
   ##分母非常大的F检验。
2.  
   pf(diff/2,df1=2,df2=1000000,lower.tail=FALSE)

我们不知道真正相关的df，但上面的总结表明df是40。 

nlmer

我想现在可以为nlmer得到正确的模型规范，但我找不到一个方便的语法来进行固定效应建模（即在这种情况下允许一些参数因组而异）--当我构建了正确的语法，nlmer无法得到答案。

基本的RE模型（没有群体效应）运行良好。

1.  
   nlmer(
2.  
   X ~ SSfpl(Day, asy, as, x, s) ~
3.  
   asy|Indi,)

根据我的理解，人们只需要构建自己的函数来封装固定效应结构；为了与nlmer一起使用，该函数还需要计算相对于固定效应参数的梯度。这有点麻烦，但可以通过修改派生函数生成的函数，使之稍微自动化。

1. 构建虚拟变量:

1.  
   mm <- model.matrix(~Group,data=d)
2.  
   grp2 <- mm[,2]

1. 构建一个函数来评估预测值及其梯度；分组结构是硬编码的。

deriv(~A+((B0+B1*grp2+B2*grp3-A)/(1+exp((x-xmid)/scale)

1. 通过插入与传递给函数的参数名称相匹配的行来查看所产生的函数，并将这些参数名称分配给梯度矩阵。

1.  
   L1 <- grep("^ +\.value +<-")
2.  
   L2 <- grep("^ +attr\(\.value",)
3.  
   eval(parse(text))

尝试一下拟合:

1.  
   nlmer(
2.  
   X ~ fpl(Day, asym, as, asymp, asR3, xmi, sca) ~
3.  
   as|Indi,
4.  
   start = list(nlpars)),data=d)

失败了（但我认为这是由于nlmer本身造成的，而不是设置有什么根本性的问题）。为了确定，我应该按照同样的思路生成一个更大的人工数据集，看看我是否能让它工作起来。

现在我们可以用稳定版（lme4.0）得到一个答案。

结果不理想

fixef(nlmerfit2)

range(predict(nlmerfit2))

我不能确定，在nlmer中是否有更简单的方法来做固定效果。

AD模型生成器

我们还可以使用AD模型生成器来解决这个问题。它可以处理更复杂的模型，比如拟合更多参数的群体效应。

部分原因是我对ADMB的熟悉程度较低，这有点费劲，最后我通过循序渐进的步骤才成功。

最小的例子

首先尝试没有随机效应、分组变量等。(即等同于上面的nls拟合）。)

1.  
   ##设置数据：调整名称，等等
2.  
   d0 <- c(list(nobs=nrow(d)),as.list(d0))
3.  
   ##起始值：调整名称，增加数值
4.  
   names(svec3) <- gsub("\.","",names(svec3)) ## 移除点
5.  
   svec3$asympR <- 0.6 ## 单一值
6.  
   ## 运行
7.  
   do_admb("algae0",
8.  
   data,
9.  
   params,
10.  
    run.opts)

结果不错

固定效应模型

现在尝试用固定效应分组，使用上面构建的虚拟变量（也可以使用if语句，或者用R[Group[i]]的for循环中的R值向量，或者（最佳选择）为R传递一个模型矩阵...）。我们必须使用elem_div而不是/来对两个向量进行元素除法。

1.  
    
2.  
   model1 <- "
3.  
   参数部分
4.  
   向量 pred(1,nobs) // 预测值
5.  
   向量Rval(1,nobs) //预测值
6.  
    
7.  
   过程部分
8.  
   pred = as+elem(Rval-asy,1.0+exp(-(Day-xmid)/scal)
9.  
   "
10.  
     

试着用模型矩阵来代替它。

1.  
    
2.  
   model1B <- "
3.  
   参数部分
4.  
   向量 pred(1,nobs) // 预测值
5.  
   向量Rval(1,nobs) //预测值
6.  
    
7.  
   过程部分
8.  
   pred = asym+ele(Rv-asy,1.0+exp(-(Da-xmi)/sc)) 。
9.  
   "
10.  
     

当然，在参数相同的情况下，也可以工作。

随机效应

现在添加随机效应。回归函数并没有完全实现随机效应模型（尽管这应该在即将到来的版本中被修复），所以我们用公式减去（n/2 log({RSS}/n)），其中RSS是残差*方和。

1.  
   model2 <- "
2.  
   参数部分
3.  
   向量 pred(1,nobs) // 预测值
4.  
   向量Rval(1,nobs) //预测值
5.  
    
6.  
   过程部分
7.  
   pred = asym+elem
8.  
   f = 0.5*no*log(norm2(X-pr)/n)+norm2(R)。
9.  
   "

由于ADMB不处理稀疏矩阵，也不惩罚循环，如果将随机效应实现为(i=1; i<=nobs; i++) Rval[i] += Rsigma*Ru[Group[i]]，效率会略高，但我是懒人/我喜欢矩阵表示的紧凑性和可扩展性.

现在我们终于可以测试R以外的参数的固定效应差异了。

1.  
    
2.  
   model3 <- "
3.  
   参数部分
4.  
   向量 prd(1,nobs) // 预测值
5.  
   向量Rl(1,nobs) // 预测值
6.  
   向量 scalal(1,nobs)
7.  
   向量xmal(1,nobs)
8.  
   sdror opr(1,nobs) //输出预测值
9.  
    
10.  
    程序部分
11.  
    Rval = XR*Rve+Rsma*（Z*Ru）。
12.  
    xmval = Xd*xdvec;....
13.  
    f = 0.5*nobs*log(norm2(X-pred)/nobs)+norm2(Ru)
14.  
    "

结果:

summary(admbfit3)

有一个非常大的AIC差异。如上文所示，对nlme拟合的似然比F测试是作为一种练习......

对于该图，最好是按组指定参数重新进行拟合，而不是按基线+对比度进行拟合。

1.  
    
2.  
   fit3B <- do_admb(,
3.  
   data,
4.  
   params,
5.  
   re,
6.  
   run.opts=run.control)

1.  
   plot2(list(cc),intercept=TRUE)
2.  
    

现在我们对标准化的问题很困扰，所以（经过一番折腾）我们可以在不同的面板上重新画出群体变化的参数。

诊断图

1.  
   ##放弃条件模式/样本-R估计值
2.  
    
3.  
   diagplot1 %+% dp2

也许这暗示了两个实验组中更大的差异？

拟合与残差

diagplot2 %+% dp2

叠加预测（虚线）:

g1 + geom_line

如果能生成*滑的预测曲线（即对中间的日值），那就更好了，但也更繁琐。

结论

• 从参数估计中得出的主要结论是，第三组下降得更早一些（xmidvec更小），同时下降得更远（Rvec更低）。

似然分析

计算一个( sigma^2_R ) 似然函数的代码并不难，但运行起来有点麻烦：它很慢，而且计算在置信度下限附*的几个点上出现了非正-无限矩阵；我运行了另一组值，试图充分覆盖这个区域。

1.  
   lapply(Rsigmavec,fitfun)
2.  
   ## 尝试填补漏洞
3.  
   lapply(Rsigmavec2,fitfun)

带有插值样条的剖面图和似然比检验分界线。 

在sigma^2_R 上的95%剖面置信区间是{0.0386,0.2169}。

我没有计算过，但转换后的剖面图（在对应于偏离度与最小偏离度的*方根偏差的 y ）上，所以二次剖面将是一个对称的V）显示，二次*似对这种情况相当糟糕 ...

1.  
   ggplot(sigma,sqrt(2*(NLL-min(NLL))+
2.  
   geom_point()

扩展

• 更多地讨论分母df问题。参数引导法/MCMC？
• 我们可以尝试在xmid和scale参数中加入随机效应。
• 在组间或作为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增长

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