关于数据挖掘:R语言nlmenlmerlme4用非线性混合模型nonlinear-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 来实现 …)

## 必须用手进行汇总
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 拟合来说成果不错，给出了正当的后果。

nlsfit1 <- nls(X ~ SSfp)
coef(nlsfit1)

能够用 gnls 来拟合组间差别（我须要指定起始值

我的第一次尝试不太胜利。

gnls(X ~ SSfpl)

但如果我只容许 asymp.R 在各组之间变动，就能运行胜利。

params=symp.R~Group

绘制预测值。

g1 + geom_line()

这些看起来很不错（如果能失去置信区间就更好了 – 须要应用 delta 法或 bootstrapping）。

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

除了 7 号样本外，没有很多证据表明个体间的变异 …… 如果咱们想疏忽个体间的变异，能够用

anova(lm(res~Individual))

大的（p\）值能够承受个体间不存在变异的无效假设 …

更个别的诊断图 – 残差与拟合，同一个体的点用线连贯。能够发现，随着平均数的减少，方差会逐步减小。

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

点击题目查阅往期内容

非线性混合效应 NLME 模型反抗哮喘药物茶碱动力学钻研

左右滑动查看更多

01

02

03

04

我不能用 nlme 来解决三个参数因组而异模型，但如果我只容许 asymp 变动，就能够运行。

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

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

退出随机效应后，参数基本就没有什么变动。

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

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

咱们能够通过 AIC 或似然比测验来比拟模型

AICtab(nlmefit1,nlmefit2,weights=TRUE)

anova(nlmefit1,nlmefit2)

能够做一个 F 测试而不是 LRT（即思考到无限大小的修改）。

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

## 分母十分大的 F 测验。pf(diff/2,df1=2,df2=1000000,lower.tail=FALSE)

咱们不晓得真正相干的 df，但下面的总结表明 df 是 40。

nlmer

我想当初能够为 nlmer 失去正确的模型标准，但我找不到一个不便的语法来进行固定效应建模（即在这种状况下容许一些参数因组而异）– 当我构建了正确的语法，nlmer 无奈失去答案。

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

 nlmer(X ~ SSfpl(Day, asy, as, x, s) ~
         asy|Indi,)

依据我的了解，人们只须要构建本人的函数来封装固定效应构造；为了与 nlmer 一起应用，该函数还须要计算绝对于固定效应参数的梯度。这有点麻烦，但能够通过批改派生函数生成的函数，使之略微自动化。

2. 构建虚构变量:

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

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

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

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

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

尝试一下拟合:

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

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

当初咱们能够用稳定版（lme4.0）失去一个答案。

后果不现实

fixef(nlmerfit2)

range(predict(nlmerfit2))

我不能确定，在 nlmer 中是否有更简略的办法来做固定成果。

AD 模型生成器

咱们还能够应用 AD 模型生成器来解决这个问题。它能够解决更简单的模型，比方拟合更多参数的群体效应。

局部起因是我对 ADMB 的相熟水平较低，这有点吃力，最初我通过循序渐进的步骤才胜利。

最小的例子

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

## 设置数据：调整名称，等等
d0 <- c(list(nobs=nrow(d)),as.list(d0))
## 起始值：调整名称，减少数值
names(svec3) <- gsub("\.","",names(svec3))  ## 移除点
svec3$asympR <- 0.6 ## 繁多值
## 运行 
do_admb("algae0",
        data,
        params,
        run.opts)

后果不错

固定效应模型

当初尝试用固定效应分组，应用下面构建的虚构变量（也能够应用 if 语句，或者用 R[Group[i]]的 for 循环中的 R 值向量，或者（最佳抉择）为 R 传递一个模型矩阵 …）。咱们必须应用 elem_div 而不是 / 来对两个向量进行元素除法。

model1 <- "
参数局部
   向量 pred(1,nobs) // 预测值
   向量 Rval(1,nobs) // 预测值

