全文鏈接：?http://tecdat.cn/?p=2596

最近我們被客戶要求撰寫關(guān)于混合效應(yīng)模型的研究報(bào)告，包括一些圖形和統(tǒng)計(jì)輸出。

在本文中，我們描述了靈活的競(jìng)爭(zhēng)風(fēng)險(xiǎn)回歸模型?；貧w模型被指定為轉(zhuǎn)移概率，也就是競(jìng)爭(zhēng)性風(fēng)險(xiǎn)設(shè)置中的累積發(fā)生率

1.混合模型是否適合您的需求？

混合模型在很多方面與線性模型相似。它估計(jì)一個(gè)或多個(gè)解釋變量對(duì)因變量的影響?；旌夏Ｐ偷妮敵鰧榻忉屩盗斜恚鼈兊男Ч笮〉墓烙?jì)值和置信區(qū)間，每種效果的p值以及至少一種模型擬合程度的度量。當(dāng)您有一個(gè)變量將數(shù)據(jù)樣本描述為可以收集的數(shù)據(jù)的子集時(shí)，應(yīng)該使用混合模型而不是簡(jiǎn)單的線性模型。

讓我們看一下正在研究的黃蜂親屬識(shí)別數(shù)據(jù)。

str(data)##?'data.frame': ? ?84?obs.?of??6?variables:##??$?Test.ID???:?int??1?2?3?4?5?6?7?8?9?10?...##??$?Observer??:?Factor?w/?4?levels?"Charles","Michelle",..:?1?4?2?4?1?3?2?2?1?2?...##??$?Relation??:?Factor?w/?2?levels?"Same","Stranger":?1?1?1?1?1?1?1?1?1?1?...##??$?Aggression:?int??4?1?15?2?1?0?2?0?3?10?...##??$?Tolerance?:?int??4?34?14?31?4?13?7?6?13?15?...##??$?Season????:?Factor?w/?2?levels?"Early","Late":?1?1?1?1?1?1?1?1?1?1?...

我感興趣的因變量是攻擊性和寬容度。侵略性是指六十分鐘內(nèi)的攻擊行為次數(shù)。寬容是指六十分鐘內(nèi)的寬容行為數(shù)量。我對(duì)關(guān)系（無(wú)論黃蜂來(lái)自相同還是不同的菌落）和季節(jié)（菌落周期的早期或晚期）對(duì)這些因變量的影響感興趣。這些影響是“固定的”，因?yàn)闊o(wú)論我在何處，如何采樣或采樣了多少只黃蜂，我在相同變量中仍將具有相同的水平：相同的菌落與不同的菌落，以及早季與晚季。

但是，還有兩個(gè)其他變量在樣本之間不會(huì)保持固定。如果我在不同的年份進(jìn)行采樣，那么觀察者的水平會(huì)有所不同。樣品之間的測(cè)試ID也會(huì)有所不同，因?yàn)槲铱偸强梢灾匦掳才拍男S蜂參與每個(gè)實(shí)驗(yàn)試驗(yàn)。每個(gè)試驗(yàn)都是我當(dāng)時(shí)收集的黃蜂的唯一子樣本。如果我能夠單獨(dú)測(cè)試黃蜂，并且如果所有觀察都對(duì)所有互動(dòng)進(jìn)行了評(píng)分，那么我將不會(huì)有任何隨機(jī)效應(yīng)。但是，相反，我的數(shù)據(jù)本來(lái)就是“塊狀”的，隨機(jī)效應(yīng)描述了這種塊狀性。

在繼續(xù)之前，您還需要考慮隨機(jī)效果的結(jié)構(gòu)。您的隨機(jī)效果是嵌套還是交叉？在我的研究中，隨機(jī)效應(yīng)是?嵌套的，因?yàn)槊總€(gè)觀察者記錄了一定數(shù)量的試，并且沒(méi)有兩個(gè)觀察者記錄了相同的試驗(yàn)，因此Test.ID嵌套在Observer中。但是說(shuō)我收集了五個(gè)不同遺傳譜系中的黃蜂。“遺傳學(xué)”的隨機(jī)效應(yīng)與觀察無(wú)關(guān)。它將與其他兩個(gè)隨機(jī)效應(yīng)_交叉_。因此，這種隨機(jī)效應(yīng)將?與其他效應(yīng)?交叉。

視頻線性混合效應(yīng)模型LMM,Linear Mixed和R語(yǔ)言實(shí)現(xiàn)

**，時(shí)長(zhǎng)12:13

2.哪種概率分布最適合您的數(shù)據(jù)？

