《R語言實戰》第四部分第十三章-廣義線性模型學習筆記

通過第八章和第九章的學習，我們可以採用線性模型預測一系列連續型和/或類別型預測變數的響應變數，但現實情況下，對於以上預測的前提條件——因變數為正態分布——常常不太合理。比如以下兩種情況：

• 結果變數可能是類別型的；

• 結構變數可能是計數型的。

本章介紹的廣義線性模型包含了非正態因變數的分析，分別為Logistic回歸（因變數為類別型）和泊松回歸（因變數為計數型），擴展了之前線性模型的框架。下面在簡單介紹廣義線性模型以及R語言中對於的glm()函數以後分別用Logistic回歸和泊松回歸解決以下兩個標準線性模型無法解決的問題：

• 預測婚姻出軌問題（二值型，出軌/未出軌）需要什麼樣的個人信息、人口統計信息和人際關係；

• 八周中所發生的癲癇次數與藥物治療有何關係。

1 廣義線性模型和glm()函數

廣義線性模型（Generalized linear model）與標準線性模型的區別主要在於響應變數Y無需滿足正態分布，只要服從指數分布族的一種分布即可，模型擬合的形式如下：

其中是條件均值的函數（稱為連接函數），連接函數以及響應變數Y所滿足的概率分布確定以後，廣義線性模型通過最大似然估計（最大似然估計）的迭代推導出各個參數。

1.1 glm()函數

R中主要由glm()函數擬合廣義線性模型，格式如下：

glm(formula, family=family(link=function), data=)n

下表列出了概率分布（family）和相應默認的連接函數（function）

以下舉例說明glm()函數對Logistic回歸和泊松回歸的擬合。

假設如下：

• 響應變數：Y

• 三個預測變數：X1，X2，X3

• 包含數據的數據框：mydata

Logistic回歸適用於二值響應變數。

假設Y服從二項分布，則線性模型的擬合形式為：

其中π=是Y的條件均值（條件均值），log(π/1-π)為連接函數，概率分布為二項分布，Logistic回歸模型擬合代碼如下：

glm(Y~X1+X2+X3, family=binomial(link="logit"), data=mydata)n

泊松回歸使用在給定時間內響應變數為事件發生數目的情形。

假設Y服從泊松分布，則線性模型的擬合形式為：

其中是Y的均值，log()為連接函數，概率分布為泊松分布，泊松回歸模型擬合代碼如下：

glm(Y~X1+X2+X3, family=position(link="log"), data=mydata)n

值得一提的是，如果令連接函數，並設定概率分布為正態（高斯）分布，那麼此時的線性回歸模型擬合代碼如下：

glm(Y~X1+X2+X3, family=gaussian(link="identity"), data=mydata)n

此時生成的結果與下列代碼的結果相同：

lm(Y~X1+X2+X3, data=mydata)n

聰明的讀者不難發現，此時的glm()函數所模擬的廣義線性模型與標準線性模型等價，也就是說，標準線性模型是廣義線性模型的一種特例。

綜上所述，這裡再簡單總結一下，廣義線性模型不直接擬合響應變數Y本身的條件均值，而不是擬合響應變數Y的條件均值的一個函數，並假設響應變數服從指數分布族中的某個分布，採用極大似然估計而非最小二乘法為推到依據，極大地擴展了線性模型的適用範圍。

1.2 連用的函數

glm()函數常用的函數與lm()函數類似。

PS：不好生意，中英文兩個版本當中表13-2都跨頁，所以只能截兩張圖上來了。

1.3 模擬擬合和回歸診斷

通常來說，可以採用第八章中描述的方法，不過也牢記如下建議：

1，繪製初始響應變數的預測值與殘差的圖形用於評價模型的適用性，代碼如下：

plot(predict(model, type="response"),n residuals(model, type="deviance"))n

其中，model為glm()函數返回的對象。

2，繪製各統計量（帽子值、學生化殘差值和Cook距離）的參考圖，找出異常大的值用於識別異常點的閾值，代碼如下：

plot(hatvalues(model))nplot(rstudent(model))nplot(cooks.distance(model))n

或者

library(car)ninfluencePlot(model)n

2 Logistic回歸

Logistic回歸在通過一系列連續型和/或類別型預測變數來預測二值型結果變數時，是一個非常有用的工具。本節我們將以AER包中的數據框Affairs為例來說明Logistic回歸的應用。在使用之前請確保已經安裝AER軟體包（install.packages("AER")）。

