自己已經(jīng)寫了好幾篇關(guān)于隨機(jī)效應(yīng)模型的文章了，今天換個(gè)角度，從傳統(tǒng)回歸和隨機(jī)效應(yīng)模型的對(duì)比中模擬出兩模型真正的差異，讓你明白加上隨機(jī)效應(yīng)到底對(duì)模型會(huì)有什么樣的改變。

回顧

在傳統(tǒng)的回歸中，我們有如下式子：

yj = βxj+ εj

就是說(shuō)x與y之間的關(guān)系是β。

但現(xiàn)在，我們有一個(gè)縱向數(shù)據(jù)，比如說(shuō)我有30個(gè)人，每個(gè)人測(cè)xy均測(cè)量15次，那么有可能會(huì)存在每個(gè)人中xy的關(guān)系不一樣的情況。

這個(gè)很好理解，畢竟人與人之間是有差異的。那么如果真的是這樣，我們對(duì)于每個(gè)人都可以擬合出一個(gè)βi，如果每個(gè)人中xy關(guān)系的基線水準(zhǔn)也不一樣，那么可能我們每個(gè)人的xy關(guān)系上會(huì)有一個(gè)特定且不一樣的截距αi。

然后，我們的式子變成了：

yij = αi + βixij + εij

上面這個(gè)式子就是最簡(jiǎn)單的隨機(jī)效應(yīng)模型。其中αi為隨機(jī)截距，βi為隨機(jī)斜率。

image

R模擬

現(xiàn)在我模擬30個(gè)人，每個(gè)人測(cè)量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)

上面的代碼產(chǎn)生了30個(gè)人，每個(gè)人有xy的測(cè)量15次，xy的關(guān)系服從以3為均值，0.2為標(biāo)準(zhǔn)差的正態(tài)分布，每個(gè)人中xy的關(guān)系都不一樣：

image

在上面的數(shù)據(jù)中，我們有兩個(gè)水平，第一個(gè)水平是xy的測(cè)量，第二個(gè)水平是人，xy是嵌套在人的水平上的。

一般線性回歸

因?yàn)槲覀冎烂總€(gè)人xy關(guān)系是不一樣，所以我們分人做個(gè)回歸，這么一個(gè)做法是一般線性回歸對(duì)多水平嵌套數(shù)據(jù)能做到的極限了：

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)

上面的代碼就可以得到每個(gè)人中xy的關(guān)系βi：

我們看一看我們用一般線性回歸估計(jì)出來(lái)的βi和我們本來(lái)模擬的有什么差異：


par(mfrow = c(2, 1))
hist(beta, main = "XY真實(shí)的斜率", 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 = "一般線性回歸估計(jì)的斜率",

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

image

可以看出來(lái)估計(jì)的斜率分布變異比真實(shí)斜率更大一點(diǎn)。此時(shí)，我們并不能說(shuō)xy的關(guān)系到底如何，因?yàn)槲覀償M合了30個(gè)β，雖然這個(gè)β的分布和真實(shí)的分布差不太多（其實(shí)變異稍大），我們無(wú)法得出真實(shí)的xy之間的關(guān)系，你說(shuō)到底30個(gè)β到底選哪個(gè)？。

再看隨機(jī)效應(yīng)模型

在R中建立隨機(jī)效應(yīng)模型需要用到lme4這個(gè)包，隨機(jī)效應(yīng)部分一般表達(dá)為：

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

在 | 的左邊你可以設(shè)定隨機(jī)截距或者隨機(jī)斜率，在右邊需要設(shè)定隨機(jī)效應(yīng)的水平；如果x有隨機(jī)截距和隨機(jī)斜率，你就可以設(shè)定左邊為“1+x”，如果x只有隨機(jī)斜率沒(méi)有隨機(jī)截距，你設(shè)定左邊為“0+x”，因?yàn)殡S機(jī)效應(yīng)都是在高水平上的變異，所以在 |右邊你就應(yīng)該將這個(gè)水平指定出來(lái)，在本例中人為高水平，對(duì)應(yīng)的變量為數(shù)據(jù)庫(kù)中的unit，所以我們將右邊設(shè)定為unit。

那么對(duì)于本例的混合效應(yīng)模型我們可以寫出代碼：

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

image

上面的代碼擬合了一個(gè)有隨機(jī)截距和隨機(jī)斜率的混合模型，此時(shí)我們得到x的固定效應(yīng)為3.055，隨機(jī)效應(yīng)為0.116，和我們?cè)仍O(shè)定的xy的關(guān)系服從以3為均值，0.2為標(biāo)準(zhǔn)差的正態(tài)分布有點(diǎn)接近了。

但是我們?cè)炔](méi)有設(shè)定人水平上的截距的變異，大家回看原來(lái)的數(shù)據(jù)公式，其實(shí)我將所有個(gè)體的截距都固定為1的，所以對(duì)數(shù)據(jù)最正確的模型應(yīng)該是隨機(jī)斜率模型，如下：

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

[圖片上傳失敗...(image-6a8cff-1612072903647)]

這次，大家再看，我們x的固定效應(yīng)為2.98，隨機(jī)效應(yīng)為0.015，基本可原先xy的關(guān)系服從以3為均值，0.2為標(biāo)準(zhǔn)差的正態(tài)分布吻合了。即通過(guò)隨機(jī)效應(yīng)模型我們正確的得到了xy的真正關(guān)系。

我們可以查看我們擬合的每個(gè)人水平上的隨機(jī)效應(yīng)：

coef(re.lm)

image

30個(gè)人，每個(gè)人都有一個(gè)相同的截距和一個(gè)基本上以3為均值，0.2為標(biāo)準(zhǔn)差的斜率。

小結(jié)

今天再寫一遍混合效應(yīng)模型，大家應(yīng)該會(huì)感覺(jué)更加清晰了，感謝大家耐心看完，自己的文章都寫的很細(xì)，代碼都在原文中，希望大家都可以自己做一做，如果對(duì)您有用請(qǐng)先收藏，再點(diǎn)贊轉(zhuǎn)發(fā)，也歡迎大家的意見(jiàn)和建議。

如果你是一個(gè)大學(xué)本科生或研究生，如果你正在因?yàn)槟愕慕y(tǒng)計(jì)作業(yè)、數(shù)據(jù)分析、論文、報(bào)告、考試等發(fā)愁，如果你在使用SPSS,R，Python，Mplus, Excel中遇到任何問(wèn)題，都可以聯(lián)系我。因?yàn)槲铱梢越o您提供最好的，最詳細(xì)和耐心的數(shù)據(jù)分析服務(wù)。

如果你對(duì)Z檢驗(yàn)，t檢驗(yàn)，方差分析，多元方差分析，回歸，卡方檢驗(yàn)，相關(guān)，多水平模型，結(jié)構(gòu)方程模型，中介調(diào)節(jié)等等統(tǒng)計(jì)技巧有任何問(wèn)題，請(qǐng)私信我，獲取最詳細(xì)和耐心的指導(dǎo)。

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

加油吧，打工人！