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()


原文发布于微信公众号 - 大数据挖掘DT数据分析(datadw)

原文发表时间:2016-10-03

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏机器人网

智能机器人语音识别技术

语音控制的基础就是语音识别技术,可以是特定人或者非特定人的。非特定人的应用更为广泛,对于用户而言不用训练,因此也更加方便。语音识别可以分为孤立词识别,连接词识别...

4375
来自专栏AI科技评论

Ian Goodfellow & ICLR-17最佳论文得主新作:验证与测试,机器学习的两大挑战

AI 科技评论按:去年年底,Ian Goodfellow与Nicolas Papernot(是的,就是ICLR 2017的最佳论文得主之一)合作开了一个博客叫c...

2726
来自专栏量子位

DeepMind的脑补AI再获新技能:看文字知场景、复杂环境、连续视频……

712
来自专栏IT派

如何在机器学习竞赛中更胜一筹?

机器学习很复杂。你可能会遇到一个令你无从下手的数据集,特别是当你处于机器学习的初期。 在这篇文章中,你将学到一些基本的关于建立机器学习模型的技巧,大多数人都从中...

3257
来自专栏AI科技评论

学界 | 传播动态学的主动监控:一种组稀疏贝叶斯学习方法

AI 科技评论按:本文作者吉林大学博士生裴红斌,本文为对他发表在 AAAI 2018 论文的独家解读稿件,未经许可不得转载。 Group Sparse Bay...

2886
来自专栏大数据文摘

是猫还是屏幕?OpenAI实力怼“神经网络多角度图像识别难被欺骗”说法

1022
来自专栏PPV课数据科学社区

AI时代就业指南:机器学习工程师求职须知

什么是机器学习? 我们来看一下机器学习是做什么的,能解决什么问题。 首先我们来看机器学习的一个类型,监督学习。 蓝色箭头部分是训练一个机器学习模型的过程。首先有...

3307
来自专栏量子位

英伟达新研究:“狗生猫,猫生万物”的多模态无监督图像转换

这是英伟达最新创造的一项技术。在最近发布的论文Multimodal Unsupervised Image-to-Image Translation中,研究人员提...

1083
来自专栏机器学习算法与Python学习

数据科学家必用的25个深度学习的开放数据集!

原文:https://www.analyticsvidhya.com/blog/2018/03/comprehensive-collection-deep-le...

46014
来自专栏携程技术中心

干货 | 携程个性化推荐算法实践

作者简介 携程基础业务研发部-数据产品和服务组,专注于个性化推荐、自然语言处理、图像识别等人工智能领域的先进技术在旅游行业的应用研究并落地产生价值。目前,团队已...

4625

扫码关注云+社区