婚外情數據一共從601個參與者身上收集了9個變數，包括一年來婚外私通的頻率以及參與者性別、年齡、婚齡、是否有小孩、宗教信仰程度(5分制,1分表示反對,5分表示非常信仰)、學歷、職業(逆向編號的戈登7種分類),還有對婚姻的自我評分(5分制,1表示非常不幸福,5表示非常幸福)。

描述性的統計信息如下：

> data(Affairs,package = "AER")n> summary(Affairs)n affairs gender age yearsmarried children n Min. : 0.000 female:315 Min. :17.50 Min. : 0.125 no :171 n 1st Qu.: 0.000 male :286 1st Qu.:27.00 1st Qu.: 4.000 yes:430 n Median : 0.000 Median :32.00 Median : 7.000 n Mean : 1.456 Mean :32.49 Mean : 8.178 n 3rd Qu.: 0.000 3rd Qu.:37.00 3rd Qu.:15.000 n Max. :12.000 Max. :57.00 Max. :15.000 n religiousness education occupation rating n Min. :1.000 Min. : 9.00 Min. :1.000 Min. :1.000 n 1st Qu.:2.000 1st Qu.:14.00 1st Qu.:3.000 1st Qu.:3.000 n Median :3.000 Median :16.00 Median :5.000 Median :4.000 n Mean :3.116 Mean :16.17 Mean :4.195 Mean :3.932 n 3rd Qu.:4.000 3rd Qu.:18.00 3rd Qu.:6.000 3rd Qu.:5.000 n Max. :5.000 Max. :20.00 Max. :7.000 Max. :5.000 n> table(Affairs$affairs)nn 0 1 2 3 7 12 n451 34 17 19 42 38 n> n

插一句，最近啟用Ubuntu系統，發現在Ubuntu下安裝一個R包可真夠費勁的，下載安裝包，用G++編譯，檢查各種頭文件，然後再配置，從輸入install.packages命令到安裝完畢要顯示好幾屏幕的提示。

比起婚外情的次數，我們更加關注二值型結果，即有或者沒有，以下代碼將affairs轉化為二值型因子ynaffair，並將此變數作為Logistic回歸的結果變數：

> Affairs$ynaffair[Affairs$affairs > 0] <- 1n> Affairs$ynaffair[Affairs$affairs == 0] <- 0n> Affairs$ynaffair <- factor(Affairs$ynaffair,n+ levels = c(0,1),n+ labels = c("NO","Yes"))n> table(Affairs$ynaffair)nn NO Yes n451 150 n> n

二值型因子轉化。

> fit.full <- glm(ynaffair ~ gender + age + yearsmarried + children + religiousness + education + occupation + rating,n+ data = Affairs, family = binomial())n> summary(fit.full)nnCall:nglm(formula = ynaffair ~ gender + age + yearsmarried + children + n religiousness + education + occupation + rating, family = binomial(), n data = Affairs)nnDeviance Residuals: n Min 1Q Median 3Q Max n-1.5713 -0.7499 -0.5690 -0.2539 2.5191 nnCoefficients:n Estimate Std. Error z value Pr(>|z|) n(Intercept) 1.37726 0.88776 1.551 0.120807 ngendermale 0.28029 0.23909 1.172 0.241083 nage -0.04426 0.01825 -2.425 0.015301 * nyearsmarried 0.09477 0.03221 2.942 0.003262 ** nchildrenyes 0.39767 0.29151 1.364 0.172508 nreligiousness -0.32472 0.08975 -3.618 0.000297 ***neducation 0.02105 0.05051 0.417 0.676851 noccupation 0.03092 0.07178 0.431 0.666630 nrating -0.46845 0.09091 -5.153 2.56e-07 ***n---nSignif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1nn(Dispersion parameter for binomial family taken to be 1)nn Null deviance: 675.38 on 600 degrees of freedomnResidual deviance: 609.51 on 592 degrees of freedomnAIC: 627.51nnNumber of Fisher Scoring iterations: 4nn> n

Logistic擬合和主要結果展示。

從結果來看，性別、是否有孩子、學歷、和職業對模型的結果貢獻不大，因此減少變數數重新擬合，並檢查新模型是否擬合得更好：

