【数据分析 R语言实战】学习笔记 第九章（下）岭回归及R实现 广义线性模型
【数据分析 R语言实战】学习笔记 第九章（下）岭回归及R实现 广义线性模型
9.4岭回归及R实现
岭回归分析是一种专用于共线性数据分析的有偏估计回归方法，实质上是一种改良的最小二乘估计法，它是通过放弃最小二乘法的无偏性，以损失部分信息、降低精度为代价获得回归系数更为符合实际、更可靠的回归方法，对病态数据的耐受性远远强于最小二乘法。
为β的岭回归估计，其中k为岭参数。显然，岭回归估计β值比最小二乘估计值稳定，当k=0时的岭回归估计就是普通最小二乘估计。
岭迹图：
根据岭迹图我们可以选择合适的k值，称为岭迹法，其一般原则是:
(1)各回归系数的岭估计基本稳定;
(2)最小二乘估计的回归系数符号不合理时，岭估计参数的符号变得合理
(3)回归系数没有不合乎实际意义的绝对值;
(4)残差平方和增加不太多。
R的核心程序包MASS中有专门用于岭回归分析的函数lm.ridge()，其调用格式为
lm.ridge(formula, data, subset, na.action,lambda = 0, model = FALSE, x = FALSE, y = FALSE, contrasts =NULL, ...)
其中，formula是回归模型公式表达形式，形如response~predictors; data指定数据的数据框;当只需要data 的一个子集参与计算时，用参数subset来设置;na.action表示遇到缺失值时应采取的行为;lambda是岭参数的标量或矢量:model, x和Y均为逻辑值，分别表示结果是否返回模型框架、设计矩阵和响应变量。
例：
首先输入数据并查看数据的多重共线性情况
> ?lm.ridge> y=c(8.4,9.6,10.4,11.4,12.2,14.2,15.8,17.9,19.6,20.8)> x1=c(82.9,88,99.9,105.3,117.7,131,148.2,161.8,174.2,184.7)> x2=c(92,93,96,94,100,101,105,112,112,112)> x3=c(17.1,21.3,25.1,29,34,40,44,49,51,53)> x4=c(94,96,97,97,100,101,104,109,111,111)> x=cbind(x1,x2,x3,x4) #将数据按列合并> xx=crossprod(x) #计算矩阵交叉积，结果为矩阵X’X> kappa(xx,exact=T) #计算条件数[1] 81630得到的条件数K=81629.86>1000，认为有严重的多重共线性。因此，不能直接用最小二乘法估计拟合回归方程。考虑用岭回归估计方法分析变量之间的关系，首先绘制岭迹图:
> library(MASS)> plot(lm.ridge(y~x1+x2+x3+x4,lambda=seq(0,0.5,0.001)))
从图中可以看出，曲线变平稳的速度很慢，很难直接得出适当的岭参数k值，而R可以通过函数select()计算出根据几个统计量得到的k值:
> select(lm.ridge(y~x1+x2+x3+x4,lambda=seq(0,0.5,0.001)))modified HKB estimator is 0.00427 modified L-W estimator is 0.0049 smallest value of GCV at 0.009 其中，H KB和L-W分别为不同方法下计算得到的岭参数估计值，据此我们选择lambda=0.0045，从而给出相应的参数估计:
> options(digits=3)> lm.ridge(y~x1+x2+x3+x4,lambda=0.0045) x1 x2 x3 x4 -17.5554 0.0887 -0.2180 0.0202 0.4073 最终得到的岭回归方程为
y=-17.5554+0.0887x1-0.218x2+0.0202x3+0.4073x4
9.5广义线性模型
9.5.1模型理论
广义线性模型(Generalized Linear Model)是一般线性模型的推广，它使因变量的总体均值通过一个非线性连接函数而依赖于线性预测值，允许响应概率分布为指数分布族中的任何一员。许多广泛应用的统计模型都属于广义线性模型，如常用于研究二元分类响应变量的Logistic回归、Poisson回归和负二项回归模型等。一个广义线性模型包含以下三个部分:
①随机成分。
②线性成分。
③连接函数g。
各种常见的指数型分布及其主要参数
典型的连接函数及对应分布
广义线性模型的参数估计一般不能用最小二乘估计，常用加权最小二乘法或最大似然法估计，各回归系数β需用迭代方法求解。
9.5.2 R语言实现
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),而二项分布默认的连接函数是logit，即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) #风险暴露数取对数> dat.glm=glm(y~type+gender,offset=logn,data=dat,family=poisson(link=log))#offset风险单位数事先已知> summary(dat.glm) #glm的输出结果Call:glm(formula = y ~ type + gender, family = poisson(link = log), data = dat, offset = logn)Deviance Residuals: 1 2 3 4 5 6 0.333 -1.498 3.792 -0.209 1.212 -1.780 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -3.797 0.236 -16.06 < 2e-16 ***type小 1.268 0.222 5.71 1.1e-08 ***type中 0.554 0.230 2.42 0.016 * gender女 1.173 0.132 8.91 < 2e-16 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1(Dispersion parameter for poisson family taken to be 1) Null deviance: 163.490 on 5 degrees of freedomResidual deviance: 21.413 on 2 degrees of freedomAIC: 61.68Number of Fisher Scoring iterations: 5估计的回归系数都是非常显著的;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) #输出结果Call:glm.nb(formula = y ~ type + gender + offset(logn), init.theta = 6.879011721, link = log)Deviance Residuals: 1 2 3 4 5 6 -0.395 -1.020 1.535 0.350 0.810 -1.646 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -3.0541 0.4172 -7.32 2.5e-13 ***type小 0.7430 0.4494 1.65 0.098 . type中 0.0208 0.4514 0.05 0.963 gender女 0.7996 0.3528 2.27 0.023 * ---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1(Dispersion parameter for Negative Binomial(6.88) family taken to be 1) Null deviance: 16.6831 on 5 degrees of freedomResidual deviance: 7.0412 on 2 degrees of freedomAIC: 60.45Number of Fisher Scoring iterations: 1 Theta: 6.88 Std. Err.: 5.24 Warning while fitting theta: alternation limit reached 2 x log-likelihood: -50.45 负二项回归拟合的模型AIC为60.45，残差Null deviance为16.6831，小于泊松i口]归拟合的残差值，说明负二项分布的广义线性模型更加稳定，但从回归系数的显著性上看，泊松回归拟合的变量系数更加显著。