R语言实现混合模型

普通的线性回归只包含两项影响因素，即固定效应（fixed-effect）和噪声（noise）。噪声是我们模型中没有考虑的随机因素。而固定效应是那些可预测因素，而且能完整的划分总体。例如模型中的性别变量，我们清楚只有两种性别，而且理解这种变量的变化对结果的影响。

那么为什么需要 Mixed-effect Model？因为有些现实的复杂数据是普通线性回归是处理不了的。例如我们对一些人群进行重复测量，此时存在两种随机因素会影响模型，一种是对某个人重复测试而形成的随机噪声，另一种是因为人和人不同而形成的随机效应（random effect）。如果将一个人的测量数据看作一个组，随机因素就包括了组内随机因素（noise）和组间随机因素（random effect）。这种嵌套的随机因素结构违反了普通线性回归的假设条件。

你可能会把人员（组间的随机效应）看作是一种分类变量放到普通线性回归模型中，但这样作是得不偿失的。有可能这个factor的level很多，可能会用去很多自由度。更重要的是，这样作没什么意义。因为人员ID和性别不一样，我们不清楚它的意义，而且它也不能完整的划分总体。也就是说样本数据中的路人甲，路人乙不能完全代表总体的人员ID。因为它是随机的，我们并不关心它的作用，只是因为它会影响到模型，所以不得不考虑它。因此对于随机效应我们只估计其方差，不估计其回归系数。

混合模型中包括了固定效应和随机效应，而随机效应有两种方式来影响模型，一种是对截距影响，一种是对某个固定效应的斜率影响。前者称为 Random intercept model，后者称为 Random Intercept and Slope Model。Random intercept model的函数结构如下

Yij = a0 + a1*Xij + bi + eij

a0: 固定截距 a1: 固定斜率 b: 随机效应（只影响截距） X: 固定效应 e: 噪声

混合线性模型有时又称为多水平线性模型或层次结构线性模型由两个部分来决定，固定效应部分+随机效应部分，

二、R语言中的线性混合模型可用包

1、nlme包

这是一个比较成熟的R包，是R语言安装时默认的包，它除了可以分析分层的线性混合模型，也可以处理非线性模型。在优势方面，个人认为它可以处理相对复杂的线性和非线性模型，可以定义方差协方差结构，可以在广义线性模型中定义几种分布函数和连接函数。

它的短板：随机效应的定义过于呆板；数据量大时速度很慢；默认情况下不能处理系谱数据；不能处理多元数据。

2、lme4包

lme4包是由Douglas Bates开发，他也是nlme包的作者之一，相对于nlme包而言，它的运行速度快一点，对于睡觉效应·随机效应的结构也可以更复杂一点，但是它的缺点也和nlme一样：不能处理协方差和相关系数结构；它可以与构建系数的包连接，比如mmpedigree包，但是结合比较脆弱。

3、ASReml-R包

ASReml-R是ASReml的R版本，它的优点：可以处理复杂的随机因子结构；可以处理多元数据；可以处理系谱数据；可以处理大批量的数据

主要的缺点：它是收费的，当然它对于不发达国家的科研机构是免费的，不过需要申请和被审核；它的用户主要是育种公司、科研机构等，它可以在各种平台上运行，包括Windows、Linux、OS X等。

二、多水平模型案例分析

案例一：

1、首先导入数据，查看一下数据的结构

数据来源：一个传统的裂区数据来说明不同软件包的用法，这个数据oats是在MASS包中，是研究大麦品种和N肥处理的裂区试验，其中品种为主区，肥料为裂区。

library(MASS)

data(oats)

names(oats) = c('block', 'variety', 'nitrogen', 'yield')

head(oats)

    block   variety nitrogen yield
1     I     Victory   0.0cwt   111
2     I     Victory   0.2cwt   130
3     I     Victory   0.4cwt   157
4     I     Victory   0.6cwt   174
5     I Golden.rain   0.0cwt   117
6     I Golden.rain   0.2cwt   114

oats$mainplot=oats$varietyoats$subplot=oats$nitrogen

2、nlme包：用这个包很简单，y-变量写在左边，然后是固定因子，然后是随机因子，注意1|block/mainplot是裂区试验残差的写法，因为里面有两个残差。代码如下：

library(lme4)

library(nlme)m1.nlme = lme(yield ~ variety*nitrogen,random = ~ 1|block/mainplot,data = oats)summary(m1.nlme)

方差分析结果为：anova(m1.nlme)                 numDF denDF   F-value p-value
(Intercept)          1    45 245.14299  <.0001
variety              2    10   1.48534  0.2724
nitrogen             3    45  37.68562  <.0001
variety:nitrogen     6    45   0.30282  0.9322

