广义线性模型

广义线性模型简介

广义线性模型是线性模型的扩展，主要是对非正态因变量的分析
广义线性拟合的核心是最大似然估计，而不是最小二乘

拟合模型如下
μy=β0+∑pj=1βjxj

其中，beta是系数，mu是优势比的对数，beta系数是对优势比的影响。通过拟合求得的就是

• 1.1 什么是非正态因变量？

  • 类别型因变量：

    1. 二值变量（是／否、活／死）
    2. 多分类因变量（差／良好／优秀）

  • 计数型因变量：一周交通事故的数目、每日酒水消耗量（通过均值和方差是否相关来判断）

我们可以通过两个例子看一下两种变量

类别型：自变量x：人口统计信息、人际关系信息、个人信息因变量y：婚姻状况计数型：自变量x：药物使用状况因变量y：八周以来累计的癫痫次数

• 1.2 glm函数介绍
  glm函数和lm函数相比，多了很多参数。

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

• formula 是结果变量～预测变量1+预测变量2
• link 后面根据拟合的不同分布函数加上不同的连接函数，也可以直接加分布函数
• data 后面加使用的数据

今天用到的两个分布

概率分布 连接函数 binomial logit poisson log

Logistic回归

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

我们主要使用AER包里的Affairs数据框为例。

• 2.1 数据集的分析
  601个参与者，收集9个变量，包括婚外情频率、性别、年龄、宗教信仰程度、学历等等

data(Affairs, package = "AER")summary(Affairs)table(Affairs$affairs)

但是其实我们需要的是二值型的预测结果，我们可以做一个转换，把affair变为ynaffair

Affairs$ynaffair[Affairs$affairs > 0] <- 1Affairs$ynaffair[Affairs$affairs == 0] <- 0Affairs$ynaffair <- factor(Affairs$ynaffair,                          levels = c(0,1),                          labels = c("NO","YES"))table(Affairs$ynaffair)

• 2.2 拟合

然后我们来做预测，把ynaffair作为结果变量进行拟合

fit.full <- glm(ynaffair ~ gender + age + yearsmarried + children +                   religiousness + education + occupation + rating,                data= Affairs, family = binomial())summary(fit.full)

从回归系数的p值可以排除几个不显著的变量，再重新拟合

fit.reduce <- glm(ynaffair ~ age + yearsmarried + religiousness +                     rating, data= Affairs, family = binomial())summary(fit.reduce)

我们可以用anova函数对两次拟合进行比较，广义线性回归应该用卡方检验

anova(fit.reduce, fit.full, test="Chisq")

两次拟合的效果显然相同

重要的是，需要回归诊断哦

先看看预测值和残差的图形

plot(predict(fit.reduce, type="response"),     residuals(fit.reduce,type= "deviance"))

遗憾的是，下面这个最方便的诊断函数用不了，只能针对lm对象

gvmodel <- gvlma(fit.reduce)summary(gvmodel)

• 2.3 解释模型

• 2.3.1 看看回归系数（优势比的变化影响因子）

coef(fit.reduce)exp(coef(fit.reduce))

其中，intercept是截距

+2.3.2 看看概率的影响（更加直观）

观察某个预测变量在各个水平对结果概率的影响（需要把其他变量设为均值）

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

使用测试数据集进行预测，可以看出不太婚姻幸福水平的人婚变的可能

testdata$prob <- predict(fit.reduce,newdata=testdata, type="response")testdata

以此类推还可以预测年龄的影响概率

关于过度离势的内容放在最后一起讲

过度离势 指的是观测变量的方差大于所属分布的方差，主要后果是会造成显著性检验不精确，导致可解释性降低

Poisson 回归

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

我们主要使用robust包里的Breslow数据框为例。

• 3.1 数据集的分析
  两个组，前后各八周，收集9个变量，包括前后发病数、年龄、治疗条件、等等

先看看数据

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)

可以看出药物的作用导致发病数和极差都变小了(判断出是泊松分布)

接下来我们开始拟合

• 3.2 拟合

然后我们来做预测，把ynaffair作为结果变量进行拟合

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

• 3.3 解释模型

• 3.3.1 看看回归系数

coef(fit)exp(coef(fit))

其中，intercept是截距

+3.3.2 也可以看看概率的影响（同上，不写了，读者请自己练习）

现在我们来看看过度离势问题

这里需要用到qcc包

我们先进行判断

qcc.overdispersion.test(breslow.dat$sumY, type="poisson")

因为p值小于0.005，所以显然存在过度离势。

这时候，你需要使用另一种family参数来拟合，那就是quasipoisson

fit.od <- glm(sumY ~ Base + Age + Trt, data= breslow.dat, family = quasipoisson())qcc.overdispersion.test(breslow.dat$sumY, type="quasipoisson")