转自个人微信公粽号【易学统计】的统计学习笔记：R语言：多水平统计模型
01 解决何种问题
- 同样是九年义务教育，凭什么别人那么优秀？显然这跟每个人，不同班级，不同学校有关系，究竟是什么样的关系呢？在临床研究中，研究成都居民和上海居民的糖尿病患病的影响因素。显然成都市民饮食偏向咸辣，上海市民饮食偏清淡，这对糖尿病的危险因素是有影响的。
- 除此以外还有上篇文章中提到的三个案例，如多次测量结局以比较两种治疗方式的治疗结果等，详细案例描述读者可参见上篇文章。
- 这些案例具有层次数据结构，两水平层次结构是第一水平在第二水平中聚集有一定相似性。如第一段不同小区为第二水平，居民为第一水平，第二段则是不同时间点测量为第一水平，受测量个体为第二水平。面对这种问题，在前两篇文章的广义估计方程和重复测量方差分析能解决，这里介绍第三种解决此类问题的模型，也是普遍应用的模型，即多水平统计模型。
02 方法说明
- 多水平统计模型可同时用于分析多层数据、多级数据、嵌套数据的方法，也叫随机效应模型。有人会问上海居民和成都居民的数据能否混合使用来研究影响因素，显然不合理，上海和成都居民存在组间差异，即第二水平差异。这两个小区的居民也存在很多相似性，不能算作独立的个体，即组内差异。
- 其常见的模型是方差成分模型和随机系数模型。它俩的区别在于随机系数模型中协变量X的系数β1j不是固定的而是随机的，即协变量对反应变量的效应在不同的第二水平间是不同的。
下面讲述该模型在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 结果解读
- 在分析之前，我们可以作图看看不同受试者的控制睡眠的天数与反应时间的关系，很明显，不同受试者控制睡眠的天数对反应时间的影响是不同的。有的呈现正相关，有的负相关，直线的截距和斜率都不一样。
- R中加载lme4包，用该包中的lmer()函数建模。
- ICC=1278/(1278+1959)=0.395，说明个体间有一定的相关性,也可以用icc()函数计算。
接下来我们加入第一水平的变量到模型中。
06R代码
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, 18
Fixed 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结果解读
- 用lmer()函数建模，其他同线性回归，一般变量表示固定效应，注意(1||Subject)写法，括号内竖线右侧的Subject表示它是一个随机效应，1表示只影响模型截距，即不同的Subject其控制睡眠天数对反应时间的效应是一样的。若1换成其他变量表示，那么表示还影响模型斜率。
- 参数REML表示限制性最大似然估计，基于全残差项，包含所有变异，为默认估计方法。
- 返回的结果中固定效应模块中，对于常数项截距，当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)
##fmod: Reaction ~ Days + (1 | Subject)
##fmod1: Reaction ~ Days + (Days | Subject)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
##fmod 4 1802.1 1814.8 -897.04 1794.1
##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结果解读
- 增加随机斜率的变量，为了判断斜率是否随着随机效应发生变化，我们可以采用方差分析anova()，加入Days后模型更优，且具有显著性差异，还可以采用ranova()函数判断随机效应项是否是显著的。表明，不同受试者之间控制睡眠的天数对反应时间的影响具有显著差异的。这个跟上面直观图片能看出来。
- 值得注意的是，如果不确定随机变量只影响截距还是会影响斜率可以采用anova()函数对比两个不同的模型，根据P值确定最优模型，且竖线之前的变量需要是数值型变量，得到最优模型后，可以按照前面的解释方法对结果进行解释。
- 还可采用之前提到的重复测量方差分析来比较，可翻看前面操作实现。
10总结
多水平统计模型既能处理连续型结局变量，又能处理分类型结局变量。又很广泛的应用，能处理常规回归模型所不能处理的非独立数据。相比于重复测量方差分析，多水平模型要求在第一水平上的观察点个数可以不等。
得到的模型，我们可以用各种泛型函数，如summary，predict，resid进行进一步处理，可考察残差齐性。
如果您觉得有用，请点赞，转发哦~
更多统计小知识，请关看 公粽号 易学统计
更多阅读
R语言：广义估计方程(GEE)
R语言|两因素重复测量方差分析
R语言|基于Cox模型pec包深度验证
R语言|中位生存时间列线图绘制
R语言|Cox模型校准度曲线绘制
R语言|中位生存时间列线图绘制
基于Lasso回归筛选变量构建Cox模型并绘制Nomogram
R语言Logistic回归模型验证及Nomogram绘制
如何进行高维变量筛选和特征选择(一)？Lasso回归
