第四部分 高級(jí)方法
第13章 廣義線性模型
廣義線性模型適用條件
假設(shè)因變量為正態(tài)分布不成立的情況，如：
?結(jié)果變量可能是類別型的二值變量（比如：是/否、通過/失敗、活著/死亡）和多分類變量（比如差/良好/優(yōu)秀）都顯然不是正態(tài)分布。
?結(jié)果變量可能是計(jì)數(shù)型的：（比如，一周交通事故的數(shù)目，每日酒水消耗的數(shù)量）。這類變量都是非負(fù)的有限值，而且它們的均值和方差通常都是相關(guān)的（正態(tài)分布變量間不是如此，而是相互獨(dú)立）。
廣義線性模型擴(kuò)展了線性模型的框架，它包含了非正態(tài)因變量的分析。
常見廣義線性模型
1、標(biāo)準(zhǔn)線性模型：是廣義線性模型的一個(gè)特例，因變量服從正態(tài)分布；
2、Logistic回歸：因變量服從二項(xiàng)分布，相應(yīng)變量是二值型（0和1）；
3、泊松回歸：因變量服從泊松分布，相應(yīng)變量是計(jì)數(shù)型。
廣義線性模型定義：
廣義線性模型(generalize linear model，GLM)是線性模型的擴(kuò)展，通過聯(lián)結(jié)函數(shù)建立響應(yīng)變量的數(shù)學(xué)期望值與線性組合的預(yù)測變量之間的關(guān)系。其特點(diǎn)是不強(qiáng)行改變數(shù)據(jù)的自然度量，數(shù)據(jù)可以具有非線性和非恒定方差結(jié)構(gòu)。是線性模型在研究響應(yīng)值的非正態(tài)分布以及非線性模型簡潔直接的線性轉(zhuǎn)化時(shí)的一種發(fā)展。
線性回歸和邏輯回歸都是廣義線性模型的一種特殊形式。
13.1 廣義線性模型和glm()函數(shù)
廣義線性模型擬合的形式為：
其中是條件均值的函數(shù)（稱為連接函數(shù)）。另外，你可放松Y為正態(tài)分布的假設(shè)，改為Y服從指數(shù)分布族中的一種分布即可。設(shè)定好連接函數(shù)和概率分布后，便可以通過最大似然估計(jì)的多次迭代推導(dǎo)出各參數(shù)值。
13.1.1 glm()函數(shù)
R中可通過glm函數(shù)（還可用其他專門的函數(shù)）擬合廣義線性模型。
glm(formula, family = gaussian, data, weights, subset,
na.action, start = NULL, etastart, mustart, offset,
control = list(...), model = TRUE, method = "glm.fit",
x = FALSE, y = TRUE, singular.ok = TRUE, contrasts = NULL, ...)
表13-1列出了概率分布（family）和相應(yīng)默認(rèn)的連接函數(shù)（function）。
假設(shè)你有一個(gè)響應(yīng)變量（Y）、三個(gè)預(yù)測變量（X1、X2、X3）和一個(gè)包含數(shù)據(jù)的數(shù)據(jù)框（mydata）。
Logistic回歸適用于二值響應(yīng)變量（0,1）。模型假設(shè)Y服從二項(xiàng)分布，線性模型的擬合形式為：
其中是Y的條件均值（即給定一系列X的值時(shí)Y = 1的概率），
為Y = 1時(shí)的優(yōu)勢比，
為對數(shù)優(yōu)勢比，或logit。本例中，
為連接函數(shù)，概率分布為二項(xiàng)分布，可用如下代碼擬合Logistic回歸模型：
glm(Y~x1+x2+x3,family=binomial(link="logit"),data=mydata)
泊松回歸適用于在給定時(shí)間內(nèi)響應(yīng)變量為事件發(fā)生數(shù)目的情形。它假設(shè)Y服從泊松分布，線性模型的擬合形式為：
其中λ是Y的均值（也等于方差）。此時(shí)，連接函數(shù)為log(λ)，概率分布為泊松分布，可用如下代碼擬合泊松回歸模型：
glm(Y~x1+x2+x3,family=poisson(link="log"),data=mydata)
值得注意的是，標(biāo)準(zhǔn)線性模型也是廣義線性模型的一個(gè)特例。如果令連接函數(shù)g(μY) =μY或恒等函數(shù)，并設(shè)定概率分布為正態(tài)（高斯）分布，那么：
glm(Y~x1+x2+x3,family=gaussian(link="identity"),data=mydata)
生成的結(jié)果與下列代碼的結(jié)果相同：
lm(Y~x1+x2+x3,data=mydata)
總之，廣義線性模型通過擬合響應(yīng)變量的條件均值的一個(gè)函數(shù)（不是響應(yīng)變量的條件均值），假設(shè)響應(yīng)變量服從指數(shù)分布族中的某個(gè)分布（并不僅限于正態(tài)分布），極大地?cái)U(kuò)展了標(biāo)準(zhǔn)線性模型。模型參數(shù)估計(jì)的推導(dǎo)依據(jù)的是極大似然估計(jì)，而非最小二乘法。
13.1.2 連用的函數(shù)
與分析標(biāo)準(zhǔn)線性模型時(shí)lm()連用的許多函數(shù)在glm()中都有對應(yīng)的形式，其中一些常見的函數(shù)見表13-2。
13.1.3 模型擬合和回歸診斷
當(dāng)評(píng)價(jià)模型的適用性時(shí)，你可以繪制初始響應(yīng)變量的預(yù)測值與殘差的圖形。例如，如下代碼可繪制一個(gè)常見的診斷圖：
plot(predict(model,type="response"),residuals(model,type="deviance"))
其中，model為glm()函數(shù)返回的對象。
R將列出帽子值（hat value）、學(xué)生化殘差值和Cook距離統(tǒng)計(jì)量的近似值。
car包繪制診斷圖
library(car)
influencePlot(model)
它可以創(chuàng)建一個(gè)綜合性的診斷圖。在后面的圖形中，橫軸代表杠桿值，縱軸代表學(xué)生化殘差值，而繪制的符號(hào)大小與Cook距離大小成正比。
13.2 Logistic回歸
當(dāng)通過一系列連續(xù)型和/或類別型預(yù)測變量來預(yù)測二值型結(jié)果變量時(shí)，Logistic回歸是一個(gè)非常有用的工具。
Logistic回歸模型的適用條件
1 因變量為二分類的分類變量或某事件的發(fā)生率，并且是數(shù)值型變量。但是需要注意，重復(fù)計(jì)數(shù)現(xiàn)象指標(biāo)不適用于Logistic回歸。
2 殘差和因變量都要服從二項(xiàng)分布。二項(xiàng)分布對應(yīng)的是分類變量，所以不是正態(tài)分布，進(jìn)而不是用最小二乘法，而是最大似然法來解決方程估計(jì)和檢驗(yàn)問題。
3 自變量和Logistic概率是線性關(guān)系。
4 各觀測對象間相互獨(dú)立。
logistic回歸和線性回歸
邏輯回歸（Logistic Regression）與線性回歸（Linear Regression）都是一種廣義線性模型（generalized linear model）。
邏輯回歸假設(shè)因變量 y 服從伯努利分布，而線性回歸假設(shè)因變量 y 服從高斯分布。
Sigmoid函數(shù)
也稱為邏輯函數(shù)（Logistic function）
從上圖可以看到sigmoid函數(shù)是一個(gè)s形的曲線，它的取值在[0,1]之間，在遠(yuǎn)離0的地方函數(shù)的值會(huì)很快接近0或者1。
《R語言實(shí)戰(zhàn)》例子，AER包Affairs數(shù)據(jù)集為例。
該數(shù)據(jù)從601個(gè)參與者身上收集了9個(gè)變量，包括一年來婚外私通的頻率以及參與者性別、年齡、婚齡、是否有小孩、宗教信仰程度（5分制，1分表示反對，5分表示非常信仰）、學(xué)歷、職業(yè)（逆向編號(hào)的戈登7種分類），還有對婚姻的自我評(píng)分（5分制，1表示非常不幸福，5表示非常幸福）。
數(shù)據(jù)的描述統(tǒng)計(jì)
data(Affairs, package="AER") # 調(diào)用AER包數(shù)據(jù)集Affairs。
head(Affairs) # 顯示數(shù)據(jù)前6行。
## affairs gender age yearsmarried children religiousness education occupation
## 4 0 male 37 10.00 no 3 18 7
## 5 0 female 27 4.00 no 4 14 6
## 11 0 female 32 15.00 yes 1 12 1
## 16 0 male 57 15.00 yes 5 18 6
## 23 0 male 22 0.75 no 2 17 6
## 29 0 female 32 1.50 no 2 17 5
## rating
## 4 4
## 5 4
## 11 4
## 16 5
## 23 3
## 29 5
summary(Affairs) # 數(shù)據(jù)描述。
## affairs gender age yearsmarried children
## Min. : 0.000 female:315 Min. :17.50 Min. : 0.125 no :171
## 1st Qu.: 0.000 male :286 1st Qu.:27.00 1st Qu.: 4.000 yes:430
## Median : 0.000 Median :32.00 Median : 7.000
## Mean : 1.456 Mean :32.49 Mean : 8.178
## 3rd Qu.: 0.000 3rd Qu.:37.00 3rd Qu.:15.000
## Max. :12.000 Max. :57.00 Max. :15.000
## religiousness education occupation rating
## Min. :1.000 Min. : 9.00 Min. :1.000 Min. :1.000
## 1st Qu.:2.000 1st Qu.:14.00 1st Qu.:3.000 1st Qu.:3.000
## Median :3.000 Median :16.00 Median :5.000 Median :4.000
## Mean :3.116 Mean :16.17 Mean :4.195 Mean :3.932
## 3rd Qu.:4.000 3rd Qu.:18.00 3rd Qu.:6.000 3rd Qu.:5.000
## Max. :5.000 Max. :20.00 Max. :7.000 Max. :5.000
table(Affairs$affairs) # 數(shù)據(jù)集Affairs中affairs列頻數(shù)統(tǒng)計(jì)。
##
## 0 1 2 3 7 12
## 451 34 17 19 42 38
將affairs轉(zhuǎn)化為二值型因子ynaffair
Affairs$ynaffair[Affairs$affairs > 0] <- 1 # 添加新列ynaffairs，原affairs列中值大于0，賦值1。
Affairs$ynaffair[Affairs$affairs == 0] <- 0 # 添加新列ynaffairs，原affairs列中值等于0，賦值0。
Affairs$ynaffair <- factor(Affairs$ynaffair,
levels=c(0,1),
labels=c("No","Yes")) # 新列轉(zhuǎn)為因子。
table(Affairs$ynaffair) # 統(tǒng)計(jì)因子頻數(shù)。
##
## No Yes
## 451 150
factor函數(shù)
factor(x, levels = sort(unique(x), na.last = TRUE),labels = levels, exclude = NA, ordered = is.ordered(x), nmax = NA)
x是要?jiǎng)?chuàng)建因子的向量，
levels用來指定因子水平值（不指定時(shí)由向量x不同值求得）；
labels用來指定水平的名字（不指定時(shí)由用水平值的對應(yīng)字符串）；exclude表示從向量x中剔除的水平值；
ordered是一個(gè)邏輯型選項(xiàng)，用來指定因子的水平是否有次序。
nmax是水平的上限數(shù)量。
對二值型因子變量進(jìn)行l(wèi)ogistic回歸
fit.full <- glm(ynaffair ~ gender + age + yearsmarried + children +
religiousness + education + occupation +rating,
data=Affairs,family=binomial()) # Logistic回歸，響應(yīng)變量為ynaffair，預(yù)測變量為gender，age，yearsmarried，children，religiousness，education，occupation，rating。
summary(fit.full) # 返回回歸結(jié)果。
##
## Call:
## glm(formula = ynaffair ~ gender + age + yearsmarried + children +
## religiousness + education + occupation + rating, family = binomial(),
## data = Affairs)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5713 -0.7499 -0.5690 -0.2539 2.5191
##
## 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
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 675.38 on 600 degrees of freedom
## Residual deviance: 609.51 on 592 degrees of freedom
## AIC: 627.51
##
## Number of Fisher Scoring iterations: 4
結(jié)果解讀：從回歸系數(shù)的p值（最后一欄）可以看到，性別、是否有孩子、學(xué)歷和職業(yè)對方程的貢獻(xiàn)都不顯著（你無法拒絕參數(shù)為0的假設(shè)）。原假設(shè)為所有預(yù)測變量對方程均有貢獻(xiàn)且相同
去除不顯著因子，重新擬合
fit.reduced <- glm(ynaffair ~ age + yearsmarried + religiousness +
rating, data=Affairs, family=binomial()) # 去除不顯著預(yù)測變量，重新擬合模型。
summary(fit.reduced) # 返回結(jié)果。
##
## Call:
## glm(formula = ynaffair ~ age + yearsmarried + religiousness +
## rating, family = binomial(), data = Affairs)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6278 -0.7550 -0.5701 -0.2624 2.3998
##
## 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 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 675.38 on 600 degrees of freedom
## Residual deviance: 615.36 on 596 degrees of freedom
## AIC: 625.36
##
## Number of Fisher Scoring iterations: 4
兩回歸模型的比較，對于廣義線性回歸，可用卡方檢驗(yàn)。
anova(fit.reduced, fit.full, test="Chisq") # 兩個(gè)擬合模型的對比。
## Analysis of Deviance Table
##
## Model 1: ynaffair ~ age + yearsmarried + religiousness + rating
## Model 2: ynaffair ~ gender + age + yearsmarried + children + religiousness +
## education + occupation + rating
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 596 615.36
## 2 592 609.51 4 5.8474 0.2108
卡方值不顯著（p=0.21），表明四個(gè)預(yù)測變量的新模型與九個(gè)完整預(yù)測變量的模型擬合程度一樣好。
13.2.1 解釋模型參數(shù)
在Logistic回歸中，響應(yīng)變量是Y=1的對數(shù)優(yōu)勢比（log）?；貧w系數(shù)含義是當(dāng)其他預(yù)測變量不變時(shí)，一單位預(yù)測變量的變化可引起的響應(yīng)變量對數(shù)優(yōu)勢比的變化。
對于二值型Logistic回歸，某預(yù)測變量n單位的變化引起的較高值上優(yōu)勢比的變化為exp(βj)^n，它反映的信息可能更為重要。
優(yōu)勢比
優(yōu)勢比(odds ratio；OR)是另外一種描述概率的方式。優(yōu)勢比將會(huì)告訴我們某種推測的概率比其反向推測的概率大多少。換句話說，優(yōu)勢比是指某種推測為真的概率與某種推測為假的概率的比值。比如下雨的概率為0.25，不下雨的概率為0.75。0.25與0.75的比值可以約分為1比3。因此，我們可以說今天將會(huì)下雨的優(yōu)勢比為1:3(或者今天不會(huì)下雨的概率比為3:1)。
coef(fit.reduced) # 查看fit.reduced模型的回歸系數(shù)。
## (Intercept) age yearsmarried religiousness rating
## 1.93083017 -0.03527112 0.10062274 -0.32902386 -0.46136144
exp(coef(fit.reduced)) # 指數(shù)化結(jié)果。
## (Intercept) age yearsmarried religiousness rating
## 6.8952321 0.9653437 1.1058594 0.7196258 0.6304248
exp(confint(fit.reduced)) # 在優(yōu)勢比尺度上得到系數(shù)95%的置信區(qū)間。
## 2.5 % 97.5 %
## (Intercept) 2.1255764 23.3506030
## age 0.9323342 0.9981470
## yearsmarried 1.0448584 1.1718250
## religiousness 0.6026782 0.8562807
## rating 0.5286586 0.7493370
13.2.2 評(píng)價(jià)預(yù)測變量對結(jié)果概率的影響
使用predict()函數(shù)，你可觀察某個(gè)預(yù)測變量在各個(gè)水平時(shí)對結(jié)果概率的影響。
創(chuàng)建一個(gè)虛擬數(shù)據(jù)集，設(shè)定年齡、婚齡和宗教信仰為它們的均值，婚姻評(píng)分的范圍為1～5。
testdata <- data.frame(rating = c(1, 2, 3, 4, 5),
age = mean(Affairs$age),
yearsmarried = mean(Affairs$yearsmarried),
religiousness = mean(Affairs$religiousness)) # 虛擬數(shù)據(jù)集testdata的創(chuàng)建。
predict函數(shù)預(yù)測數(shù)據(jù)集testdata中rating變量的概率
testdata$prob <- predict(fit.reduced, newdata=testdata, type="response") # 數(shù)據(jù)集testdata中rating變量的概率預(yù)測。
testdata # 顯示數(shù)據(jù)集。
## 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
predict函數(shù)預(yù)測數(shù)據(jù)集testdata中age變量的概率
testdata <- data.frame(rating = mean(Affairs$rating),
age = seq(17, 57, 10),
yearsmarried = mean(Affairs$yearsmarried),
religiousness = mean(Affairs$religiousness)) # age變量定義。
testdata$prob <- predict(fit.reduced, newdata=testdata, type="response") # 數(shù)據(jù)集testdata中age變量的概率預(yù)測。
testdata # 顯示結(jié)果。
## 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
13.2.3 過度離勢
抽樣于二項(xiàng)分布的數(shù)據(jù)的期望方差是σ2 = nπ(1-π)，n為觀測數(shù)，π為屬于Y = 1組的概率。
所謂過度離勢，即觀測到的響應(yīng)變量的方差大于期望的二項(xiàng)分布的方差。過度離勢會(huì)導(dǎo)致奇異的標(biāo)準(zhǔn)誤檢驗(yàn)和不精確的顯著性檢驗(yàn)。
當(dāng)出現(xiàn)過度離勢時(shí)，仍可使用glm()函數(shù)擬合Logistic回歸，但此時(shí)需要將二項(xiàng)分布改為類二項(xiàng)分布（quasibinomial distribution）
檢測過度離勢的一種方法是比較二項(xiàng)分布模型的殘差偏差與殘差自由度。
當(dāng)比值比1大很多，你便可認(rèn)為存在過度離勢。
你還可以對過度離勢進(jìn)行檢驗(yàn)。為此，你需要擬合模型兩次，第一次使用 family =
"binomial"，第二次使用family = "quasibinomial"。假設(shè)第一次glm()返回對象記為fit，
第二次返回對象記為fit.od，那么：
pchisq(summary(fit.od) dispersion * fit df.residual, fit df.residual, lower = F))
提供的p值即可對零假設(shè)H0：Φ = 1與備擇假設(shè)H1：Φ≠ 1進(jìn)行檢驗(yàn)。若p很小（小于0.05），你便可拒絕零假設(shè)。
對Affairs數(shù)據(jù)集過度離勢的判斷
fit.aff <- glm(ynaffair ~ age + yearsmarried + religiousness +
rating, family = binomial(), data = Affairs) # 模型fit.aff構(gòu)建。
fit.od <- glm(ynaffair ~ age + yearsmarried + religiousness +
rating, family = quasibinomial(), data = Affairs) # 模型fit.od構(gòu)建。
pchisq(summary(fit.od)$dispersion * fit.aff$df.residual,
fit.aff$df.residual, lower = F) # 過度離勢檢驗(yàn)。
## [1] 0.340122
13.2.4 擴(kuò)展
13.3 泊松回歸
當(dāng)通過一系列連續(xù)型和/或類別型預(yù)測變量來預(yù)測計(jì)數(shù)型結(jié)果變量時(shí)，泊松回歸是一個(gè)非常有用的工具。
泊松回歸（英語：Poisson regression）是用來為計(jì)數(shù)資料和列聯(lián)表建模的一種回歸分析。泊松回歸假設(shè)反應(yīng)變量Y是泊松分布，并假設(shè)它期望值的對數(shù)可被未知參數(shù)的線性組合建模。泊松回歸模型有時(shí)（特別是當(dāng)用作列聯(lián)表模型時(shí)）又被稱作對數(shù)-線性模型。
泊松回歸常用來研究單位時(shí)間/單位面積/單位空間內(nèi)某事件的發(fā)生數(shù)的影響因素。
泊松回歸要求響應(yīng)變量滿足泊松分布。Poisson分布的概率密度函數(shù)：
P表示單位時(shí)間/單位面積/單位空間內(nèi)某事件發(fā)生k次的概率。λ是這個(gè)分布唯一的參數(shù)，λ是該分布的均值，也是該分布的方差，即均值等于方差，λ越大Poisson分布越逼近正態(tài)分布。
泊松回歸的模型：
β0表示各個(gè)“自變量取值為0”時(shí)觀測頻數(shù)的自然對數(shù)值；βi則表示其他自變量取值不變時(shí)，自變量xi每改變一個(gè)單位，所引起的觀測頻數(shù)λ的自然對數(shù)值的改變量，βi為正表示自變量每增加一個(gè)單位觀測頻數(shù)λ會(huì)增加，βi為負(fù)值表示自變量每增加一個(gè)單位觀測頻數(shù)λ會(huì)減少。
《R語言實(shí)戰(zhàn)》例子
遭受輕微或嚴(yán)重間歇性癲癇的病人的年齡和癲癇發(fā)病數(shù)收集了數(shù)據(jù)，包含病人被隨機(jī)分配到藥物組或者安慰劑組前八周和隨機(jī)分配后八周兩種情況。響應(yīng)變量為sumY（隨機(jī)化后八周內(nèi)癲癇發(fā)病數(shù)），預(yù)測變量為治療條件（Trt）、年齡（Age）和前八周內(nèi)的基礎(chǔ)癲癇發(fā)病數(shù)（Base）。之所以包含基礎(chǔ)癲癇發(fā)病數(shù)和年齡，是因?yàn)樗鼈儗憫?yīng)變量有潛在影響。在解釋這些協(xié)變量后，我們感興趣的是藥物治療是否能減少癲癇發(fā)病數(shù)。
數(shù)據(jù)統(tǒng)計(jì)匯總信息查詢：
data(breslow.dat, package="robust") # 調(diào)用數(shù)據(jù)集breslow.dat。
names(breslow.dat) # 查看數(shù)據(jù)集breslow.dat的列名。
## [1] "ID" "Y1" "Y2" "Y3" "Y4" "Base" "Age" "Trt" "Ysum"
## [10] "sumY" "Age10" "Base4"
summary(breslow.dat[c(6, 7, 8, 10)]) # 6，7，8，10列描述性統(tǒng)計(jì)。
## Base Age Trt sumY
## Min. : 6.00 Min. :18.00 placebo :28 Min. : 0.00
## 1st Qu.: 12.00 1st Qu.:23.00 progabide:31 1st Qu.: 11.50
## Median : 22.00 Median :28.00 Median : 16.00
## Mean : 31.22 Mean :28.34 Mean : 33.05
## 3rd Qu.: 41.00 3rd Qu.:32.00 3rd Qu.: 36.00
## Max. :151.00 Max. :42.00 Max. :302.00
響應(yīng)變量圖形考察。
opar <- par(no.readonly=TRUE) # 更改當(dāng)前變量環(huán)境。
par(mfrow=c(1, 2)) # 圖形參數(shù)設(shè)定。
attach(breslow.dat) # 數(shù)據(jù)集加入搜索引擎。
hist(sumY, breaks=20, xlab="Seizure Count",
main="Distribution of Seizures") # 直方圖繪制變量sumY。
boxplot(sumY ~ Trt, xlab="Treatment", main="Group Comparisons") # 箱線圖考察變量sumY。
par(opar) # 還原默認(rèn)變量環(huán)境。
泊松回歸模型擬合
fit24 <- glm(sumY ~ Base + Age + Trt, data=breslow.dat, family=poisson()) # 泊松回歸模型擬合。
summary(fit24) # 返回?cái)M合結(jié)果。
##
## Call:
## glm(formula = sumY ~ Base + Age + Trt, family = poisson(), data = breslow.dat)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -6.0569 -2.0433 -0.9397 0.7929 11.0061
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.9488259 0.1356191 14.370 < 2e-16 ***
## Base 0.0226517 0.0005093 44.476 < 2e-16 ***
## Age 0.0227401 0.0040240 5.651 1.59e-08 ***
## Trtprogabide -0.1527009 0.0478051 -3.194 0.0014 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (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
13.3.1 解釋模型參數(shù)
coef(fit24) # 獲取模型系數(shù)。
## (Intercept) Base Age Trtprogabide
## 1.94882593 0.02265174 0.02274013 -0.15270095
在泊松回歸中，因變量以條件均值的對數(shù)形式ln(λ)來建模。年齡的回歸參數(shù)為0.0227，表明保持其他預(yù)測變量不變，年齡增加一歲，癲癇發(fā)病數(shù)的對數(shù)均值將相應(yīng)增加0.03。截距項(xiàng)即當(dāng)預(yù)測變量都為0時(shí)，癲癇發(fā)病數(shù)的對數(shù)均值。由于不可能為0歲，且調(diào)查對象的基礎(chǔ)癲癇發(fā)病數(shù)均不為0，因此本例中截距項(xiàng)沒有意義。
exp(coef(fit24)) # 指數(shù)化回歸系數(shù)。
## (Intercept) Base Age Trtprogabide
## 7.0204403 1.0229102 1.0230007 0.8583864
現(xiàn)在可以看到，保持其他變量不變，年齡增加一歲，期望的癲癇發(fā)病數(shù)將乘以1.023。這意味著年齡的增加與較高的癲癇發(fā)病數(shù)相關(guān)聯(lián)。更為重要的是，一單位Trt的變化（即從安慰劑到治療組），期望的癲癇發(fā)病數(shù)將乘以0.86，也就是說，保持基礎(chǔ)癲癇發(fā)病數(shù)和年齡不變，服藥組相對于安慰劑組癲癇發(fā)病數(shù)降低了20%。
另外需要牢記的是，與Logistic回歸中的指數(shù)化參數(shù)相似，泊松模型中的指數(shù)化參數(shù)對響應(yīng)變量的影響都是成倍增加的，而不是線性相加。
13.3.2 過度離勢
泊松分布的方差和均值相等。當(dāng)響應(yīng)變量觀測的方差比依據(jù)泊松分布預(yù)測的方差大時(shí)，泊松回歸可能發(fā)生過度離勢。
可能發(fā)生過度離勢的原因有如下幾個(gè)（Coxe et al.，2009）。
?遺漏了某個(gè)重要的預(yù)測變量。
?可能因?yàn)槭录嚓P(guān)。在泊松分布的觀測中，計(jì)數(shù)中每次事件都被認(rèn)為是獨(dú)立發(fā)生的。以癲癇數(shù)據(jù)為例，這意味著對于任何病人，每次癲癇發(fā)病的概率與其他癲癇發(fā)病的概率相互獨(dú)立。但是這個(gè)假設(shè)通常都無法滿足。對于某個(gè)的病人，在已知他已經(jīng)發(fā)生了39次癲癇時(shí)，第一次發(fā)生癲癇的概率不可能與第40次發(fā)生癲癇的概率相同。
?在縱向數(shù)據(jù)分析中，重復(fù)測量的數(shù)據(jù)由于內(nèi)在群聚特性可導(dǎo)致過度離勢。此處并不討論縱向泊松模型。
與Logistic回歸類似，此處如果殘差偏差與殘差自由度的比例遠(yuǎn)遠(yuǎn)大于1，那么表明存在過度離勢。
殘差偏差與殘差自由度的比例核算
deviance(fit24)/df.residual(fit24) # 殘差偏差與殘差自由度的比例。
## [1] 10.1717
泊松模型過度離勢的檢驗(yàn)
library(qcc) # 調(diào)用qcc包。
qcc.overdispersion.test(breslow.dat$sumY, type="poisson") # 泊松模型過度離勢的檢驗(yàn)。
##
## Overdispersion test Obs.Var/Theor.Var Statistic p-value
## poisson data 62.87013 3646.468 0
通過用family = "quasipoisson"替換family = "poisson"，你仍然可以使用glm()
函數(shù)對該數(shù)據(jù)進(jìn)行擬合。
泊松回歸模型構(gòu)建
fit.od1 <- glm(sumY ~ Base + Age + Trt, data=breslow.dat,
family=quasipoisson()) # 泊松回歸模型。
summary(fit.od1) # 返回模型結(jié)果。
##
## Call:
## glm(formula = sumY ~ Base + Age + Trt, family = quasipoisson(),
## data = breslow.dat)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -6.0569 -2.0433 -0.9397 0.7929 11.0061
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.948826 0.465091 4.190 0.000102 ***
## Base 0.022652 0.001747 12.969 < 2e-16 ***
## Age 0.022740 0.013800 1.648 0.105085
## Trtprogabide -0.152701 0.163943 -0.931 0.355702
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 11.76075)
##
## Null deviance: 2122.73 on 58 degrees of freedom
## Residual deviance: 559.44 on 55 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
13.3.3 擴(kuò)展
13.4 小結(jié)
參考資料：
- 《R語言實(shí)戰(zhàn)》（中文版），人民郵電出版社，2013.
- 廣義線性模型-百度，https://baike.baidu.com/item/%E5%B9%BF%E4%B9%89%E7%BA%BF%E6%80%A7%E6%A8%A1%E5%9E%8B
- 廣義線性模型一（Generalized Linear Models，GLM），http://www.itdecent.cn/p/97368848f727
- 泊松回歸-百度，https://baike.baidu.com/item/%E6%B3%8A%E6%9D%BE%E5%9B%9E%E5%BD%92
- 泊松分布，https://zhuanlan.zhihu.com/p/279330443
- 優(yōu)勢比，百度，https://baike.baidu.com/item/%E4%BC%98%E5%8A%BF%E6%AF%94