假設(shè)您已決定要運(yùn)行混合模型。接下來(lái)要做的是找到最適合您數(shù)據(jù)的概率分布。有很多測(cè)試方法。請(qǐng)注意，負(fù)二項(xiàng)式和伽馬分布只能處理正數(shù)，而泊松分布只能處理正整數(shù)。二項(xiàng)分布和泊松分布與其他分布不同，因?yàn)樗鼈兪请x散的而不是連續(xù)的，這意味著它們可以量化不同的，可數(shù)的事件或這些事件的概率?，F(xiàn)在讓我們?yōu)槲业腁ggression變量找到一個(gè)合適的分布。

require(car)##?正在加載所需的包:?carrequire(MASS)#?必須為非零的分布qqp(Aggression,?"norm")

#?lnorm?表示對(duì)數(shù)正態(tài)qqp(Aggression,?"lnorm")

# qqp需要估計(jì)負(fù)二項(xiàng)式，泊松和伽瑪分布的參數(shù)。您可以使用fitdistr函數(shù)生成估算值。保存輸出并提取每個(gè)參數(shù)的估計(jì)值，如下所示。fitdistr(rAggression,?"Negative?Binomial")

qqp(Aggressio,?"pois",?estimate)

fitdistr(Aggression.t,?"gamma")

查看我使用qqp生成的圖。y軸表示觀測(cè)值，x軸表示通過(guò)分布建模的分位數(shù)。紅色實(shí)線表示理想分布擬合，紅色虛線表示理想分布擬合的置信區(qū)間。您想選擇最大的觀測(cè)值落在虛線之間的分布。在這種情況下，這就是對(duì)數(shù)正態(tài)分布，其中只有一個(gè)觀測(cè)值落在虛線之外。現(xiàn)在，我可以嘗試擬合模型。

點(diǎn)擊標(biāo)題查閱往期內(nèi)容

R語(yǔ)言用Rshiny探索lme4廣義線性混合模型（GLMM）和線性混合模型（LMM）

左右滑動(dòng)查看更多

01

02

03

04

3.如何將混合模型擬合到您的數(shù)據(jù)

3a.如果您的數(shù)據(jù)是正態(tài)分布的

首先，請(qǐng)注意：如果您的數(shù)據(jù)最適合對(duì)數(shù)正態(tài)分布，?請(qǐng)不要對(duì)其進(jìn)行_變換_。?由于變換使模型結(jié)果的解釋更加困難。

如果數(shù)據(jù)呈正態(tài)分布，則可以使用線性混合模型（LMM）。該函數(shù)的第一個(gè)參數(shù)是一個(gè)公式，形式為y?x1 + x2 ...等，其中y是因變量，而x1，x2等是解釋變量。交叉隨機(jī)效應(yīng)的形式為（1 | r1）+（1 | r2）...，而嵌套隨機(jī)效應(yīng)的形式為（1 | r1 / r2）。

在這里，您可以指定混合模型將使用最大似然還是受限最大似然來(lái)估計(jì)參數(shù)。如果您的隨機(jī)效應(yīng)是嵌套的，或者只有一個(gè)隨機(jī)效應(yīng)，并且您的數(shù)據(jù)是平衡的（即，每個(gè)因子組中的樣本量相似），則將REML設(shè)置為FALSE，因?yàn)槟梢允褂米畲笏迫宦?。如果交叉了隨機(jī)效果，請(qǐng)不要設(shè)置REML參數(shù)，因?yàn)闊o(wú)論如何它默認(rèn)為TRUE。

為了避免這一切看起來(lái)太抽象，讓我們嘗試一些數(shù)據(jù)。我們將有關(guān)八哥歌曲研究的一些數(shù)據(jù)。在這項(xiàng)研究中，我們對(duì)雄性和雌性八哥歌曲之間的差異以及社會(huì)地位，不同的鳥類的歌唱是否不同感興趣。我們的隨機(jī)效應(yīng)是社會(huì)群體。歌曲的平均音高符合正態(tài)概率分布。