3、lme4包

lme4包的语法也相似，随机效应有着和nlme相同的语法，不同的是lme4包它的结果给出了随机效应的标准差，而不是方差。library(lme4)m1.lme4 = lmer(yield ~ variety*nitrogen + (1|block/mainplot),data = oats)summary(m1.lme4)anova(m1.lme4) 

Analysis of Variance Table of type III  with  Satterthwaite 
approximation for degrees of freedom

                  Sum Sq Mean Sq NumDF DenDF F.value    Pr(>F)    
variety            526.1   263.0     2    10   1.485    0.2724    
nitrogen         20020.5  6673.5     3    45  37.686 2.458e-12 ***
variety:nitrogen   321.8    53.6     6    45   0.303    0.9322    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

案例二：利用nlme包

下面我们以faraway包中的psid数据来举例，该数据集是对美国人的收入情况进行调查所得到的，其中包括了年龄、教育、性别、时间和个体ID这几个变量，我们希望了解这些因素对收入的影响。

age educ sex income year person 1 31 12 M 6000 68 1 2 31 12 M 5300 69 1 3 31 12 M 5200 70 1 4 31 12 M 6900 71 1 5 31 12 M 7500 72 1 6 31 12 M 8000 73 1

如果假设认为这些调查对象是同质的，也就是个体间没有差异性，那么可以将数据完全汇集（complete pooling）到一起，直接利用lm函数进行回归。但这个混合效应模型的同质假设往往不成立，数据汇集导致过度简化。另一种思路是假设研究的异质性，将不同的个体分别进行回归，从而得到针对特定个体的估计值，这称为不汇集（no pooling）。但这种方法导致每个回归所用到的样本减少，从而难以估计统计量的标准差。

多层回归模型的思路是前两者的折中，所以又称为部分汇集（partial pooling）。在R语言中我们使用mgcv包中的lmer函数来完成这项工作。首先载入faraway包以便读取psid数据集，然后加载mgcv包，再将年份数据中心化以方便解释模型，最后用lmer函数进行建模。

install.packages("faraway")library(faraway)library(nlme)  age educ sex income year person cyear
1  31   12   M   6000   68      1   -10
2  31   12   M   5300   69      1    -9
3  31   12   M   5200   70      1    -8
4  31   12   M   6900   71      1    -7
5  31   12   M   7500   72      1    -6
6  31   12   M   8000   73      1    -5

library(mgcv)psid$cyear <- psid$year-78head(psid)model1=lmer(log(income) ~ cyear*sex +age+educ+(cyear|person),psid) 

lmer函数使用和lm是类似的，一般变量表示固定效应，括号内竖线右侧的person表示它是一个随机效应，它与模型中其它变量相加，而且与年份cyear变量相乘，影响其斜率。这就是一个随机效应模型。如果认为随机效应只影响模型截距，那么固定效应回归模型可以用下面的公式

model2=lmer(log(income) ~ cyear*sex +age+educ+(1|person),psid)


为了判断哪一个模型更为合适，可以使用anova函数，从结果中观察到P值很小，判断应当使用model1 

anova(model1,model2)

Data: psid
Models:
..1: log(income) ~ cyear * sex + age + educ + (1 | person)
object: log(income) ~ cyear * sex + age + educ + (cyear | person)
       Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)    
..1     8 3943.8 3987.1 -1963.9   3927.8                            
object 10 3805.5 3859.6 -1892.7   3785.5 142.3      2  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

得到的模型结果还可以用各种泛型函数如summary、predict、resid进行进一步处理。当响应变量不符合正态分布假设时，还可以用广义多层回归进行(glmer)建模

案例三：

1、Generate a longitudinal dataset and convert it into long format

library(MASS) dat.tx.a <- mvrnorm(n=250, mu=c(30, 20, 28), Sigma=matrix(c(25.0, 17.5, 12.3, 17.5, 25.0, 17.5, 12.3, 17.5, 25.0), nrow=3, byrow=TRUE)) dat.tx.b <- mvrnorm(n=250, mu=c(30, 20, 22), Sigma=matrix(c(25.0, 17.5, 12.3, 17.5, 25.0, 17.5, 12.3, 17.5, 25.0), nrow=3, byrow=TRUE))

dat <- data.frame(rbind(dat.tx.a, dat.tx.b)) names(dat) <- c('measure.1', 'measure.2', 'measure.3')

