多层回归模型简介

多层回归模型简介

多层回归模型（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