str(starlings)##?'data.frame': ? ?28?obs.?of??5?variables:##??$?Individual?:?Factor?w/?28?levels?"B-40917","B-41205",..:?4?5?6?15?3?16?8?13?20?14?...##??$?Sex????????:?Factor?w/?2?levels?"F","M":?2?2?2?2?2?1?1?1?1?2?...##??$?Group??????:?Factor?w/?5?levels?"DRT1","MRC1",..:?2?5?5?4?4?4?4?4?4?4?...##??$?Social.Rank:?Factor?w/?2?levels?"Breeder","Helper":?2?1?1?1?2?2?2?2?1?2?...##??$?Mean.Pitch?:?num??2911?2978?3313?3268?3312?...summary(lmm)##?Linear?mixed?model?fit?by?maximum?likelihood??['lmerMod']##?Formula:?Mean.Pitch?~?Sex?+?Social.Rank?+?(1?|?Group)##????Data:?starlings##?##??????AIC??????BIC???logLik?deviance?df.resid?##????389.3????396.0???-189.7????379.3???????23?##?##?Scaled?residuals:?##?????Min??????1Q??Median??????3Q?????Max?##?-2.0559?-0.6272??0.0402??0.5801??2.0110?##?##?Random?effects:##??Groups???Name????????Variance?Std.Dev.##??Group????(Intercept)?????0??????0?????##??Residual?????????????44804????212?????##?Number?of?obs:?28,?groups:?Group,?5##?##?Fixed?effects:##???????????????????Estimate?Std.?Error?t?value##?(Intercept)?????????3099.0???????82.2????37.7##?SexM??????????????????51.7???????81.3?????0.6##?Social.RankHelper????-45.0???????82.4????-0.5##?##?Correlation?of?Fixed?Effects:##?????????????(Intr)?SexM??##?SexM????????-0.630???????##?Scl.RnkHlpr?-0.668??0.106

讓我們看看結(jié)果。首先，我們獲得一些模型擬合的度量，包括AIC，BIC，對(duì)數(shù)似然度和偏差。然后我們得到由隨機(jī)效應(yīng)解釋的方差估計(jì)。這個(gè)數(shù)字很重要，因?yàn)槿绻c零沒(méi)有區(qū)別，那么您的隨機(jī)效應(yīng)可能并不重要，您可以繼續(xù)進(jìn)行常規(guī)的線性模型建立。接下來(lái)，我們將對(duì)固定效應(yīng)進(jìn)行估算，帶有標(biāo)準(zhǔn)誤差。這些信息可能足以滿足您的需求。一些期刊將這些模型的結(jié)果報(bào)告為帶有置信區(qū)間的效應(yīng)大小。當(dāng)然，當(dāng)我查看固定效應(yīng)估算值時(shí)，我已經(jīng)可以看出，性別和社會(huì)地位之間的平均音高沒(méi)有差異。但是有些期刊希望您報(bào)告p值。

如果您想要一些p值，則需要使用Anova函數(shù)。

##?Analysis?of?Deviance?Table?(Type?II?Wald?chisquare?tests)##?##?Response:?Mean.Pitch##?????????????Chisq?Df?Pr(>Chisq)##?Sex???????????0.4??1???????0.52##?Social.Rank???0.3??1???????0.58

Anova函數(shù)進(jìn)行了Wald檢驗(yàn)，該檢驗(yàn)告訴我們我們對(duì)性別和社會(huì)地位對(duì)音高的影響的估計(jì)p值。

擬合線性混合模型時(shí)，可能會(huì)遇到一種復(fù)雜情況。R可能會(huì)有“無(wú)法收斂”錯(cuò)誤，通常將其表述為“沒(méi)有收斂就達(dá)到了迭代限制”。這意味著您的模型有太多因素，樣本量不夠大，無(wú)法擬合。然后，您應(yīng)該做的是從模型中刪除固定效果和隨機(jī)效果，然后進(jìn)行比較以找出最合適的效果。一次刪除固定效果和隨機(jī)效果。保持固定效果不變，并一次刪除一個(gè)隨機(jī)效果，然后找出最合適的效果。然后保持隨機(jī)效果不變，并一次刪除固定效果。在這里，我只有一種隨機(jī)效果，

anova(noranklmm,?nosexlmm,?nofixedlmm)##?Data:?starlings##?Models:##?nofixedlmm:?Mean.Pitch?~?1?+?(1?|?Group)##?noranklmm:?Mean.Pitch?~?Sex?+?(1?|?Group)##?nosexlmm:?Mean.Pitch?~?Social.Rank?+?(1?|?Group)##????????????Df?AIC?BIC?logLik?deviance?Chisq?Chi?Df?Pr(>Chisq)##?nofixedlmm??3?386?390???-190??????380????????????????????????##?noranklmm???4?388?393???-190??????380??0.48??????1???????0.49##?nosexlmm????4?388?393???-190??????380??0.00??????0???????1.00

