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

原文出處：拓端數(shù)據(jù)部落公眾號

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

混合模型在很多方面與線性模型相似。它估計一個或多個解釋變量對因變量的影響。混合模型的輸出將為解釋值列表，它們的效果大小的估計值和置信區(qū)間，每種效果的p值以及至少一種模型擬合程度的度量。當(dāng)您有一個變量將數(shù)據(jù)樣本描述為可以收集的數(shù)據(jù)的子集時，應(yīng)該使用混合模型而不是簡單的線性模型。

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

1. str(data)

2. ## 'data.frame': 84 obs. of ?6 variables:

3. ## ?$ Test.ID ? : int ?1 2 3 4 5 6 7 8 9 10 ...

4. ## ?$ Observer ?: Factor w/ 4 levels "Charles","Michelle",..: 1 4 2 4 1 3 2 2 1 2 ...

5. ## ?$ Relation ?: Factor w/ 2 levels "Same","Stranger": 1 1 1 1 1 1 1 1 1 1 ...

6. ## ?$ Aggression: int ?4 1 15 2 1 0 2 0 3 10 ...

7. ## ?$ Tolerance : int ?4 34 14 31 4 13 7 6 13 15 ...

8. ## ?$ Season ? ?: Factor w/ 2 levels "Early","Late": 1 1 1 1 1 1 1 1 1 1 ...

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

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

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

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

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

1. require(car)

2. ## 正在加載所需的包: car

3. require(MASS)

4. # 必須為非零的分布

5. qqp(Aggression, "norm")

1. # lnorm 表示對數(shù)正態(tài)

2. 


3. 


4. qqp(Aggression, "lnorm")

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

2. fitdistr(rAggression, "Negative Binomial")

qqp(Aggressio, "pois", estimate)

fitdistr(Aggression.t, "gamma")

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

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

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

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

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

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

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

1. str(starlings)

2. ## 'data.frame': 28 obs. of ?5 variables:

3. ## ?$ Individual : Factor w/ 28 levels "B-40917","B-41205",..: 4 5 6 15 3 16 8 13 20 14 ...

4. ## ?$ Sex ? ? ? ?: Factor w/ 2 levels "F","M": 2 2 2 2 2 1 1 1 1 2 ...

5. ## ?$ Group ? ? ?: Factor w/ 5 levels "DRT1","MRC1",..: 2 5 5 4 4 4 4 4 4 4 ...

6. ## ?$ Social.Rank: Factor w/ 2 levels "Breeder","Helper": 2 1 1 1 2 2 2 2 1 2 ...

7. ## ?$ Mean.Pitch : num ?2911 2978 3313 3268 3312 ...

8. 


9. summary(lmm)

10. ## Linear mixed model fit by maximum likelihood ?['lmerMod']

11. ## Formula: Mean.Pitch ~ Sex + Social.Rank + (1 | Group)

12. ## ? ?Data: starlings

13. ##

14. ## ? ? ?AIC ? ? ?BIC ? logLik deviance df.resid

15. ## ? ?389.3 ? ?396.0 ? -189.7 ? ?379.3 ? ? ? 23

16. ##

17. ## Scaled residuals:

18. ## ? ? Min ? ? ?1Q ?Median ? ? ?3Q ? ? Max

19. ## -2.0559 -0.6272 ?0.0402 ?0.5801 ?2.0110

20. ##

21. ## Random effects:

22. ## ?Groups ? Name ? ? ? ?Variance Std.Dev.

23. ## ?Group ? ?(Intercept) ? ? 0 ? ? ?0

24. ## ?Residual ? ? ? ? ? ? 44804 ? ?212

25. ## Number of obs: 28, groups: Group, 5

26. ##

27. ## Fixed effects:

28. ## ? ? ? ? ? ? ? ? ? Estimate Std. Error t value

29. ## (Intercept) ? ? ? ? 3099.0 ? ? ? 82.2 ? ?37.7

30. ## SexM ? ? ? ? ? ? ? ? ?51.7 ? ? ? 81.3 ? ? 0.6

31. ## Social.RankHelper ? ?-45.0 ? ? ? 82.4 ? ?-0.5

32. ##

33. ## Correlation of Fixed Effects:

34. ## ? ? ? ? ? ? (Intr) SexM

35. ## SexM ? ? ? ?-0.630

36. ## Scl.RnkHlpr -0.668 ?0.106

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

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

2. ## Analysis of Deviance Table (Type II Wald chisquare tests)

3. ##

4. ## Response: Mean.Pitch

5. ## ? ? ? ? ? ? Chisq Df Pr(>Chisq)

