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

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

混合效應(yīng)邏輯回歸用于建立二元結(jié)果變量的模型，其中，當(dāng)數(shù)據(jù)被分組或同時(shí)存在固定和隨機(jī)效應(yīng)時(shí)，結(jié)果的對(duì)數(shù)幾率被建模為預(yù)測(cè)變量的線性組合。

混合效應(yīng)邏輯回歸的例子

例1：一個(gè)研究人員對(duì)40所不同大學(xué)的申請(qǐng)進(jìn)行抽樣調(diào)查，以研究預(yù)測(cè)大學(xué)錄取的因素。預(yù)測(cè)因素包括學(xué)生的高中GPA、課外活動(dòng)和SAT分?jǐn)?shù)。一些學(xué)校的選擇性較多或較少，所以每所學(xué)校的基準(zhǔn)錄取概率是不同的。學(xué)校層面的預(yù)測(cè)因素包括學(xué)校是公立還是私立，目前學(xué)生與教師的比例，以及學(xué)校的排名。

例2：一家大型HMO想知道哪些病人和醫(yī)生的因素與病人的肺癌在治療后是否得到緩解最相關(guān)，這是一項(xiàng)關(guān)于肺癌病人的治療效果和生活質(zhì)量的研究的一部分。

例3：一家電視臺(tái)想知道時(shí)間和廣告活動(dòng)如何影響人們是否觀看電視節(jié)目。他們對(duì)四個(gè)城市的人進(jìn)行了為期六個(gè)月的抽樣調(diào)查。每個(gè)月，他們都會(huì)詢問(wèn)人們?cè)谶^(guò)去一周是否觀看了某個(gè)節(jié)目。三個(gè)月后，他們?cè)谒膫€(gè)城市中的兩個(gè)城市推出了一個(gè)新的廣告活動(dòng)，并繼續(xù)監(jiān)測(cè)人們是否觀看了該節(jié)目。

數(shù)據(jù)描述

在這個(gè)例子中，我們將使用一個(gè)數(shù)據(jù)集來(lái)探討關(guān)于肺癌的例子。我們收集了病人的各種結(jié)果，他們被包含在醫(yī)生身上，而醫(yī)生又被包含在醫(yī)院里。還有一些醫(yī)生層面的變量，比如我們將在例子中使用的 "醫(yī)生經(jīng)驗(yàn)"。

現(xiàn)在我們要對(duì)我們的連續(xù)預(yù)測(cè)變量進(jìn)行繪圖。數(shù)據(jù)的可視化可以幫助我們理解分布情況，發(fā)現(xiàn)編碼錯(cuò)誤（例如，我們知道一個(gè)變量的取值范圍是0到7，但我們?cè)趫D中看到了999），并讓我們了解變量之間的關(guān)系。例如，我們可能看到兩個(gè)預(yù)測(cè)因子高度相關(guān)，于是決定只在模型中包括一個(gè)，或者我們可能注意到兩個(gè)變量之間有曲線關(guān)系。數(shù)據(jù)可視化是一種快速、直觀的方式，可以一次性檢查所有這些情況。如果你的大多數(shù)預(yù)測(cè)因子看起來(lái)都是相互獨(dú)立的，數(shù)據(jù)很好。例如，如果它們是獨(dú)立的，當(dāng)你輸入另一個(gè)預(yù)測(cè)因子時(shí)，一個(gè)預(yù)測(cè)因子的估計(jì)值不應(yīng)該有太大變化（盡管標(biāo)準(zhǔn)誤差和顯著性檢驗(yàn)可能會(huì)有）。我們可以通過(guò)簡(jiǎn)單地查看數(shù)據(jù)來(lái)了解所有這些信息以及判斷如何建模。

我們的連續(xù)預(yù)測(cè)因子之間似乎沒(méi)有強(qiáng)的線性關(guān)系。讓我們看看我們的變量在癌癥階段中的分布情況。因?yàn)樽≡簳r(shí)間是以天為單位的，我們可以用氣泡圖來(lái)研究癌癥階段與它的關(guān)系。每個(gè)氣泡的面積與具有這些數(shù)值的觀察值的數(shù)量成正比。對(duì)于連續(xù)的預(yù)測(cè)因子，我們使用小提琴圖。所有的原始數(shù)據(jù)都按癌癥階段分開(kāi)顯示。

我們?cè)黾恿诵√崆賵D。小提琴圖只是圍繞繪圖軸反映的核密度圖。我們將小提琴圖繪制在具有透明度的抖動(dòng)點(diǎn)之上，這樣就可以看到原始數(shù)據(jù)。