請(qǐng)注意，該方差分析函數(shù)與我們用來(lái)評(píng)估模型中固定效果的重要性的方差分析函數(shù)不同。方差分析函數(shù)用于比較模型。p值表明模型之間沒(méi)有明顯的重要差異。我們還可以比較AIC值，請(qǐng)注意，具有最低AIC值的模型是完全沒(méi)有固定影響的模型，這符合我們的理解，即性別和社會(huì)地位對(duì)歌曲的音調(diào)沒(méi)有影響。無(wú)論采用哪種方法，請(qǐng)務(wù)必在稿件中報(bào)告用于選擇最佳模型的p值或AIC值。

3b.如果您的數(shù)據(jù)不是正態(tài)分布的

您會(huì)看到，用于估計(jì)模型中影響大小的REML和最大似然法做出了不適用于數(shù)據(jù)的正態(tài)假設(shè)，因此您必須使用其他方法進(jìn)行參數(shù)估計(jì)。問(wèn)題在于，存在許多替代的估算方法，每種估算方法都使用不同的R包運(yùn)行，并且很難確定哪種方法合適。

首先，我們需要測(cè)試是否可以使用懲罰擬似然（PQL）。PQL是一種靈活的技術(shù)，可以處理非正常數(shù)據(jù)，不平衡設(shè)計(jì)和交叉隨機(jī)效應(yīng)。但是，如果您的因變量符合離散計(jì)數(shù)分布（例如泊松或二項(xiàng)式）且均值小于5，或者您的因變量為二元變量，則會(huì)產(chǎn)生偏差估計(jì)。

Aggression變量適合對(duì)數(shù)正態(tài)分布，該分布不是離散分布。這意味著我們可以繼續(xù)使用PQL方法。但是在繼續(xù)之前，讓我們回到轉(zhuǎn)變?yōu)檎龖B(tài)的問(wèn)題。

將分布設(shè)置為對(duì)數(shù)正態(tài)，我們將族設(shè)置為高斯，并將鏈接設(shè)置為log。

##?????lmListsummary(PQL)##?Linear?mixed-effects?model?fit?by?maximum?likelihood##??Data:?recog?##???AIC?BIC?logLik##????NA??NA?????NA##?##?Random?effects:##??Formula:?~1?|?Observer##?????????(Intercept)##?StdDev:??????0.3312##?##??Formula:?~1?|?Test.ID?%in%?Observer##?????????(Intercept)?Residual##?StdDev:??????0.5295????7.128##?##?Variance?function:##??Structure:?fixed?weights##??Formula:?~invwt?##?Fixed?effects:?Aggression.t?~?Relation?+?Season?##???????????????????Value?Std.Error?DF?t-value?p-value##?(Intercept)???????1.033????0.5233?55???1.974??0.0535##?RelationStranger??1.210????0.4674?55???2.589??0.0123##?SeasonLate???????-1.333????0.5983?23??-2.228??0.0359##??Correlation:?##??????????????????(Intr)?RltnSt##?RelationStranger?-0.855???????##?SeasonLate???????-0.123??0.000##?##?Standardized?Within-Group?Residuals:##??????Min???????Q1??????Med???????Q3??????Max?##?-4.86916?-0.29958?-0.08012??0.14280??5.93336?##?##?Number?of?Observations:?84##?Number?of?Groups:?##??????????????Observer?Test.ID?%in%?Observer?##?????????????????????4????????????????????28

該模型表明季節(jié)對(duì)侵略性有影響，也就是說(shuō)，在菌落周期后期收集的黃蜂比早期收集的黃蜂侵略性小。這也表明黃蜂之間的關(guān)系有影響。他們相對(duì)于巢友更可能對(duì)陌生人有侵略性。我將這些統(tǒng)計(jì)數(shù)據(jù)與估計(jì)值，標(biāo)準(zhǔn)誤，t值和p值一起報(bào)告。

那么，如果您的因變量的平均值小于5，或者您有一個(gè)二元因變量，而您不能使用PQL，該怎么辦？在這里，您可以使用兩種選擇：拉普拉斯（Laplace）近似?和馬爾可夫鏈蒙特卡羅算法（MCMC）。拉普拉斯近似值最多可以處理3個(gè)隨機(jī)效果。除此之外，您還必須使用MCMC。

讓我們從一個(gè)可以使用拉普拉斯逼近的例子開始。我們將使用學(xué)生在學(xué)校的學(xué)習(xí)情況的數(shù)據(jù)。出于本示例的目的，我將僅將數(shù)據(jù)子集化為幾個(gè)感興趣的變量，并將“ repeatgr”變量簡(jiǎn)化為二元因變量。