6. ## Sex ? ? ? ? ? 0.4 ?1 ? ? ? 0.52

7. ## Social.Rank ? 0.3 ?1 ? ? ? 0.58

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

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

1. anova(noranklmm, nosexlmm, nofixedlmm)

2. ## Data: starlings

3. ## Models:

4. ## nofixedlmm: Mean.Pitch ~ 1 + (1 | Group)

5. ## noranklmm: Mean.Pitch ~ Sex + (1 | Group)

6. ## nosexlmm: Mean.Pitch ~ Social.Rank + (1 | Group)

7. ## ? ? ? ? ? ?Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)

8. ## nofixedlmm ?3 386 390 ? -190 ? ? ?380

9. ## noranklmm ? 4 388 393 ? -190 ? ? ?380 ?0.48 ? ? ?1 ? ? ? 0.49

10. ## nosexlmm ? ?4 388 393 ? -190 ? ? ?380 ?0.00 ? ? ?0 ? ? ? 1.00

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

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

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

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

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

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

2. ## ? ? lmList

3. summary(PQL)

4. ## Linear mixed-effects model fit by maximum likelihood

5. ## ?Data: recog

6. ## ? AIC BIC logLik

7. ## ? ?NA ?NA ? ? NA

8. ##

9. ## Random effects:

10. ## ?Formula: ~1 | Observer

11. ## ? ? ? ? (Intercept)

12. ## StdDev: ? ? ?0.3312

13. ##

14. ## ?Formula: ~1 | Test.ID %in% Observer

15. ## ? ? ? ? (Intercept) Residual

16. ## StdDev: ? ? ?0.5295 ? ?7.128

17. ##

18. ## Variance function:

19. ## ?Structure: fixed weights

20. ## ?Formula: ~invwt

21. ## Fixed effects: Aggression.t ~ Relation + Season

22. ## ? ? ? ? ? ? ? ? ? Value Std.Error DF t-value p-value

23. ## (Intercept) ? ? ? 1.033 ? ?0.5233 55 ? 1.974 ?0.0535

24. ## RelationStranger ?1.210 ? ?0.4674 55 ? 2.589 ?0.0123

25. ## SeasonLate ? ? ? -1.333 ? ?0.5983 23 ?-2.228 ?0.0359

26. ## ?Correlation:

27. ## ? ? ? ? ? ? ? ? ?(Intr) RltnSt

28. ## RelationStranger -0.855

29. ## SeasonLate ? ? ? -0.123 ?0.000

30. ##

31. ## Standardized Within-Group Residuals:

32. ## ? ? ?Min ? ? ? Q1 ? ? ?Med ? ? ? Q3 ? ? ?Max

33. ## -4.86916 -0.29958 -0.08012 ?0.14280 ?5.93336

34. ##

35. ## Number of Observations: 84

36. ## Number of Groups:

37. ## ? ? ? ? ? ? ?Observer Test.ID %in% Observer

38. ## ? ? ? ? ? ? ? ? ? ? 4 ? ? ? ? ? ? ? ? ? ?28

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

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

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

1. 


2. str(bdf)

3. ## 'data.frame': 2287 obs. of ?4 variables:

4. ## ?$ schoolNR: Factor w/ 131 levels "1","2","10","12",..: 1 1 1 1 1 1 1 1 1 1 ...

5. ## ?$ Minority: Factor w/ 2 levels "N","Y": 1 2 1 1 1 2 2 1 2 2 ...

6. ## ?$ ses ? ? : num ?23 10 15 23 10 10 23 10 13 15 ...

7. ## ?$ repeatgr: Factor w/ 3 levels "0","1","2": 1 1 1 1 1 1 1 1 2 1 ...

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

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

1. 


2. summary(GHQ)

3. ## Generalized linear mixed model fit by maximum likelihood (Adaptive

4. ## ? Gauss-Hermite Quadrature, nAGQ = 25) [glmerMod]

5. ## ?Family: binomial ( logit )

6. ## Formula: repeatgr ~ Minority + ses + ses * Minority + (1 | schoolNR)

7. ## ? ?Data: bdf

8. ##

9. ## ? ? ?AIC ? ? ?BIC ? logLik deviance df.resid

10. ## ? 1672.8 ? 1701.4 ? -831.4 ? 1662.8 ? ? 2282

11. ##

12. ## Scaled residuals:

13. ## ? ?Min ? ? 1Q Median ? ? 3Q ? ?Max

14. ## -0.964 -0.407 -0.314 -0.221 ?5.962

15. ##

16. ## Random effects:

17. ## ?Groups ? Name ? ? ? ?Variance Std.Dev.