因?yàn)镮L6和CRP都有偏斜分布的傾向，所以我們?cè)赮軸上使用了平方根刻度。分布看起來(lái)相當(dāng)正常和對(duì)稱，你仍然可以看到長(zhǎng)的右尾，即使使用了平方根刻度（注意，只有刻度被轉(zhuǎn)移，數(shù)值本身沒(méi)有被轉(zhuǎn)換，這很重要，因?yàn)檫@讓你看到并解釋實(shí)際的分?jǐn)?shù)，而不是分?jǐn)?shù)的平方根）。

因?yàn)楹茈y看到二元變量在連續(xù)變量的水平上如何變化，我們可以反過(guò)來(lái)看看二元結(jié)果的每個(gè)水平上的連續(xù)變量的分布。

分析方法

下面是一個(gè)分析方法的列表：

• 混合效應(yīng)邏輯回歸，是本頁(yè)面的重點(diǎn)。

• 混合效應(yīng)probit回歸與混合效應(yīng)logistic回歸非常相似，但它使用的是正態(tài)CDF而不是logistic CDF。兩者都對(duì)二元結(jié)果進(jìn)行建模，可以包括固定和隨機(jī)效應(yīng)。

• 固定效應(yīng)邏輯回歸在這種情況下是有限的，因?yàn)樗赡芎雎粤吮匾碾S機(jī)效應(yīng)和/或數(shù)據(jù)中的非獨(dú)立性。

• 固定效應(yīng)的probit回歸在這種情況下是有限的，因?yàn)樗赡芎雎粤吮匾碾S機(jī)效應(yīng)或數(shù)據(jù)中的非獨(dú)立性。

• 有聚類穩(wěn)健標(biāo)準(zhǔn)差的Logistic回歸。這些可以調(diào)整非獨(dú)立性，但不允許有隨機(jī)效應(yīng)。

• 有聚類穩(wěn)健標(biāo)準(zhǔn)差的Probit回歸。這些可以調(diào)整非獨(dú)立性，但不允許有隨機(jī)效應(yīng)。

混合效應(yīng)邏輯回歸

下面我們使用glmer命令估計(jì)混合效應(yīng)邏輯回歸模型，Il6、CRP和住院時(shí)間為患者水平的連續(xù)預(yù)測(cè)因素，癌癥階段為患者水平的分類預(yù)測(cè)因素（I、II、III或IV），經(jīng)驗(yàn)為醫(yī)生水平的連續(xù)預(yù)測(cè)因素，還有DID的隨機(jī)截距，醫(yī)生ID。

第一部分告訴我們，估計(jì)值是基于自適應(yīng)高斯-赫米特的似然性近似。為了避免出現(xiàn)不收斂的警告，我們用參數(shù)control=glmerControl(optimizer="bobyqa")指定不同的優(yōu)化器。

下一節(jié)給我們提供了可用于比較模型的基本信息，接著是隨機(jī)效應(yīng)估計(jì)值。這表示對(duì)數(shù)尺度上截距的估計(jì)變化。如果有其他隨機(jī)效應(yīng)，比如隨機(jī)斜率，它們也會(huì)出現(xiàn)在這里。最上面的部分最后是觀察值的總數(shù)和第2級(jí)觀察值的數(shù)量。在我們的案例中，這包括病人（8,525）和醫(yī)生（407）的總數(shù)。

最后一節(jié)是固定效應(yīng)估計(jì)值的表格。這些估計(jì)值代表回歸系數(shù)。這些是未標(biāo)準(zhǔn)化的，而且是在對(duì)數(shù)尺度上。估計(jì)值后面是它們的標(biāo)準(zhǔn)誤差（SE）。系數(shù)估計(jì)的近似值可能比SEs的近似值穩(wěn)定得更快。Wald檢驗(yàn)，（frac{Estimate}{SE}），依賴于漸進(jìn)理論，這里指的是當(dāng)最高級(jí)別的單位大小收斂到無(wú)窮大時(shí)，這些檢驗(yàn)將呈正態(tài)分布，并由此得出p值（鑒于真實(shí)估計(jì)值為0，獲得觀察估計(jì)值或更極端的概率）。

獲得置信區(qū)間（CI）。我們可以使用SE來(lái)獲得粗略的區(qū)間估計(jì)。

如果我們需要比值比而不是對(duì)數(shù)刻度上的系數(shù)，則可以對(duì)估計(jì)值和CI求冪。

多層bootstrapping（自助法）

