R语言：多水平统计模型
01
解决何种问题
同样是九年义务教育，凭什么别人那么优秀？显然这跟每个人，不同班级，不同学校有关系，究竟是什么样的关系呢？在临床研究中，研究成都居民和上海居民的糖尿病患病的影响因素。显然成都市民饮食偏向咸辣，上海市民饮食偏清淡，这对糖尿病的危险因素是有影响的。
除此以外还有上篇文章中提到的三个案例，如多次测量结局以比较两种治疗方式的治疗结果等，详细案例描述读者可参见上篇文章。
这些案例具有层次数据结构，两水平层次结构是第一水平在第二水平中聚集有一定相似性。如第一段不同小区为第二水平，居民为第一水平，第二段则是不同时间点测量为第一水平，受测量个体为第二水平。面对这种问题，在前两篇文章的广义估计方程和重复测量方差分析能解决，这里介绍第三种解决此类问题的模型，也是普遍应用的模型，即多水平统计模型。
02
方法说明
多水平统计模型可同时用于分析多层数据、多级数据、嵌套数据的方法，也叫随机效应模型。有人会问上海居民和成都居民的数据能否混合使用来研究影响因素，显然不合理，上海和成都居民存在组间差异，即第二水平差异。这两个小区的居民也存在很多相似性，不能算作独立的个体，即组内差异。
其常见的模型是方差成分模型和随机系数模型。它俩的区别在于随机系数模型中协变量X的系数β1j不是固定的而是随机的，即协变量对反应变量的效应在不同的第二水平间是不同的。
方差成分模型基本形式：
注意一个指标组内相关ICC，它反应了第二水平内个体间的相关程度，即水平1在水平2中的聚集度或者相似性。若该值为0则说明数据不具备层次结构，可忽略第二水平，简化为传统单水平模型，反之则不能忽略。
随机系数模型基本形式：
完整模型(水平1和水平2上都有解释变量)：
其中W1j为第二层的解释变量(可以有多个)，可以在零模型(不含任何协变量)和完整模型之间，根据研究目的，设置不同的随机成分和固定成分，构建一系列模型。该模型涉及到的参数解释，读者可参考孙振球主编的医学统计学书籍。
该模型的估计过程，分两步。
第一：分别在每个第二水平中进行第一水平变量的回归运算，将得到j个第一水平模型的截距和系数。
第二：将第一步所得到的系数作为因变量，将第二水平的某些变量作为自变量进行回归拟合，将得到第二水平模型。
具体的建模分析过程：
第一步拟合常数项模型，计算ICC。
第二步将第二水平的变量加入模型。
第三步将第一水平的变量加入。
第四步确定第一水平变量的斜率是固定的还是随机的。
第五步计入交互作用项。
下面讲述该模型在R中的操作。
03
加载数据
数据集说明：这是R自带的控制睡眠时间通过一系列测试研究反应时间的数据集，Reaction是平均反应时间(ms)，Days是控制睡眠的天数，Subject是受试者编号。
04
R代码
library(lme4) #这个是模型包，先加载library(ggplot2)data("sleepstudy") ##调用数据集ggplot(data=sleepstudy,aes(x=Days,y=Reaction,color=Subject,group = as.factor(Subject)))+geom_point()+geom_smooth(method=lm,se=FALSE, col="black", size=.5, alpha=.8)
mod <- lmer(Reaction~1+(1|Subject),data=sleepstudy)summary(mod)library(sjstats)icc(mod) ##0.395##返回结果##Random effects:## Groups Name Variance Std.Dev.## Subject (Intercept) 1278 35.75## Residual 1959 44.26##Number of obs: 180, groups: Subject, 18##Fixed effects:## Estimate Std. Error df t value Pr(>|t|)##(Intercept) 298.51 9.05 17.00 32.98 <2e-16 ***
05
结果解读
1.在分析之前，我们可以作图看看不同受试者的控制睡眠的天数与反应时间的关系，很明显，不同受试者控制睡眠的天数对反应时间的影响是不同的。有的呈现正相关，有的负相关，直线的截距和斜率都不一样。
2.R中加载lme4包，用该包中的lmer()函数建模。
3..ICC=1278/(1278+1959)=0.395，说明个体间有一定的相关性,也可以用icc()函数计算。
接下来我们加入第一水平的变量到模型中。
06
R代码
fmod <- lmer(Reaction~Days+(1|Subject),sleepstudy,REML = T)summary(fmod) ##REML 采用限制性最大似然法,Random effects: 随机效应## Groups Name Variance Std.Dev.## Subject (Intercept) 1378.2 37.12## Residual 960.5 30.99##Number of obs: 180, groups: Subject, 18Fixed effects: 固定效应##Estimate Std. Error df t value Pr(>|t|)##(Intercept) 251.4051 9.7467 22.8102 25.79 <2e-16 ***##Days 10.4673 0.8042 161.0000 13.02 <2e-16 ***
07
结果解读
1.用lmer()函数建模，其他同线性回归，一般变量表示固定效应，注意(1||Subject)写法，括号内竖线右侧的Subject表示它是一个随机效应，1表示只影响模型截距，即不同的Subject其控制睡眠天数对反应时间的效应是一样的。若1换成其他变量表示，那么表示还影响模型斜率。
2.参数REML表示限制性最大似然估计，基于全残差项，包含所有变异，为默认估计方法。
3.返回的结果中固定效应模块中，对于常数项截距，当Days为0，其反应时间为251.4051，控制睡眠的天数每增加1天，反应时间增加10.4673ms，控制睡眠天数与反应时间显著相关。
08
R代码
fmod1 <- lmer(Reaction~Days+(Days|Subject),sleepstudy,REML = T)summary(fmod1)##返回结果Random effects:## Groups Name Variance Std.Dev. Corr## Subject (Intercept) 611.90 24.737## Days 35.08 5.923 0.07## Residual 654.94 25.592##Number of obs: 180, groups: Subject, 18##Fixed effects:## Estimate Std. Error df t value Pr(>|t|)##(Intercept) 251.405 6.824 17.005 36.843 < 2e-16 ***##Days 10.467 1.546 16.995 6.771 3.27e-06 ***##模型比较anova(fmod,fmod1)###fmod1: Reaction ~ Days + (Days | Subject)## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)###fmod1 6 1763.9 1783.1 -875.97 1751.9 42.139 2 7.072e-10 ***
##检验变量d 随机效应是否显著ranova(fmod1)##ANOVA-like table for random-effects: Single term deletions##Model:##Reaction ~ Days + (Days | Subject)## npar logLik AIC LRT Df Pr(>Chisq)##<none> 6 -871.81 1755.6##Days in (Days | Subject) 4 -893.23 1794.5 42.837 2 4.99e-10 ***
09
结果解读
1.增加随机斜率的变量，为了判断斜率是否随着随机效应发生变化，我们可以采用方差分析anova()，加入Days后模型更优，且具有显著性差异，还可以采用ranova()函数判断随机效应项是否是显著的。表明，不同受试者之间控制睡眠的天数对反应时间的影响具有显著差异的。这个跟上面直观图片能看出来。
2.值得注意的是，如果不确定随机变量只影响截距还是会影响斜率可以采用anova()函数对比两个不同的模型，根据P值确定最优模型，且竖线之前的变量需要是数值型变量，得到最优模型后，可以按照前面的解释方法对结果进行解释。
3.还可采用之前提到的重复测量方差分析来比较，可翻看前面操作实现。
10
总结
多水平统计模型既能处理连续型结局变量，又能处理分类型结局变量。又很广泛的应用，能处理常规回归模型所不能处理的非独立数据。相比于重复测量方差分析，多水平模型要求在第一水平上的观察点个数可以不等。
得到的模型，我们可以用各种泛型函数，如summary，predict，resid进行进一步处理，可考察残差齐性。
11
参考文献
孙振球等. 医学统计学. 人民卫生出版社。
余松林.混合线性模型的应用[J].中国医院统计,2006,13(1):70-75
更多阅读：