原文链接：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一起应用，该函数还须要计算绝对于固定效应参数的梯度。这有点麻烦，但能够通过批改派生函数生成的函数，使之略微自动化。
- 构建虚构变量:
mm <- model.matrix(~Group,data=d)grp2 <- mm[,2]- 构建一个函数来评估预测值及其梯度；分组构造是硬编码的。
deriv(~A+((B0+B1*grp2+B2*grp3-A)/(1+exp((x-xmid)/scale)- 通过插入与传递给函数的参数名称相匹配的行来查看所产生的函数，并将这些参数名称调配给梯度矩阵。
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多层（档次）线性模型模型