回归分析
1.概述
- 选择回归模型和回归方程
- 一般线性模型被解释变量是一个服从正态分布的连续型数值变量。若研究它如何受多个数值型解释变量的影响，则选择的回归模型是一元或多元回归模型；若研究它如何受到离散型数值变量以及分类型变量的影响，则选择带有虚拟变量的回归模型
- 广义线性模型的被解释变量是0-1变量。若研究它如何受多个解释变量的影响，则建立logistic模型。若被解释变量是计数变量：最近一个月内的购物次数...，采用泊松分布
- 回归方程的参数估计与检验
- 最小二/极大似然
- 回归方程的显著性、回归系数的显著性
- 回归诊断
- 回归模型对数据是有假设的，若无法满足这些假设，则方程可能因为存在很大偏差而没有意义
- 回归方程可能受到异常点的影响。是否需要在排除其影响下重新建立方程
- 多重共线性的诊断
- 模型验证
- 回归方程的主要目的是预测，预测的精度或误差是模型验证的重要内容。
2.建立回归模型
2.1模型
- 由解释变量引起的Y的线性变化
- 由误差项引起的
2.2假设
alpson服从正态分布
2.3参数估计
2.4 R程序
lm(R公式，data=数据框名)
coefficients(回归分析结果对象名)
1.画图
CarData<-read.table(file="CarData.txt",header=TRUE) CarData$ModelYear<-as.factor(CarData$ModelYear)
#查看相关关系 pairs(~MPG+weight+displacement+horsepower,data=CarData)
2.建立模型
##########建立线性回归经验方程 Fit<-lm(MPG~weight+displacement+horsepower,data=CarData) #############浏览回归分析结果 >coefficients(Fit) (Intercept) weight displacement horsepower 44.855935695 -0.005351593 -0.005768819 -0.041674144
方程显著；回归系数:displacement不显著
2.4.回归诊断
#方程显著性
>summary(Fit) Call: lm(formula = MPG ~ weight + displacement + horsepower, data = CarData) Residuals: Min 1Q Median 3Q Max -11.3347 -2.8028 -0.3402 2.2037 16.2409 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 44.8559357 1.1959200 37.507 < 2e-16 *** weight -0.0053516 0.0007124 -7.513 4.04e-13 *** displacement -0.0057688 0.0065819 -0.876 0.38132 horsepower -0.0416741 0.0128139 -3.252 0.00125 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 4.241 on 388 degrees of freedom (6 observations deleted due to missingness) Multiple R-squared: 0.707, Adjusted R-squared: 0.7047 F-statistic: 312 on 3 and 388 DF, p-value: < 2.2e-16 >confint(Fit) 2.5 % 97.5 % (Intercept) 42.50464109 47.207230298 weight -0.00675215 -0.003951036 displacement -0.01870945 0.007171815 horsepower -0.06686744 -0.016480849
- F：方程显著性
- t:系数显著性
因为方程中包含系数不显著的解释变量，所以模型需要改进
Fit<-lm(MPG~weight+horsepower,data=CarData)
summary(Fit)
Call:
lm(formula = MPG ~ weight + horsepower, data = CarData)
Residuals:
Min 1Q Median 3Q Max
-11.0762 -2.7340 -0.3312 2.1752 16.2601
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 45.6402108 0.7931958 57.540 < 2e-16 ***
weight -0.0057942 0.0005023 -11.535 < 2e-16 ***
horsepower -0.0473029 0.0110851 -4.267 2.49e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.24 on 389 degrees of freedom
(6 observations deleted due to missingness)
Multiple R-squared: 0.7064, Adjusted R-squared: 0.7049
F-statistic: 467.9 on 2 and 389 DF, p-value: < 2.2e-16
2.5预测
FitMPG<-predict(Fit,CarData,type="response") ##预测
plot(CarData$weight,CarData$MPG,pch=1,xlab="自重",ylab="MPG")
points(CarData$weight,FitMPG,pch=10,col=2)
legend("topright",c("实际值","拟合值"),pch=c(1,10),col=c(1,2))
3.回归诊断
3.1误差项是否符合高斯-马尔科夫定理
3.1.1残差与误差
模型中的误差项是那些与解释变量无关的因素导致的偏差，误差无法衡量，所以用残差来进行研究
3.1.2高斯-马尔科夫定理
- 独立性：第i个观测的误差与第j个误差不相关，在截面数据中容易满足
- 正态性：第二与第三个可从残差图形中直观判断。
- 等方差性
1.独立性检验
library("car")
>durbinWatsonTest(Fit) ##独立性检验
lag Autocorrelation D-W Statistic p-value
1 0.5575821 0.880059 0
Alternative hypothesis: rho != 0
2.正态性与等差性
############拟合值与残差值 fitted(Fit) residuals(Fit) #######绘制残差图 par(mfrow=c(2,2)) plot(Fit)
- 考察的是 正态分布假定
- 考察队是残差项等方差假定。本例中，残差随着拟合值的增大而有扩散的趋势，不满足等方差假定。
解决异方差办法：
- 采用加权最小二乘法。普通最小二乘法的离差平方和中个观测变量对平方和的贡献相等，而方差偏大的贡献更大，所以回归线偏向方差大的观测。加权最小二乘的权值估计很重要，若能找到异方差与解释变量的函数形式，就能方便的估计各个观测的权值，否则应采用异方差-稳健t统计量进行显著性检验。
2.对被解释变量做BOX-COX变换后在建立模型，最常见的是取log。
##########重新建立对数-水平的二元回归经验方程 Fit<-lm(log(MPG)~weight+horsepower,data=CarData) summary(Fit) par(mfrow=c(2,2)) plot(Fit)
- 对数-水平模型，x变动一单位引起的y变化的百分比semi-elasticity
- 图中残差的到了很好的改善，基本满足了正态分布与等差性要求
- 进一步思考：如果所建立的模型没有一楼对被解释变量有重要影响的解释变量，那么残差不应与拟合值有关系
不满足正态性的处理
>summary(powerTransform(CarData$MPG)) #非正态性处理
bcPower Transformation to Normality
Est Power Rounded Pwr Wald Lwr bnd Wald Upr Bnd
CarData$MPG 0.1974 0 -0.0907 0.4854
Likelihood ratio tests about transformation parameters
LRT df pval
LR test, lambda = (0) 1.809026 1 1.786251e-01
LR test, lambda = (1) 29.356835 1 6.020380e-08
等方差性的处理
>spreadLevelPlot(Fit) #等方差性假定图形 Suggested power transformation: 0.4920781 >ncvTest(Fit) #等方差性假定检验 Non-constant Variance Score Test Variance formula: ~ fitted.values Chisquare = 3.720084 Df = 1 p = 0.05376168 >crPlots(Fit) ####被解释变量和解释变量的线性相关性判断
解释变量与被解释变量的线性关系
crPlots(Fit)
线性关系较明显，基本满足高斯条件。
3.2诊断数据中的异常点
3.2.1高杠杆值点
在解释变量方向上取异常值，在被解释变量上正常的点。
hatvalues(回归分析结果对象名)
LeveragePlot<-function(fit){
Np<-length(coefficients(fit))-1
N<-length(fitted(fit))
plot(hatvalues(fit),main="观测点的杠杆值序列图",ylab="杠杆值",xlab="观测编号")
abline(2*(Np+1)/N,0,col="red",lty=2)
abline(3*(Np+1)/N,0,col="red",lty=2)
identify(1:N,hatvalues(fit),names(hatvalues(fit)))
}
LeveragePlot(Fit)
3.2.2离群点
rstudent()#学生化残差 outlierTest()#离群点
###############探测离群点
library("car")
Fit<-lm(log(MPG)~weight+horsepower,data=CarData)
rstudent(Fit)
outlierTest(Fit)
Fit<-lm(log(MPG)~weight+horsepower,data=CarData[-388,])
outlierTest(Fit)
> outlierTest(Fit)
No Studentized residuals with Bonferonni p < 0.05
Largest |rstudent|:
rstudent unadjusted p-value Bonferonni p
388 3.224743 0.0013677 0.53613
> Fit<-lm(log(MPG)~weight+horsepower,data=CarData[-388,])
> outlierTest(Fit)
No Studentized residuals with Bonferonni p < 0.05
Largest |rstudent|:
rstudent unadjusted p-value Bonferonni p
112 -2.999241 0.0028815 NA
第一个测出388是离群点，p-value=0.54,不能拒绝原假设，所以Fit<-lm(log(MPG)~weight+horsepower,data=CarData[-388,])
3.2.3强影响点
包含或者剔除改点会使得回归方程的截距或者斜率有较大变化。
cooks.distance()
利用图形展示各观测点的库克距离
#######################探测强影响点 Fit<-lm(log(MPG)~weight+horsepower,data=CarData) par(mfrow=c(2,1)) plot(cooks.distance(Fit),main="Cook's distance",cex=0.5) #获得Cook距离
#用identify函数做互动可以获得高库克距离的观测点 Np<-length(coefficients(Fit))-1#系数长度 N<-length(fitted(Fit))#样本量 CutLevel<-4/(N-Np-1)#公式，把判断标准画出来
plot(Fit,which=4)#自动绘制6幅图，第四幅更显著，就只显示第四幅 abline(CutLevel,0,lty=2,col="red")
、
分别查看weight与horsepower对斜率的影响
library("car")
avPlots(Fit,ask=FALSE,onepage=TRUE,id.method="identify")
##############异常观测点的综合展示 Fit<-lm(log(MPG)~weight+horsepower,data=CarData) influencePlot(Fit,id.method="identify",main="异常观测点的可视化")
3.3多重共线性的诊断
1.容忍度
2.方差膨胀因子
> library("car")
> vif(Fit)
weight horsepower
3.959228 3.959228
有一定的共线性，但是不强
4.回归建模策略
1.拟合优度最高
r-square越接近1，拟合优度越高
Fit1<-lm(log(MPG)~weight+horsepower,data=CarData) Fit2<-lm(log(MPG)~weight+horsepower+displacement,data=CarData) summary(Fit1) summary(Fit2) anova(Fit1,Fit2) Analysis of Variance Table Model 1: log(MPG) ~ weight + horsepower Model 2: log(MPG) ~ weight + horsepower + displacement Res.Df RSS Df Sum of Sq F Pr(>F) 1 389 9.5871 2 388 9.5360 1 0.051059 2.0775 0.1503
增加变量会使误差平方和减小：9.58-9.53
这里的F指的是偏F统计量
2.带惩罚的拟合优度最高点模型。
AIC&BIC
AIC=－2ln(模型中极大似然函数值)+2(模型中未知参数个数)
SBC=－2ln(模型中极大似然函数值)+ln(n)(模型中未知参数个数)
越小越好
> Fit1<-lm(log(MPG)~weight+horsepower,data=CarData)
> Fit2<-lm(log(MPG)~weight+horsepower+displacement,data=CarData)
> AIC(Fit1,Fit2)
df AIC
Fit1 4 -334.2030
Fit2 5 -334.2963
> BIC(Fit1,Fit2)
df BIC
Fit1 4 -318.3179
Fit2 5 -314.4399
选第一个模型：BIC更小
3.解释变量的筛选
- 向前筛选
- 向后筛选
- 逐步筛选
#######不同解释变量筛选策略下的AIC
Fit<-lm(log(MPG)~weight+horsepower+displacement+cylinders+acceleration,data=CarData)
library("MASS")
stepAIC(Fit,direction="backward")
Start: AIC=-1449.07
log(MPG) ~ weight + horsepower + displacement + cylinders + acceleration
Df Sum of Sq RSS AIC
- displacement 1 0.00007 9.4316 -1451.1#剔除会让AIC减少的最多
- acceleration 1 0.03725 9.4687 -1449.5
<none> 9.4315 -1449.1
- cylinders 1 0.07211 9.5036 -1448.1
- horsepower 1 0.49772 9.9292 -1430.9
- weight 1 1.09267 10.5242 -1408.1
Step: AIC=-1451.06
log(MPG) ~ weight + horsepower + cylinders + acceleration
Df Sum of Sq RSS AIC
- acceleration 1 0.03733 9.4689 -1451.5
<none> 9.4316 -1451.1
- cylinders 1 0.13744 9.5690 -1447.4
- horsepower 1 0.53766 9.9692 -1431.3
- weight 1 1.35833 10.7899 -1400.3
Step: AIC=-1451.51
log(MPG) ~ weight + horsepower + cylinders
Df Sum of Sq RSS AIC
<none> 9.4689 -1451.5
- cylinders 1 0.11820 9.5871 -1448.7
- horsepower 1 0.68621 10.1551 -1426.1
- weight 1 2.06375 11.5326 -1376.2
Call:
lm(formula = log(MPG) ~ weight + horsepower + cylinders, data = CarData)
Coefficients:
(Intercept) weight horsepower cylinders
4.1169705 -0.0002178 -0.0022721 -0.0242519
全子集回归模型
