【原】R数据分析：用lme4包拟合线性和非线性混合效应模型

快一个月没更新文章啦，今天收到好几个粉丝的催更私信，好的吧，实在对不住大家期待的眼神，看样子不能再拖啦，想想写啥好呢，大家咨询比较多的，混合模型算一个，今天就继续给大家写写混合模型如何做吧。

混合模型一般都可以用lme4这个包解决，lme4既可以做线性混合模型，也可以做广义线性混合模型还可以做非线性混合模型，大家有需要可以只研究这一个包就行。

所谓混合模型就是既有固定效应又有随机效应的模型：

“mixedeffects”, denotes a model that incorporates both fixed- and random-effects terms in a linear predictor expression from which the conditional mean of the response can be evaluated

第一部分 线性混合模型

直接上例子，数据是来自一篇研究睡眠剥夺的文献，整个数据大概长下图这样，其中我们的受试者在day0的时候可以睡到自然醒，在之后的日子里所有的受试者就只能睡3个小时了，我们的响应变量是Reaction，就是对受试者做的测验的响应时间，我现在关心睡眠剥夺后，响应时间的变化情况：

对于这么一个纵向数据，

我们来捋一捋：我们只有18个人受试者，每个受试者随访10次，我们需要明白的是，此时我们的每一次测量是嵌套在人的水平上的，我们可以认为，不同人自己的10次测量是有强烈的相关性的，而不同人之间的这种关系又不一定是相同的。

直观一点，我们可以画出来每一天所有人响应时间和睡眠剥夺的变化，画出来就是下图：

可以看到我们上面的这个大图是由很多个小图组成的，每一个小图中横轴就是睡眠剥夺的时间，纵轴是反应时间。每个小图就代表着我们要研究的睡眠剥夺和反应时间的关系（具体到人），但是我们也应该注意到这种关系在不同的人上是不同的，体现在：关系的斜率不同和截距不同。（这个关系的不同可以很明显的在图中看出来）

所以我们就可以拟合一个带有随机效应的混合模型：

fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)

运行代码后得到下面的结果：

结果中有随机效应的标准差和固定效应的β估计，我得到的截距是251.4，斜率是10.5，这两个系数就是我们研究的总体关系的表示，通常需要在文献中汇报，就意味着睡眠不剥夺的时候人的反应时间是251.4，而睡眠每剥夺多一天反应时间增长10.5。

上面这个是最简单的混合模型。我们继续看：

lme4包高水平设置介绍

• 混合模型公式

对于一个常见的混合模型，我们可以在lme4包中写出来如下差不多的混合模型公式：

resp ~ FEexpr + (REexpr1 | factor1) + (REexpr2 | factor2) + ...

这个公式中FEexpr就是固定效应，(REexpr1 | factor1)and(REexpr2 | factor2)都是随机效应，理论上你可以弄很多个随机效应但是实际操作中我们不会关心那么多。

• 理解混合模型公式

我们看到每一个随机效应在公式中的表达都是(expr | factor)这样的。竖杠前面的expr就是一个常规的回归公式，竖杠后面的factor就是一个常规的因子，你可以把竖杠想象成回归公式和因子的交互：

One way to think about the vertical bar operator is as a special kind of interaction betweenthe model matrix and the grouping factor。This interaction ensures that the columns of themodel matrix have different effects for each level of the grouping factor.

这种交互的意思就是在因子的不同水平，我们的回归是不一样的，这也正好和我们前面的解释相对应，就是在不同的人的水平睡眠剥夺和响应时间的关系不一样。

写到这，希望大家能记住下面这张表：

这个表就给我们展示了常见的随机效应的设置，比如(1 | g)，就是说在因子g的不同水平，我们响应变量的截距都不一样。表中的第二行有个offset，表示没有固定效应。如果我们的数据是一个三层嵌套数据，我们可以用第三行的设定来表示随机截距；如果你的数据没有直接嵌套但是在g1和g2的不同水平上存在相关，那么可以用第四行的设定，这个在项目反应理论中比较常见。