從GLMMs進(jìn)行推斷是很復(fù)雜的。除了在每個(gè)層次（尤其是最高層次）有很多觀測(cè)值的情況下，假設(shè)（frac{Estimate}{SE}）是正態(tài)分布可能不準(zhǔn)確。人們提出了各種替代方法，包括蒙特卡洛模擬、貝葉斯估計(jì)和bootstrapping。每種方法的實(shí)施都可能很復(fù)雜。我們將重點(diǎn)討論一個(gè)小的bootstrapping例子。

Bootstrapping是一種重抽樣方法，就是利用有限的樣本資料經(jīng)由多次重復(fù)抽樣，重新建立起足以代表母體樣本分布的新樣本。它決不是完美的，但它在概念上是直接易懂的，而且容易在代碼中實(shí)現(xiàn)。一個(gè)缺點(diǎn)是，它對(duì)計(jì)算要求很高。對(duì)于大型數(shù)據(jù)集或復(fù)雜的模型，每個(gè)模型的運(yùn)行需要幾分鐘，在成千上萬(wàn)的樣本上進(jìn)行估計(jì)，很容易需要幾個(gè)小時(shí)或幾天。在本頁(yè)的例子中，我們使用了非常少的樣本，但在實(shí)踐中你會(huì)使用更多的樣本。

對(duì)于單層次模型，我們可以實(shí)現(xiàn)簡(jiǎn)單的隨機(jī)抽樣，并進(jìn)行替換，以進(jìn)行bootstrapping。對(duì)于多層次數(shù)據(jù)，我們希望以與數(shù)據(jù)生成機(jī)制相同的方式重新取樣。我們從最高級(jí)別開(kāi)始重新取樣，然后逐級(jí)向下。在我們的案例中，我們首先將從醫(yī)生那里取樣，然后在每個(gè)取樣的醫(yī)生中，我們將從他們的病人那里取樣。要做到這一點(diǎn)，我們首先需要寫(xiě)一個(gè)函數(shù)，在每個(gè)層次上重新取樣。

現(xiàn)在，我們將重新對(duì)數(shù)據(jù)進(jìn)行取樣，并采取100次重復(fù)。同樣在實(shí)踐中，你可能會(huì)采取數(shù)千次。我們?cè)O(shè)置種子，以便我們的結(jié)果可以重復(fù)。你也很可能需要比你最終想要的更多的重復(fù)樣本，因?yàn)樵S多樣本可能不收斂，所以你不能從它們那里得到估計(jì)。

接下來(lái)，我們?cè)谥匦氯拥臄?shù)據(jù)上重新擬合模型。首先，我們存儲(chǔ)原始模型的估計(jì)值，我們將用它作為自助模型的起始值。然后，我們建立一個(gè)有4個(gè)節(jié)點(diǎn)的本地集群。接下來(lái)，我們導(dǎo)出數(shù)據(jù)并在集群上加載。最后，我們寫(xiě)一個(gè)函數(shù)來(lái)擬合模型并返回估計(jì)值。對(duì)glmer()的調(diào)用被封裝在try中，因?yàn)椴皇撬械哪Ｐ投寄茉谥匦虏蓸拥臄?shù)據(jù)上收斂。這樣可以捕捉到錯(cuò)誤并返回，而不是停止處理。

現(xiàn)在我們已經(jīng)有了數(shù)據(jù)、本地集群和擬合函數(shù)的設(shè)置，我們準(zhǔn)備實(shí)際進(jìn)行bootstrapping了。來(lái)自所有節(jié)點(diǎn)的結(jié)果被匯總回一個(gè)單一的列表，存儲(chǔ)在對(duì)象res中。一旦完成，我們就可以關(guān)閉本地集群，終止額外的R實(shí)例并釋放了內(nèi)存。

現(xiàn)在我們有了bootstrapping法的結(jié)果，我們可以對(duì)其進(jìn)行總結(jié)。首先，我們計(jì)算成功收斂的模型的數(shù)量。我們可以計(jì)算成功的平均數(shù)，以看到收斂的比例。

接下來(lái)我們把引導(dǎo)結(jié)果列表轉(zhuǎn)換成矩陣，然后計(jì)算每個(gè)參數(shù)的2.5和97.5百分位。最后，我們可以將結(jié)果制成表格，包括原始估計(jì)值和標(biāo)準(zhǔn)誤差、平均引導(dǎo)估計(jì)值，以及bootstrap 的置信區(qū)間。

預(yù)測(cè)的概率和繪圖