str(bdf) ##?'data.frame': ? ?2287?obs.?of??4?variables: ##??$?schoolNR:?Factor?w/?131?levels?"1","2","10","12",..:?1?1?1?1?1?1?1?1?1?1?... ##??$?Minority:?Factor?w/?2?levels?"N","Y":?1?2?1?1?1?2?2?1?2?2?... ##??$?ses?????:?num??23?10?15?23?10?10?23?10?13?15?... ##??$?repeatgr:?Factor?w/?3?levels?"0","1","2":?1?1?1?1?1?1?1?1?2?1?...

假設(shè)我們要找出是否屬于少數(shù)民族和社會(huì)經(jīng)濟(jì)地位會(huì)影響學(xué)生復(fù)讀成績(jī)的可能性。我們的因變量是“ repeatgr”，指示學(xué)生是否重復(fù)了成績(jī)。少數(shù)族身份是二元Y / N類別，社會(huì)經(jīng)濟(jì)地位由“ ses”表示，其數(shù)字范圍為10至50，其中50是最富有的。我們的隨機(jī)因素是“ schoolNR”，它代表從中采樣學(xué)生的學(xué)校。因?yàn)橐蜃兞渴嵌模晕覀冃枰哂卸?xiàng)式分布的廣義線性混合模型，并且由于我們的隨機(jī)效應(yīng)少于五個(gè)，因此可以使用Laplace近似?。

嚴(yán)格來(lái)說(shuō)，Laplace近似??是一種稱為Gauss-Hermite正交（GHQ）的參數(shù)估算方法的特例，需要進(jìn)行一次迭代。GHQ由于重復(fù)迭代而比Laplace更為準(zhǔn)確，但在第一次迭代后變得不那么靈活，因此您只能將它用于一個(gè)隨機(jī)效果。我們唯一的隨機(jī)效應(yīng)是“ schoolNR”。

summary(GHQ)##?Generalized?linear?mixed?model?fit?by?maximum?likelihood?(Adaptive##???Gauss-Hermite?Quadrature,?nAGQ?=?25)?[glmerMod]##??Family:?binomial?(?logit?)##?Formula:?repeatgr?~?Minority?+?ses?+?ses?*?Minority?+?(1?|?schoolNR)##????Data:?bdf##?##??????AIC??????BIC???logLik?deviance?df.resid?##???1672.8???1701.4???-831.4???1662.8?????2282?##?##?Scaled?residuals:?##????Min?????1Q?Median?????3Q????Max?##?-0.964?-0.407?-0.314?-0.221??5.962?##?##?Random?effects:##??Groups???Name????????Variance?Std.Dev.##??schoolNR?(Intercept)?0.264????0.514???##?Number?of?obs:?2287,?groups:?schoolNR,?131##?##?Fixed?effects:##???????????????Estimate?Std.?Error?z?value?Pr(>|z|)????##?(Intercept)???-0.45194????0.20763???-2.18?????0.03?*??##?MinorityY??????0.47957????0.47345????1.01?????0.31????##?ses???????????-0.06205????0.00798???-7.78??7.5e-15?***##?MinorityY:ses??0.01196????0.02289????0.52?????0.60????##?---##?Signif.?codes:??0?'***'?0.001?'**'?0.01?'*'?0.05?'.'?0.1?'?'?1##?##?Correlation?of?Fixed?Effects:##?????????????(Intr)?MnrtyY?ses???##?MinorityY???-0.400??????????????##?ses?????????-0.905??0.369???????##?MinortyY:ss??0.299?-0.866?-0.321Laplace?<-?glmer(repeatgr?~?Minority?+?ses?+?ses?*?Minority?+?(1?|?schoolNR), ????data?=?bdf,?family?=?binomial(link?=?"logit"))??#?Contrast?to?the?Laplace?approximation##?Warning:?Model?failed?to?converge?with?max|grad|?=?0.0458484?(tol?=?0.001)summary(Laplace)##?Generalized?linear?mixed?model?fit?by?maximum?likelihood?(Laplace##???Approximation)?[glmerMod]##??Family:?binomial?(?logit?)##?Formula:?repeatgr?~?Minority?+?ses?+?ses?*?Minority?+?(1?|?schoolNR)##????Data:?bdf##?##??????AIC??????BIC???logLik?deviance?df.resid?##???1672.8???1701.5???-831.4???1662.8?????2282?##?##?Scaled?residuals:?##????Min?????1Q?Median?????3Q????Max?##?-0.960?-0.407?-0.316?-0.222??5.950?##?##?Random?effects:##??Groups???Name????????Variance?Std.Dev.##??schoolNR?(Intercept)?0.258????0.508???##?Number?of?obs:?2287,?groups:?schoolNR,?131##?##?Fixed?effects:##???????????????Estimate?Std.?Error?z?value?Pr(>|z|)????##?(Intercept)???-0.45456????0.20611???-2.21????0.027?*??##?MinorityY??????0.48005????0.47121????1.02????0.308????##?ses???????????-0.06191????0.00791???-7.83??4.9e-15?***##?MinorityY:ses??0.01194????0.02274????0.53????0.600????##?---##?Signif.?codes:??0?'***'?0.001?'**'?0.01?'*'?0.05?'.'?0.1?'?'?1##?##?Correlation?of?Fixed?Effects:##?????????????(Intr)?MnrtyY?ses???##?MinorityY???-0.400??????????????##?ses?????????-0.906??0.369???????##?MinortyY:ss??0.299?-0.866?-0.321

