【原】R数据分析：嵌套数据分析为什么要用加随机效应？终于解释清楚了

自己已经写了好几篇关于随机效应模型的文章了，今天换个角度，从传统回归和随机效应模型的对比中模拟出两模型真正的差异，让你明白加上随机效应到底对模型会有什么样的改变。

回顾

在传统的回归中，我们有如下式子：

yj = βxj+ εj

就是说x与y之间的关系是β。

但现在，我们有一个纵向数据，比如说我有30个人，每个人测xy均测量15次，那么有可能会存在每个人中xy的关系不一样的情况。

这个很好理解，毕竟人与人之间是有差异的。那么如果真的是这样，我们对于每个人都可以拟合出一个βi，如果每个人中xy关系的基线水准也不一样，那么可能我们每个人的xy关系上会有一个特定且不一样的截距αi。

然后，我们的式子变成了：

yij = αi + βixij + εij

上面这个式子就是最简单的随机效应模型。其中αi为随机截距，βi为随机斜率。

R模拟

现在我模拟30个人，每个人测量xy15次：

J <- 15
N <- 30

test.df <- data.frame( unit = sort(rep(c(1:N),J)),

J = rep(c(1:J),N) , x = rnorm(n = J*N) )
beta <- 3 + .2*rnorm(N)

test.df$beta <- beta[test.df$unit]
test.df$y <- 1 + test.df$x * test.df$beta + .75*rnorm(n = J*N)
head(test.df, 18)

上面的代码产生了30个人，每个人有xy的测量15次，xy的关系服从以3为均值，0.2为标准差的正态分布，每个人中xy的关系都不一样：

在上面的数据中，我们有两个水平，第一个水平是xy的测量，第二个水平是人，xy是嵌套在人的水平上的。

一般线性回归

因为我们知道每个人xy关系是不一样，所以我们分人做个回归，这么一个做法是一般线性回归对多水平嵌套数据能做到的极限了：

beta.hat <- list()
for(i in 1:N){
unit.lm <- lm(y ~ x, data = subset(test.df, unit == i) )
beta.hat[i] <- coef(unit.lm)[2]
}
beta.hat <- as.numeric(beta.hat)

上面的代码就可以得到每个人中xy的关系βi：

我们看一看我们用一般线性回归估计出来的βi和我们本来模拟的有什么差异：



par(mfrow = c(2, 1))
hist(beta, main = "XY真实的斜率", col = "blue",

xlab = expression(beta[i]), cex.axis = 1.5, cex.lab = 1.5,

breaks = seq(from = 2.4, to = 3.6, by = .1) )

hist(as.numeric(beta.hat), main = "一般线性回归估计的斜率",

xlab = expression(hat(beta)[i]), col = "blue", cex.axis = 1.5,

cex.lab = 1.5, breaks = seq(from = 2.4, to = 3.6, by = .1) )

可以看出来估计的斜率分布变异比真实斜率更大一点。此时，我们并不能说xy的关系到底如何，因为我们拟合了30个β，虽然这个β的分布和真实的分布差不太多（其实变异稍大），我们无法得出真实的xy之间的关系，你说到底30个β到底选哪个？。

再看随机效应模型

在R中建立随机效应模型需要用到lme4这个包，随机效应部分一般表达为：

(formula for random terms | unit for which these terms apply).

在 | 的左边你可以设定随机截距或者随机斜率，在右边需要设定随机效应的水平；如果x有随机截距和随机斜率，你就可以设定左边为“1+x”，如果x只有随机斜率没有随机截距，你设定左边为“0+x”，因为随机效应都是在高水平上的变异，所以在 |右边你就应该将这个水平指定出来，在本例中人为高水平，对应的变量为数据库中的unit，所以我们将右边设定为unit。

那么对于本例的混合效应模型我们可以写出代码：

library(lme4)
re.lm <- lmer(y ~ x + (1+x|unit), data = test.df)
summary(re.lm)

上面的代码拟合了一个有随机截距和随机斜率的混合模型，此时我们得到x的固定效应为3.055，随机效应为0.116，和我们原先设定的xy的关系服从以3为均值，0.2为标准差的正态分布有点接近了。

但是我们原先并没有设定人水平上的截距的变异，大家回看原来的数据公式，其实我将所有个体的截距都固定为1的，所以对数据最正确的模型应该是随机斜率模型，如下：

re.lm <- lmer(y ~ x + (0+x|unit), data = test.df)
summary(re.lm)

这次，大家再看，我们x的固定效应为2.98，随机效应为0.015，基本可原先xy的关系服从以3为均值，0.2为标准差的正态分布吻合了。即通过随机效应模型我们正确的得到了xy的真正关系。

我们可以查看我们拟合的每个人水平上的随机效应：

coef(re.lm)

30个人，每个人都有一个相同的截距和一个基本上以3为均值，0.2为标准差的斜率。

小结

今天再写一遍混合效应模型，大家应该会感觉更加清晰了，感谢大家耐心看完，自己的文章都写的很细，代码都在原文中，希望大家都可以自己做一做，如果对您有用请先收藏，再点赞转发，也欢迎大家的意见和建议。

如果你是一个大学本科生或研究生，如果你正在因为你的统计作业、数据分析、论文、报告、考试等发愁，如果你在使用SPSS,R，Python，Mplus, Excel中遇到任何问题，都可以联系我。因为我可以给您提供最好的，最详细和耐心的数据分析服务。

如果你对Z检验，t检验，方差分析，多元方差分析，回归，卡方检验，相关，多水平模型，结构方程模型，中介调节等等统计技巧有任何问题，请私信我，获取最详细和耐心的指导。

If you are a student and you are worried about you statistical #Assignments, #Data #Analysis, #Thesis, #reports, #composing, #Quizzes, Exams.. And if you are facing problem in #SPSS, #R-Programming, #Excel, Mplus, then contact me. Because I could provide you the best services for your Data Analysis.

Are you confused with statistical Techniques like z-test, t-test, ANOVA, MANOVA, Regression, Logistic Regression, Chi-Square, Correlation, Association, SEM, multilevel model, mediation and moderation etc. for your Data Analysis...??

Then Contact Me. I will solve your Problem...

加油吧，打工人！