之前学过标准回归分析，用预测变量预测响应变量，响应变量要符合正态分布。当结果变量为二值型或者计数型可以分别用Logistic回归，泊松回归进行阐释建模；可均通过基础glm()函数进行分析。它们统属于广义线性模型概念（而标准线性模型也是广义线性模型的一个特例）。概念与原理详见p280

要点一、Logistic回归

通过一系列连续型/类别型变量来预测二值型结果变量。

• 实验示例：研究一些因素对是否会出轨（会/不会，0/1）的影响。数据包Affairs在AER包内。

1、首先创建二值型因子变量

data(Affairs, package="AER")
#Affairs为有601个观测，9个变量的数据表。
summary(Affairs)
table(Affairs$affairs)
#大致瞅一下不同出轨次数的人数
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"))
#进一步转化为二值因子
table(Affairs$ynaffair)
str(Affairs$ynaffair)

2、然后glm()Logistic回归建模

fit.full <- glm(ynaffair ~ gender + age + yearsmarried + children + 
                  religiousness + education + occupation +rating,
                data=Affairs,family=binomial())
#格式基本与lm()相同，重要的是famly= 的参数选择
summary(fit.full)
#由结果看出仅有4个变量有显著意义。

summary(fit.full)

• 由上述返回数据，可以尝试进一步单独选择这四个变量进行分析。

fit.reduced <- glm(ynaffair ~ age + yearsmarried + religiousness + 
                     rating, data=Affairs, family=binomial())
summary(fit.reduced)
#结果p值表明拟合的更好，不信？那再对比检查下

• anova() 可用于模型拟合优度的比较，在第七节曾使用过。

anova(fit.reduced, fit.full, test="Chisq")
#结果p值无显著意义，表明二者拟合程度一样好，于是选择4个变量的简单模型更符合要求。

summary(fit.reduced)

3、解释模型参数

• 主要是对回归系数的解释；由于连接函数为对数优势比，所以解释意义比较复杂

coef(fit.reduced)

此时回归系数表示当其他预测变量不变时，以单位预测变量的变化可引起的响应变量对数优势比的变化。

优势比=π/(1-π) ，其中π为二值（0/1）分布中，1（出轨）出现的概率。

exp(coef(fit.reduced))

指数转换后，比如控制其它变量不变，婚外情的优势比将乘以1.106。

指数转换

值得注意的是：与标准回归不同，泊松回归、Logistic回归对响应变量的影响是成倍增加的，不是线性增加的。即年龄增加10岁，就是乘以1.106的10次方

4、预测结果概率

即根据回归模型结果，在某一变量不同下水平下，预测可能的出轨概率（0/1中1的概率）
（1）首先创建一个新表（dataframe格式，变量名相同），控制其它变量为平均水平时，某一变量为所有可能的水平。

testdata <- data.frame(rating = c(1, 2, 3, 4, 5),
                       age = mean(Affairs$age),
                       yearsmarried = mean(Affairs$yearsmarried),
                       religiousness = mean(Affairs$religiousness))

（2）根据回归结果，用predict()函数预测不同水平的概率

testdata$prob <- predict(fit.reduced, newdata=testdata, type="response")
testdata  #查看结果

testdata

• 换个年龄变量再看看

testdata <- data.frame(rating = mean(Affairs$rating),
                       age = seq(17, 57, 10), 
                       yearsmarried = mean(Affairs$yearsmarried),
                       religiousness = mean(Affairs$religiousness))
testdata$prob <- predict(fit.reduced, newdata=testdata, type="response")
testdata

5、检测过度离势

• 定义：观测到的响应变量的方差大于期望的二巷分布的方差。
• 意义：过度离势会导致奇异的标准误检验和不精确的显著性检验。
• 思路：检验φ=残差偏差/残差自由度；若值接近1，则没有过度离势。

#法1
deviance(fit.reduced)/df.residual(fit.reduced)
#返回结果为1.032

#法2  假设检验，零假设为φ=1，若p>0.05，则无过度离势
fit <- glm(ynaffair ~ age + yearsmarried + religiousness +
             rating, family = binomial(), data = Affairs)
#第一次使用family = binomial()拟合；
fit.od <- glm(ynaffair ~ age + yearsmarried + religiousness +
                rating, family = quasibinomial(), data = Affairs)
#第二次使用family = binomial()拟合；
pchisq(summary(fit.od)$dispersion * fit$df.residual,  
       fit$df.residual, lower = F)
# 返回结果为0.34

若出现过度离势，仍可使用glm() 函数拟合Logistic回归，但是需要将二项分布改为类二项分布（quasibinomial distribution）.

要点二、泊松回归

通过一系列连续型/类别型变量来预测计数型结果变量。

• 实验示例：研究癫痫发病次数（sumY）与用药（Trt）的关系，此外有两个协变量。数据包breslow.dat在robust包内。

1、确定目标数据，并观察计数型变量分布

data(breslow.dat, package="robust")
names(breslow.dat)
summary(breslow.dat[c(6, 7, 8, 10)])
#查看目标数据分布
opar <- par(no.readonly=TRUE)
#绘图观察响应变量分布
par(mfrow=c(1, 2))
attach(breslow.dat)
hist(sumY, breaks=20, xlab="Seizure Count", 
     main="Distribution of Seizures")
boxplot(sumY ~ Trt, xlab="Treatment", main="Group Comparisons")
par(opar)

计数型变量分布

左图反映了因变量（y）的偏倚特性；右图中经药物治疗的发病数似乎减少且方差也减小了（泊松分布中较小的方差伴随着较小的均值）

与标准最小二乘回归不同，泊松回归不关注方差异质性（即协变量的影响）。

2、拟合泊松回归

fit <- glm(sumY ~ Base + Age + Trt, data=breslow.dat, family=poisson())
#协变量还是要放在表达式的前面
summary(fit)

泊松回归

3、解释模型参数

coef(fit)
#同样是对数作为连接函数，参数经过log转换过
exp(coef(fit))
#如上进行指数化系数，便于解释意义

image.png

如上图结果，表明保持其他变量不变，年龄增加一岁，发病数会乘以1.023；更为重要的是保持协变量（Base、Age）不变，服药组相对安慰组癫痫发病数降低了近20%。

• 预测新模型看看

test <- data.frame(Base = 50,
                   Age = c(20:30),
                   Trt = "placebo")
test$sumY <- round(predict(fit, test, type = "response"), 0)

4、检验过度离势

定义，意义等基本同前，不过有两点要注意：

• 泊松分布的方差和均值相等；
• 处理计数型数据经常发生过度离势。

法1
deviance(fit.reduced)/df.residual(fit.reduced)
#返回结果为10，远大于1
法2：假设检验 qcc包
library(qcc)
qcc.overdispersion.test(breslow.dat$sumY, type="poisson")
#返回p值<0.05,进一步表明存在过度离势

• 若存在过度离势，可使用family="quasipoisson"替换family="poisson"

fit.od <- glm(sumY ~ Base + Age + Trt, data=breslow.dat,
              family=quasipoisson())
summary(fit.od)

summary(fit.od)

但尴尬的是，Trt的p值大于0.05；即考虑过度离势，并控制协变量，并没有重组的证据表明药物相对于未用药组能够显著降低发病次数。仅是拿本数据举例阐释泊松模型~

关于Logistic回归与泊松回归的用法，教材还有一些扩展，分别在p289与p294就不深入学习了。
参考教材《R语言实战（第2版）》