您可以看到在這種情況下，Laplace和GHQ之間沒(méi)有重要區(qū)別。兩者都表明，社會(huì)經(jīng)濟(jì)課對(duì)學(xué)生復(fù)讀成績(jī)的可能性有非常顯著的影響，盡管即使采用logit變換，我們也可以看到影響量很小。但是，使用此方法時(shí)，還需要考慮其他一些因素。當(dāng)用于過(guò)度分散的數(shù)據(jù)時(shí)，即合并的殘差遠(yuǎn)大于殘差自由度時(shí)，它變得不準(zhǔn)確。使用此方法時(shí)，應(yīng)檢查模型以確保數(shù)據(jù)不會(huì)過(guò)度分散。

????##?n×n方差-協(xié)方差矩陣中方差參數(shù)的數(shù)量????vpars?<-?function(m)?{????????nrow(m)?*?(nrow(m)?+?1)/2????}????#?接下來(lái)計(jì)算剩余自由度????rdf?<-?nrow(model.frame(model))?-?model.df ????#?提取皮爾遜殘差????prat?<-?Pearson.chisq/rdf ????#?生成一個(gè)p值。如果小于0.05，則數(shù)據(jù)過(guò)于分散。

這些學(xué)校數(shù)據(jù)并不分散。如果是，可以將過(guò)度分散建模為隨機(jī)效應(yīng)，每個(gè)觀察值具有一個(gè)隨機(jī)效應(yīng)水平。在這種情況下，我將使用學(xué)生ID作為隨機(jī)效果變量。

我們繼續(xù)討論泊松數(shù)據(jù)的均值太小或因變量是分類的情況，并且我們有五個(gè)或五個(gè)以上隨機(jī)效應(yīng)。考慮有關(guān)大麥農(nóng)民的這些數(shù)據(jù)。想象一下，要使大麥?zhǔn)粘僧a(chǎn)生利潤(rùn)，大麥?zhǔn)粘傻氖杖氡仨毚笥?40。

farmers?<-?farmers[c(1:5,?8)] str(farmers)##?'data.frame': ? ?102?obs.?of??6?variables:##??$?year????:?int??1901?1901?1901?1901?1902?1902?1902?1902?1902?1902?...##??$?farmer??:?Factor?w/?18?levels?"Allardyce","Dooley",..:?10?5?3?18?10?5?18?17?4?12?...##??$?place???:?Factor?w/?16?levels?"Arnestown","Bagnalstown",..:?3?16?14?11?3?16?11?4?8?6?...##??$?district:?Factor?w/?4?levels?"CentralPlain",..:?2?2?1?1?2?2?1?1?4?4?...##??$?gen?????:?Factor?w/?2?levels?"Archer","Goldthorpe":?1?1?1?1?1?1?1?1?1?1?...##??$?profit??:?num??1?1?1?1?1?1?1?1?1?1?...

在這種情況下，假設(shè)我們根本沒(méi)有解釋變量。我們對(duì)哪些隨機(jī)效應(yīng)會(huì)影響利潤(rùn)感興趣。畢竟，隨機(jī)效應(yīng)是改變因變量方差的因素。為此，我們將處理隨機(jī)效應(yīng)，而且還為隨機(jī)效應(yīng)提供置信區(qū)間。

所有模型都對(duì)數(shù)據(jù)中方差的分布進(jìn)行假設(shè)，但是在貝葉斯方法中，這些假設(shè)是明確的，因此我們需要指定這些假設(shè)的分布。在貝葉斯統(tǒng)計(jì)中，我們稱這些?先驗(yàn)。

