1. 前言 這篇文檔,是為那些想了解混合線性模型的人準(zhǔn)備的。這里面很多部分,可以在很多領(lǐng)域中使用 。我們假定大家對(duì)一些矩陣和線性回歸的理論有所了解,但是更高級(jí)的知識(shí)只有模糊的認(rèn)識(shí),希望對(duì)你有所幫助。
混合線性模型,不同的學(xué)科有不同的名稱,使用不同的術(shù)語(yǔ)去描述同樣的東西。 這里試圖用一種比較簡(jiǎn)明的方法進(jìn)行闡述,我希望這不會(huì)讓你更迷惑。
2. 介紹混合模型是用于群集數(shù)據(jù)情況的極其有用的建模工具。擁有重復(fù)測(cè)量觀測(cè)數(shù)據(jù)的數(shù)據(jù)或以其他方式聚集觀測(cè)數(shù)據(jù)的數(shù)據(jù)(例如學(xué)校內(nèi)的學(xué)生,地理區(qū)域內(nèi)的城市)非常普遍?;旌夏P涂梢砸远喾N方式處理此類數(shù)據(jù),但是對(duì)于剛開始使用的術(shù)語(yǔ),尤其是跨學(xué)科的術(shù)語(yǔ),可能有點(diǎn)令人生畏。
關(guān)于混合模型,您可能會(huì)遇到一些術(shù)語(yǔ):
Random intercepts and slopes 隨機(jī)截距和斜率
Random effects 隨機(jī)效應(yīng)
Random coefficients 隨機(jī)系數(shù)
Varying coefficients 系數(shù)變化
Intercepts and slopes-as-outcomes 攔截和結(jié)果傾斜
Hierarchical linear models 分層線性模型
Growth curve models (possibly Latent GCM) 增長(zhǎng)曲線模型(可能是潛在的GCM)
Mixed effects models 混合效果模型
所有描述混合模型的名稱, 有些可能更具歷史性,有些則更多地出現(xiàn)在特定學(xué)科中,有些則可能引用某種數(shù)據(jù)結(jié)構(gòu)(例如多級(jí)群集),而另一些則是特殊情況。
混合效應(yīng)或簡(jiǎn)單混合模型通常是指固定效應(yīng)和隨機(jī)效應(yīng)的混合。我更喜歡混合模型一詞,因?yàn)樗芎?jiǎn)單并且沒有暗示特定的結(jié)構(gòu)。
3. 標(biāo)準(zhǔn)線性模型首先,讓我們從標(biāo)準(zhǔn)線性模型開始,以熟悉該表示法。為了使事情盡可能簡(jiǎn)單,同時(shí)又可以推廣到常見數(shù)據(jù)情況,我將假設(shè)一些感興趣的變量y和一個(gè)連續(xù)/數(shù)字協(xié)變量。
這里:
e 是殘差,它滿足平均數(shù)為0,方差為Ve的正態(tài)分布 如果寫為矩陣的形式:
μ = X %*% β y = rnorm(n, μ, σ2)
在嘗試達(dá)到平衡時(shí),我懷疑這種方法可能會(huì)在不同程度上成功或失敗,但是在符號(hào)和代碼之間(很多教科書演示中都缺少這種東西),我希望事情會(huì)很清楚。
4. 應(yīng)用實(shí)例讓我們看一些數(shù)據(jù),開始考慮混合模型。我將使用lme4軟件包中的sleepstudy數(shù)據(jù)。以下描述來自相應(yīng)的幫助文件。
? 睡眠剝奪研究對(duì)象每天的平均反應(yīng)時(shí)間。在第0天,受試者具有正常的睡眠量。從那天晚上開始,他們每晚只能睡3個(gè)小時(shí)。觀察結(jié)果代表每天對(duì)每個(gè)受試者進(jìn)行的一系列測(cè)試的平均反應(yīng)時(shí)間(以毫秒為單位)。
? 讓我們使用標(biāo)準(zhǔn)的線性模型來探討持續(xù)睡眠剝奪對(duì)反應(yīng)時(shí)間的影響。
這里用線性回歸模型,Days為x變量,Reaction為y變量(還有人和我一樣,對(duì)因變量和自變量摸不著頭腦的么,用x變量和y變量更容易理解有沒有?。?/p># > data(sleepstudy, package='lme4') # > slim = lm(Reaction ~ Days, data=sleepstudy) # > summary(slim) # # Call: # lm(formula = Reaction ~ Days, data = sleepstudy) # # Residuals: # Min 1Q Median 3Q Max # -110.848 -27.483 1.546 26.142 139.953 # # Coefficients: # Estimate Std. Error t value Pr(>|t|) # (Intercept) 251.405 6.610 38.033 < 2e-16 *** # Days 10.467 1.238 8.454 9.89e-15 *** # --- # Signif. codes: 0 '***’ 0.001 '**’ 0.01 '*’ 0.05 '.’ 0.1 ' ’ 1 # # Residual standard error: 47.71 on 178 degrees of freedom # Multiple R-squared: 0.2865,Adjusted R-squared: 0.2825 # F-statistic: 71.46 on 1 and 178 DF, p-value: 9.894e-15
在天數(shù)為正的情況下,我們看到更多的睡眠剝奪會(huì)導(dǎo)致反應(yīng)時(shí)間增加/變慢。但是讓我們繪制數(shù)據(jù)。 這告訴我們什么?黑線是我們當(dāng)前模型的建議,即假設(shè)每個(gè)人的出發(fā)點(diǎn)和軌跡都相同。但是,我們看到對(duì)象的起點(diǎn)可能相差100毫秒之多。此外,雖然斜率通常為正,但有些斜率隨時(shí)間變化很小甚至沒有變化。換句話說,個(gè)人的截距和坡度都有明顯的變化。我們將在最后討論這些數(shù)據(jù)。
5. 所有可能的混線性模型分析這個(gè)數(shù)據(jù)因此,我們要考慮數(shù)據(jù)的集群性質(zhì)。與其像上面的SLiM中那樣忽略聚類,不如考慮為每個(gè)人運(yùn)行完全獨(dú)立的回歸。但是,這些模型通常只需要很少的數(shù)據(jù)就可以運(yùn)行,并且會(huì)被過度上下文化。正如我們將看到的,混合模型將允許每個(gè)人的隨機(jī)截距和斜率,并在不因個(gè)人而異的情況下考慮聚類。
如何描述這個(gè)模型?事實(shí)證明,它可以并且以多種方式顯示,具體取決于您正在查看的文本或文章。以下內(nèi)容受Gelman&Hill(2007)的啟發(fā),他們展示了編寫混合模型的五種方法。為簡(jiǎn)單起見,我們通常只關(guān)注隨機(jī)截距模型,但有時(shí)會(huì)超出該范圍。前幾對(duì)公式僅需要了解標(biāo)準(zhǔn)回歸模型,但其他模型描述則需要更多知識(shí)。
5.1 Mixed Model 1a: Allowing coefficients to vary across groups 在上面,c簇內(nèi)的每個(gè)觀測(cè)i都有一個(gè)截距α,這取決于它所屬的c簇。αc假設(shè)為正態(tài)分布,平均μα和方差τ2。μα是我們?cè)赟LiM方法中看到的總截距。e是SLiM中描述的正態(tài)分布的平均零。對(duì)于每一個(gè)模型描述,我將注意到一個(gè)主要的參考,在那里人們可以看到它的形式與特定的文本或文章幾乎相同。它將不是唯一一個(gè)這樣做的引用,但至少它應(yīng)該是一個(gè)提供一些額外視角的引用。
5.2 Mixed Model 1b: Multilevel model 5.3 Mixed Model 2: Combining separate local regressions 在這里插入圖片描述 5.4 Mixed Model 3a: Design matrix for random component 5.5 Mixed Model 3b: Design matrix again 5.6 Mixed Model 3c: General notation
5.7 Mixed Model 4a: Regression with multiple error terms 5.8 Mixed Model 4b: Conditional vs. marginal model
5.9 Mixed Model 5b: Multivariate normal model 5.10 Mixed Model 6: Penalized regression
5.11 Mixed Model 7: Bayesian mixed model 6 模擬數(shù)據(jù)運(yùn)行混合模型這里,設(shè)置:Va = 0.50.5 = 0.25
Ve = 1 1 =1
隨機(jī)因子blup:gamma_
截距:3
固定因子blue:0.75
# setup set.seed(1234) nclus = 100 # number of groups clus = factor(rep(1:nclus, each=10)) # cluster variable n = length(clus) # total n # parameters sigma = 1 # residual sd tau = .5 # re sd gamma_ = rnorm(nclus, mean=0, sd=tau) # random effects e = rnorm(n, mean=0, sd=sigma) # residual error intercept = 3 # fixed effects b1 = .75 # data x = rnorm(n) # covariate y = intercept + b1*x + gamma_[clus] + e # see model 1 d = data.frame(x, y, clus=clus) head(d) str(d)
「使用lme4包運(yùn)行」
library(lme4) lmeMod = lmer(y ~ x + (1|clus), data=d) summary(lmeMod)
估算的結(jié)果可以看出:Va = 0.224,和0.25比較接近
Ve = 0.97,和1比較接近
blue:0.799,和0.75接近
截距:2.9,和3接近
提取blup值:
「asreml代碼」
> library(asreml) > mod1 = asreml(y ~ x, random = ~ clus, residual = ~ idv(units),data=d) Model fitted using the sigma parameterization. ASReml 4.1.0 Wed Apr 5 16:34:50 2020 LogLik Sigma2 DF wall cpu 1 -3817.282 1.0 998 16:34:50 0.0 2 -2862.495 1.0 998 16:34:50 0.0 3 -1811.528 1.0 998 16:34:50 0.0 4 -1082.178 1.0 998 16:34:50 0.0 5 -688.386 1.0 998 16:34:50 0.0 6 -572.668 1.0 998 16:34:50 0.0 7 -553.615 1.0 998 16:34:50 0.0 8 -552.690 1.0 998 16:34:50 0.0 9 -552.687 1.0 998 16:34:50 0.0 > summary(mod1)$varcomp component std.error z.ratio bound %ch clus 0.2247008 0.04605034 4.879461 P 0 units!units 0.9755909 0.04601438 21.201871 P 0 units!R 1.0000000 NA NA F 0 > coef(mod1)$fixed effect x 0.7994379 (Intercept) 2.9008683
結(jié)果是一致的
比較設(shè)定的blup值和計(jì)算的blup值的相關(guān)系數(shù):
> cor(gamma_,coef(mod1)$random) effect [1,] 0.838686
7 sleepstudy數(shù)據(jù)運(yùn)行混合模型sleepMod = lmer(Reaction ~ Days + (Days|Subject), data=sleepstudy) summary(sleepMod)
> # sleepstudy 數(shù)據(jù) > sleepMod = lmer(Reaction ~ Days + (Days|Subject), data=sleepstudy) > summary(sleepMod) Linear mixed model fit by REML ['lmerMod'] Formula: Reaction ~ Days + (Days | Subject) Data: sleepstudy REML criterion at convergence: 1743.6 Scaled residuals: Min 1Q Median 3Q Max -3.9536 -0.4634 0.0231 0.4633 5.1793 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 t value (Intercept) 251.405 6.824 36.843 Days 10.467 1.546 6.771 Correlation of Fixed Effects: (Intr) Days -0.138
查看每個(gè)subject的效應(yīng)值以及截距:
> ranef(sleepMod) $Subject (Intercept) Days 308 2.2575329 9.1992737 309 -40.3942719 -8.6205161 310 -38.9563542 -5.4495796 330 23.6888704 -4.8141448 331 22.2585409 -3.0696766 332 9.0387625 -0.2720535 333 16.8389833 -0.2233978 334 -7.2320462 1.0745075 335 -0.3326901 -10.7524799 337 34.8865253 8.6290208 349 -25.2080191 1.1730997 350 -13.0694180 6.6142185 351 4.5777099 -3.0152825 352 20.8614523 3.5364062 369 3.2750882 0.8722876 370 -25.6110745 4.8222518 371 0.8070591 -0.9881730 372 12.3133491 1.2842380 with conditional variances for “Subject”
如果我們對(duì)其進(jìn)行作圖:
我們可以看到混合模型的好處,因?yàn)槲覀儠?huì)有結(jié)合了個(gè)體特定影響的預(yù)測(cè),預(yù)測(cè)的更準(zhǔn)確。
8 其它主題我將簡(jiǎn)要提及其他一些主題,但這些主題不會(huì)改變到目前為止討論的一般方法。
增加分組的協(xié)變量(Cluster level covariates ) 注意隨機(jī)因子是鑲嵌結(jié)構(gòu),還是交互結(jié)構(gòu) 你可能注意lme4包中沒有給出p-value值,軟件不會(huì)直接給出(除非用的是貝葉斯框架),其它軟件包給出p-value,不一定說明他就是正確的。 隨機(jī)因子有關(guān)系矩陣?響應(yīng)變量是二分類?還有很多問題需要考慮,有些數(shù)據(jù)不適合用混合模型去分析 9. 匯總在模型描述和代碼之間,希望您對(duì)標(biāo)準(zhǔn)的混合模型框架有了很好的了解?;旌夏P褪菍?duì)標(biāo)準(zhǔn)glm的非常靈活的擴(kuò)展,可以直接與加性模型,空間模型和其他模型建立聯(lián)系,因此可以將它們帶到很遠(yuǎn)。我可以說在lme4,mgcv和brms之間,將有很多很多方法可以以多種方式瀏覽其數(shù)據(jù)。祝您研究順利!
10. 參考文獻(xiàn)? Bates, Douglas, Martin M?chler, Benjamin Bolker, and Steven Walker. 2015. “Fitting Linear Mixed-Effects Models Using Lme4.”
Fahrmeir, Ludwig, Thomas Kneib, Stefan Lang, and Brian Marx. 2013. “Regression.”
Gelman, Andrew, and Jennifer Hill. 2007. “Data Analysis Using Regression and Multilevel/Hierarchical Models.”
Gelman, Andrew, John Carlin, Hal Stern, David Dunson, Aki Vehtari, and Donald Rubin. 2013. “Bayesian Data Analysis.”
Wood, Simon. 2006. “Generalized Additive Models.”
? 英文原文:https://m-clark./docs/mixedModels/mixedModels.html#standard_linear_model
「?jìng)€(gè)人公眾號(hào):育種數(shù)據(jù)分析之放飛自我」