> fit.reduced <- glm(ynaffair ~ age + yearsmarried + religiousness + rating,n+ data = Affairs, family = binomial())n> summary(fit.reduced)nnCall:nglm(formula = ynaffair ~ age + yearsmarried + religiousness + n rating, family = binomial(), data = Affairs)nnDeviance Residuals: n Min 1Q Median 3Q Max n-1.6278 -0.7550 -0.5701 -0.2624 2.3998 nnCoefficients:n Estimate Std. Error z value Pr(>|z|) n(Intercept) 1.93083 0.61032 3.164 0.001558 ** nage -0.03527 0.01736 -2.032 0.042127 * nyearsmarried 0.10062 0.02921 3.445 0.000571 ***nreligiousness -0.32902 0.08945 -3.678 0.000235 ***nrating -0.46136 0.08884 -5.193 2.06e-07 ***n---nSignif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1nn(Dispersion parameter for binomial family taken to be 1)nn Null deviance: 675.38 on 600 degrees of freedomnResidual deviance: 615.36 on 596 degrees of freedomnAIC: 625.36nnNumber of Fisher Scoring iterations: 4nn> n

新模型每個回歸係數都非常顯著，對於兩個嵌套模型可以使用anova()函數對它們進行比較，對於廣義線性回歸可以採用卡方檢驗。

> anova(fit.full,fit.reduced,test = "Chisq")nAnalysis of Deviance TablennModel 1: ynaffair ~ gender + age + yearsmarried + children + religiousness + n education + occupation + ratingnModel 2: ynaffair ~ age + yearsmarried + religiousness + ratingn Resid. Df Resid. Dev Df Deviance Pr(>Chi)n1 592 609.51 n2 596 615.36 -4 -5.8474 0.2108n> n

卡方值不顯著，表明四個預測變數的新模型與九個完整預測變數的模型擬合程度一樣好，因此可以簡化模型。

2.1 解釋模型參數

在Logistic回歸當中，響應變數是Y=1的對數優勢比（log），因此為了更好的解釋性，需要對結果進行指數化。

> coef(fit.reduced)n (Intercept) age yearsmarried religiousness rating n 1.93083017 -0.03527112 0.10062274 -0.32902386 -0.46136144 n> exp(coef(fit.reduced))n (Intercept) age yearsmarried religiousness rating n 6.8952321 0.9653437 1.1058594 0.7196258 0.6304248 n

可以看到婚齡增加一年，婚外情的優勢比將乘以1.106(保持年齡、宗教信仰和婚姻評定不變)；相反，年齡增加一歲，婚外情的的優勢比則乘以0.965。因此，隨著婚齡的增加和年齡、宗教信仰與婚姻評分的降低，婚外情優勢比將上升。類似的，考慮到預測變數不能為0，截距項在此處沒有什麼特定含義。

如有必要，可以使用confint()函數獲取係數的置信區間，也可以進行指數化。

> confint(fit.reduced)nWaiting for profiling to be done...n 2.5 % 97.5 %n(Intercept) 0.75404303 3.150622807nage -0.07006400 -0.001854759nyearsmarried 0.04388142 0.158562400nreligiousness -0.50637196 -0.155156981nrating -0.63741235 -0.288566411n> exp(confint(fit.reduced))nWaiting for profiling to be done...n 2.5 % 97.5 %n(Intercept) 2.1255764 23.3506030nage 0.9323342 0.9981470nyearsmarried 1.0448584 1.1718250nreligiousness 0.6026782 0.8562807nrating 0.5286586 0.7493370n> n

2.2 評價預測變數對結果概率的影響

通常，以概率的方式思考比使用優勢比更加的直觀。predict()函數可以觀察某個預測變數在各個水平時對結果概率的影響。我們需要創建一個包含感興趣預測變數值的虛擬數據集，然後用predict()函數作用與上述虛擬數據集，預測這些值的結果概率。從本質上講，就是對某一個預測變數值進行敏感性分析。下面以婚姻評分這一預測變數為例。

#創建虛擬數據集ntestdata <- data.frame(rating=c(1, 2, 3, 4, 5), age=mean(Affairs$age),nyearsmarried=mean(Affairs$yearsmarried),nreligiousness=mean(Affairs$religiousness))n#使用測試數據集預測相應的概率ntestdata$prob <- predict(fit.reduced, newdata=testdata, type="response")ntestdadan

結果如下：