##??Iterations?=?3001:12991##??Thinning?interval??=?10##??Sample?size??=?1000?##?##??DIC:?76.61?##?##??G-structure:??~year##?##??????post.mean?l-95%?CI?u-95%?CI?eff.samp##?year??????18.2????0.266?????68.6?????18.1##?##????????????????~farmer##?##????????post.mean?l-95%?CI?u-95%?CI?eff.samp##?farmer??????1.93???0.0789?????6.68?????99.7##?##????????????????~place##?##???????post.mean?l-95%?CI?u-95%?CI?eff.samp##?place??????1.54???0.0711?????4.84??????173##?##????????????????~gen##?##?????post.mean?l-95%?CI?u-95%?CI?eff.samp##?gen??????12.1???0.0874?????28.6?????1000##?##????????????????~district##?##??????????post.mean?l-95%?CI?u-95%?CI?eff.samp##?district??????1.78???0.0639?????5.63?????1000##?##??R-structure:??~units##?##???????post.mean?l-95%?CI?u-95%?CI?eff.samp##?units?????????1????????1????????1????????0##?##??Location?effects:?profit?~?1?##?##?????????????post.mean?l-95%?CI?u-95%?CI?eff.samp?pMCMC##?(Intercept)??????3.47????-1.16????10.75??????220??0.14

如果我們看一下均值和置信區(qū)間，我們可以看到，真正影響利潤(rùn)可變性的唯一隨機(jī)效應(yīng)是年份和基因型。也就是說(shuō)，大麥?zhǔn)粘色@利的可能性每年之間以及基因型之間變化很大，但是在農(nóng)民或地方之間變化不大。

結(jié)束：了解您的數(shù)據(jù)

除非您熟悉這些分析，否則您將無(wú)法真正知道哪種分析對(duì)您的數(shù)據(jù)適用，而熟悉它們的最佳方法是繪制它們。通常，我的第一步是繪制我感興趣的變量的密度圖。

ggplot(recog,?aes(x?=?Aggression)+?geom_density()?+

在這里，我按季節(jié)和關(guān)系（兩個(gè)固定效應(yīng)）劃分了數(shù)據(jù)。我們可以立即看到數(shù)據(jù)集包含一個(gè)極端正的異常值；大多數(shù)觀測(cè)值都介于0到20之間。我們還可以看到，后期觀測(cè)值的很大一部分等于零。

繪圖對(duì)于評(píng)估模型擬合也很重要。您可以通過(guò)各種方式繪制擬合值來(lái)判斷適合的模型最能描述數(shù)據(jù)。一個(gè)簡(jiǎn)單的應(yīng)用是繪制模型的殘差。如果您將模型想象為通過(guò)數(shù)據(jù)散點(diǎn)圖的最佳擬合線，則殘差為散點(diǎn)圖中各點(diǎn)與最佳擬合線之間的距離。如果模型適合，則將殘差與擬合值作圖，則應(yīng)該看到隨機(jī)散布。如果散布不是隨機(jī)的，則意味著還有其他隨機(jī)或固定的影響。

讓我們嘗試?yán)L制擬合八哥的歌曲音高的混合模型的殘差。此圖所做的是創(chuàng)建一條表示零的水平虛線：與最佳擬合線的零偏差平均值。它還創(chuàng)建了一條實(shí)線，代表與最佳擬合線的殘差。我希望實(shí)線覆蓋虛線。

結(jié)果很好：與最佳擬合線的偏差趨于零。

讓我們嘗試一種以圖形方式比較MCMC模型的技術(shù)。我們將回到大麥農(nóng)戶的示例。

Trace(MCMC,?log?=?TRUE)

這些隨機(jī)效果看起來(lái)不像白噪聲。因此，讓我們嘗試使用更多迭代來(lái)重新擬合模型。這需要更多的計(jì)算量，但會(huì)產(chǎn)生更準(zhǔn)確的結(jié)果。

Trace(MCMC2,?log?=?TRUE)

現(xiàn)在，看起來(lái)都更接近直線周圍的白噪聲，這表明模型更好?，F(xiàn)在，讓我們比較兩個(gè)模型之間隨機(jī)效應(yīng)方差的置信區(qū)間。

#?將兩個(gè)模型的估計(jì)值和置信區(qū)間放在一起rbind?(covariances,?Gcovariances)#?創(chuàng)建一個(gè)數(shù)據(jù)框架，其中包含模型和隨機(jī)效應(yīng)的因素data.frame(coint,?model?=?factor(el1",?"model2"),?eac), ????random.effect?=?dimnames(i) ##?Warning:?some?row.names?duplicated:?6,7,8,9,10?-->?row.names?NOT?used #?將估計(jì)和置信區(qū)間繪制為按模型分組的箱形圖 ggplot(conf.int+?geom_crossbar(aes(y.95..CI, ????y.95..CI=?model=?"dodge")