在lme4中，默认认为同一个模型的截距和斜率是存在相关的，如果你确定截距和斜率无关那么设定随机效应的时候就可以用两个竖杠，或者把截距和斜率分开来写，就是说(x || g)和x +(1 | g) + (0 + x | g)表达的随机效应都是一样的。

比如如果我认为睡眠剥夺和反应时间随机效应的截距和斜率无关，我便可以做如下设定：

fm2 <- lmer(Reaction ~ Days + (Days || Subject), sleepstudy)#截距和斜率无关的设定

有时候我们拟合一个后又想尝试对模型进行改变，但又不想重写，此时就可以直接对相似的模型基础上进行更新：

• 模型的更新

比如我想在fm1的基础上去掉随机斜率只留随机截距，我就可以用updata写出如下代码：

fm3 <- update(fm1, . ~ . - (Days | Subject) + (1 | Subject))#模型的更新

到底哪一个模型更好呢？

可以用anova方法进行模型间的比较：

anova(fm1, fm2, fm3)

运行代码会输出比较的结果：

其中，从模型比较的结果可以看出，给模型增加一个截距和斜率无关的随机效应相比会使得模型的deviance变小，进一步将随机效应设定为相关，并不能够显著地减小deviance，从而我们就可以知道fm2才是对数据拟合最好的模型。

第二部分 非线性混合模型

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

The fixed-effects parameters describe the general patterns of the data and random-effects parameters describe specific clusters. If the model is nonlinear in the parameters,it is called a nonlinear mixed-effects model (Davidian &Giltinan, 2003)

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

看下面这个图，这个图描绘了不同人用了茶碱过后的反应，时间是横轴，残留是纵轴，和开篇线性模型中睡眠剥夺和反应时间的例子一样，我们把每个人的关系都做了图出来，不过从图中可以明显看出这种关系并不是简单线性的。

其实这种不是线性的关系存在的情况很多。

比如渐进回归：

再比如逻辑增长：

此时我们要注意到像这两非线性关系模型的参数都不是简单的一个斜率加个截距了。都有φ1，φ2，φ3三个额外参数。

这儿先给大家写一个逻辑增长的实际例子：我现在有一个关于树木周径的数据集，每棵树随访了7次，每次随访记录树的年龄age，和周径，我现在想研究在所有树木中时间和周径的关系。

很自然，我们可以想到不同的树这个关系应该是不一样的，我们想探求的一定是考虑了树水平的变异之后的总体关系，所以不妨先画出来每个树的关系：

从图中可以看到我们总共有5棵树，基本关系是一致的，但存在些许变异相关（所以考虑混合模型），而且这个关系并不是线性的（时间越大周径基本不改变了），所以我们应该考虑非线性的混合模型。

具体地，我们可以用nlmer方法来拟合非线性混合模型，方法参数包括3部分：首先是响应变量，然后是非线性函数，然后是混合效应公式：

The formula argument fornlmeris in three parts: the response, the nonlinear model function depending on covariates and a set of nonlinear model (nm) parameters, and the mixed-effects formula.

比如对我们的数据我就可以写出如下SSlogis方法的代码：

print(nm1 <- nlmer(circumference ~ SSlogis(age, Asym, xmid, scal) ~ Asym | Tree,
Orange, start = c(Asym = 200, xmid = 770, scal = 120)), corr = FALSE)

此时我们选择的非线性函数是逻辑增长函数SSlogis，刚刚给大家解释了这个函数是有3个参数的，在上面的代码中，age是我们的预测变量，Asym, xmid, scal分别是额外的三个参数（之前的逻辑增长的式子和Asym/(1+exp((xmid-input)/scal))等同）：

进一步，拟合逻辑增长是要我们给出这些参数的初始值的，然后从初始值通过梯度下降寻找各个参数的最优解：

SSlogis has an attribute called "initial", which is a function that nls can call to compute reasonable starting values for fitting a logistic function to the input data.

所以我们看到代码中都给出了响应参数的初始值。

运行上面代码后输出如下结果：

我们可以看到结果中的固定效应里面有Asym，Xmid，scal参数的估计结果。

那么这些参数如何解释呢？

留个悬念我们下期再更。关注关注关注，嘿嘿。

小结