---
title: GLM
tags: 算法
---
前言
之前我們使用R進行線性分析時，我們得到了滿足正態(tài)分布下的線性擬合模型，但在許多實際情況中，假定滿足正態(tài)分布的假設(shè)并不合理
結(jié)果變量為類別型：包括二值變量（是/否、通過/失敗、存活/死亡）和多類別變量（優(yōu)/良/可/差）等。
結(jié)果變量為計數(shù)型的（一周交通事故的數(shù)目、每日酒水消耗的數(shù)量），這類變量都是非負(fù)的有限值，它們的均值和方差通常都是相關(guān)的（而正態(tài)分布變量間則是相互獨立）。
為此，我們在一般線性的基礎(chǔ)上進行深入學(xué)習(xí)，來研究包含了非正態(tài)因變量的廣義線性模型。
一、GLM函數(shù)概況
在R語言中，我們使用glm()函數(shù)構(gòu)建廣義線性模型，調(diào)用語法如下：
# family為擬合所屬的函數(shù)族
# function為對應(yīng)的連接函數(shù)
glm(formula,family=family(link = function),data=)
常用的分布族和連接函數(shù)見下表：
1. logistic回歸適用于二值響應(yīng)變量，連接函數(shù)為logit函數(shù)，概率分布為二項分布：
glm(Y ~ X1 + X2 + X3, family=binomial(link="logit"), data=data)
2. 泊松回歸適用于在給定時間內(nèi)響應(yīng)變量為事件發(fā)生數(shù)目的情形，其假設(shè)Y服從泊松分布，連接函數(shù)為log(λ)，概率分布為泊松分布：
glm(Y ~ X1 + X2 + X3, family=poisson(link="log"), data=data)
除基本公式外，廣義線性模型還會用到其他輔助函數(shù)，包括summary()，coef()等，這些函數(shù)的使用和線性模型中l(wèi)m()的配套輔助函數(shù)一致，具體可參閱：R Language Learning：線性回歸（五）
總之，廣義線性模型通過擬合響應(yīng)變量的條件均值的一個函數(shù)（不是響應(yīng)變量的條件均值），并假設(shè)響應(yīng)變量服從指數(shù)分布族中的某個分布（不限于正態(tài)分布），從而極大地擴展了標(biāo)準(zhǔn)線性模型。模型參數(shù)估計的推導(dǎo)依據(jù)是極大似然估計，而非最小二乘法。
二、Logistics回歸分析
當(dāng)因變量為二分類或多分類時，Logistics回歸是非常重要的模型。由于Logistics回歸對資料的正態(tài)性和方差齊性不做要求、對自變量類型也不做要求，使得Logistic回歸模型在眾多領(lǐng)域被廣泛使用。
例如，想探討胃癌發(fā)生的危險因素，可以選擇兩組人群，一組是胃癌組，一組是非胃癌組，兩組人群肯定有不同的體征和生活方式等。這里的因變量為是否胃癌，即“是”或“否”，為兩分類變量，自變量就可以包括很多了，例如年齡、性別、飲食習(xí)慣等。自變量既可以是連續(xù)的，也可以是分類的。通過logistic回歸分析，就可以大致了解到底哪些因素是胃癌的危險因素。
估計的Logistics回歸方程為：
[公式] or [公式]
Logistics回歸模型的表達(dá)形式為：[公式]
模型解讀：[公式]為暴露于某種狀態(tài)下的結(jié)局概率，[公式]是一種變量變換方式，表示對[公式]進行[公式]變換，[公式]為偏回歸系數(shù)，表示在其他自變量不變的條件下，每變化一個單位[公式]的估計值。
Logistics回歸是通過最大似然估計求解常數(shù)項和偏回歸系數(shù)，基本思想時當(dāng)從總體中隨機抽取n個樣本后，最合理的參數(shù)估計量應(yīng)該使得這n個樣本觀測值的概率最大。最大似然法的基本思想是先建立似然函數(shù)與對數(shù)似然函數(shù)，再通過使對數(shù)似然函數(shù)最大求解相應(yīng)的參數(shù)值，所得到的估計值稱為參數(shù)的最大似然估計值。
1. 數(shù)據(jù)準(zhǔn)備（類別型變量進行0/1量化）
首先，我們選用AER包中的Affairs數(shù)據(jù)集來構(gòu)建Logistics回歸模型，這個數(shù)據(jù)集記錄了一組婚外情數(shù)據(jù)，其中包括參與者性別、年齡、婚齡、是否有小孩、宗教信仰程度（5分制，1表示反對，5表示非常信仰）、學(xué)歷、職業(yè)和婚姻的自我評分（5分制，1表示非常不幸福，5表示非常幸福）。
在使用數(shù)據(jù)集之前，載入AER包
> library(AER)
> data(Affairs,package="AER")
對于這個數(shù)據(jù)集，我們關(guān)注是否出軌，即這是一個二值型結(jié)果（出軌過/從未出軌）。因此，我們接下來將'affaris'特征轉(zhuǎn)化為二值型因子'ynaffair'，該二值型因子即可以作為Logistic回歸的結(jié)果變量。
> Affairs$ynaffair[Affairs$affairs > 0] <- 1
> Affairs$ynaffair[Affairs$affairs== 0] <- 0
> Affairs$ynaffair <-factor(Affairs$ynaffair,levels=c(0,1),labels=c("No","Yes"))
2. Logistics模型構(gòu)建：
> myfit <- glm(ynaffair ~ gender + age + yearsmarried + children + religiousness + education + occupation + rating, data=Affairs, family=binomial())
> summary(myfit)
Coefficients:
? ? ? ? ? ? ? Estimate Std. Error z value Pr(>|z|)? ?
(Intercept)? ? 1.37726? ? 0.88776? 1.551 0.120807? ?
gendermale? ? 0.28029? ? 0.23909? 1.172 0.241083? ?
age? ? ? ? ? -0.04426? ? 0.01825? -2.425 0.015301 *?
yearsmarried? 0.09477? ? 0.03221? 2.942 0.003262 **
childrenyes? ? 0.39767? ? 0.29151? 1.364 0.172508? ?
religiousness -0.32472? ? 0.08975? -3.618 0.000297 ***
education? ? ? 0.02105? ? 0.05051? 0.417 0.676851? ?
occupation? ? 0.03092? ? 0.07178? 0.431 0.666630? ?
rating? ? ? ? -0.46845? ? 0.09091? -5.153 2.56e-07 ***
-----------------------------------------------
Signif. codes:? 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
從回歸系數(shù)的p值（最后一欄）可以看到，性別、是否有孩子、學(xué)歷和職業(yè)對方程的貢獻(xiàn)都不顯著（p值較大）。去除這些變量重新擬合模型，檢驗新模型是否擬合得更好。
> myfit_reduced <- glm(ynaffair ~ age + yearsmarried + religiousness + rating, data=Affairs, family=binomial())
> summary(myfit_reduced)
Coefficients:
? ? ? ? ? ? ? Estimate Std. Error z value Pr(>|z|)? ?
(Intercept)? ? 1.93083? ? 0.61032? 3.164 0.001558 **
age? ? ? ? ? -0.03527? ? 0.01736? -2.032 0.042127 *?
yearsmarried? 0.10062? ? 0.02921? 3.445 0.000571 ***
religiousness -0.32902? ? 0.08945? -3.678 0.000235 ***
rating? ? ? ? -0.46136? ? 0.08884? -5.193 2.06e-07 ***
從最后一列可以看到，新模型的每個回歸系數(shù)都非常顯著（p<0.05）
3. 多模型比較
由于這個兩模型是嵌套關(guān)系，所以我們可以使用anova()函數(shù)對它們進行比較。而對于廣義線性回歸，可用卡方檢驗進行判斷。
> anova(myfit,myfit_reduced,test="Chisq")
Analysis of Deviance Table
Model 1: ynaffair ~ gender + age + yearsmarried + children + religiousness +
? ? education + occupation + rating
Model 2: ynaffair ~ age + yearsmarried + religiousness + rating
? Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1? ? ? 592? ? 609.51? ? ? ? ? ? ? ? ? ?
2? ? ? 596? ? 615.36 -4? -5.8474? 0.2108
檢驗結(jié)果的卡方值不顯著（p=0.21），這表明，四個預(yù)測變量的新模型和九個預(yù)測變量的舊模型模擬程度其實差不多。一般情況下我們更傾向于選擇四個預(yù)測變量的新模型。
4. 評價預(yù)測變量對結(jié)果概率的影響（predict函數(shù)）
既然已經(jīng)有了模型，那么給定一個用戶行為，預(yù)測其出軌的概率。可以使用predict()函數(shù)進行預(yù)測，以下代碼生成了一些虛擬數(shù)據(jù)并進行了預(yù)測，主要評測婚姻評分對出軌概率的影響。
# 構(gòu)建dataframe
> testdata <- data.frame(rating=c(1, 2, 3, 4, 5), age=mean(Affairs$age), yearsmarried=mean(Affairs$yearsmarried), religiousness=mean(Affairs$religiousness))
# 添加一列為概率預(yù)測值
> testdata$prob <- predict(myfit_reduced, newdata=testdata, type="response")
> testdata
? rating? ? ? age yearsmarried religiousness? ? ? prob
1? ? ? 1 32.48752? ? 8.177696? ? ? 3.116473 0.5302296
2? ? ? 2 32.48752? ? 8.177696? ? ? 3.116473 0.4157377
3? ? ? 3 32.48752? ? 8.177696? ? ? 3.116473 0.3096712
4? ? ? 4 32.48752? ? 8.177696? ? ? 3.116473 0.2204547
5? ? ? 5 32.48752? ? 8.177696? ? ? 3.116473 0.1513079
從結(jié)果可以看到，當(dāng)婚姻評分從1分（很不幸福）增加到5分（非常幸福）時，出軌概率從0.53降低至0.15。
> testdata <- data.frame(rating=mean(Affairs$rating), age=seq(17, 57, 10), yearsmarried=mean(Affairs$yearsmarried), religiousness=mean(Affairs$religiousness))
> testdata$prob <- predict(myfit_reduced, newdata=testdata, type="response")
> testdata
? rating age yearsmarried religiousness? ? ? prob
1 3.93178? 17? ? 8.177696? ? ? 3.116473 0.3350834
2 3.93178? 27? ? 8.177696? ? ? 3.116473 0.2615373
3 3.93178? 37? ? 8.177696? ? ? 3.116473 0.1992953
4 3.93178? 47? ? 8.177696? ? ? 3.116473 0.1488796
5 3.93178? 57? ? 8.177696? ? ? 3.116473 0.1094738
在年齡的影響方面，可以看出，當(dāng)其他變量不變、年齡從17增加到57時，出軌的概率從0.34降低至0.11。
5. 過度離勢
過度離勢指的是觀測到的響應(yīng)變量方差大于期望的二項分布方差，過度離勢會導(dǎo)致奇異的標(biāo)準(zhǔn)誤檢驗和不精確的顯著性檢驗。
當(dāng)出現(xiàn)過度離勢時，仍可使用glm()函數(shù)擬合Logistic回歸，但此時需要將二項分布改為類二項分布。
檢測過度離勢的一種方法是比較二項分布模型的殘差偏差與殘差自由度（在四變量模型中分別是615.36和596），如果兩者的比值比1大很多，則可以認(rèn)為存在過度離勢。
6. Logistic回歸模型的進階學(xué)習(xí)
擴展的Logistic回歸和變種如下所示：
穩(wěn)健Logistic回歸：當(dāng)擬合Logistic回歸模型數(shù)據(jù)出現(xiàn)離群點和強影響點時，robust包中的glmRob()函數(shù)可用來擬合穩(wěn)健的廣義線性模型
多項分布回歸：當(dāng)響應(yīng)變量包含兩個以上的無序類別（例如已婚、寡居、離婚）時，可使用mlogit包中的mlogit()函數(shù)擬合多項Logistic回歸
序數(shù)Logistic回歸：當(dāng)響應(yīng)變量是一組有序的類別（例如信用為好/良/差）時，可使用rms包中的lrm()函數(shù)擬合序數(shù)Logistic回歸。
三、泊松回歸
泊松回歸的自變量是連續(xù)性或類別型變量，因變量是計數(shù)型的變量。泊松回歸因變量通常局限在一個固定長度時間段內(nèi)進行測量（如過去一年交通事故數(shù)），整個觀測集中時間長度都是不變的。Poisson回歸主要有兩個假設(shè)，首先，具有相同特征和同時的不同對象的人時風(fēng)險是同質(zhì)的，其次，當(dāng)樣本量越來越大時，頻數(shù)的均數(shù)趨近于方差。
調(diào)用模型的公式如下：
> myfit <- glm(y~x1+x2+...+xn,data=,family = poisson)
1. 數(shù)據(jù)準(zhǔn)備
robust包中Breslow癲癇數(shù)據(jù)記錄了治療初期八周內(nèi)，抗癲癇藥物對癲癇發(fā)病數(shù)的影響。響應(yīng)變量為sumY（隨機化后八周內(nèi)癲癇發(fā)病數(shù)），預(yù)測變量為治療條件（Trt）、年齡（Age）和前八周內(nèi)的基礎(chǔ)癲癇發(fā)病數(shù)（Base），在這個數(shù)據(jù)集中，我們感興趣的是藥物治療能否減少癲癇發(fā)病數(shù)。
> data(breslow.dat, package="robust")
接下來我們考察響應(yīng)變量，基礎(chǔ)和隨機化后的癲癇發(fā)病數(shù)都有很高的偏度，從圖中可以清楚地看到因變量的偏倚特性以及可能的離群點。
> library(ggplot2)
> p1 <- ggplot(breslow.dat,aes(sumY)) + geom_histogram(color='green')
> p2 <- ggplot(breslow.dat,aes(Trt,sumY)) + geom_bar()
# 組合圖形
> library(gridExtra)
> grid.arrange(p1,p2,nrow=1)
從圖中可以清楚的看到因變量的偏移特性及可能的離群點。藥物治療下癲癇的發(fā)病數(shù)似乎變小，且方差也變小了。
2. 構(gòu)建模型
> myfit <- glm(sumY ~ Base + Age + Trt, data=breslow.dat, family=poisson())
> summary(myfit)
3. 過度離勢
在泊松分布中同樣可能存在過度離勢，泊松分布的方差和均值相等，當(dāng)響應(yīng)變量觀測的方差比依據(jù)泊松分布預(yù)測的方差大時，可能存在過度離勢。
> summary(myfit)
(Dispersion parameter for poisson family taken to be 1)
? ? Null deviance: 2122.73? on 58? degrees of freedom
Residual deviance:? 559.44? on 55? degrees of freedom
AIC: 850.71
Number of Fisher Scoring iterations: 5
與Logistic回歸類似，但此例中殘差偏差（559.44）與殘差自由度（55）的比例為10.17，這一值遠(yuǎn)大于1，故存在過度離勢。在這種情況下，可以考慮使用類泊松回歸替代泊松回歸（family=quasipoisson()）。
除此之外，我們還可以使用qcc包來檢驗泊松模型是否存在過度離勢。
> library(qcc)
> qcc.overdispersion.test(breslow.dat$sumY, type="poisson")
4. 泊松回歸模型的進階學(xué)習(xí)：
時間段變化的泊松回歸：如果每個觀測的時間段長短不同，則因變量需要從計數(shù)更改為比率，即發(fā)生次數(shù)除以觀測時間，使用glm()的offset選項即可，offset=log(time)，其中time是每名病人的觀測時間
零膨脹的泊松回歸：在一個數(shù)據(jù)集中，0計數(shù)的數(shù)目時常比用泊松模型預(yù)測得多。在Affairs數(shù)據(jù)集中，很有可能有一群從未出軌的對象，他們稱為數(shù)據(jù)集的結(jié)構(gòu)零值?？梢杂昧闩蛎浀牟此苫貧w來分析這種情況，它將同時擬合兩個模型，一個用來預(yù)測哪些人又會出軌，另一個用來預(yù)測排除了結(jié)構(gòu)零值之后的調(diào)查對象會出軌多少次。可以將其看作是Logistic回歸（預(yù)測結(jié)構(gòu)零值）和泊松回歸（預(yù)測無結(jié)構(gòu)零值觀測的計數(shù)）的組合，pscl包的zeroinfl()函數(shù)可做零膨脹泊松回歸
穩(wěn)健泊松回歸：robust包中的glmRob()函數(shù)可以擬合穩(wěn)健廣義線性模型，包括穩(wěn)健泊松回歸。
References：
R語言實戰(zhàn)（第2版） (豆瓣)
R學(xué)習(xí)筆記－11 廣義線型模型
Generalized Linear Models
RPubs - 廣義線性模型
What does a generalized linear model do?