這些結(jié)果很適合放在表格中或研究文本中；但是，數(shù)字的解釋可能很麻煩。圖形展示有助于解釋，也有助于演講。
在一個(gè)邏輯模型中，結(jié)果通常是

• 對(duì)數(shù)幾率（也叫對(duì)數(shù)），這是線性化

• 指數(shù)化的對(duì)數(shù)幾率，不在線性尺度上

• 概率

對(duì)于表格來(lái)說(shuō)，人們經(jīng)常呈現(xiàn)的是幾率比。對(duì)于可視化來(lái)說(shuō)，對(duì)數(shù)或概率比例是最常見(jiàn)的。每種方法都有一些優(yōu)點(diǎn)和缺點(diǎn)。對(duì)數(shù)表很方便，因?yàn)樗蔷€性化的，這意味著一個(gè)預(yù)測(cè)因素增加1個(gè)單位，結(jié)果就會(huì)增加一個(gè)系數(shù)單位，而且無(wú)論其他預(yù)測(cè)因素的水平如何。缺點(diǎn)是這個(gè)量表的可解釋性不強(qiáng)。讀者很難對(duì)對(duì)數(shù)有一個(gè)直觀的理解。相反，概率是一個(gè)很好的尺度，可以直觀地理解結(jié)果；但是，它們不是線性的。這意味著預(yù)測(cè)因子增加一個(gè)單位，不等于概率的恒定增加--概率的變化取決于為其他預(yù)測(cè)因子選擇的值。在普通邏輯回歸中，你可以保持所有預(yù)測(cè)因子不變，只改變你感興趣的預(yù)測(cè)因子。然而，在混合效應(yīng)邏輯模型中，隨機(jī)效應(yīng)也對(duì)結(jié)果產(chǎn)生影響。因此，如果你保持一切不變，那么只有當(dāng)所有協(xié)變量保持不變，并且你在同一組或具有相同隨機(jī)效應(yīng)的一組時(shí)，結(jié)果的概率變化才是真的。

我們將探討一個(gè)平均邊際概率的例子。這比條件概率需要更多的工作，因?yàn)槟惚仨殲槊恳唤M計(jì)算單獨(dú)的條件概率，然后將其平均化。

首先，讓我們使用這里的符號(hào)來(lái)定義一般程序。我們通過(guò)獲取?

并將感興趣的特定預(yù)測(cè)因子，比如說(shuō)在j列，設(shè)置為常數(shù)來(lái)創(chuàng)建?

。如果我們只關(guān)心預(yù)測(cè)器的一個(gè)值，那就是

。然而，更常見(jiàn)的是，我們希望預(yù)測(cè)因子有一定的取值范圍，以便繪制預(yù)測(cè)概率在其范圍內(nèi)的變化情況。我們可以通過(guò)獲取預(yù)測(cè)模型的觀察范圍，并在該范圍內(nèi)均勻地抽取k個(gè)樣本。例如，假設(shè)我們的預(yù)測(cè)模型的范圍是5到10，我們想要6個(gè)樣本，

，所以每個(gè)樣本將與前一個(gè)樣本相隔1，它們將是

. 然后我們創(chuàng)建不同的k個(gè)不同的Xi，其中

，在每種情況下，第j列被設(shè)置為某個(gè)常數(shù)。然后我們計(jì)算：

這些是所有不同的線性預(yù)測(cè)因子。最后，我們采取

，這就得到?

，這是原始尺度上的條件期望，在我們的例子中是概率。然后我們可以取每個(gè)

的期望值，并將其與我們感興趣的預(yù)測(cè)因子的值作對(duì)比。我們還可以繪制圖表，不僅顯示平均邊際預(yù)測(cè)概率，而且還顯示預(yù)測(cè)概率的分布。

你可能已經(jīng)注意到，這些估計(jì)值中有很多變數(shù)。我們?cè)谑褂?

時(shí)，只將我們感興趣的預(yù)測(cè)因子保持在一個(gè)常數(shù)，這使得所有其他預(yù)測(cè)因子都能在原始數(shù)據(jù)中取值。另外，我們把?

留在我們的樣本中，這意味著有些組的代表性比其他組要高或低。如果我們想的話，我們可以對(duì)所有的群體進(jìn)行重新加權(quán)，使其具有同等的權(quán)重。在這個(gè)例子中，我們選擇讓所有這些東西保持原樣，是基于這樣的假設(shè)：我們的樣本確實(shí)是我們感興趣的人群的良好代表。我們沒(méi)有試圖挑選有意義的值來(lái)保持協(xié)變量（，而是使用了我們樣本的值。這也表明，如果我們的樣本能很好地代表總體，那么平均邊際預(yù)測(cè)概率就能很好地代表我們總體中新的隨機(jī)樣本的概率。