18. ## ?schoolNR (Intercept) 0.264 ? ?0.514

19. ## Number of obs: 2287, groups: schoolNR, 131

20. ##

21. ## Fixed effects:

22. ## ? ? ? ? ? ? ? Estimate Std. Error z value Pr(>|z|)

23. ## (Intercept) ? -0.45194 ? ?0.20763 ? -2.18 ? ? 0.03 *

24. ## MinorityY ? ? ?0.47957 ? ?0.47345 ? ?1.01 ? ? 0.31

25. ## ses ? ? ? ? ? -0.06205 ? ?0.00798 ? -7.78 ?7.5e-15 ***

26. ## MinorityY:ses ?0.01196 ? ?0.02289 ? ?0.52 ? ? 0.60

27. ## ---

28. ## Signif. codes: ?0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

29. ##

30. ## Correlation of Fixed Effects:

31. ## ? ? ? ? ? ? (Intr) MnrtyY ses

32. ## MinorityY ? -0.400

33. ## ses ? ? ? ? -0.905 ?0.369

34. ## MinortyY:ss ?0.299 -0.866 -0.321

35. Laplace <- glmer(repeatgr ~ Minority + ses + ses * Minority + (1 | schoolNR),

36. data = bdf, family = binomial(link = "logit")) ?# Contrast to the Laplace approximation

37. ## Warning: Model failed to converge with max|grad| = 0.0458484 (tol = 0.001)

38. summary(Laplace)

39. ## Generalized linear mixed model fit by maximum likelihood (Laplace

40. ## ? Approximation) [glmerMod]

41. ## ?Family: binomial ( logit )

42. ## Formula: repeatgr ~ Minority + ses + ses * Minority + (1 | schoolNR)

43. ## ? ?Data: bdf

44. ##

45. ## ? ? ?AIC ? ? ?BIC ? logLik deviance df.resid

46. ## ? 1672.8 ? 1701.5 ? -831.4 ? 1662.8 ? ? 2282

47. ##

48. ## Scaled residuals:

49. ## ? ?Min ? ? 1Q Median ? ? 3Q ? ?Max

50. ## -0.960 -0.407 -0.316 -0.222 ?5.950

51. ##

52. ## Random effects:

53. ## ?Groups ? Name ? ? ? ?Variance Std.Dev.

54. ## ?schoolNR (Intercept) 0.258 ? ?0.508

55. ## Number of obs: 2287, groups: schoolNR, 131

56. ##

57. ## Fixed effects:

58. ## ? ? ? ? ? ? ? Estimate Std. Error z value Pr(>|z|)

59. ## (Intercept) ? -0.45456 ? ?0.20611 ? -2.21 ? ?0.027 *

60. ## MinorityY ? ? ?0.48005 ? ?0.47121 ? ?1.02 ? ?0.308

61. ## ses ? ? ? ? ? -0.06191 ? ?0.00791 ? -7.83 ?4.9e-15 ***

62. ## MinorityY:ses ?0.01194 ? ?0.02274 ? ?0.53 ? ?0.600

63. ## ---

64. ## Signif. codes: ?0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

65. ##

66. ## Correlation of Fixed Effects:

67. ## ? ? ? ? ? ? (Intr) MnrtyY ses

68. ## MinorityY ? -0.400

69. ## ses ? ? ? ? -0.906 ?0.369

70. ## MinortyY:ss ?0.299 -0.866 -0.321

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

1. 


2. ## n×n方差-協(xié)方差矩陣中方差參數(shù)的數(shù)量

3. 


4. 


5. vpars <- function(m) {

6. nrow(m) * (nrow(m) + 1)/2

7. }

8. # 接下來計算剩余自由度

9. 


10. rdf <- nrow(model.frame(model)) - model.df

11. # 提取皮爾遜殘差

12. 


13. 


14. prat <- Pearson.chisq/rdf

15. # 生成一個p值。如果小于0.05，則數(shù)據(jù)過于分散。

16. 


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

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

1. 


2. farmers <- farmers[c(1:5, 8)]

3. str(farmers)

4. ## 'data.frame': 102 obs. of ?6 variables:

5. ## ?$ year ? ?: int ?1901 1901 1901 1901 1902 1902 1902 1902 1902 1902 ...

6. ## ?$ farmer ?: Factor w/ 18 levels "Allardyce","Dooley",..: 10 5 3 18 10 5 18 17 4 12 ...

7. ## ?$ place ? : Factor w/ 16 levels "Arnestown","Bagnalstown",..: 3 16 14 11 3 16 11 4 8 6 ...

