R之回归分析广义线性模型(Generalized Linear Model)glm

1. 介绍

广义线性模型(Generalized  Linear Model)是一般线性模型的推广，它使因变量的总体均值通过一个非线性连接函数而依赖于线性预测值，允许响应概率分布为指数分布族中的任何一员。许多广泛应用的统计模型都属于广义线性模型，如常用于研究二元分类响应变量的Logistic回归、Poisson回归和负二项回归模型等。一个广义线性模型包含以下三个部分:

①随机成分。

②线性成分。

③连接函数g。

各种常见的指数型分布及其主要参数

典型的连接函数及对应分布

广义线性模型的参数估计一般不能用最小二乘估计，常用加权最小二乘法或最大似然法估计，各回归系数β需用迭代方法求解。

2. 实现

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, contrasts = NULL, ...)

其中，

formula为拟合公式，与函数lm()中的参数formula用法相同;

family用于指定分布族，包括正态分布(gaussian)、二项分布(binomial)、泊松分布(poisson)和伪伽马分布(Gamma)；

分布族还可以通过选项link来指定连接函数，默认值为family=gaussian (link=identity),二项分布默认值为family=binomial(link=logit);

data指定数据集;

offset指定线性函数的常数部分，通常反映已知信息;

control用于对待估参数的范围进行设置。

例：

车险保单索赔次数分组数据

已知索赔次数服从泊松分布，相应的连接函数常用对数连接函数，模型可以写为

下面用R实现，首先建立数据集，分类变量直接输入定性的取值即可，glm()分析时会自动转换成矩阵X,注意参数family的写法。

> dat=data.frame(

y=c(42, 37, 10, 101, 73, 14),

n=c(500, 1200, 100, 400, 500, 300),

type=rep(c('小','中','大'),2),

gender=rep(c('男','女'),each=3)

)

> dat$logn=log(dat$n)  #风险暴露数取对数

#offset风险单位数事先已知

> dat.glm=glm(y~type+gender,offset=logn,data=dat,family=poisson(link=log))

> summary(dat.glm)  #glm的输出结果

估计的回归系数都是非常显著的;Null deviance可以认为是模型的残差，它的值越小说明模型拟合效果越好;模型的AIC统计量为61.68，它和deviance一起可以用来作为判断标准，选取合适的分布族和链接函数。

下面通过作图来观察模型拟合的效果，首先提取模型的预测值，注意函数predict()提取的是线性部分的拟合值，在对数连接函数下，要得到Y的拟合值，应当再做一次指数变换。以实际观测值为横坐标，模型拟合值为纵坐标作图，散点越接近直线y=x，说明模型的拟合效果越好。

> dat.pre=predict(dat.glm)

> layout(1)  #取消绘图区域分割

> plot(y,exp(dat.pre),xlab='观测值',ylab='拟合值',main="索赔次数的拟合效果",pch="*")

> abline(0,1)  #添加直线y=x，截距为0，斜率为1

若假设上例中的索赔次数服从负二项分布，在R中应输入指令:

> library(MASS)

> attach(dat)

> dat.glmnb=glm.nb(y~type+gender+offset(logn))  #负二项回归

> summary(dat.glmnb)  #输出结果

负二项回归拟合的模型AIC为60.45，残差Null deviance为16.6831，小于泊松回归拟合的残差值，说明负二项分布的广义线性模型更加稳定，但从回归系数的显著性上看，泊松回归拟合的变量系数更加显著。