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...

加油吧,打工人!

(0)

相关推荐