《R语言实战》自学笔记70-广义线性模型

第四部分 高级方法

第13章 广义线性模型

广义线性模型适用条件

假设因变量为正态分布不成立的情况，如：

结果变量可能是类别型的二值变量（比如：是/否、通过/失败、活着/死亡）和多分类变量（比如差/良好/优秀）都显然不是正态分布。

结果变量可能是计数型的：（比如，一周交通事故的数目，每日酒水消耗的数量）。这类变量都是非负的有限值，而且它们的均值和方差通常都是相关的（正态分布变量间不是如此，而是相互独立）。

广义线性模型扩展了线性模型的框架，它包含了非正态因变量的分析。

常见广义线性模型

1、标准线性模型：是广义线性模型的一个特例，因变量服从正态分布；
2、Logistic回归：因变量服从二项分布，相应变量是二值型（0和1）；
3、泊松回归：因变量服从泊松分布，相应变量是计数型。

广义线性模型定义：

广义线性模型(generalize linear model，GLM)是线性模型的扩展，通过联结函数建立响应变量的数学期望值与线性组合的预测变量之间的关系。其特点是不强行改变数据的自然度量，数据可以具有非线性和非恒定方差结构。是线性模型在研究响应值的非正态分布以及非线性模型简洁直接的线性转化时的一种发展。

线性回归和逻辑回归都是广义线性模型的一种特殊形式。

13.1 广义线性模型和glm()函数

广义线性模型拟合的形式为：

{\rm{g(}}{\mu _Y}{\rm{)}} = {\beta _0} + {\sum\nolimits_{j = 1}^p \beta _j}{X_j}
其中g(μ_Y)是条件均值的函数（称为连接函数）。另外，你可放松Y为正态分布的假设，改为Y服从指数分布族中的一种分布即可。设定好连接函数和概率分布后，便可以通过最大似然估计的多次迭代推导出各参数值。

13.1.1 glm()函数

R中可通过glm函数（还可用其他专门的函数）拟合广义线性模型。

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）和相应默认的连接函数（function）。

image.png

假设你有一个响应变量（Y）、三个预测变量（X1、X2、X3）和一个包含数据的数据框（mydata）。
Logistic回归适用于二值响应变量（0,1）。模型假设Y服从二项分布，线性模型的拟合形式为：
ln(\frac {\pi}{1-\pi})=\beta_0 + \sum\nolimits_{j = 1}^p{\beta_jX_j}

其中\pi = \mu_Y是Y的条件均值（即给定一系列X的值时Y = 1的概率），\frac {\pi}{1-\pi}为Y = 1时的优势比，log(\frac {\pi}{1-\pi})为对数优势比，或logit。本例中，log(\frac {\pi}{1-\pi})为连接函数，概率分布为二项分布，可用如下代码拟合Logistic回归模型：

glm(Y~x1+x2+x3,family=binomial(link="logit"),data=mydata)

泊松回归适用于在给定时间内响应变量为事件发生数目的情形。它假设Y服从泊松分布，线性模型的拟合形式为：

ln(\lambda)=\beta_0 + \sum\nolimits_{j = 1}^p{\beta_jX_j}

其中λ是Y的均值（也等于方差）。此时，连接函数为log(λ)，概率分布为泊松分布，可用如下代码拟合泊松回归模型：
glm(Y~x1+x2+x3,family=poisson(link="log"),data=mydata)

值得注意的是，标准线性模型也是广义线性模型的一个特例。如果令连接函数g(μY) =μY或恒等函数，并设定概率分布为正态（高斯）分布，那么：
glm(Y~x1+x2+x3,family=gaussian(link="identity"),data=mydata)
生成的结果与下列代码的结果相同：
lm(Y~x1+x2+x3,data=mydata)
总之，广义线性模型通过拟合响应变量的条件均值的一个函数（不是响应变量的条件均值），假设响应变量服从指数分布族中的某个分布（并不仅限于正态分布），极大地扩展了标准线性模型。模型参数估计的推导依据的是极大似然估计，而非最小二乘法。

