原文鏈接:http:///?p=23717Logistic回歸,也稱為L(zhǎng)ogit模型,用于對(duì)二元結(jié)果變量進(jìn)行建模。在Logit模型中,結(jié)果的對(duì)數(shù)概率被建模為預(yù)測(cè)變量的線性組合。 例子例1. 假設(shè)我們對(duì)影響一個(gè)政治候選人是否贏得選舉的因素感興趣。結(jié)果(因)變量是二元的(0/1);贏或輸。我們感興趣的預(yù)測(cè)變量是花在競(jìng)選上的錢,花在競(jìng)選上的時(shí)間,以及候選人是否是現(xiàn)任者。 例2. 一個(gè)研究者對(duì)GRE(研究生入學(xué)考試成績(jī))、GPA(平均分)和本科院校的聲望等變量如何影響研究生院的錄取感興趣。因變量,錄取/不錄取,是一個(gè)二元變量。 數(shù)據(jù)的描述對(duì)于我們下面的數(shù)據(jù)分析,我們將在例2的基礎(chǔ)上展開(kāi)關(guān)于進(jìn)入研究生院的分析。我們生成了假設(shè)的數(shù)據(jù),這些數(shù)據(jù)可以在R中從我們的網(wǎng)站上獲得。請(qǐng)注意,R在指定文件位置時(shí)需要正斜杠(/)而不是反斜杠(),該文件在你的硬盤上。 ##查看數(shù)據(jù)的前幾行head(mydata) 這個(gè)數(shù)據(jù)集有一個(gè)二元因(結(jié)果,因果)變量,叫做錄取。有三個(gè)預(yù)測(cè)變量:gre、gpa和rank。我們將把gre和gpa這兩個(gè)變量視為連續(xù)變量。變量rank的值為1到4。排名為1的院校有最高的聲望,而排名為4的院校有最低的聲望。我們可以通過(guò)使用總結(jié)來(lái)獲得整個(gè)數(shù)據(jù)集的基本描述。為了得到標(biāo)準(zhǔn)差,我們使用sapply對(duì)數(shù)據(jù)集中的每個(gè)變量應(yīng)用sd函數(shù)。 你可能考慮的分析方法以下是你可能遇到過(guò)的一些分析方法的清單。所列的一些方法是相當(dāng)合理的,而其他的方法可能有局限性。
使用logit模型下面的代碼使用glm(廣義線性模型)函數(shù)估計(jì)一個(gè)邏輯回歸模型。首先,我們將等級(jí)轉(zhuǎn)換為一個(gè)因子變量,以表明等級(jí)應(yīng)被視為一個(gè)分類變量。 rank <- factor(rank)由于我們給我們的模型起了個(gè)名字(mylogit),R不會(huì)從我們的回歸中產(chǎn)生任何輸出。為了得到結(jié)果,我們使用summary命令。
點(diǎn)擊標(biāo)題查閱往期內(nèi)容 左右滑動(dòng)查看更多 我們可以使用confint函數(shù)來(lái)獲得系數(shù)估計(jì)值的置信區(qū)間。注意,對(duì)于logistic模型,置信區(qū)間是基于剖析的對(duì)數(shù)似然函數(shù)。我們也可以通過(guò)使用默認(rèn)的方法,只根據(jù)標(biāo)準(zhǔn)誤差來(lái)獲得CI。 我們可以用wald.test函數(shù)來(lái)檢驗(yàn)等級(jí)的整體效應(yīng)。系數(shù)表中系數(shù)的順序與模型中項(xiàng)的順序相同。這一點(diǎn)很重要,因?yàn)閣ald.test函數(shù)是按照系數(shù)在模型中的順序來(lái)參考的。我們使用wald.test函數(shù)。b提供了系數(shù),而Sigma提供了誤差項(xiàng)的方差協(xié)方差矩陣,最后Terms告訴R模型中哪些項(xiàng)要被測(cè)試,在本例中,4、5、6項(xiàng)是等級(jí)水平的三個(gè)項(xiàng)。 卡方檢驗(yàn)統(tǒng)計(jì)量為20.9,有三個(gè)自由度,P值為0.00011,表明等級(jí)的總體影響在統(tǒng)計(jì)上是顯著的。 我們還可以檢驗(yàn)關(guān)于不同等級(jí)的系數(shù)差異的其他假設(shè)。下面我們測(cè)試等級(jí)=2的系數(shù)是否等于等級(jí)=3的系數(shù)。下面的第一行代碼創(chuàng)建了一個(gè)向量l,定義了我們要執(zhí)行的測(cè)試。在這種情況下,我們要測(cè)試等級(jí)=2的項(xiàng)和等級(jí)=3的項(xiàng)(即模型中的第4和第5項(xiàng))的差異(減法)。為了對(duì)比這兩個(gè)項(xiàng),我們把其中一個(gè)項(xiàng)乘以1,另一個(gè)項(xiàng)乘以-1。下面的第二行代碼使用L=l來(lái)告訴R,我們希望以向量l為基礎(chǔ)進(jìn)行測(cè)試(而不是像上面那樣使用Terms選項(xiàng))。 wald.test(b , Sigma , L = l)1個(gè)自由度的卡方檢驗(yàn)統(tǒng)計(jì)量為5.5,P值為0.019,表明等級(jí)=2的系數(shù)和等級(jí)=3的系數(shù)之間的差異具有統(tǒng)計(jì)學(xué)意義。 你也可以對(duì)系數(shù)進(jìn)行指數(shù)化,并將其解釋為概率。為了得到指數(shù)化的系數(shù),你要告訴R你要進(jìn)行指數(shù)化(exp),你要指數(shù)化的對(duì)象叫做coefficients,它是mylogit的一部分(coef(mylogit))。我們可以使用同樣的邏輯,通過(guò)對(duì)之前的置信區(qū)間進(jìn)行指數(shù)化,得到概率及其置信區(qū)間。為了把這些都放在一個(gè)表中,我們用cbind把系數(shù)和置信區(qū)間按列綁定起來(lái)。 ## 概率比##概率和95%CI 現(xiàn)在我們可以說(shuō),gpa增加一個(gè)單位,被研究生院錄?。ㄅc未被錄?。┑膸茁示蜁?huì)增加2.23倍。請(qǐng)注意,截距的幾率一般不會(huì)被解釋。 你也可以使用預(yù)測(cè)概率來(lái)幫助你理解模型。預(yù)測(cè)概率可以針對(duì)分類和連續(xù)預(yù)測(cè)變量進(jìn)行計(jì)算。為了創(chuàng)建預(yù)測(cè)的概率,我們首先需要?jiǎng)?chuàng)建一個(gè)新的數(shù)據(jù)框架,其中包含我們希望自變量采取的數(shù)值,來(lái)創(chuàng)建我們的預(yù)測(cè)。 我們將首先計(jì)算每個(gè)等級(jí)值的預(yù)測(cè)錄取概率,保持gre和gpa的平均值。首先,我們創(chuàng)建并查看數(shù)據(jù)框架。 data.frame(mean(gre), mean(gpa), factor(1:4))## 查看數(shù)據(jù)框 這些對(duì)象的名稱必須與上述邏輯回歸中的變量相同(例如,在本例中,gre的平均值必須被命名為gre)。現(xiàn)在我們有了要用來(lái)計(jì)算預(yù)測(cè)概率的數(shù)據(jù)框,我們可以告訴R來(lái)創(chuàng)建預(yù)測(cè)概率。下面的第一行代碼非常緊湊,我們將把它拆開(kāi)來(lái)討論各個(gè)部分的作用。newdata1$rankP告訴R,我們要在數(shù)據(jù)集(數(shù)據(jù)框)newdata1中創(chuàng)建一個(gè)名為rankP的新變量,命令的其余部分告訴R,rankP的值應(yīng)該是使用predict( )函數(shù)進(jìn)行的預(yù)測(cè)。括號(hào)內(nèi)的選項(xiàng)告訴R,預(yù)測(cè)應(yīng)該基于mylogit分析,預(yù)測(cè)變量的值來(lái)自newdata1,預(yù)測(cè)的類型是預(yù)測(cè)的概率(type="response")。代碼的第二行列出數(shù)據(jù)框newdata1中的值。這是預(yù)測(cè)概率的表格。 predict(mylogit, newdata, type)在上面的輸出中,我們看到,在保持gre和gpa的平均值的情況下,來(lái)自最高聲望的本科院校(排名=1)的學(xué)生被研究生課程錄取的預(yù)測(cè)概率為0.52,而來(lái)自排名最低的院校(排名=4)的學(xué)生為0.18。我們可以做一些非常類似的事情,創(chuàng)建一個(gè)預(yù)測(cè)概率的表格,改變gre和排名的值。我們將繪制這些圖表,因此我們將在每個(gè)等級(jí)值(即1、2、3和4)上創(chuàng)建100個(gè)200至800的gre值。 gre = rep(seq(from = 200, to = 800, length.out = 100),4), mean(gpa), factor(rep(1:4, each = 100)) 生成預(yù)測(cè)概率的代碼(下面第一行)與之前的相同,只是我們還要提供標(biāo)準(zhǔn)誤差,這樣我們就可以繪制一個(gè)置信區(qū)間。我們?cè)阪溄訕?biāo)度上得到估計(jì)值,并將預(yù)測(cè)值和置信區(qū)間都反過(guò)來(lái)轉(zhuǎn)化為概率。 PredictedProbLL <- plogis(fit - (1.96 * se.fit)) UL <- plogis(fit + (1.96 * se.fit)) ##查看最終數(shù)據(jù)集的前幾行 使用預(yù)測(cè)概率的圖表來(lái)理解和/或展示模型也是有幫助的。我們將使用ggplot2軟件包來(lái)繪制圖表。下面我們用預(yù)測(cè)的概率和95%的置信區(qū)間做一個(gè)圖。 ggplot( aes(x = gre, y = Predicted))我們也可能希望看到我們的模型擬合程度的方法。在比較相互比較的模型時(shí),這可能特別有用。summary(mylogit)產(chǎn)生的輸出包括擬合指數(shù)(顯示在系數(shù)下面),包括無(wú)效和偏差殘差以及AIC。衡量模型擬合度的一個(gè)指標(biāo)是整個(gè)模型的顯著性。這個(gè)測(cè)試問(wèn)的是有預(yù)測(cè)因子的模型是否比只有截距的模型(即空模型)明顯更適合。檢驗(yàn)統(tǒng)計(jì)量是帶有預(yù)測(cè)因子的模型與無(wú)效模型的殘差。檢驗(yàn)統(tǒng)計(jì)量是分布式的卡方,自由度等于當(dāng)前模型和無(wú)效模型之間的自由度差異(即模型中預(yù)測(cè)變量的數(shù)量)。為了找到兩個(gè)模型的偏差差異(即檢驗(yàn)統(tǒng)計(jì)量),我們可以使用以下命令。 with(mylogit, null.deviance - deviance)## \[1\] 41.5 兩個(gè)模型之間差異的自由度等于模式中預(yù)測(cè)變量的數(shù)量,可以用以下方法得到。 with(mylogit, df.null - df.residual)## \[1\] 5 最后,可以用P值得到。 ## \[1\] 7.58e-085個(gè)自由度的卡方為41.46,相關(guān)的P值小于0.001,這告訴我們,我們的模型作為一個(gè)整體的擬合度明顯好于一個(gè)空模型。這有時(shí)被稱為似然比檢驗(yàn)(偏差殘差為-2*對(duì)數(shù)似然)。要查看模型的對(duì)數(shù)似然,我們可以輸入。 logLik(mylogit)## 'log Lik.' -229 (df=6) 需要考慮的事項(xiàng)
參考文獻(xiàn)Hosmer, D. & Lemeshow, S. (2000). Applied Logistic Regression (Second Edition). New York: John Wiley & Sons, Inc. Long, J. Scott (1997). Regression Models for Categorical and Limited Dependent Variables. Thousand Oaks, CA: Sage Publications. |
|
來(lái)自: 拓端數(shù)據(jù) > 《待分類》