,其中I为N×N单位矩阵。
用最小二乘法原理,令残差平方和
最小,得到
为β的最优线性无偏估计量〔高斯-马尔可夫定理〕。
2.𝜎2的估计和T检验
选取𝜎2的估计量:
如此
假设t值的绝对值相当大,就可以在适当选定的置信水平上否认原假设,参数的1-α置信区间可由下式得出:
其中tα/2为与α%显著水平有关的t分布临界值。
3.R2和F检验
假设因变量不具有0平均值,如此必须对R2做如下改良:
随着模型中增添新的变量,R2的值必定会增大,为了去掉这种增大的干扰,还需要对R2进展修正〔校正拟合优度对自由度的依赖关系〕:
做假设检验:
H0:
𝛽1=…=𝛽N=0;H1:
𝛽1…,𝛽N至少有一个≠0;
使用F统计量做检验,
假设F值较大,如此否认原假设。
4.回归诊断
〔1〕残差图分析
残差图就是以残差
为纵坐标,某一个适宜的自变量为横坐标的散点图。
回归模型中总是假定误差项是独立的正态分布随机变量,且均值为零和方差相等为𝜎2.如果模型适合于观察到的数据,那么残差作为误差的无偏估计,应根本反映误差的假设特征。
即残差图应该在零点附近对称地密布,越远离零点的地方就疏散〔在形象上似有正态趋势〕,如此认为模型与数据拟合得很好。
假设残差图呈现如图〔a〕所示的形式,如此认为建立的回归模型正确,更进一步再诊断“学生化残差〞是否具有正态性:
图〔b〕明确数据有异常点,应处理掉它重新做回归分析〔在SAS的REG回归过程步中用来度量异常点影响大小的统计量是COOKD统计量〕;
图〔c〕残差随x的增大而增大,图〔d〕残差随x的增大而先增后减,都属于异方差。
此时应该考虑在回归之前对数据y或x进展变换,实现方差稳定后再拟合回归模型。
原如此上,当误差方差变化不太快时取变换
;当误差方差变化较快时取变换logy或lny;当误差方差变化很快时取变换1/y;还有其他变换,如著名的Box-Cox幂变换
.
图〔e〕〔f〕表示选用回归模型是错误的。
〔2〕共线性
回归分析中很容易发生模型中两个或两个以上的自变量高度相关,从而引起最小二乘估计可能很不准确〔称为共线性问题〕。
在实际中最常见的问题是一些重要的自变量很可能由于在假设检验中t值不显著而被不恰当地剔除了。
共线性诊断问题就是要找出哪些变量间存在共线性关系。
〔3〕误差的独立性
回归分析之前,要检验误差的独立性。
假设误差项不独立,那么回归模型的许多处理,包括误差项估计、假设检验等都将没有推导依据。
由于残差是误差的合理估计,因此检验统计量通常是建立在残差的根底上。
检验误差独立性的最常用方法,是对残差的一阶自相关性进展Durbin-Watson检验。
H0:
误差项是相互独立的;H1:
误差项是相关的
检验统计量:
DW接近于0,表示残差中存在正自相关;如果DW接近于4,表示残差中存在负自相关;如果DW接近于2,表示残差独立性。
二、R语言实现
还是用lm()函数实现,不同是需要设置更复杂的formula格式:
y~x1+x2——只考虑自变量的主效应〔y=k1x1+k2x2〕,y~.表示全部自变量的主效应;
y~x1+x2+x1:
x2——考虑主效应和交互效应〔y=k1x1+k2x2+k3x1x2〕;
y~x1*x2——考虑全部主效应和交互效应的简写〔效果同上〕;
y~(x1+x2+x3)^2——考虑主效应以与至2阶以下的交互效应,相当于x1+x2+x3+x1:
x2+x2:
x3+x1:
x3
y~x1%in%x2——x1含于x2,相当于x2+x2:
x1
y~(x1+x2)^2-x1:
x2——表示从(x1+x2)^2中去掉x1:
x2
y~x1+I((x2+x3)^2)——使用I()函数,相当于用(x2+x3)^2计算出新变量h,然后y~x1+h
function——在表达式中使用数学函数,例如log(y)~x1+x2
三、实例
例2现有1990~2009年财政收入的数据revenue.txt:
各变量分别表示:
y:
财政收入〔亿元〕
x1:
第一产业国内生产总值〔亿元〕
x2:
第二产业国内生产总值〔亿元〕
x3:
第三产业国内生产总值〔亿元〕
x4:
人口数〔万人〕
x5:
社会消费品零售总额〔亿元〕
x6:
受灾面积〔万公顷〕
做多元线性回归分析。
setwd("E:
/办公资料/R语言/R语言学习系列/codes")
revenue=read.table("revenue.txt",header=TRUE)
lm.reg=lm(y~x1+x2+x3+x4+x5+x6,revenue)
summary(lm.reg)
Residuals:
Min1QMedian3QMax
-295.71-173.5226.5990.16370.01
Coefficients:
EstimateStd.ErrortvaluePr(>|t|)
(Intercept)6.046e+043.211e+0318.8298.12e-11***
x1-1.171e-018.638e-02-1.3560.19828
x23.427e-023.322e-021.0320.32107
x36.182e-014.103e-0215.0671.31e-09***
x4-5.152e-012.930e-02-17.5851.91e-10***
x5-1.104e-012.878e-02-3.8370.00206**
x6-1.864e-021.023e-02-1.8230.09143.
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:
234.8on13degreesoffreedom
MultipleR-squared:
0.9999,AdjustedR-squared:
0.9999
说明:
拟合优度R2=0.9999,效果非常好。
但是多元回归时,自变量个数越多,拟合优度必然越好,所以还要看检验效果和回归系数是否显著。
结果解释、回归方程、回归预测与前文类似〔略〕。
结合显著性代码可看出:
x1和x2不显著,x6只在0.1显著水平下显著,故应考虑剔除x1和x2.
R语言中提供了update()函数,用来在原模型的根底上进展修正,还可以对变量进展运算,其根本格式为:
update(object,formula.,...,evaluate=TRUE)
其中,object为前面拟合好的原模型对象;
formula指定模型的格式,原模型不变的局部用“.〞表示,只写出需要修正的地方即可,例如
update(lm.reg,.~.+x7)表示添加一个新的变量
update(lm.reg,sqrt(.)~.)表示对因变量y开方,再重新拟合回归模型
lm.reg2<-update(lm.reg,.~.-x1-x2)#剔除自变量x1,x2
summary(lm.reg2)
Residuals:
Min1QMedian3QMax
-325.62-147.5414.07108.28427.42
Coefficients:
EstimateStd.ErrortvaluePr(>|t|)
(Intercept)6.339e+042.346e+0327.0203.89e-14***
x36.584e-011.548e-0242.523<2e-16***
x4-5.438e-011.981e-02-27.4453.09e-14***
x5-1.392e-011.918e-02-7.2562.80e-06***
x6-1.803e-029.788e-03-1.8420.0854.
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:
233.6on15degreesoffreedom
MultipleR-squared:
0.9999,AdjustedR-squared:
0.9999
〔三〕逐步回归
多元线性回归模型中,并不是所有的自变量都与因变量有显著关系,有时有些自变量的作用可以忽略。
这就需要考虑怎样从所有可能有关的自变量中挑选出对因变量有显著影响的局部自变量。
逐步回归的根本思想是,将变量一个一个地引入或剔出,引入或剔出变量的条件是“偏相关系数〞经检验是显著的,同时每引入或剔出一个变量后,对已选入模型的变量要进展逐个检验,将不显著变量剔除或将显著的变量引入,这样保证最后选入的所有自变量都是显著的。
逐步回归每一步只有一个变量引入或从当前的回归模型中剔除,当没有回归因子能够引入或剔出模型时,该过程停止。
R语言中,用step()函数进展逐步回归,以AIC信息准如此作为选入和剔除变量的判别条件。
AIC是日本统计学家赤池弘次,在熵概念的根底上建立的:
AIC=2(p+1)-2ln(L)
其中,p为回归模型的自变量个数,L是似然函数。
注:
AIC值越小越被优先选入。
根本格式:
step(object,direction=,steps=,k=2,...)
其中,object是线性模型或广义线性模型的返回结果;
direction确定逐步回归的方法,默认“both〞综合向前向后法,“backward〞向后法〔先把全部自变量参加模型,假设无统计学意义如此剔出模型〕,“forward〞向前法〔先将局部自变量参加模型,再逐个添加其它自变量,假设有统计学意义如此选入模型〕;
steps表示回归的最大步数,默认1000;
k默认=2,输出为AIC值,=log(n)有时输出BIC或SBC值。
另外,有时还需要借助使用drop1(object)和add1(object)函数,其中object为逐步回归的返回结果,判断剔除或选入一个自变量,AIC值的变化情况,以筛选选入模型的自变量。
lm.step<-step(lm.reg)
summary(lm.step)
Call:
lm(formula=y~x3+x4+x5+x6,data=revenue)
Residuals:
Min1QMedian3QMax
-325.62-147.5414.07108.28427.42
Coefficients:
EstimateStd.ErrortvaluePr(>|t|)
(Intercept)6.339e+042.346e+0327.0203.89e-14***
x36.584e-011.548e-0242.523<2e-16***
x4-5.438e-011.981e-02-27.4453.09e-14***
x5-1.392e-011.918e-02-7.2562.80e-06***
x6-1.803e-029.788e-03-1.8420.0854.
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:
233.6on15degreesoffreedom
MultipleR-squared:
0.9999,AdjustedR-squared:
0.9999
说明:
默认输出每步的结果〔略〕,进展了3步回归,逐步剔除了自变量x1和x2,AIC值逐步减小,最终得到最优的模型。
drop1(lm.step)
Singletermdeletions
Model:
y~x3+x4+x5+x6
DfSumofSqRSSAIC
3<-lm(y~x3+x4+x5,revenue)
3)
Call:
lm(formula=y~x3+x4+x5,data=revenue)
Residuals:
Min1QMedian3QMax
-336.34-186.821.5289.46437.84
Coefficients:
EstimateStd.ErrortvaluePr(>|t|)
(Intercept)6.284e+042.494e+0325.1912.66e-14***
x36.614e-011.651e-0240.066<2e-16***
x4-5.467e-012.118e-02-25.8131.81e-14***
x5-1.412e-012.053e-02-6.8773.72e-06***
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:
250.5on16degreesoffreedom
MultipleR-squared:
0.9999,AdjustedR-squared:
0.9998
说明:
使用drop1()函数考察分别剔除每个自变量,AIC值变化的情况,可以看出不剔除x6与剔除x6,AIC值只从222.40变大到224.47,相对其它自变量变化很小。
所以,可以考虑剔除掉x6,重新做多元线性回归。
〔四〕回归诊断
回归分析之后,还需要从残差的随机性、强影响分析、共线性方面进展诊断。
一、残差诊断
(1)残差
y.res<-lm.reg3$residuals#回归模型的残差
y.fit<-predict(lm.reg3)#回归模型的预测值
plot(y.res~y.fit,main="残差图")#绘制残差图,以预测值作为横坐标
说明:
从图形看,残差分布比拟均匀,大致满足随机性。
shapiro.test(y.res)#残差的正态性检验
Shapiro-Wilknormalitytest
说明:
p值=0.2622>0.05,承受原假设,即残差服从正态分布。
(2)标准化残差
残差与数据的数量级有关,除以标准误差后得到标准化残差。
理想的标准化残差服从N(0,1).
rs<-rstandard(lm.reg3)#得到标准化残差
plot(rs~y.fit,main="标准残差图")
shapiro.test(rs)#标准化残差的正态性检验
Shapiro-Wilknormalitytest
data:
rs
(3)学生化残差
为了回避标准化残差的方差齐性假设,使用学生化残差。
rst<-rstudent(