本节书摘来自华章计算机《数学建模：基于R》一书中的第2章，第2.1节,作者：薛 毅 更多章节内容可以访问云栖社区“华章计算机”公众号查看。
2.1 回归分析
在许多实际问题中,经常会遇到需要同时考虑几个变量的情况.例如,在电路中,会遇到电压、电流和电阻之间的关系；在炼钢过程中,会遇到钢水中的碳含量和钢材的物理性能（如强度、延伸率等）之间的关系；在医学上,经常测量人的身高、体重,研究人的血压与年龄的关系等,这些变量之间是相互制约的.
通常,变量间的关系有两大类.一类是变量间有完全确定的关系,可用函数关系式来表示，如电路中的欧姆定律I=U/R,其中I表示电流,U表示电压,R表示电阻.
另一类是变量间有一定的关系,但由于情况错综复杂,无法精确研究,或由于存在不可避免的误差等原因,以致它们的关系无法用函数形式表示出来.研究这类变量之间的关系就需要通过大量试验或观测来获取数据,用统计方法去寻找它们之间的关系,这种关系反映了变量间的统计规律.研究这类统计规律的方法之一便是回归分析.
在回归分析中,把变量分成两类.一类是因变量,它们通常是实际问题中所关心的一些指标,通常用Y表示,而影响因变量取值的另一类变量称为自变量,用X1,X2,…,Xp来表示.
2.1.1 线性回归模型
1.模型
设X1,X2,…,Xp为自变量,Y为因变量.它们之间的相关关系可用如下线性关系来描述：Y=β0+β1X1+…+βpXp+ε（2.1）其中ε～N(0,σ2),β0,β1,…,βp和σ2是未知参数.若p=1,则式（2.1）称为一元回归模型，否则称为多元线性回归模型.
设(xi1,xi2,…,xip,yi)（i=1,2,…,n）是(X1,X2,…,Xp,Y)的n次独立观测值,则多元线性模型(2.1)可表示为yi=β0+β1xi1+…+βpxip+εi, i=1,2,…,n（2.2）其中εi～N(0,σ2),且独立同分布.
为书写方便,常采用矩阵形式,令Y=y1
y2
yn, β=β0
β1
βp, X=1x11…x1p
1x21…x2p
1xn1…xnp, ε=ε1
ε2
εn则多元线性模型(2.2)可表示为Y=Xβ+ε（2.3）其中Y为由响应变量构成的n维向量,X为n×(p+1)阶设计矩阵,β为p+1维参数向量,ε是n维误差向量,并且满足E(ε)=0,var(ε)=σ2In.
求参数β的估计值本质上是最小二乘估计,即求函数Q(β)=(Y-Xβ)T(Y-Xβ)（2.4）的最小值β.由多元函数极值的必要条件和充分条件,得到β=(XTX)-1XTY（2.5）从而可得经验回归方程为Y=Xβ,称ε=Y－Y为残差向量.通常取σ2=εTε/(n－p－1)（2.6）为σ2的估计,也称为σ2的最小二乘估计.可以证明：Eσ2=σ2.
2.检验
β的估计值β是否适用,一个重要的识别方法就是进行检验.这里讨论三种检验方法.
(1) 回归系数的检验.回归系数的显著性检验是检验变量Xj的系数是否为0,即原假设为 由于软件可以提供β0的检验情况,所以这里j从0开始,下同.Hj0:βj=0, Hj1:βj≠0, j=0,1,…,p当Hj0成立时,统计量Tj=βjσcjj~t(n－p－1), j=0,1,…,p（2.7）其中cjj是C=(XTX)-1的对角线上第j个元素.如果t统计量的P值<α(通常取0.05),则拒绝原假设,认为βj≠0.
(2) 回归方程的检验.回归方程的显著性检验是检验是否可用线性方程来处理数据,也就是说,方程的系数是否全为0,即假设检验为H0:β0=β1=…=βp=0, H1:β0,β1,…,βp不全为0当H0成立时,统计量F=SSR/pSSE/(n－p－1)~F(p,n－p－1)（2.8）其中SSR=∑ni=1(yi－y)2, SSE=∑ni=1(yi－yi)2
y=1n∑ni=1yi, yi=β0+β1xi1+…+βpxip通常称SSR为回归平方和,称SSE为残差平方和.如果F统计量的P值<α,则拒绝原假设,即可以用线性方程来处理问题.
(3) 相关性检验.相关系数的平方定义为R2=SSRSST(2.9)其中SST=∑ni=1(yi－y)2,称为总离差平方和.通常用R2来衡量Y与X1,X2,…,Xp之间相关的密切程度.当R2接近于0,可以认为Y与X1,X2,…,Xp之间不相关,接近于1表示相关.因此,可以使用R2作为衡量自变量与因变量是否相关的重要指标.
3.计算
在R中,lm()函数可以完成回归方程的计算工作,其使用格式为lm(formula, data, subset, weights, na.action, method = "qr",
model = TRUE, x = FALSE, y = FALSE, qr = TRUE,
singular.ok = TRUE, contrasts = NULL, offset, ...)参数formula为模型公式,形如y ~ 1 + x1 + x2的形式,表示常数项、X1和X2的系数.如果去掉公式中的1,其意义不变.如果需要拟合成齐次线性模型,其公式改为y ~ 0 + x1 + x2,或者改为y ~ -1 + x1 + x2.
data为数据框,由样本数据构成,每一行为一个样本,每一列为一项指标.
subset为可选项,表示所使用的样本子集.weights为可选向量,表示样本所对应的权重.na.action为函数,表示当出现缺失数据(NA)时的处理方法.
method为估计回归系数的计算方法,默认值为"qr",即QR分解方法.
model、x、y和qr为逻辑变量,如果取TRUE,函数的返回值将给出模型的框架、模型矩阵、响应变量及QR分解.
singular.ok为逻辑变量,取FALSE表示奇异值拟合是错误的.
contrasts为可选列表.offset为NULL或者数值向量....为附加参数.
例2.1 根据经验,在人的身高相等的情况下,血压的收缩压Y（千帕）与体重X1（千克）、年龄X2（岁数）有关.现收集了13个男子的数据,如表2.1所示.试建立Y关于X1和X2的线性回归方程.
表2.1 数据表体重/kg年龄/岁数收缩压/kPa176.050120291.520141385.520124482.530126（续）
体重/kg年龄/岁数收缩压/kPa579.030117680.550125774.560123879.050125985.0401321076.5551231182.0401321295.0401551392.520147
解 将数据以表格形式存入文件(exam0201.data),使用read.table()函数将数据读出,再用lm()计算.> rt <- read.table("exam0201.data")
lm.sol <- lm(Y ~ 1 + X1 + X2, data = rt)使用summary()函数可显示全部的计算结果,如> summary(lm.sol)这里列出部分(重要)结果,并予以解释.Coefficients:
EstimateStd.Errort valuePr(>|t|)
(Intercept)-62.9633616.99976-3.7040.004083 **
X12.136560.1753412.1852.53e-07 *
X20.400220.083214.8100.000713 *
Residual standard error: 2.854 on 10 degrees of freedom
Multiple R-squared: 0.9461, Adjusted R-squared: 0.9354
F-statistic: 87.84 on 2 and 10 DF, p-value: 4.531e-07 所列结果的第一段表示系数的估计值、标准差、t统计量和t统计量对应的P值以及显著性标记.
所列结果的第二段分别给出了:残差的标准差,即式(2.6)中的σ.自由度,即n-p-1.相关系数的平方,即R2.修正相关系数的平方,这个值会小于R2,其目的是不要轻易作出自变量与因变量相关的判断.F统计量,F统计量的自由度,即(p,n-p-1),以及F统计量对应的P值.
从计算结果可以看出,回归方程通过系数的检验(t统计量的P值<0.05),方程的检验(F统计量的P值<0.05)和相关性检验(R2=0.9461).因此,回归方程为Y=－62.96+2.136X1+0.4002X2回归系数的显著性检验只是问题的一个方面,有时还可以进一步作回归系数β的区间估计.由式(2.7)可得到βj(j=0,1,…,p)的置信区间［βi－σcjjtα/2(n－p－1), βi+σcjjtα/2(n－p－1)］（2.10）在R中,可以用confint()函数列出回归系数的置信区间,其使用格式为confint(object, parm, level = 0.95,...)参数object为lm生成的对象.parm为需要作区间估计的参数,默认值为全部参数.level为置信水平,默认值为0.95.
例如,对于例2.1,下面的命令> confint(lm.sol)列出β的95％的置信区间2.5 %97.5 %
(Intercept) -100.8411862-25.0855320
X11.74587092.5272454
X20.21480770.5856246上述区间不包含零,同样说明系数是显著的.
4.预测
经过检验,如果回归效果显著,就可以利用回归方程进行预测.所谓预测,就是对给定的回归自变量的值,预测对应的回归因变量所有可能的取值范围,因此,这是一个区间估计问题.
给定x0=(x01,x02,…,x0p)T,回归方程的真实值为y0=β0+β1x01+…+βpx0p+ε0但其值是未知的,只能得到相应的估计值y0=β0+β1x01+…+βpx0p（2.11）设置信水平为1－α,则y0的预测区间为［y0tα/2(n－p－1)σ1+xT0(XTX)－1x0］（2.12）其中X为设计矩阵,x0=(1,x01,x02,…,x0p)T.
E(y0)的置信区间为［y0tα/2(n－p－1)σxT0(XTX)－1x0］（2.13）其中X与x0的意义与式(2.12)的意义相同.
在R中,可用predict()函数计算y0的估计值、y0的预测区间,以及E(y0)的置信区间,其使用格式为predict(object, newdata, se.fit = FALSE,
scale = NULL, df = Inf,
interval = c("none", "confidence", "prediction"),
level = 0.95, type = c("response", "terms"),
terms = NULL, na.action = na.pass,
pred.var = res.var/weights, weights = 1, ...)参数object为lm()函数生成的对象.
newdata为数据框,由预测点构成,如果该值缺省,将计算已知数据的回归值.
se.fit为逻辑变量,表示是否输出预测值的标准差、自由度和残差尺度信息,默认值为FALSE.
scale为计算标准差的尺度参数.df为尺度参数的自由度.
interval为字符串,表示作预测的区间类型,取"none"(默认值)表示不作计算;取"confidence"表示计算E(y0)置信区间;取"prediction"表示计算y0预测区间.
level为置信水平,默认值为0.95,在计算预测区间或置信区间时用到.
type为预测类型,或者取"response"(默认值),或者取"terms".
terms为选择项,只能取NULL,1,2等,当type="terms"时为全部选择项.
na.action为函数,表示在newdata中出现缺失数据(NA)时的处理方法,默认的预测为NA.
pred.var为预测区间的方差.weights为数值向量或单侧模型公式,用于预测时方差的权....为附加参数.
例2.2（继例2.1） 设x0=(80,40)T,求y0的估计值y0,以及y0的预测区间和E(y0)的置信区间（取置信水平为0.95）.
解 注意,数据的输入必须是数据框的形式,下面是R程序和计算结果.> newdata <- data.frame(X1 = 80, X2 = 40)
predict(lm.sol, newdata, interval="prediction")
fitlwrupr
1 123.9699117.2889130.6509
predict(lm.sol, newdata, interval="confidence")
fitlwrupr
1 123.9699121.9183126.02152.1.2 回归诊断
所谓回归诊断的问题,其主要内容有以下几个方面:
(1) 关于误差项是否满足独立性、等方差性(也称为方差齐性)和正态性;
(2) 选择的线性模型是否合适;
(3) 是否存在异常样本点;
(4) 自变量之间是否存在高度相关, 即是否有多重共线性现象.
1.残差的检验
关于残差的检验是检验模型的误差是否满足正态性和方差齐性.残差检验中最简单且直观的方法是画出模型的残差图.以残差εi为纵坐标,以拟合值yi,或者是对应的数据观测序号i,或者是数据观测时间为横坐标作散点图,这类图形统称为残差图.
为检验建立的多元线性回归模型是否合适,可以通过回归值Y与残差的散点图来检验.其方法是画出回归值Y与普通残差的散点图((Yi,εi),i=1,2,…,n),或者画出回归值Y与标准残差的散点图((Yi,ri),i=1,2,…,n),其图形可能会出现下面三种情况(如图2.1所示).
图2.1 回归值Y与残差的散点图
对于图2.1a的情况,不论回归值Y的大小,而残差εi(或标准化残差ri)具有相同的分布,并满足模型的各假设条件;对于图2.1b的情况,表示回归值Y的大小与残差的波动大小有关系,即等方差的假设有问题;对于图2.1c,表示线性模型不合适,应考虑非线性模型.
对于图2.1a,如果大部分点都落在中间部分,而只有少数几个点落在外边,则这些点对应的样本,可能有异常值存在.
如何判断哪些是异常的呢?通常的方法是画出标准化残差图.由3σ原则可知,在残差中,大约会有95％的点落在2σ的范围内.如果某个点的残差落在2σ(对于标准化残差就是2)范围之外,就有理由认为该点可能是异常值点.
最后,还需要对残差作正态性检验.较为简单的方法是画出残差的正态Q-Q散点图,若这些点位于一条直线上,则说明残差服从正态分布;否则,认为不服从正态分布.
在R中,residuals()函数(或者resid()函数)计算回归模型的残差,rstandard()函数计算回归模型的标准化(也称为内学生化)残差,rstudent()函数计算回归模型的外学生化残差.所谓外学生化残差,就是删除第i个样本数据后得到的标准化残差.这些函数的使用格式为residuals(object, ...)
resid(object, ...)
rstandard(model, infl=lm.influence(model, do.coef=FALSE),
sd = sqrt(deviance(model)/df.residual(model)), ...)
rstudent(model, infl=lm.influence(model, do.coef=FALSE),
res = infl$wt.res, ...)参数object或model为lm生成的对象.infl为lm.influence返回值得到的影响结构.sd为模型的标准差.res为模型残差.
例2.3（继例2.1） 对例2.1得到的模型作残差分析,画出残差图.
解 画出残差图,其中横坐标是回归值(由fitted.values()函数计算),纵坐标是标准化残差.pre <- fitted.values(lm.sol); rst <- rstandard(lm.sol)
plot(pre, rst, xlab = "Fitted Values",
ylab = "Standardized Residuals", xlim = c(116, 158))
text(pre, rst, adj = c(-0.5, 0.5))所绘图形如图2.2所示.从图形可以看出,残差基本满足齐方差性.
图2.2 例2.1模型的残差散点图
仔细分析,3号样本、6号样本和13号样本残差的绝对值超过1.5,也就是说,这三个样本点的残差超过2σ,这说明它们可能是异常值点.
在R中,plot()函数提供了非常方便的残差分析图,其使用格式为plot(x, which = c(1:3, 5),
caption = list("Residuals vs Fitted", "Normal Q-Q",
"Scale-Location", "Cook's distance",
"Residuals vs Leverage",
expression("Cook's dist vs Leverage " * h［ii］ / (1 - h［ii］))),
panel = if(add.smooth) panel.smooth else points,
sub.caption = NULL, main = "",
ask = prod(par("mfcol")) <length(which) && dev.interactive(),
...,
id.n = 3, labels.id = names(residuals(x)),
cex.id = 0.75, qqline = TRUE,
cook.levels = c(0.5, 1.0),
add.smooth = getOption("add.smooth"),
label.pos = c(4, 2), cex.caption = 1)参数x是由lm生成的对象.which是开关变量,其值为1～6,默认值为子集{1,2,3,5},即绘出第1、2、3和5号散点图.
该函数可以画出6张诊断图,其中,第1张是残差与预测值的残点图,第2张是残差的正态Q-Q图,第3张是标准化残差的平方根与预测值的散点图,第4张是Cook距离图,第5张是残差与高杠杆值(也称为帽子值)的散点图,第6张是Cook距离与高杠杆值的散点图.
例如,对于例2.1,可用plot(lm.sol, 1)绘出残差图,plot(lm.sol, 2)绘出正态Q-Q图,用以检查残差的方差齐性和正态性.
- Box-Cox变换
在作回归分析时,通常假设回归方程的残差εi具有齐性,即等方差,也就是图2.1a所示的正常情况.如果残差不满足齐性,则其计算结果可能会出现问题.一种解决问题的方法是Box-Cox变换.
在出现异方差(图2.1b)的情况下,Box-Cox变换可以使回归方程的残差满足齐性要求,它对Y作如下变换Y(λ)=Yλ－1λ,λ≠0
lnY,λ=0（2.14）其中λ为待定参数.
Box-Cox变换主要有两项工作.第一项是作变换,这一点很容易由式(2.14)得到.第二项是确定参数λ的值,这项工作较为复杂,需要用极大似然估计的方法才能确定出λ的值.
在R软件中,boxcox()函数可以绘出不同参数下对数似然函数的目标值,可以通过图形来选择参数λ的值,其使用格式为boxcox(object, lambda = seq(-2, 2, 1/10), plotit = TRUE,
interp, eps = 1/50, xlab = expression(lambda),
ylab = "log-Likelihood", ...)参数object为lm生成的对象.
lambda为参数λ,默认值为-2到2,间隔值是0.1.
plotit为逻辑变量,表示是否画出图形,默认值为TRUE.
interp为逻辑变量,表示在计算时是否使用三次样条插值.
eps为控制精度,默认值为0.02.
xlab为横轴的标记,默认值为"lambda".ylab为纵轴的标记, 默认值为"log-Likelihood".
...为附加参数.
注意:在调用函数boxcox()之前,需要先加载MASS程序包,或使用命令library(MASS).
例2.4 某公司为了研究产品的营销策略,对产品的销售情况进行了调查.设Y为某地区该产品的家庭人均购买量（单位:元）,X为家庭人均收入（单位:元）.表2.2给出了53个家庭的数据.试通过这些数据建立Y与X的关系式.
解 将数据存入数据文件(exam0204.data)中,该文件是纯文本文件,其中前6行为自变量X的数据,后6行为因变量Y的数据.
读取数据,作线性回归.X <- scan("exam0204.data", nlines = 6)
Y <- scan("exam0204.data", skip = 6)
lm.sol <- lm(Y ~ X); summary(lm.sol)从计算结果(略)可以看出, 回归方程通过系数的检验和方程的检验.但这个模型是否就适合于这组数据呢?再作进一步的分析.library(MASS) #% 加载MASS程序包
作图, 共4张
op <- par(mfrow = c(2, 2), mar = .4 + c(4, 4, 1, 1),
oma = c(0, 0, 2, 0))
第1张, 残差与预测散点图
plot(fitted(lm.sol), resid(lm.sol),
cex = 1.2, pch = 21, col = "red", bg = "orange",
xlab = "Fitted Value", ylab = "Residuals")
第2张, 确定参数lambda
boxcox(lm.sol, lambda = seq(0, 1, by = 0.1))
##%% Box-Cox变换后, 作回归分析
lambda <- 0.55; Ylam <- (Y^lambda-1)/lambda
lm.lam <- lm(Ylam ~ X); summary(lm.lam)
第3张, 变换后残差与预测散点图
plot(fitted(lm.lam), resid(lm.lam),
cex = 1.2, pch=21, col = "red", bg = "orange",
xlab = "Fitted Value", ylab = "Residuals")
第4张, 回归曲线和相应的散点
beta0<-lm.lam$coefficients［1］
beta1<-lm.lam$coefficients［2］
curve((1+lambda*(beta0+beta1*x))^(1/lambda),
from = min(X), to = max(X), col = "blue", lwd = 2,
xlab = "X", ylab = "Y")
points(X, Y, pch = 21, cex = 1.2, col = "red", bg = "orange")
mtext("Box-Cox Transformations", outer = TRUE, cex = 1.5)
par(op)程序给出的图形(共4张)如图2.3所示.
图2.3 Box-Cox变换前与变换后的情况
从第1张图可以看出,由原始数据得到的残差图呈喇叭口形状,属于异方差情况,这样的数据需要作Box-Cox变换.在变换前先确定参数λ(调用函数boxcox),得到第2张图.从第2张图中看到,当λ=0.55时,对数似然函数达到最大值,因此,选择参数λ=0.55.作Box-Cox变换(Y(λ)=Yλ-1λ),变换后再作回归分析,然后画出残差的散点图(第3张图),从第3张图看出,喇叭口形状有很大改善.第4张图是给出曲线Y=(1+λβ0+λβ1X)1/λ和相应的散点图.
3.多重共线性现象
当自变量彼此相关时,回归模型的结果可能非常令人困惑.估计的效应会由于模型中的其他自变量而改变数值,甚至是改变符号.故在分析时,了解自变量间的关系的影响是很重要的.这一复杂问题常称为共线性或多重共线性.
记x(1),x(2),…,x(p)为自变量X1,X2,…,Xp经过标准化得到的向量,X=(x(1),x(2),…,x(p)).设λ为XTX的一个特征值,φ为对应的特征向量,其模长为1,即φTφ=1.若λ≈0,则XTXφ=λφ≈0用φT左乘上式,得到φTXTXφ=λφTφ=λ≈0(2.15)所以Xφ≈0,即向量x(1),x(2),…,x(p)之间有近似的线性关系,也就是说,自变量之间存在多重共线性.
度量多重共线性严重程度的一个重要指标是方矩XTX的条件数,即κ(XTX)=XTX·(XTX)-1=λmax(XTX)λmin(XTX)其中X是数据标准化得到的设计矩阵,λmax(XTX)和λmin(XTX)分别表示方矩XTX的最大和最小特征值.
直观上,条件数刻画了XTX的特征值差异的大小.从实际应用的经验角度,一般若κ<100,则认为多重共线性的程度很小;若100≤κ≤1000,则认为存在中等程度或较强的多重共线性;若κ>1000,则认为存在严重的多重共线性.
例2.5(法国经济分析数据) 考虑进口总额Y与三个自变量:国内总产值X1,存储量X2,总消费量X3（单位：10亿法郎）之间的关系.现收集了1949年至1959年共11年的数据,如表2.3所示.试对此数据作回归分析.
表2.3 法国经济分析数据X1X2X3Y1149.34.2108.115.92161.24.1114.816.43171.53.1123.219.04175.53.1126.919.15180.81.1132.118.86190.72.2137.720.47202.12.1146.022.78212.45.6154.126.59226.15.0162.328.110231.95.1164.327.611239.00.7167.626.3
解 将数据以及表格形式存放在文件(exam0205.data)中,读取数据,作回归分析：france <- read.table("exam0205.data")
lm.sol <- lm(y ~ 1 + x1 + x2 + x3, data = france)
lm.sum <- summary(lm.sol)
coef(lm.sum)
EstimateStd.Errort valuePr(>|t|)
(Intercept)-10.127988161.21215996-8.35532316.899183e-05
x1-0.051396160.07027999-0.73130584.883443e-01
x20.586949040.09461842 6.20332744.438135e-04
x30.286848680.10220811 2.80651572.627710e-02虽然方程的系数均通过检验,但仔细研究,发现它并不合理.回到问题本身,Y代表进口量,X1表示国内总产值,而对应系数的符号为负,也就是说,国内的总产值越高,其进口量越少,这与实际情况显然是不相符的.问其原因,三个变量存在着多重共线性,其计算如下:> X <- as.matrix(france［1:3］); min(eigen(cor(X))$values)
［1］ 0.002690889λmin=0.0027≈0.
4.岭估计
当设计矩阵存在着复共线性关系时,最小二乘估计的性质不够理想,有时甚至很坏.在这种情况下,需要一些新的估计方法.在这些方法中,岭估计是最有影响且应用较为广泛的估计方法.对于线性模型y=Xβ+ε岭估计的回归系数定义为β(k)=(XTX+kI)－1XTy(2.16)其中k>0为可选择参数,称为岭参数或偏参数.
当k取不同值时,得到不同的估计,因此,岭估计β(k)为一类估计.当k=0时,β(0)=(XTX)-1XTy就是通常的最小二乘估计.从严格意义上讲,最小二乘估计是岭估计类中的一个估计.
在R中,lm.ridge()函数(需要加载MASS程序包)计算岭估计,其使用格式为lm.ridge(formula, data, subset, na.action,
lambda = 0, model = FALSE,
x = FALSE, y = FALSE, contrasts = NULL, ...)参数formula为模型公式.data为数据框.
subset为可选向量,表示观察值的子集.na.action为函数,表示当数据中出现缺失数据时的处理方法，默认值为NA.
lambda为岭参数,即岭估计中的参数k,默认值为0.
例2.6 对例2.5的法国经济分析数据作岭回归分析.
解 在岭回归中,最困难的是确定岭参数k*.
最常用的一种方法是岭迹法,以横坐标为k,纵坐标为βi(k),画出岭参数的图形,即岭迹.原则上应该选取βi(k)稳定的最小k值为k.也就是说,当k≥k时,βi(k)不再改变符号.画图(程序名:exam0206.R).library(MASS)
lm.rid <- lm.ridge(y ~ 1 + x1 + x2 + x3, data = france,
lambda = seq(0, 0.2, length = 50))
plot(lm.rid)绘的岭迹图形如图2.4所示.从图形可以看出,取k*=0.04即能满足岭迹选择k的条件.计算结果如下:> lm.ridge(y ~ 1 + x1 + x2 + x3, data = france, lambda = 0.04)
x1 x2 x3
-9.4956573 0.0198393 0.5976872 0.1828705 图2.4 法国经济分析数据的岭迹曲线
三个系数均为正数,所以此模型更合理一些.
2.1.3 逐步回归
1.“最优”回归方程的选择
在实际问题中,影响因变量Y的因素很多,人们可以从中挑选若干个变量建立回归方程,这便涉及变量选择的问题.
一般来说,如果在一个回归方程中忽略了对Y有显著影响的自变量,那么所建立的方程必与实际有较大的偏离,但变量选得过多,使用就不方便,特别地,当方程中含有对Y影响不大的变量时,可能因为误差平方和的自由度的减小而使σ2的估计增大,从而影响使用回归方程作预测的精度.因此,适当地选择变量以建立一个“最优”的回归方程是十分重要的.
什么是“最优”回归方程呢？对于这个问题有许多不同的准则,在不同的准则下“最优”回归方程也可能不同.这里讲的“最优”是指从可供选择的所有变量中选出对Y有显著影响的变量建立方程,并且在方程中不含对Y无显著影响的变量.
在上述意义下,可以有多种方法来获得“最优”回归方程,如“一切子集回归法”“前进法”“后退法”“逐步回归法”等,其中“逐步回归法”由于计算机程序简便,因而使用较为普遍.
2.逐步回归的计算
R提供了较为方便的“逐步回归”计算函数——step()函数,它是以AIC AIC是Akaike Information Criterion的缩写,是由日本统计学家赤池弘次提出的.信息统计量为准则,通过选择最小的AIC信息统计量,来达到删除或增加变量的目的.step()函数的使用格式为step(object, scope, scale = 0,
direction = c("both", "backward", "forward"),
trace = 1, keep = NULL, steps = 1000, k = 2, ...)参数object为lm或glm得到的回归模型.
scope为逐步回归的搜索区域,或者由单一公式定义,或者由包含upper和lower两个公式构成的列表定义,默认值为全部变量.
scale用于AIC统计量的确定,默认值为0,表示scale是由极大似然估计得到.
direction为逐步回归的搜索方向,取"both"(默认值)表示“一切子集回归法”(即可增加变量也可减少变量),取"backward"表示“后退法”(只减少变量),取"forward"表示“前进法”(只增加变量).当scope使用默认值时,则direction只使用"backward".
trace为数值,如果取正值,则列出逐步回归的过程,默认值为1.keep为过滤函数.steps为正整数,表示逐步回归的最大步数,默认值为1000.
k为正数,表示自由度数目的倍数,只有当k = 2时,其返回值中的AIC才是真正的AIC值....为附加参数.
例2.7(Hald水泥问题) 某种水泥在凝固时放出的热量Y（卡/克）与水泥中四种化学成分X1（3CaO·Al2O3含量的百分比）,X2（3CaO·SiO2含量的百分比）,X3（4CaO·Al2O3·Fe2O3含量的百分比）和X4（2CaO·SiO2含量的百分比）有关,现测得13组数据（见表2.4）.希望从中选出主要的变量,建立Y关于它们的线性回归方程.
表2.4 Hald水泥问题数据X1X2X3X4YX1X2X3X4Y172666078.58131224472.52129155274.39254182293.131156820104.3102147426115.94113184787.611140233483.8575263395.9121166912113.361155922109.2131068812109.47371176102.7
解 将数据以表格形式存入文本文件(exam0207.data),读取数据,作回归分析：cement <- read.table("exam0207.data")
lm.sol <- lm(Y ~ X1 + X2 + X3 + X4, data = cement)
summary(lm.sol)运算后(结果略)可以看到:回归系数没有一项通过检验(取α=0.05),效果不好.
下面用函数step()作逐步回归.> lm.step <- step(lm.sol)
Start: AIC=26.94
Y ~ X1 + X2 + X3 + X4
Df Sum of Sq RSS AIC
- X310.109147.97324.974
- X410.247048.11125.011
- X212.972550.83625.728
47.86426.944 - X1125.950973.81530.576
Step: AIC=24.97
Y ~ X1 + X2 + X4
Df Sum of Sq RSS AIC
47.9724.974
- X419.9357.9025.420
- X2126.7974.7628.742
- X11820.91868.8860.629从程序的运行结果可以看到,用全部变量作回归方程时,AIC值为26.94.接下来显示的数据告诉我们,如果去掉变量X3,则相应的AIC值为24.97;如果去掉变量X4,则相应的AIC值为25.42.后面的类推.由于去掉变量X3可以使AIC达到最小,因此,R软件自动去掉变量X3,进行下一轮计算.
在下一轮计算中,无论去掉哪一个变量,AIC值均会升高.因此,R软件终止计算,得到“最优”的回归方程.
下面分析一下计算结果.用函数summary()提取相关信息(只列系数部分).> su`javascript
mmary(lm.step)
Coefficients:
EstimateStd.Errort valuePr(>|t|)
(Intercept)71.648314.14245.0660.000675 *
X11.45190.117012.4105.78e-07 *
X20.41610.18562.2420.051687
.`
X4-0.23650.1733-1.3650.205395由显示结果看到:回归系数检验的显著性水平有很大提高,但变量X2,X4系数检验的显著性水平仍不理想.下面如何处理呢?
从step函数的最后运算结果来看,如果去掉变量X4,则AIC值会从24.97增加到25.42,是增加最少的.另外,除AIC准则外,残差的平方和也是逐步回归的重要指标之一,从直观来看,拟合越好的方程,残差的平方和应越小.去掉变量X4,残差的平方和上升9.93,也是最少的.因此,从这两项指标来看,应该再去掉变量X4.用人工的方法去掉X4,看一下计算结果.> lm.opt<-lm(Y ~ X1+X2, data = ce`javascript
ment); summary(lm.opt)
Coefficients:
EstimateStd.Errort valuePr(>|t|)
(Intercept)52.577352.2861723.005.46e-10 *
X11.468310.1213012.112.69e-07 *
X20.662250.0458514.445.03e-08 *
Residual standard error: 2.406 on 10 degrees of freedom
Multiple R-Squared: 0.9787, Adjusted R-squared: 0.9744
F-statistic: 229.5 on 2 and 10 DF, p-value: 4.407e-09这个结果应该还是满意的,因为所有的检验均是显著的.最后得到“最优”的回归方程为Y=52.58+1.4683X1+0.6623X2改变step()函数中的某些参数,可能得到不同的结果.例如,> step(lm.sol, trace = 0, k = 3)
Call:
lm(formula = Y ~ X1 + X2, data = cement)
Coefficients:
(Intercept)X1X2
52.57731.46830.6623直接去掉变量X3和X4.这里取trace = 0,所以不列出中间的计算过程.
从增加变量的角度考虑逐步回归,例如> lm0 <- lm(Y ~ 1, data = cement)
> lm.ste <- step(lm0, scope = ~X1+X2+X3+X4, k = 4)(计算结果略).由于这里取k = 4,最后还是剩下变量X1和X2.