結(jié)果很好，因?yàn)閮蓚€(gè)模型之間的估算值非常相似，但是在第二個(gè)模型中，對(duì)年的置信區(qū)間明顯較小，說(shuō)明這個(gè)估計(jì)更好。圖中可以證明第二種模型的推論，即基因型和年份是變異的主要因素。

本文摘選?《?基于R語(yǔ)言混合效應(yīng)模型（mixed model）案例研究?》?，點(diǎn)擊“閱讀原文”獲取全文完整資料。

點(diǎn)擊標(biāo)題查閱往期內(nèi)容

R語(yǔ)言建立和可視化混合效應(yīng)模型mixed effect model
【視頻】線性混合效應(yīng)模型(LMM,Linear Mixed Models)和R語(yǔ)言實(shí)現(xiàn)案例
R語(yǔ)言 線性混合效應(yīng)模型實(shí)戰(zhàn)案例
R語(yǔ)言用潛類別混合效應(yīng)模型(Latent Class Mixed Model ,LCMM)分析老年癡呆年齡數(shù)據(jù)
R語(yǔ)言貝葉斯廣義線性混合（多層次/水平/嵌套）模型GLMM、邏輯回歸分析教育留級(jí)影響因素?cái)?shù)據(jù)R語(yǔ)言估計(jì)多元標(biāo)記的潛過(guò)程混合效應(yīng)模型（lcmm）分析心理測(cè)試的認(rèn)知過(guò)程
R語(yǔ)言因子實(shí)驗(yàn)設(shè)計(jì)nlme擬合非線性混合模型分析有機(jī)農(nóng)業(yè)施氮水平
R語(yǔ)言非線性混合效應(yīng) NLME模型(固定效應(yīng)&隨機(jī)效應(yīng))對(duì)抗哮喘藥物茶堿動(dòng)力學(xué)研究
R語(yǔ)言用線性混合效應(yīng)（多水平/層次/嵌套）模型分析聲調(diào)高低與禮貌態(tài)度的關(guān)系
R語(yǔ)言LME4混合效應(yīng)模型研究教師的受歡迎程度R語(yǔ)言nlme、nlmer、lme4用（非）線性混合模型non-linear mixed model分析藻類數(shù)據(jù)實(shí)例
R語(yǔ)言混合線性模型、多層次模型、回歸模型分析學(xué)生平均成績(jī)GPA和可視化
R語(yǔ)言線性混合效應(yīng)模型（固定效應(yīng)&隨機(jī)效應(yīng)）和交互可視化3案例
R語(yǔ)言用lme4多層次（混合效應(yīng)）廣義線性模型（GLM），邏輯回歸分析教育留級(jí)調(diào)查數(shù)據(jù)R語(yǔ)言 線性混合效應(yīng)模型實(shí)戰(zhàn)案例
R語(yǔ)言混合效應(yīng)邏輯回歸（mixed effects logistic）模型分析肺癌數(shù)據(jù)
R語(yǔ)言如何用潛類別混合效應(yīng)模型（LCMM）分析抑郁癥狀
R語(yǔ)言基于copula的貝葉斯分層混合模型的診斷準(zhǔn)確性研究
R語(yǔ)言建立和可視化混合效應(yīng)模型mixed effect model
R語(yǔ)言LME4混合效應(yīng)模型研究教師的受歡迎程度
R語(yǔ)言 線性混合效應(yīng)模型實(shí)戰(zhàn)案例
R語(yǔ)言用Rshiny探索lme4廣義線性混合模型（GLMM）和線性混合模型（LMM）
R語(yǔ)言基于copula的貝葉斯分層混合模型的診斷準(zhǔn)確性研究
R語(yǔ)言如何解決線性混合模型中畸形擬合(Singular fit)的問(wèn)題
基于R語(yǔ)言的lmer混合線性回歸模型
R語(yǔ)言用WinBUGS 軟件對(duì)學(xué)術(shù)能力測(cè)驗(yàn)建立層次（分層）貝葉斯模型
R語(yǔ)言分層線性模型案例
R語(yǔ)言用WinBUGS 軟件對(duì)學(xué)術(shù)能力測(cè)驗(yàn)（SAT）建立分層模型
使用SAS，Stata，HLM，R，SPSS和Mplus的分層線性模型HLM
R語(yǔ)言用WinBUGS 軟件對(duì)學(xué)術(shù)能力測(cè)驗(yàn)建立層次（分層）貝葉斯模型
SPSS中的多層（等級(jí)）線性模型Multilevel linear models研究整容手術(shù)數(shù)據(jù)
用SPSS估計(jì)HLM多層（層次）線性模型模型