8. ## ?$ district: Factor w/ 4 levels "CentralPlain",..: 2 2 1 1 2 2 1 1 4 4 ...

9. ## ?$ gen ? ? : Factor w/ 2 levels "Archer","Goldthorpe": 1 1 1 1 1 1 1 1 1 1 ...

10. ## ?$ profit ?: num ?1 1 1 1 1 1 1 1 1 1 ...

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

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

2. ## ?Iterations = 3001:12991

3. ## ?Thinning interval ?= 10

4. ## ?Sample size ?= 1000

5. ##

6. ## ?DIC: 76.61

7. ##

8. ## ?G-structure: ?~year

9. ##

10. ## ? ? ?post.mean l-95% CI u-95% CI eff.samp

11. ## year ? ? ?18.2 ? ?0.266 ? ? 68.6 ? ? 18.1

12. ##

13. ## ? ? ? ? ? ? ? ?~farmer

14. ##

15. ## ? ? ? ?post.mean l-95% CI u-95% CI eff.samp

16. ## farmer ? ? ?1.93 ? 0.0789 ? ? 6.68 ? ? 99.7

17. ##

18. ## ? ? ? ? ? ? ? ?~place

19. ##

20. ## ? ? ? post.mean l-95% CI u-95% CI eff.samp

21. ## place ? ? ?1.54 ? 0.0711 ? ? 4.84 ? ? ?173

22. ##

23. ## ? ? ? ? ? ? ? ?~gen

24. ##

25. ## ? ? post.mean l-95% CI u-95% CI eff.samp

26. ## gen ? ? ?12.1 ? 0.0874 ? ? 28.6 ? ? 1000

27. ##

28. ## ? ? ? ? ? ? ? ?~district

29. ##

30. ## ? ? ? ? ?post.mean l-95% CI u-95% CI eff.samp

31. ## district ? ? ?1.78 ? 0.0639 ? ? 5.63 ? ? 1000

32. ##

33. ## ?R-structure: ?~units

34. ##

35. ## ? ? ? post.mean l-95% CI u-95% CI eff.samp

36. ## units ? ? ? ? 1 ? ? ? ?1 ? ? ? ?1 ? ? ? ?0

37. ##

38. ## ?Location effects: profit ~ 1

39. ##

40. ## ? ? ? ? ? ? post.mean l-95% CI u-95% CI eff.samp pMCMC

41. ## (Intercept) ? ? ?3.47 ? ?-1.16 ? ?10.75 ? ? ?220 ?0.14

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

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

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

1. 


2. ggplot(recog, aes(x = Aggression)+ geom_density() +

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

繪圖對于評估模型擬合也很重要。您可以通過各種方式繪制擬合值來判斷適合的模型最能描述數(shù)據(jù)。一個簡單的應(yīng)用是繪制模型的殘差。如果您將模型想象為通過數(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ī)效果看起來不像白噪聲。因此，讓我們嘗試使用更多迭代來重新擬合模型。這需要更多的計算量，但會產(chǎn)生更準(zhǔn)確的結(jié)果。

Trace(MCMC2, log = TRUE)

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

1. # 將兩個模型的估計值和置信區(qū)間放在一起

2. 


3. 


4. rbind (covariances, Gcovariances)

5. # 創(chuàng)建一個數(shù)據(jù)框架，其中包含模型和隨機(jī)效應(yīng)的因素

6. 


7. data.frame(coint, model = factor(el1", "model2"), eac),

8. random.effect = dimnames(i)

9. ## Warning: some row.names duplicated: 6,7,8,9,10 --> row.names NOT used

10. # 將估計和置信區(qū)間繪制為按模型分組的箱形圖

11. 


12. 


13. ggplot(conf.int+ geom_crossbar(aes(y.95..CI,

14. y.95..CI= model= "dodge")

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

最受歡迎的見解

1.基于R語言的lmer混合線性回歸模型

2.R語言用Rshiny探索lme4廣義線性混合模型（GLMM）和線性混合模型（LMM）

3.R語言線性混合效應(yīng)模型實(shí)戰(zhàn)案例

4.R語言線性混合效應(yīng)模型實(shí)戰(zhàn)案例2

5.R語言線性混合效應(yīng)模型實(shí)戰(zhàn)案例

6.線性混合效應(yīng)模型Linear Mixed-Effects Models的部分折疊Gibbs采樣

7.R語言LME4混合效應(yīng)模型研究教師的受歡迎程度

8.R語言中基于混合數(shù)據(jù)抽樣(MIDAS)回歸的HAR-RV模型預(yù)測GDP增長

9.使用SAS，Stata，HLM，R，SPSS和Mplus的分層線性模型HLM