13.1.2 连用的函数

与分析标准线性模型时lm()连用的许多函数在glm()中都有对应的形式，其中一些常见的函数见表13-2。

image.png

13.1.3 模型拟合和回归诊断

当评价模型的适用性时，你可以绘制初始响应变量的预测值与残差的图形。例如，如下代码可绘制一个常见的诊断图：
plot(predict(model,type="response"),residuals(model,type="deviance"))
其中，model为glm()函数返回的对象。
R将列出帽子值（hat value）、学生化残差值和Cook距离统计量的近似值。

car包绘制诊断图
library(car)
influencePlot(model)
它可以创建一个综合性的诊断图。在后面的图形中，横轴代表杠杆值，纵轴代表学生化残差值，而绘制的符号大小与Cook距离大小成正比。

13.2 Logistic回归

当通过一系列连续型和/或类别型预测变量来预测二值型结果变量时，Logistic回归是一个非常有用的工具。

Logistic回归模型的适用条件
1 因变量为二分类的分类变量或某事件的发生率，并且是数值型变量。但是需要注意，重复计数现象指标不适用于Logistic回归。
2 残差和因变量都要服从二项分布。二项分布对应的是分类变量，所以不是正态分布，进而不是用最小二乘法，而是最大似然法来解决方程估计和检验问题。
3 自变量和Logistic概率是线性关系。
4 各观测对象间相互独立。

logistic回归和线性回归

逻辑回归（Logistic Regression）与线性回归（Linear Regression）都是一种广义线性模型（generalized linear model）。
逻辑回归假设因变量 y 服从伯努利分布，而线性回归假设因变量 y 服从高斯分布。

Sigmoid函数
也称为逻辑函数（Logistic function）
g(z)=\frac{1}{1+e^{-z}}

image.png

从上图可以看到sigmoid函数是一个s形的曲线，它的取值在[0,1]之间，在远离0的地方函数的值会很快接近0或者1。

《R语言实战》例子，AER包Affairs数据集为例。

该数据从601个参与者身上收集了9个变量，包括一年来婚外私通的频率以及参与者性别、年龄、婚龄、是否有小孩、宗教信仰程度（5分制，1分表示反对，5分表示非常信仰）、学历、职业（逆向编号的戈登7种分类），还有对婚姻的自我评分（5分制，1表示非常不幸福，5表示非常幸福）。

数据的描述统计

data(Affairs, package="AER") # 调用AER包数据集Affairs。
head(Affairs) # 显示数据前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) # 数据描述。
##     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) # 数据集Affairs中affairs列频数统计。
## 
##   0   1   2   3   7  12 
## 451  34  17  19  42  38

将affairs转化为二值型因子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")) # 新列转为因子。
table(Affairs$ynaffair) # 统计因子频数。
## 
##  No Yes 
## 451 150

factor函数
factor(x, levels = sort(unique(x), na.last = TRUE),labels = levels, exclude = NA, ordered = is.ordered(x), nmax = NA)
x是要创建因子的向量，
levels用来指定因子水平值（不指定时由向量x不同值求得）；
labels用来指定水平的名字（不指定时由用水平值的对应字符串）；exclude表示从向量x中剔除的水平值；
ordered是一个逻辑型选项，用来指定因子的水平是否有次序。
nmax是水平的上限数量。

对二值型因子变量进行logistic回归

fit.full <- glm(ynaffair ~ gender + age + yearsmarried + children + 
                  religiousness + education + occupation +rating,
                data=Affairs,family=binomial()) # Logistic回归，响应变量为ynaffair，预测变量为gender，age，yearsmarried，children，religiousness，education，occupation，rating。
summary(fit.full) # 返回回归结果。
## 
## 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

结果解读：从回归系数的p值（最后一栏）可以看到，性别、是否有孩子、学历和职业对方程的贡献都不显著（你无法拒绝参数为0的假设）。原假设为所有预测变量对方程均有贡献且相同