根据指定的均值和协方差生成多元正态数据：MASS包中的mvrnorm()函数 mvrnorm(n,mean,sigma) measure.1 measure.2 measure.3 1 25.31761 20.89468 34.65525 2 19.30890 22.82886 33.24505 3 25.53250 16.27554 24.61767 4 27.52445 23.07668 32.62275 5 35.89934 28.24790 36.07344 6 25.71556 13.36878 26.69632

dat <- data.frame(subject.id = factor(1:500), tx = rep(c('A', 'B'), each = 250), dat)

subject.id tx measure.1 measure.2 measure.3
1          1  A  25.31761  20.89468  34.65525
2          2  A  19.30890  22.82886  33.24505
3          3  A  25.53250  16.27554  24.61767
4          4  A  27.52445  23.07668  32.62275
5          5  A  35.89934  28.24790  36.07344
6          6  A  25.71556  13.36878  26.69632

rm(dat.tx.a, dat.tx.b) dat <- reshape(dat,varying = c('measure.1','measure.2','measure.3'),idvar = 'subject.id', direction = 'long') subject.id tx time measure 1.1 1 A 1 25.31761 2.1 2 A 1 19.30890 3.1 3 A 1 25.53250 4.1 4 A 1 27.52445 5.1 5 A 1 35.89934 6.1 6 A 1 25.71556 subject.id tx time measure 495.3 495 B 3 29.14236 496.3 496 B 3 17.03742 497.3 497 B 3 18.75601 498.3 498 B 3 20.80353 499.3 499 B 3 31.94328 500.3 500 B 3 25.78684

2、Multilevel growth models 有很多R包可以做多水平分析，其中 lme4对于一般模型（二分类及离散输出）比较适合，

另外一个nlme包 比较适合连续输出变量（正态或高斯分布）

install.packages('lme4') library(Matrix) library(lme4)

m <- lmer(measure ~ time + (1 | subject.id), data = dat)

summary(m)

注：符号 DV ~ IV + (1 | rand.int) ，其中 DV 为输出变量，IV 为自变量, 1 为自变量的系数或斜率,

rand.int 为随机截距变量

Likewise, a random slopes model is specified using the syntax DV ~ IV + (rand.slope | rand.int).

结果显示如下：

Linear mixed model fit by REML ['lmerMod'] Formula: measure ~ time + (1 | subject.id) Data: dat REML criterion at convergence: 9750.6 Scaled residuals: Min 1Q Median 3Q Max -2.41363 -0.69275 0.04628 0.69361 2.47612 Random effects:

Groups Name Variance Std.Dev. subject.id (Intercept) 8.136 2.852 Residual 32.268 5.681 Number of obs: 1500, groups: subject.id, 500 Fixed effects: Estimate Std. Error t value (Intercept) 29.2715 0.4085 71.66 time -2.2832 0.1796 -12.71 Correlation of Fixed Effects: (Intr) time -0.880

3、 Multilevel growth models with approximate p-values

install.packages('lmerTest') library(lmerTest) m <- lmer(measure ~ time + (1 | subject.id), data = dat) summary(m)

结果是一样的，不过多了P值

4、Calculating 95% CI and PI

dat.new <- data.frame(time = 1:3) dat.new$measure <- predict(m, dat.new, re.form = NA) m.mat <- model.matrix(terms(m), dat.new) dat.new$var <- diag(m.mat %*% vcov(m) %*% t(m.mat)) + VarCorr(m)$subject.id[1] dat.new$pvar <- dat.new$var + sigma(m)^2 dat.new$ci.lb <- with(dat.new, measure - 1.96 * sqrt(var)) dat.new$ci.ub <- with(dat.new, measure + 1.96 * sqrt(var)) dat.new$pi.lb <- with(dat.new, measure - 1.96 * sqrt(pvar)) dat.new$pi.ub <- with(dat.new, measure + 1.96 * sqrt(pvar))

5、Plot the average values

library(ggplot2) ggplot(data = dat.new, aes(x = time, y = measure)) + geom_line(data = dat, alpha = .02, aes(group = subject.id)) + geom_errorbar(width = .02, colour = 'red', aes(x = time - .02, ymax = ci.ub, ymin = ci.lb)) + geom_line(colour = 'red', linetype = 'dashed', aes(x = time-.02)) + geom_point(size = 3.5, colour = 'red', fill = 'white', aes(x = time - .02)) + geom_errorbar(width = .02, colour = 'blue', aes(x = time + .02, ymax = pi.ub, ymin = pi.lb)) + geom_line(colour = 'blue', linetype = 'dashed', aes(x = time + .02)) + geom_point(size = 3.5, colour = 'blue', fill = 'white', aes(x = time + .02)) + theme_bw()