> testdatan rating age yearsmarried religiousness probn1 1 32.48752 8.177696 3.116473 0.5302296n2 2 32.48752 8.177696 3.116473 0.4157377n3 3 32.48752 8.177696 3.116473 0.3096712n4 4 32.48752 8.177696 3.116473 0.2204547n5 5 32.48752 8.177696 3.116473 0.1513079n

從結果上看，當婚姻評分從1(很不幸福)變為5(非常幸福)時，婚外情概率從0.53降低到了0.15(假定年齡、婚齡和宗教信仰不變)。類似的方法我們可以探究每一個預測變數對結果概率的影響。

2.3 過度離勢

抽樣於二項分布（二項分布）的數據的期望方差是,其中n為觀測數，Pi為屬於Y=1組的概率。所謂過度離勢，就是觀測到的響應變數的方差大於期望的二項分布本身的方差。過度離勢會導致奇異的標準誤檢驗和不精確的顯著性檢驗。此時仍然可以使用glm()函數來擬合Logistic回歸，不過要用類二項分布來代替二項分布。

檢測是否過度離勢的一種方法是比較二項分布模型的殘差偏差與殘差自由度，如果比值=殘差偏差/殘差自由度比1大很多，則可以認為存在過度離勢。

> deviance(fit.reduced)/df.residual(fit.reduced)n[1] 1.03248n

婚外情的例子比值非常接近於1，表明沒有過度離勢。

此外，通過分布取用不同的分布（二項分布和類二項分布）對模型擬合兩次，然後用卡方檢驗對過度離勢進行檢驗。

> fit <- glm(ynaffair ~ age + yearsmarried + religiousness + rating, family = binomial(),data = Affairs)n> fit.od <- glm(ynaffair ~ age + yearsmarried + religiousness + rating, family = quasibinomial(),data = Affairs)n> pchisq(summary(fit.od)$dispersion * fit$df.residual,n+ fit$df.residual,lower =F)n[1] 0.340122n> n

此時p值不顯著，這進一步增加了我們認為不存在過度離勢的信心。

2.3 擴展

R中還有不少擴展的Logistic回歸和變種：

• 穩健Logistic回歸 robust 包中的 glmRob() 函數可用來擬合穩健的廣義線性模型，包括穩健Logistic回歸。
• 多項分布回歸 若響應變數包含兩個以上的無序類別(比如,已婚/寡居/離婚)，便可使用 mlogit 包中的 mlogit() 函數擬合多項Logistic回歸。
• 序數Logistic回歸 若響應變數是一組有序的類別(比如,信用風險為差/良/好)，便可使用 rms 包中的 lrm() 函數擬合序數Logistic回歸。

如果我們將興趣從是否有婚外情轉移到過去一年中婚外情的次數，便可以直接泊松回歸對計數型數據進行分析。

3 泊松回歸

泊松回歸（泊松回歸）通常用於通過一系列連續型和/或類別型預測變數來預測計數型結果變數。以下將採用robust包中的Breslow數據闡述泊松回歸的應用。

Breslow數據的響應變數為 sumY (隨機化後八周內癲癇發病數)，預測變數為治療條件( Trt )、年齡( Age )和前八周內的基礎癲癇發病數（base）。包含基礎癲癇發病數和年齡的原因在於他們對響應變數有潛在影響。我們最終感興趣的是藥物治療是否能減少癲癇發病數。

首先，安裝robust包：

install.packages("robust")n

然後查看數據集的統計匯總信息：

> data(breslow.dat, package = "robust")n> names(breslow.dat)n [1] "ID" "Y1" "Y2" "Y3" "Y4" "Base" "Age" "Trt" "Ysum" n[10] "sumY" "Age10" "Base4"n> summary(breslow.dat[c(6,7,8,10)])n Base Age Trt sumY n Min. : 6.00 Min. :18.00 placebo :28 Min. : 0.00 n 1st Qu.: 12.00 1st Qu.:23.00 progabide:31 1st Qu.: 11.50 n Median : 22.00 Median :28.00 Median : 16.00 n Mean : 31.22 Mean :28.34 Mean : 33.05 n 3rd Qu.: 41.00 3rd Qu.:32.00 3rd Qu.: 36.00 n Max. :151.00 Max. :42.00 Max. :302.00 n> n

雖然數據集一共有12個變數，但是我們只關心之前描述的四個變數。對於響應變數，我們用如下代碼生成下列圖形：