去除不显著因子，重新拟合

fit.reduced <- glm(ynaffair ~ age + yearsmarried + religiousness + 
                     rating, data=Affairs, family=binomial()) # 去除不显著预测变量，重新拟合模型。
summary(fit.reduced) # 返回结果。
## 
## 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

两回归模型的比较，对于广义线性回归，可用卡方检验。

anova(fit.reduced, fit.full, test="Chisq") # 两个拟合模型的对比。
## 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），表明四个预测变量的新模型与九个完整预测变量的模型拟合程度一样好。

13.2.1 解释模型参数

在Logistic回归中，响应变量是Y=1的对数优势比（log）。回归系数含义是当其他预测变量不变时，一单位预测变量的变化可引起的响应变量对数优势比的变化。
对于二值型Logistic回归，某预测变量n单位的变化引起的较高值上优势比的变化为exp(βj)^n，它反映的信息可能更为重要。

优势比
优势比(odds ratio；OR)是另外一种描述概率的方式。优势比将会告诉我们某种推测的概率比其反向推测的概率大多少。换句话说，优势比是指某种推测为真的概率与某种推测为假的概率的比值。比如下雨的概率为0.25，不下雨的概率为0.75。0.25与0.75的比值可以约分为1比3。因此，我们可以说今天将会下雨的优势比为1:3(或者今天不会下雨的概率比为3:1)。

coef(fit.reduced) # 查看fit.reduced模型的回归系数。
##   (Intercept)           age  yearsmarried religiousness        rating 
##    1.93083017   -0.03527112    0.10062274   -0.32902386   -0.46136144
exp(coef(fit.reduced)) # 指数化结果。
##   (Intercept)           age  yearsmarried religiousness        rating 
##     6.8952321     0.9653437     1.1058594     0.7196258     0.6304248
exp(confint(fit.reduced)) # 在优势比尺度上得到系数95%的置信区间。
##                   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 评价预测变量对结果概率的影响

使用predict()函数，你可观察某个预测变量在各个水平时对结果概率的影响。

创建一个虚拟数据集，设定年龄、婚龄和宗教信仰为它们的均值，婚姻评分的范围为1～5。

testdata <- data.frame(rating = c(1, 2, 3, 4, 5),
                       age = mean(Affairs$age),
                       yearsmarried = mean(Affairs$yearsmarried),
                       religiousness = mean(Affairs$religiousness)) # 虚拟数据集testdata的创建。

predict函数预测数据集testdata中rating变量的概率

testdata$prob <- predict(fit.reduced, newdata=testdata, type="response") # 数据集testdata中rating变量的概率预测。
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

predict函数预测数据集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") # 数据集testdata中age变量的概率预测。
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

13.2.3 过度离势

抽样于二项分布的数据的期望方差是σ2 = nπ(1-π)，n为观测数，π为属于Y = 1组的概率。
所谓过度离势，即观测到的响应变量的方差大于期望的二项分布的方差。过度离势会导致奇异的标准误检验和不精确的显著性检验。
当出现过度离势时，仍可使用glm()函数拟合Logistic回归，但此时需要将二项分布改为类二项分布（quasibinomial distribution）
检测过度离势的一种方法是比较二项分布模型的残差偏差与残差自由度。
\phi = \frac {残差偏差}{残差自由度}
当比值比1大很多，你便可认为存在过度离势。

你还可以对过度离势进行检验。为此，你需要拟合模型两次，第一次使用 family =
"binomial"，第二次使用family = "quasibinomial"。假设第一次glm()返回对象记为fit，
第二次返回对象记为fit.od，那么：
pchisq(summary(fit.od) dispersion * fit df.residual, fit df.residual, lower = F))
提供的p值即可对零假设H0：Φ = 1与备择假设H1：Φ≠ 1进行检验。若p很小（小于0.05），你便可拒绝零假设。

对Affairs数据集过度离势的判断