現(xiàn)在我們有了一些背景和理論，我們看看如何實(shí)際去計(jì)算這些東西。我們得到一個(gè)住院時(shí)間（我們感興趣的預(yù)測(cè)因子）的摘要，然后在其范圍內(nèi)得到100個(gè)值，用于預(yù)測(cè)。我們復(fù)制一份數(shù)據(jù)，這樣我們就可以固定其中一個(gè)預(yù)測(cè)因子的值，然后使用預(yù)測(cè)函數(shù)來(lái)計(jì)算預(yù)測(cè)值。默認(rèn)情況下，所有的隨機(jī)效應(yīng)都被包括在內(nèi)。?

現(xiàn)在我們有了所有的預(yù)測(cè)概率， 可以可視化它們。例如，我們可以看一下少數(shù)不同停留時(shí)間的平均邊際預(yù)測(cè)概率。

我們還可以加上下限和四分位數(shù)?？梢钥吹?0%的預(yù)測(cè)概率所處的范圍。

除了改變住院時(shí)間之外，我們還可以對(duì)癌癥階段的每一級(jí)做同樣的平均邊際預(yù)測(cè)概率。

對(duì)于一個(gè)住院10天的第四期肺癌患者，其癌癥得到緩解的機(jī)會(huì)看起來(lái)相當(dāng)小。看起來(lái)分布也是偏斜的。我們可以檢查一下僅針對(duì)該組的預(yù)測(cè)概率分布。

即使使用平方根尺度，將較低的數(shù)值拉長(zhǎng)，它仍然是極其偏斜的。據(jù)估計(jì)，絕大多數(shù)人的病情緩解的概率不到0.1。

三層混合效應(yīng)邏輯回歸

我們已經(jīng)深入研究了一個(gè)帶有隨機(jī)截距的兩級(jí)邏輯模型。這是最簡(jiǎn)單的混合效應(yīng)邏輯模型?，F(xiàn)在我們要簡(jiǎn)要地看一下如何增加第三層次和隨機(jī)斜率效應(yīng)以及隨機(jī)截距。

下面我們估計(jì)一個(gè)三層邏輯模型，醫(yī)生有一個(gè)隨機(jī)截距，醫(yī)院有一個(gè)隨機(jī)截距。在這個(gè)例子中，醫(yī)生被嵌套在醫(yī)院內(nèi)，也就是說(shuō)，每個(gè)醫(yī)生屬于一家而且只有一家醫(yī)院。另一種情況有時(shí)被稱為 "交叉分類"，意思是一個(gè)醫(yī)生可能屬于多家醫(yī)院，比如該醫(yī)生的一些病人來(lái)自A醫(yī)院，另一些來(lái)自B醫(yī)院。在glmer中，你不需要指定組是嵌套還是交叉分類，R可以根據(jù)數(shù)據(jù)計(jì)算出來(lái)。

輸出告訴我們族（二元結(jié)果的二項(xiàng)式）和鏈接函數(shù)（logit）。接著是通常的擬合指數(shù)和隨機(jī)效應(yīng)的方差。在這種情況下，醫(yī)生之間和醫(yī)院之間的截距（在對(duì)數(shù)賠率尺度上）的變化。還顯示了標(biāo)準(zhǔn)差（只是方差的平方根，而不是估計(jì)方差的標(biāo)準(zhǔn)誤差）。我們還得到了每個(gè)層次上的單位的數(shù)量。最后是固定效應(yīng)。

看一下條件模型的分布也是很有用的，下面我們用“毛毛蟲(chóng)圖”來(lái)做。藍(lán)點(diǎn)是帶有誤差條的條件模型。我們對(duì)醫(yī)生和醫(yī)院都是這樣做的。例如，對(duì)于醫(yī)生來(lái)說(shuō)，我們可以看到一個(gè)有點(diǎn)長(zhǎng)的右尾，即極端的正值比負(fù)值多。

我們也可以在模型中加入隨機(jī)斜率。我們只是要為 "住院時(shí)間 "增加一個(gè)隨機(jī)斜率，這個(gè)斜率在不同的醫(yī)生之間變化。就像在常規(guī)的R公式中一樣，我們使用+運(yùn)算符來(lái) "添加 "一個(gè)效應(yīng)。

最受歡迎的見(jiàn)解

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

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

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

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

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

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

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

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

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