> opar <- par(no.readonly = TRUE)n> par(mfow=c(1,2))nWarning message:nIn par(mfow = c(1, 2)) : "mfow" is not a graphical parametern> par(mfrow=c(1,2))n> attach(breslow.dat)n> hist(sumY, breaks = 20, xlab = "Seizure Count",n+ main = "Distribution of Seizures")n> boxplot(sumY~Trt, xlab="Treatment",main = "Group Comparisons")n> par(opar)n> n

上圖清楚地展示了因變數的偏倚特性和可能存在的離群點，而且在藥物治療下癲癇發病數似乎變小了，且方差也變小了。下面用泊松回歸進行擬合：

> fit <- glm(sumY ~ Base + Age + Trt, data = breslow.dat, family = poisson())n> summary(fit)nnCall:nglm(formula = sumY ~ Base + Age + Trt, family = poisson(), data = breslow.dat)nnDeviance Residuals: n Min 1Q Median 3Q Max n-6.0569 -2.0433 -0.9397 0.7929 11.0061 nnCoefficients:n Estimate Std. Error z value Pr(>|z|) n(Intercept) 1.9488259 0.1356191 14.370 < 2e-16 ***nBase 0.0226517 0.0005093 44.476 < 2e-16 ***nAge 0.0227401 0.0040240 5.651 1.59e-08 ***nTrtprogabide -0.1527009 0.0478051 -3.194 0.0014 ** n---nSignif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1nn(Dispersion parameter for poisson family taken to be 1)nn Null deviance: 2122.73 on 58 degrees of freedomnResidual deviance: 559.44 on 55 degrees of freedomnAIC: 850.71nnNumber of Fisher Scoring iterations: 5nn> n

輸出結果中的Coefficients結果列出了偏差、回歸參數、標準誤差和參數為0的檢驗。顯然，所有的p值都小於0.05，結果都非常顯著。

3.1 解釋模型參數

Logistic回歸是對數優勢比，在泊松回歸中，因變數以條件均值的對數形式來建模，響應的初始模型參數為對數均值，同樣也可以進行指數化以探求因變數的初始尺度上解釋回歸係數。以下代碼給出了著兩者的參數：

> coef(fit)n (Intercept) Base Age Trtprogabide n 1.94882593 0.02265174 0.02274013 -0.15270095 n> exp(coef(fit))n (Intercept) Base Age Trtprogabide n 7.0204403 1.0229102 1.0230007 0.8583864 n> n

年齡的回歸參數為0.0227，表明保持其他預測變數不變,年齡增加一歲，癲癇發病數的對數均值將相應增加0.03。截距項為預測變數都為0時，發病數的對數均值，顯然年齡不可能為0，而且調查對象的基礎發病數也都不為0，因此截距項沒有實際意義。

從指數化的係數可以看出，保持其他變數不變,年齡增加一歲，期望的癲癇發病數將乘以1.023。這意味著年齡的增加與較高的癲癇發病數相關聯。更為重要的是,一單位 Trt 的變化(即從安慰劑到治療組),期望的癲癇發病數將乘以0.86，換句話說，保持基礎癲癇發病數和年齡不變，服藥組相對於安慰劑組癲癇發病數降低了20%。

需要注意的是，與Logistic回歸中的指數化參數相似,泊松模型中的指數化參數對響應變數的影響都是成倍增加的，而不是線性相加。

3.2 過度離勢

泊松分布（泊松分布）的期望和方差相等，當響應變數觀測的方差比依據泊松分布預測的方差大時，泊松回歸可能發生過度離勢，而且發生的概率很大。可能發生過度離勢的原因有如下幾個：

• 遺漏了某個重要的預測變數；

• 可能因為事件相關，在在泊松分布的觀測中,計數中每次事件都被認為是獨立發生的。

• 在縱向數據分析中，重複測量的數據由於內在群聚特性可導致過度離勢。

如果過度離勢發生了，在模型中你無法進行解釋，那麼有可能你會發現並不真實存在的效應。

判斷過度離勢是否發生的準則依然是殘差偏差與殘差自由度的比例，如果遠遠大於1則表明存在過度離勢。對於上述癲癇發病數數據，它的比例是：

> deviance(fit)/df.residual(fit)n[1] 10.1717n> n

很顯然，比例遠大於1，存在過度離勢。

qcc包提供了一個對泊松模型過度離勢的檢驗方法。