fit.aff <- glm(ynaffair ~ age + yearsmarried + religiousness +
             rating, family = binomial(), data = Affairs) # 模型fit.aff构建。
fit.od <- glm(ynaffair ~ age + yearsmarried + religiousness +
                rating, family = quasibinomial(), data = Affairs) # 模型fit.od构建。
pchisq(summary(fit.od)$dispersion * fit.aff$df.residual,  
       fit.aff$df.residual, lower = F) # 过度离势检验。
## [1] 0.340122

13.2.4 扩展

13.3 泊松回归

当通过一系列连续型和/或类别型预测变量来预测计数型结果变量时，泊松回归是一个非常有用的工具。
泊松回归（英语：Poisson regression）是用来为计数资料和列联表建模的一种回归分析。泊松回归假设反应变量Y是泊松分布，并假设它期望值的对数可被未知参数的线性组合建模。泊松回归模型有时（特别是当用作列联表模型时）又被称作对数-线性模型。
泊松回归常用来研究单位时间/单位面积/单位空间内某事件的发生数的影响因素。
泊松回归要求响应变量满足泊松分布。Poisson分布的概率密度函数：
P(y|x=k)=\frac{\lambda^k}{k!}e^{-\lambda}
P表示单位时间/单位面积/单位空间内某事件发生k次的概率。λ是这个分布唯一的参数，λ是该分布的均值，也是该分布的方差，即均值等于方差，λ越大Poisson分布越逼近正态分布。

泊松回归的模型：

ln\lambda=\beta_0+\beta_1x_1+\beta_2x_2+....\beta_ix_i

β0表示各个“自变量取值为0”时观测频数的自然对数值；βi则表示其他自变量取值不变时，自变量xi每改变一个单位，所引起的观测频数λ的自然对数值的改变量，βi为正表示自变量每增加一个单位观测频数λ会增加，βi为负值表示自变量每增加一个单位观测频数λ会减少。

《R语言实战》例子
遭受轻微或严重间歇性癫痫的病人的年龄和癫痫发病数收集了数据，包含病人被随机分配到药物组或者安慰剂组前八周和随机分配后八周两种情况。响应变量为sumY（随机化后八周内癫痫发病数），预测变量为治疗条件（Trt）、年龄（Age）和前八周内的基础癫痫发病数（Base）。之所以包含基础癫痫发病数和年龄，是因为它们对响应变量有潜在影响。在解释这些协变量后，我们感兴趣的是药物治疗是否能减少癫痫发病数。

数据统计汇总信息查询：

data(breslow.dat, package="robust") # 调用数据集breslow.dat。
names(breslow.dat) # 查看数据集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列描述性统计。
##       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

响应变量图形考察。

opar <- par(no.readonly=TRUE) # 更改当前变量环境。
par(mfrow=c(1, 2)) # 图形参数设定。
attach(breslow.dat) # 数据集加入搜索引擎。
hist(sumY, breaks=20, xlab="Seizure Count", 
     main="Distribution of Seizures") # 直方图绘制变量sumY。
boxplot(sumY ~ Trt, xlab="Treatment", main="Group Comparisons") # 箱线图考察变量sumY。
par(opar) # 还原默认变量环境。
image.png

泊松回归模型拟合

fit24 <- glm(sumY ~ Base + Age + Trt, data=breslow.dat, family=poisson()) # 泊松回归模型拟合。
summary(fit24) # 返回拟合结果。
## 
## 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 解释模型参数

coef(fit24) # 获取模型系数。
##  (Intercept)         Base          Age Trtprogabide 
##   1.94882593   0.02265174   0.02274013  -0.15270095

在泊松回归中，因变量以条件均值的对数形式ln(λ)来建模。年龄的回归参数为0.0227，表明保持其他预测变量不变，年龄增加一岁，癫痫发病数的对数均值将相应增加0.03。截距项即当预测变量都为0时，癫痫发病数的对数均值。由于不可能为0岁，且调查对象的基础癫痫发病数均不为0，因此本例中截距项没有意义。

