星期四, 十二月 08, 2011

多层回归模型简介

多层回归模型(Multi-level model)中有很多容易混淆的概念,因为很多概念是来源于不同的专业背景。首先让我们先罗列这些名词进行区分,再来R语言来举例。

多层回归模型通常涉及到对同一个体进行反复测量,这样得到的数据就不再相互独立而是存在某种相关性,所以普通线性回归不再适用。当这种反复测量是在不同时点上进行时,这就称为面板数据分析(panel data analysis)或者纵向数据分析(longitudinal data analysis)。


Fixed Effect:固定效应,也就是普通线性回归中处理的预测变量,它在模型中对响应变量的期望造成影响,其因子水平包含明确的信息,通常是在实验中加以控制的因素。

Random Effect:随机效应,也就是一个随机变量,估计它是没有意义的。它在模型中对响应变量的方差造成影响,其因子水平包含信息不清楚,通常在实验中无法控制的因素,例如在面板数据分析中,不同的个体差异就是随机效应。

Mix effect model:混合效应模型,模型中包含了fixed effect和random effect,根据random effect的影响,又区别为对模型截距的影响(random intercept)和对模型截距与斜率的影响(random intercept and slope)。

假设一个面板数据中个体之间存在差异,那么它是一个random effect。它的影响分两种情况,一种是只影响截距,即是它与模型中其它预测变量相加,如果对不同个体分别作线性回归,那么得到回归线截距会不同,但回归线平行,此时又称为固定效应回归模型(Fixed Effects Regression Model)。另一种同时影响模型截距与斜率,是指它还与其它变量相乘,那么分别得到的回归线截距和斜率都不同,此时又称为随机效应模型(Random Effects Regression Model)

如果模型中不考虑random effect,也就是认为个体和时间都没有显著性差异,此时模型退化为混合估计模型(Pooled Regression Model),此时可以直接把面板数据混合在一起用普通线性回归估计参数。

下面我们以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语言中我们使用lme4包中的lmer函数来完成这项工作。首先载入faraway包以便读取psid数据集,然后加载mgcv包,再将年份数据中心化以方便解释模型,最后用lmer函数进行建模。
-----------------
library(faraway)
library(lme4)
psid$cyear <- psid$year-78 
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:
model2: log(income) ~ cyear * sex + age + educ + (1 | person)
model1: log(income) ~ cyear * sex + age + educ + (cyear | person)
       Df    AIC    BIC  logLik  Chisq Chi Df
model2  8 3943.9 3987.2 -1963.9              
model1 10 3805.6 3859.7 -1892.8 142.27      2
       Pr(>Chisq)    
model2               
model1  < 2.2e-16 ***
-----------------

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

参考资料:
环境与生态统计--R语言的应用
The R Book
A_Handbook_of_Statistical_Analyses_Using_R
Extending_the_Linear_Model_with_R__Generalized_Linear__Mixed_Effects_and_Nonparametric_Regression_Model

1 条评论:

  1. “另一种同时影响模型截距与斜率,是指它还与其它变量相乘,那么分别得到的回归线截距和斜率都不同,此时又称为随机效应模型(Random Effects Regression Model)。”这里有点疑问。随机效应应该是指individual和regressor不相关。

    回复删除