过程局部
   pred = as+elem(Rval-asy,1.0+exp(-(Day-xmid)/scal) 
"

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

model1B <- "
参数局部
   向量 pred(1,nobs) // 预测值
   向量 Rval(1,nobs) // 预测值

过程局部
   pred = asym+ele(Rv-asy,1.0+exp(-(Da-xmi)/sc))。"

当然，在参数雷同的状况下，也能够工作。

随机效应

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

model2 <- "
参数局部
   向量 pred(1,nobs) // 预测值
   向量 Rval(1,nobs) // 预测值

过程局部
   pred = asym+elem
   f = 0.5*no*log(norm2(X-pr)/n)+norm2(R)。"

因为 ADMB 不解决稠密矩阵，也不惩办循环，如果将随机效应实现为(i=1; i<=nobs; i++) Rval[i] += Rsigma*Ru[Group[i]]，效率会略高，但我是懒人 / 我喜爱矩阵示意的紧凑性和可扩展性.

当初咱们终于能够测试 R 以外的参数的固定效应差别了。

model3 <- "
参数局部
   向量 prd(1,nobs) // 预测值
   向量 Rl(1,nobs) // 预测值
   向量 scalal(1,nobs)
   向量 xmal(1,nobs)
   sdror opr(1,nobs) // 输入预测值

程序局部
   Rval = XR*Rve+Rsma*（Z*Ru）。xmval = Xd*xdvec;....
   f = 0.5*nobs*log(norm2(X-pred)/nobs)+norm2(Ru)
"

后果:

summary(admbfit3)

有一个十分大的 AIC 差别。如上文所示，对 nlme 拟合的似然比 F 测试是作为一种练习 ……

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

fit3B <- do_admb(,
        data,
        params,
        re,
        run.opts=run.control)

plot2(list(cc),intercept=TRUE)

当初咱们对标准化的问题很困扰，所以（通过一番折腾）咱们能够在不同的面板上从新画出群体变动的参数。

诊断图

## 放弃条件模式 / 样本 - R 估计值

diagplot1 %+% dp2

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

拟合与残差

diagplot2 %+% dp2

叠加预测（虚线）:

g1 + geom_line

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

论断

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

似然剖析

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

lapply(Rsigmavec,fitfun)
## 尝试填补破绽
lapply(Rsigmavec2,fitfun)

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

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

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

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

扩大

• 更多地探讨分母 df 问题。参数疏导法 /MCMC？
• 咱们能够尝试在 xmid 和 scale 参数中退出随机效应。
• 在组间或作为 X 的函数的方差（无论是残差还是个体间的方差）中可能有额定的模式。

点击文末 “浏览原文”

获取全文残缺代码数据资料。

本文选自《R 语言 nlme、nlmer、lme4 用（非）线性混合模型 non-linear mixed model 剖析藻类数据实例》。

点击题目查阅往期内容

R 语言混合线性模型、多层次模型、回归模型剖析学生均匀问题 GPA 和可视化
R 语言线性混合效应模型（固定效应 & 随机效应）和交互可视化 3 案例
R 语言用 lme4 多层次（混合效应）狭义线性模型（GLM），逻辑回归剖析教育留级考察数据 R 语言 线性混合效应模型实战案例
R 语言混合效应逻辑回归（mixed effects logistic）模型剖析肺癌数据
R 语言如何用潜类别混合效应模型（LCMM）剖析抑郁症状
R 语言基于 copula 的贝叶斯分层混合模型的诊断准确性钻研
R 语言建设和可视化混合效应模型 mixed effect model
R 语言 LME4 混合效应模型钻研老师的受欢迎水平
R 语言 线性混合效应模型实战案例
R 语言用 Rshiny 摸索 lme4 狭义线性混合模型（GLMM）和线性混合模型（LMM）
R 语言基于 copula 的贝叶斯分层混合模型的诊断准确性钻研
R 语言如何解决线性混合模型中畸形拟合 (Singular fit) 的问题
基于 R 语言的 lmer 混合线性回归模型
R 语言用 WinBUGS 软件对学术能力测验建设档次（分层）贝叶斯模型
R 语言分层线性模型案例
R 语言用 WinBUGS 软件对学术能力测验（SAT）建设分层模型
应用 SAS，Stata，HLM，R，SPSS 和 Mplus 的分层线性模型 HLM
R 语言用 WinBUGS 软件对学术能力测验建设档次（分层）贝叶斯模型
SPSS 中的多层（等级）线性模型 Multilevel linear models 钻研整容手术数据
用 SPSS 预计 HLM 多层（档次）线性模型模型