exp(coef(fit24)) # 指数化回归系数。
##  (Intercept)         Base          Age Trtprogabide 
##    7.0204403    1.0229102    1.0230007    0.8583864

现在可以看到，保持其他变量不变，年龄增加一岁，期望的癫痫发病数将乘以1.023。这意味着年龄的增加与较高的癫痫发病数相关联。更为重要的是，一单位Trt的变化（即从安慰剂到治疗组），期望的癫痫发病数将乘以0.86，也就是说，保持基础癫痫发病数和年龄不变，服药组相对于安慰剂组癫痫发病数降低了20%。
另外需要牢记的是，与Logistic回归中的指数化参数相似，泊松模型中的指数化参数对响应变量的影响都是成倍增加的，而不是线性相加。

13.3.2 过度离势

泊松分布的方差和均值相等。当响应变量观测的方差比依据泊松分布预测的方差大时，泊松回归可能发生过度离势。
可能发生过度离势的原因有如下几个（Coxe et al.，2009）。
遗漏了某个重要的预测变量。
可能因为事件相关。在泊松分布的观测中，计数中每次事件都被认为是独立发生的。以癫痫数据为例，这意味着对于任何病人，每次癫痫发病的概率与其他癫痫发病的概率相互独立。但是这个假设通常都无法满足。对于某个的病人，在已知他已经发生了39次癫痫时，第一次发生癫痫的概率不可能与第40次发生癫痫的概率相同。
在纵向数据分析中，重复测量的数据由于内在群聚特性可导致过度离势。此处并不讨论纵向泊松模型。
与Logistic回归类似，此处如果残差偏差与残差自由度的比例远远大于1，那么表明存在过度离势。

残差偏差与残差自由度的比例核算

deviance(fit24)/df.residual(fit24) # 残差偏差与残差自由度的比例。
## [1] 10.1717

泊松模型过度离势的检验

library(qcc) # 调用qcc包。
qcc.overdispersion.test(breslow.dat$sumY, type="poisson") # 泊松模型过度离势的检验。
##                    
## Overdispersion test Obs.Var/Theor.Var Statistic p-value
##        poisson data          62.87013  3646.468       0

通过用family = "quasipoisson"替换family = "poisson"，你仍然可以使用glm()
函数对该数据进行拟合。

泊松回归模型构建

fit.od1 <- glm(sumY ~ Base + Age + Trt, data=breslow.dat,
              family=quasipoisson()) # 泊松回归模型。
summary(fit.od1) # 返回模型结果。
## 
## 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 扩展

13.4 小结

参考资料：

  1. 《R语言实战》（中文版），人民邮电出版社，2013.
  2. 广义线性模型-百度，https://baike.baidu.com/item/%E5%B9%BF%E4%B9%89%E7%BA%BF%E6%80%A7%E6%A8%A1%E5%9E%8B
  3. 广义线性模型一（Generalized Linear Models，GLM），//www.greatytc.com/p/97368848f727
  4. 泊松回归-百度，https://baike.baidu.com/item/%E6%B3%8A%E6%9D%BE%E5%9B%9E%E5%BD%92
  5. 泊松分布，https://zhuanlan.zhihu.com/p/279330443
  6. 优势比，百度，https://baike.baidu.com/item/%E4%BC%98%E5%8A%BF%E6%AF%94
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 211,884评论 6赞 492
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 90,347评论 3赞 385
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 157,435评论 0赞 348
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 56,509评论 1赞 284
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 65,611评论 6赞 386
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 49,837评论 1赞 290
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 38,987评论 3赞 408
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 37,730评论 0赞 267
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 44,194评论 1赞 303
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,525评论 2赞 327
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,664评论 1赞 340
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 34,334评论 4赞 330
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 39,944评论 3赞 313
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,764评论 0赞 21
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,997评论 1赞 266
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 46,389评论 2赞 360
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 43,554评论 2赞 349

推荐阅读更多精彩内容