拓端tecdat|R語言混合線性模型、多層次模型、回歸模型分析學生平均成績GPA和可視化
原文鏈接：http://tecdat.cn/?p=23159?
原文出處：拓端數據部落公眾號
混合模型在統(tǒng)計學領域已經存在了很長時間。例如，標準的方差分析方法可以被看作是混合模型的特殊情況。最近，混合模型有多種應用和擴展，使其能夠涵蓋各種不同的數據情況。
術語
對于不熟悉的人來說，圍繞混合模型的術語，特別是跨學科的術語，可能有點令人困惑。你可能遇到的關于這些類型的模型的一些術語包括。
方差分量
隨機截距和斜率
隨機效應
隨機系數
變化的系數
截距和/或斜率作為結果
分層線性模型
多層次模型（意味著多層次的分層聚類數據
增長曲線模型（可能是Latent GCM）。
混合效應模型
都描述了混合模型的類型。 混合效應，或簡稱混合，模型一般指固定效應和隨機效應的混合。對于一般的模型，我更喜歡用 "混合模型 "或 "隨機效應模型"，因為它們是更簡單的術語，沒有暗示具體的結構，而且后者還可以適用于許多人在使用其他術語時不會想到的擴展情況。關于混合效應，固定效應指的是人們在線性回歸模型中看到的典型主效應，即混合模型的非隨機部分，在某些情況下，它們被稱為總體平均效應。無論如何定義，隨機效應只是那些特定于觀察樣本的效應。本文所概述的方法主要涉及觀察樣本是某種分組因素水平的情況。
聚類的種類
數據可能有一個或多個聚類來源，而且聚類可能是分層的，比如聚類是嵌套在其他聚類中。一個例子是對學生進行多次的學術能力測試（重復觀察嵌套在學生中，學生嵌套在學校中，學校嵌套在地區(qū)中）。在其他情況下，不存在嵌套結構。一個例子是反應時間實驗，參與者執(zhí)行同一組任務。雖然觀察結果是在個人內部嵌套的，但觀察結果也是根據任務類型進行聚類的。有人用嵌套和交叉這兩個術語來區(qū)分這些情況。此外，聚類可能是平衡的或不平衡的。
在下面的內容中，我們將看到所有這些數據情況下的混合效應模型。一般來說，我們的方法將是相同的，因為這種聚類實際上更多的是數據的屬性而不是模型的屬性。然而，了解混合模型在處理各種數據情況時的靈活性是很重要的。
隨機截距模型
下面我們展示混合模型的最簡單和最常見的情況，即我們有一個單一的分組/群組結構的隨機效應。這通常被稱為隨機截距模型。
例子：學生大學平均成績GPA
下面我們將評估預測大學平均成績（GPA）的因素。200名學生中的每個人都被評估了六次（前三年的每個學期），因此我們在學生中進行了觀察分組。我們還有其他變量，如狀態(tài)、性別和高中GPA。有些會以標簽和數字的形式出現。
標準回歸模型
現在說說基礎模型。我們可以用幾種不同的方式來展示它。首先，我們從一個標準回歸開始，來確定我們的方向。?
我們對截距和時間的影響有系數（b）。誤差（?）被假定為正態(tài)分布，平均值為0，標準差為σ。?
另一種寫法是強調gpa的基本數據生成過程的模型，如下。?
更嚴格地說，GPA和μ變量有一個隱含的下標來表示每個觀察值，但你也可以把它看成是單個個體在單個時間點的模型。
混合模型
描述
現在我們展示一種描述混合模型的方法，其中包括每個學生的獨特效應。考慮一下下面這個單一學生的模型。這表明，學生的特定效應，即GPA的偏差只是因為該學生是誰，可以被看作是一個額外的方差來源。
?
我們（通常）會對學生效應做如下假設。
因此，學生效應是隨機的，具體來說是正態(tài)分布，均值為零，有一定的估計標準差（τ）。換句話說，從概念上講，這個混合模型和標準回歸之間的唯一區(qū)別是學生效應，平均而言，學生效應是沒有影響的，但通常會因學生的不同而有一定的變化，平均而言是τ。
如果我們重新排列，我們反而可以關注模型的系數，而不是作為額外的誤差來源。
或者更簡潔的說
這樣一來，我們就會有針對學生的截距，因為每個人都會有自己獨特的影響加到總體截距上，導致每個人的截距不同。
現在我們看到截距是正態(tài)分布的，有總體截距的平均值和一些標準差。因此，這通常被稱為隨機截距模型。?
多層次模型
另一種顯示混合模型的方式常見于多層次模型的文獻中。它被更明確地顯示為一個兩部分的回歸模型，一個在樣本層面，一個在學生層面。
然而，在將第二層次的部分 "插入 "第一層次后，它與前者是相同的。
請注意，我們并沒有一個針對學生情景的效應。在這種情況下，情景被說成是一個固定效應，而沒有隨機成分。但情況絕對不是這樣的，我們后面會看到。
應用
可視化
在這里，我們繪制GPA與情景（即學期）的關系，來了解起點和趨勢的變化。
plot(occasion, gpa,smooth(method = 'lm')
?
所有學生的路徑都以淡色路徑顯示，10個樣本以粗體顯示。我們稍后要做的回歸所估計的總體趨勢顯示為紅色。有兩件事很突出。一是學生在開始時有很大的變數。第二，雖然GPA的總體趨勢是隨著時間的推移而上升但個別學生在這個軌跡上可能會有所不同。
標準回歸
因此，讓我們開始吧。首先，我們來看看回歸，只把時間指標作為協變量，我們把它當作數字。
lm(gpa ~ occasion)## summary(lm)
?
上面的數據告訴我們，開始時，即學期為零時，平均GPA，用截距表示，是2.6。此外，隨著我們從一個學期到另一個學期，我們可以預期GPA會增加大約0.11分。這將是很好的，除了我們忽略了聚類。這樣做的一個副作用是，我們的標準誤差是有偏差的，因此基于標準誤差的統(tǒng)計學意義的主張會有偏差。更重要的是，我們無法探索學生效應，但學生效應很有意義。
分組回歸
另一種方法是對每個學生單獨進行回歸。然而，這種方法有很多缺點--當有很多組的時候，它不容易被總結，通常每個組內的數據很少，無法做到這一點（如在這個案例中），而且模型是過度的背景化，意味著他們忽略了學生的共同點。我們將在后面比較這樣的方法和混合模型。
運行混合模型
接下來我們運行一個混合模型，它將允許學生的特定效應。在下文中，代碼看起來就像你用lm做的回歸一樣，但有一個額外的部分來指定組，即學生的效應。(1|student)意味著我們允許截距（用1表示）因學生而異。使用混合模型，我們可以得到與回歸相同的結果，但將有更多的內容可以討論。
library(lme4)gpa_mixed = lmer(gpa ~ occasion + (1 | student), data = gpa)## summary(gpa_mixed)
?
首先我們看到，截距和時間的系數，即在這里可以稱為固定效應，與我們在標準回歸中看到的相同，其解釋也相同。另一方面，這里的標準誤差是不同的，盡管最后我們的結論在統(tǒng)計意義上是相同的。請?zhí)貏e注意，截距的標準誤差已經增加。從概念上講，你可以認為允許每個人的隨機截距允許我們獲得關于個人的信息，同時認識到關于總體平均數的不確定性。
雖然我們有系數和標準誤差，但你可能已經注意到，lme4并沒有提供p值，這有幾個原因。這有幾個原因，即對于混合模型，我們基本上是在處理不同的樣本量，群組內的Nc，可能因群組而異（甚至是單個觀測值），以及N個總觀測值，這使我們在參考分布、自由度以及如何近似 "最佳 "解決方案方面處于一種模糊的狀態(tài)。其他程序自動提供p值，好像沒有什么問題，而且沒有告訴你他們使用哪種方法來計算p值（有幾種）。此外，這些近似值在某些情況下可能非常差，或者做出可能不適合情況的假設。
然而，得到置信區(qū)間是比較直接的，如下所示7。
confint(gpa)
?
方差成分
與標準回歸輸出相比，有一點是新的，那就是學生效應的估計方差/標準差（在我們之前的公式描述中是ττ）。這告訴我們，平均而言，當我們從一個學生轉移到另一個學生時，GPA會有多大的變動。換句話說，即使在根據時間點進行預測后，每個學生都有自己獨特的偏差，而這個值（就標準偏差而言）是整個學生的估計平均偏差。
另一種解釋方差輸出的方法是注意學生方差在總數中的百分比，或0.064 / 0.122 = 52%。這也被稱為類內相關，因為它也是對群內相關的估計，我們將在后面看到。
隨機效應的估計
運行模型后，我們實際上可以得到學生效應的估計值。我為前五個學生展示了兩種方法，既可以作為隨機效應，也可以作為隨機截距（即截距+隨機效應）。
ef(mixed)$student %>% head(5)
?
coef
?
請注意，我們不允許學期的變化，所以它對所有的學生都是一個恒定的，也就是固定的效果。
通常情況下，我們對這些效應非常感興趣，并希望對它們有一些不確定性的感覺?？梢酝ㄟ^預測區(qū)間來實現?；蛘呖梢灾苯尤タ此鼈兊膱D。
Interval(mixed) ? # 用于各種模型的預測，可能使用新的數據sim(mixed) ? ? ? ? ? ? # ?隨機效應估計值的平均值、中位數和SD值plot(mixed)) ?# 繪制區(qū)間估計值
下面的圖是每個學生的估計隨機效應及其區(qū)間估計。隨機效應是正態(tài)分布，其平均值為零，由水平線表示。不包括零的區(qū)間用粗體表示。
預測
現在讓我們來看看標準預測與特定群組預測的對比。與大多數R模型一樣，我們可以在模型上使用預測函數。
predict(mixed, re.form=NA) %>% head()
在上面的代碼中，我們指定不使用隨機效應re.form=NA，因此，我們對觀測值的預測和我們從標準線性模型中得到的預測差不多。
predict(mixed, re.form=NA)predict(lm)
但是每個人都有自己獨特的截距，所以讓我們看看當我們加入這些信息時，預測結果有什么不同。
predict(mixed)
?
根據估計的學生效應，學生開始于高或低于所有學生的估計截距。下面是無條件預測與包含了前兩個學生的隨機截距的條件預測的直觀對比。
plot(x = occasion,y = gpa, color = student,prediction, group = student,y = prediction)
我們可以看到，由于截距不同，混合模型的預測結果發(fā)生了偏移。對于這些學生來說，這種轉變反映了他們相對較差的起點。
聚類層面的協變量
請注意我們把混合模型描述為一個多層次模型。
?
如果我們在模型中加入學生層面的協變量，例如性別，我們就會有以下結果。
其中，在插入后，我們仍然有與之前相同的模型，只是增加了一個預測因素。
因此，增加群組級協變量對我們思考模型的方式沒有任何不尋常的影響。我們只是把它們添加到我們的預測變量集合中。還要注意的是，我們可以將聚類層面的協變量創(chuàng)建為平均值或其他一些觀察層面變量的總結。當聚類代表地理單位，而觀察對象是人時，這一點尤其常見。例如，我們可以將收入作為個人層面的協變量，并使用中位數來代表地理區(qū)域的整體財富。
混合模型基礎知識總結
混合模型使我們能夠考慮到數據中的聚類。我們更好地理解了目標變量的變異性來源。我們還得到了模型中參數的具體組別估計，使我們能夠準確地了解各組之間的差異。此外，這反過來又允許我們進行特定群體的預測，從而進行更準確的預測，假設由于聚類而存在明顯的變異。簡而言之，即使在最簡單的情況下，混合模型也有很多好處。
練習
睡眠
在這個練習中，我們將使用lme4軟件包中的睡眠研究數據。以下是對它的描述。
睡眠限制研究中的受試者每天的平均反應時間。在第0天，受試者有正常的睡眠量。從那天晚上開始，他們被限制在每晚3小時的睡眠時間。觀察結果代表了每天給每個受試者的一系列測試的平均反應時間（以毫秒計）。
在加載軟件包后，可以按以下方式加載數據。我展示了最初的幾個觀察結果。
運行一個以因變量為目標變量，以天為預測變量的回歸。
運行一個混合模型，用隨機截距表示Subject。
解釋方差分量和固定效應。
添加聚類級協變量
用GPA數據重新運行混合模型，添加聚類級協變量性別或高中GPA（highgpa），或兩者。解釋結果的所有方面。
在模型中加入聚類層面的協變量后，學生的方差發(fā)生了什么變化？
模擬混合模型
下面代表了一種模擬隨機截距模型的簡單方法。
set.seed(1234) ?# 復制你的結果N = Ngroups * NperGroupy = 2 + .5 * x + u[groups] + e
以上哪些代表了固定效應和隨機效應？現在運行混合模型。
結果與你所期望的一致嗎？
在下面的內容中，我們將改變數據的各個方面，然后在每次改變后重新運行模型，然后像以前一樣進行總結并得到置信區(qū)間。對于每一項，都要特別注意結果中至少有一點變化。
?首先計算類內相關系數。
.
此外，創(chuàng)建隨機效應的密度圖。
改變隨機效應方差/SD和/或殘差/SD，注意你對ICC的新估計，并像以前一樣繪制隨機效應。
將數值重置為原始值。將Ngroups改為50。你在置信區(qū)間估計中看到了什么不同？
將Ngroups重新設置為100?，F在將NperGroup改為10，并再次注意到置信區(qū)間與基本條件有什么不同。
最受歡迎的見解
1.基于R語言的lmer混合線性回歸模型
2.R語言用Rshiny探索lme4廣義線性混合模型（GLMM）和線性混合模型（LMM）
3.R語言線性混合效應模型實戰(zhàn)案例
4.R語言線性混合效應模型實戰(zhàn)案例2
5.R語言線性混合效應模型實戰(zhàn)案例
6.線性混合效應模型Linear Mixed-Effects Models的部分折疊Gibbs采樣
7.R語言LME4混合效應模型研究教師的受歡迎程度
8.R語言中基于混合數據抽樣(MIDAS)回歸的HAR-RV模型預測GDP增長
9.使用SAS，Stata，HLM，R，SPSS和Mplus的分層線性模型HLM
