自己已經(jīng)寫了好幾篇關(guān)于隨機(jī)效應(yīng)模型的文章了曼玩，今天換個角度鳞骤，從傳統(tǒng)回歸和隨機(jī)效應(yīng)模型的對比中模擬出兩模型真正的差異，讓你明白加上隨機(jī)效應(yīng)到底對模型會有什么樣的改變黍判。
回顧
在傳統(tǒng)的回歸中豫尽，我們有如下式子：
yj = βxj+ εj
就是說x與y之間的關(guān)系是β。
但現(xiàn)在顷帖，我們有一個縱向數(shù)據(jù)美旧，比如說我有30個人，每個人測xy均測量15次贬墩，那么有可能會存在每個人中xy的關(guān)系不一樣的情況榴嗅。
這個很好理解，畢竟人與人之間是有差異的陶舞。那么如果真的是這樣嗽测，我們對于每個人都可以擬合出一個βi，如果每個人中xy關(guān)系的基線水準(zhǔn)也不一樣肿孵，那么可能我們每個人的xy關(guān)系上會有一個特定且不一樣的截距αi唠粥。
然后，我們的式子變成了：
yij = αi + βixij + εij
上面這個式子就是最簡單的隨機(jī)效應(yīng)模型停做。其中αi為隨機(jī)截距厅贪，βi為隨機(jī)斜率。
R模擬
現(xiàn)在我模擬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)
上面的代碼產(chǎn)生了30個人养涮，每個人有xy的測量15次，xy的關(guān)系服從以3為均值眉抬，0.2為標(biāo)準(zhǔn)差的正態(tài)分布贯吓，每個人中xy的關(guān)系都不一樣：
在上面的數(shù)據(jù)中，我們有兩個水平蜀变，第一個水平是xy的測量悄谐，第二個水平是人，xy是嵌套在人的水平上的库北。
一般線性回歸
因為我們知道每個人xy關(guān)系是不一樣爬舰，所以我們分人做個回歸，這么一個做法是一般線性回歸對多水平嵌套數(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)
上面的代碼就可以得到每個人中xy的關(guān)系β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的關(guān)系到底如何，因為我們擬合了30個β杂腰，雖然這個β的分布和真實的分布差不太多（其實變異稍大）垃你，我們無法得出真實的xy之間的關(guān)系，你說到底30個β到底選哪個？惜颇。
再看隨機(jī)效應(yīng)模型
在R中建立隨機(jī)效應(yīng)模型需要用到lme4這個包皆刺，隨機(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ī)斜率沒有隨機(jī)截距锨亏，你設(shè)定左邊為“0+x”林说，因為隨機(jī)效應(yīng)都是在高水平上的變異，所以在 |右邊你就應(yīng)該將這個水平指定出來屯伞，在本例中人為高水平腿箩，對應(yīng)的變量為數(shù)據(jù)庫中的unit，所以我們將右邊設(shè)定為unit劣摇。
那么對于本例的混合效應(yīng)模型我們可以寫出代碼：
library(lme4)
re.lm <- lmer(y ~ x + (1+x|unit), data = test.df)
summary(re.lm)
上面的代碼擬合了一個有隨機(jī)截距和隨機(jī)斜率的混合模型珠移，此時我們得到x的固定效應(yīng)為3.055，隨機(jī)效應(yīng)為0.116末融，和我們原先設(shè)定的xy的關(guān)系服從以3為均值钧惧，0.2為標(biāo)準(zhǔn)差的正態(tài)分布有點接近了。
但是我們原先并沒有設(shè)定人水平上的截距的變異勾习，大家回看原來的數(shù)據(jù)公式浓瞪，其實我將所有個體的截距都固定為1的，所以對數(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)分布吻合了湿右。即通過隨機(jī)效應(yīng)模型我們正確的得到了xy的真正關(guān)系诅妹。
我們可以查看我們擬合的每個人水平上的隨機(jī)效應(yīng)：
coef(re.lm)
30個人，每個人都有一個相同的截距和一個基本上以3為均值毅人，0.2為標(biāo)準(zhǔn)差的斜率吭狡。
小結(jié)
今天再寫一遍混合效應(yīng)模型，大家應(yīng)該會感覺更加清晰了丈莺，感謝大家耐心看完划煮，自己的文章都寫的很細(xì)，代碼都在原文中场刑，希望大家都可以自己做一做般此，如果對您有用請先收藏蚪战，再點贊轉(zhuǎn)發(fā)牵现，也歡迎大家的意見和建議铐懊。
如果你是一個大學(xué)本科生或研究生，如果你正在因為你的統(tǒng)計作業(yè)瞎疼、數(shù)據(jù)分析科乎、論文、報告贼急、考試等發(fā)愁茅茂，如果你在使用SPSS,R，Python太抓，Mplus, Excel中遇到任何問題空闲，都可以聯(lián)系我。因為我可以給您提供最好的走敌，最詳細(xì)和耐心的數(shù)據(jù)分析服務(wù)碴倾。
如果你對Z檢驗，t檢驗掉丽，方差分析跌榔，多元方差分析，回歸捶障，卡方檢驗僧须，相關(guān)，多水平模型项炼，結(jié)構(gòu)方程模型担平，中介調(diào)節(jié)等等統(tǒ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...
加油吧，打工人空免！
往期內(nèi)容：
從“我丑到我自己了”說起——混合效應(yīng)模型續(xù)
重復(fù)測量數(shù)據(jù)分析系列：混合效應(yīng)模型基礎(chǔ)
