R語(yǔ)言混合效應(yīng)邏輯回歸（mixed effects logistic）模型分析肺癌數(shù)據(jù)|附代碼數(shù)據(jù)
下載全文：http://tecdat.cn/?p=22302
在混合效應(yīng)邏輯回歸用于建立二元結(jié)果變量的模型，其中，當(dāng)數(shù)據(jù)被分組或同時(shí)存在固定和隨機(jī)效應(yīng)時(shí)，結(jié)果的對(duì)數(shù)幾率被建模為預(yù)測(cè)變量的線性組合（點(diǎn)擊文末“閱讀原文”獲取完整代碼數(shù)據(jù)）。
混合效應(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)"。
within(hp,?{??已婚?<-?factor(已婚,?levels?=?0:1,?labels?=?c("否",?"是"))??DID?<-?factor(DID)??HID?<-?factor(HID)??癌癥階段?<-?factor(癌癥階段)
現(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)了解所有這些信息以及判斷如何建模。
ggpairs(hp\[,?c("IL6",?"CRP",?"住院時(shí)間",?"醫(yī)生經(jīng)驗(yàn)")\],
我們的連續(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)顯示。
點(diǎn)擊標(biāo)題查閱往期內(nèi)容
R語(yǔ)言建立和可視化混合效應(yīng)模型mixed effect model
轉(zhuǎn)存失敗重新上傳取消
左右滑動(dòng)查看更多
轉(zhuǎn)存失敗重新上傳取消
01
02
03
04
我們?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ù)的平方根）。
ggplot(????stat_sum(aes(size?=?..n..,?group?=?1))?+??scale\_size\_area(max_size=10)
ggplot(tp,??+??geom_jitter(alpha?=?.1)?+??geom_violin(alpha?=?.75)?+
因?yàn)楹茈y看到二元變量在連續(xù)變量的水平上如何變化，我們可以反過(guò)來(lái)看看二元結(jié)果的每個(gè)水平上的連續(xù)變量的分布。
ggplot(tmp,???geom_boxplot()?+??facet\_wrap(~variable,?scales="free\_y")
分析方法
下面是一個(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ì)模型并將結(jié)果存儲(chǔ)在m中#?輸出結(jié)果，固定效果之間不相關(guān)print(m,?corr?=?FALSE)
第一部分告訴我們，估計(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ì)。
#?帶有95％CI的估計(jì)表?cbind(Est?=?fixef(m),?LL?=?fixef(m)?-?1.96?*?se,?UL?=?fixef(m)?+?1.96?*????se))
如果我們需要比值比而不是對(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è)層次上重新取樣。
sample?<-?function(dat,?cluvar)?{????id?<-?sample(cid,?size?=??cid?*?reps,?replace?=?TRUE)????if?(replace)?{???????lapply(seq_along(cid),?function(i)?{????????????cbind(?ID?=?i,?RowID?=?sample(?(dat\[,?cluvar\]?==?recid\[i\]),????????????????size?=?length(?(dat\[,?cluvar\]?==?recid\[i\])?)????????})????}?else?{???????lapply(seq_along(cid),?function(i)?{????????????cbind(?ID?=?i,?RowID?=??(dat\[,?cluvar\]?==?recid\[i\]))???????)????}?)
現(xiàn)在，我們將重新對(duì)數(shù)據(jù)進(jìn)行取樣，并采取100次重復(fù)。同樣在實(shí)踐中，你可能會(huì)采取數(shù)千次。我們?cè)O(shè)置種子，以便我們的結(jié)果可以重復(fù)。你也很可能需要比你最終想要的更多的重復(fù)樣本，因?yàn)樵S多樣本可能不收斂，所以你不能從它們那里得到估計(jì)。
tmp?<-?sample(hdp,?"DID",?reps?=?100)
接下來(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ò)誤并返回，而不是停止處理。
makeCluster(4)clusterEvalQ(cl,?require(lme4))
boot?<-?function(i)?{
????object?<-?try(glmer(緩解?~?IL6?+?CRP?+?癌癥階段?+?住院時(shí)間??subset?=?Replicate?==?i,?family?=?binomial,
?????if?(class(object)?==?"try-error")
????????return(object)
}
現(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)存。
?parLapplyLB(cl,?X?=?levels(data$Replicate),?fun?=?boot)#?關(guān)閉集群stopCluster(cl)
現(xiàn)在我們有了bootstrapping法的結(jié)果，我們可以對(duì)其進(jìn)行總結(jié)。首先，我們計(jì)算成功收斂的模型的數(shù)量。我們可以計(jì)算成功的平均數(shù)，以看到收斂的比例。
#?計(jì)算成功收斂的模型的比例succ??<-?sapply(res,?is.numeric)mean(succ??)
##?\[1\]?1
接下來(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ū)間。
#?合并成功的結(jié)果?do.call(cbind,?res\[success\])#?計(jì)算出95%CI的2.5和97.5百分位數(shù)。?apply(?res,?1，quantile,?probs?=?c(0.025,?0.975)
##???????????????????????2.5%?????97.5%
##?(Intercept)???????-3.61982?-0.985404
##?IL6???????????????-0.08812?-0.029664
##?CRP???????????????-0.04897??0.006824
##?癌癥階段II?????-0.60754?-0.228019
##?癌癥階段III????-1.30217?-0.754609
##?癌癥階段IV?????-2.91414?-2.002643
##?住院時(shí)間??????-0.21596?-0.046420
##?醫(yī)生經(jīng)驗(yàn)?????????0.06819??0.207223
##?NewID.(Intercept)??2.03868??2.476366
#?所有結(jié)果取整輸出
round(?able,?3)
##????????????????????Est????SE?BootMean???2.5%??97.5%
##?(Intercept)?????-2.053?0.531???-2.205?-3.620?-0.985
##?IL6?????????????-0.057?0.012???-0.059?-0.088?-0.030
##?CRP?????????????-0.021?0.010???-0.022?-0.049??0.007
##?癌癥階段II???-0.414?0.076???-0.417?-0.608?-0.228
##?癌癥階段III??-1.003?0.098???-1.043?-1.302?-0.755
##?癌癥階段IV???-2.337?0.158???-2.460?-2.914?-2.003
##?住院時(shí)間????-0.121?0.034???-0.142?-0.216?-0.046
##?醫(yī)生經(jīng)驗(yàn)???????0.120?0.027????0.128??0.068??0.207
##?DID.(Intercept)??2.015????NA????2.263??2.039??2.476
預(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)。?
?#?計(jì)算預(yù)測(cè)的概率并存儲(chǔ)在列表中?lapply(jvalues,?function(j)?{????predict(m,?newdata?=?tmpdat,?type?=?"response")
現(xiàn)在我們有了所有的預(yù)測(cè)概率， 可以可視化它們。例如，我們可以看一下少數(shù)不同停留時(shí)間的平均邊際預(yù)測(cè)概率。
#幾種不同時(shí)間的平均邊際預(yù)測(cè)概率
?t(sapply(pp,?function(x)?{????c(M?=?mean(x),?quantile(x,?c(0.25,?0.75)))#?加入時(shí)間的數(shù)值并轉(zhuǎn)換為數(shù)據(jù)框??as.data.frame(cbind(plotdat,??values))#?顯示前幾行head(plotdat
#?繪制平均邊際預(yù)測(cè)概率ggplot(plotdat)?+?geom_line()?+
我們還可以加上下限和四分位數(shù)?？梢钥吹?0%的預(yù)測(cè)概率所處的范圍。
ggplot(plotdat)?+?geom\_linerange(?)?+?geom\_line(size?=?2)
除了改變住院時(shí)間之外，我們還可以對(duì)癌癥階段的每一級(jí)做同樣的平均邊際預(yù)測(cè)概率。
#?計(jì)算預(yù)測(cè)的概率并存儲(chǔ)在一個(gè)列表中?lapply(?癌癥階段,?function(stage)?{????predict(m,?newdata?=?tmpdat,?type?=?"response")??})#?獲得每個(gè)級(jí)別的癌癥階段的所有j值的平均值和四分位數(shù)?lapply(?probs,?function(X)?{????c(M=mean(x),?quantile(x,?c(.25,?.75)))#?放到一個(gè)數(shù)據(jù)框??do.call(rbind,?plotdat2)#?添加癌癥階段?factor(rep(levels(癌癥階段?))#?顯示前幾行數(shù)據(jù)head(plotdat2)
#?繪制ggplot(plotdat2?)?+??geom_ribbon(aes(ymin?=?Lower,?ymax?=?Upper,?fill?=?癌癥階段),?alpha?=?.15)?+
對(duì)于一個(gè)住院10天的第四期肺癌患者，其癌癥得到緩解的機(jī)會(huì)看起來(lái)相當(dāng)小?？雌饋?lái)分布也是偏斜的。我們可以檢查一下僅針對(duì)該組的預(yù)測(cè)概率分布。
ggplot(??aes(Probs))?+?geom_histogram()?+
即使使用平方根尺度，將較低的數(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)。
#?輸出沒(méi)有固定效應(yīng)之間相關(guān)性的mod結(jié)果print(m3a,?corr=FALSE)
輸出告訴我們族（二元結(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ù)值多。
lattice::dotplot(?(m)?))##?$DID
##?$HID
我們也可以在模型中加入隨機(jī)斜率。我們只是要為 "住院時(shí)間 "增加一個(gè)隨機(jī)斜率，這個(gè)斜率在不同的醫(yī)生之間變化。就像在常規(guī)的R公式中一樣，我們使用+運(yùn)算符來(lái) "添加 "一個(gè)效應(yīng)。
#?輸出沒(méi)有固定效應(yīng)之間相關(guān)性的mod結(jié)果print(m3b,?corr?=?FALSE)
##?$DID
##?$HID
點(diǎn)擊文末“閱讀原文”
獲取全文完整資料。
本文選自《R語(yǔ)言混合效應(yīng)邏輯回歸（mixed effects logistic）模型分析肺癌數(shù)據(jù)》。
點(diǎn)擊標(biāo)題查閱往期內(nèi)容
多水平模型、分層線性模型HLM、混合效應(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多層（層次）線性模型模型