> install.packages("qcc")#安裝qcc包nInstalling package into 『/home/jijiwhywhy/R/x86_64-pc-linux-gnu-library/3.3』n(as 『lib』 is unspecified)ntrying URL https://mirrors.tuna.tsinghua.edu.cn/CRAN/src/contrib/qcc_2.6.tar.gznContent type application/octet-stream length 159062 bytes (155 KB)n==================================================ndownloaded 155 KBnn* installing *source* package 『qcc』 ...n** package 『qcc』 successfully unpacked and MD5 sums checkedn** Rn** datan** demon** instn** byte-compile and prepare package for lazy loadingn** helpn*** installing help indicesn** building package indicesn** testing if installed package can be loadedn* DONE (qcc)nnThe downloaded source packages are innt『/tmp/Rtmpxk2no3/downloaded_packages』n> library(qcc)#導入qcc包nPackage qcc, version 2.6nType citation("qcc") for citing this R package in publications.n> qcc.overdispersion.test(breslow.dat$sumY, type = "poisson")#檢驗n nOverdispersion test Obs.Var/Theor.Var Statistic p-valuen poisson data 62.87013 3646.468 0n> n

顯著性檢驗的p值小於0.05，進一步確認存在過度離勢。

同樣的，可以用類泊松分布代替泊松分布來嘗試解決這個問題。

> fit.od <- glm(sumY ~ Base + Age + Trt, data=breslow.dat,n+ family=quasipoisson())n> summary(fit.od)nnCall:nglm(formula = sumY ~ Base + Age + Trt, family = quasipoisson(), n data = breslow.dat)nnDeviance Residuals: n Min 1Q Median 3Q Max n-6.0569 -2.0433 -0.9397 0.7929 11.0061 nnCoefficients:n Estimate Std. Error t value Pr(>|t|) n(Intercept) 1.948826 0.465091 4.190 0.000102 ***nBase 0.022652 0.001747 12.969 < 2e-16 ***nAge 0.022740 0.013800 1.648 0.105085 nTrtprogabide -0.152701 0.163943 -0.931 0.355702 n---nSignif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1nn(Dispersion parameter for quasipoisson family taken to be 11.76075)nn Null deviance: 2122.73 on 58 degrees of freedomnResidual deviance: 559.44 on 55 degrees of freedomnAIC: NAnnNumber of Fisher Scoring iterations: 5nn> n

凡事都有兩面性，使用類泊松(quasi-Poisson)方法所得的參數估計與泊松方法相同，但標準誤變大了許多，並且標準誤越大將會導致 Trt (和 Age )的p值越大於0.05。當考慮過度離勢,並控制基礎癲癇數和年齡時，並沒有充足的證據表明藥物治療相對於使用安慰劑能顯著降低癲癇發病次數。

3.3 擴展

R中擴展的泊松模型變種主要有如下三種：

• 時間段變化的泊松回歸 需要引入一個記錄每個觀測的時間長度的變數，並將模型從 修改為： 同時，為了擬合新模型，需要使用glm()函數當中的offset選項。以癲癇數據為例：

  fit <- glm(sumY ~ Base + Age + Trt, data=breslow.dat,noffset= log(time), family=poisson)n

• 零膨脹的泊松回歸 零膨脹的泊松回歸主要解決0計數的數目比泊松模型預測的數目多的問題，以婚外情數據為例，它將同時擬合兩個模型：一個用來預測哪些人又會發生婚外情，另外一個用來預測排除了婚姻忠誠者後的調查對象會發生多少次婚外情。pscl包中的zeroinfl()函數可以做零膨脹泊松回歸。
• 穩健泊松回歸 robust包中的glmRob()函數可以擬合穩健廣義線性模型，包括穩健泊松回歸，主要應對存在離群點和強影響點的問題。

4 小結

本章主要介紹廣義線性模型用於分析非正態的響應變數，包括分類數據和離散的計數型數據。並探討了對應的解決上述問題的兩個模型：Logistic回歸模型和泊松回歸模型，具體包括擬合後參數的解釋以及過度離勢的判別和診斷，以及一些針對特殊情況的變種。

到目前為止，我們介紹的每種統計方法都是直接處理可觀測、可記錄的變數。下一章我們將了解一些處理潛變數的統計模型。

ps：清明假期第一天，早起跑了一個小時，不禁套用信樂團的一句話：我不抽煙，我不喝酒，我愛學習和跑步。

最後祝大家假期愉快，希望我寫的東西能夠對大家有用，歡迎留言討論交流。

推薦閱讀：