R语言做线性回归

1.回归的多面性

回归类型	用途
简单线性	个量化的解释变量来预测一个量化的响应变量（一个因变量、一个自变量）
多项式	一个量化的解释变量预测一个量化的响应变量，模型的关系是 n阶多项式（一个预测变量，但同时包含变量的幂）
多元线性	用两个或多个量化的解释变量预测一个量化的响应变量（不止一个预测变量）
多变量	用一个或多个解释变量预测多个响应变量
Logistic	用一个或多个解释变量预测一个类别型变量
泊松	用一个或多个解释变量预测一个代表频数的响应变量
Cox比例风险	用一个或多个解释变量预测一个事件（死亡、失败或旧病复发）发生的时间 时间序列对误差项相关的时间序列数据建模
非线性	用一个或多个量化的解释变量预测一个量化的响应变量，不过模型是非线性的
非参数	用一个或多个量化的解释变量预测一个量化的响应变量，模型的形式源 自数据形式，不事先设定
稳健	用一个或多个量化的解释变量预测一个量化的响应变量，能抵御强影响点的干扰

2.OLS回归

OLS回归是通过预测变量的加权和来预测量化的因变量，其中权重是通过数据估计而得以的参数。

使残差平方和最小

 

为能够恰当地解释OLS模型的系数，数据必须满足以下统计假设：

（1）      正态性对于固定的自变量，因变量值成正态分布

（2）      独立性 Yi值之间相互独立

（3）      线性   因变量与自变量之间为线性相关

（4）      同方差性因变量的方差不随自变量的水平不同而变化，即不变方差或同方差性

3. 用lm（）拟合回归模型

拟合线性模型最基本的函数就是lm（），格式为：

myfit<-lm(formula,data)

formula指要拟合的模型形式，data是一个数据框，包含了用于拟合模型的数据

formula形式如下：Y~X1+X2+……+Xk （~左边为响应变量，右边为各个预测变量，预测变量之间用+符号分隔）

R表达式中常用的符号
符号	用途
~	分隔符号，左边为响应变量，右边为解释变量，eg：要通过x、z和w预测y，代码为y~x+z+w
+	分隔预测变量
：	表示预测变量的交互项  eg：要通过x、z及x与z的交互项预测y，代码为y~x+z+x:z
*	表示所有可能交互项的简洁方式，代码y~x*z*w可展开为y~x+z+w+x:z+x:w+z:w+x:z:w
^	表示交互项达到某个次数，代码y~(x+z+w)^2可展开为y~x+z+w+x:z+x:w+z:w
.	表示包含除因变量外的所有变量，eg：若一个数据框包含变量x、y、z和w，代码y~.可展开为y~x+z+w
-	减号，表示从等式中移除某个变量，eg：y~(x+z+w)^2-x:w可展开为y~x+z+w+x:z+z:w
-1	删除截距项，eg：表示y~x-1拟合y在x上的回归，并强制直线通过原点
I（）	从算术的角度来解释括号中的元素。Eg：y~x+(z+w)^2将展开为y~x+z+w+z:w。相反，代码y~x+I((z+w)^2)将展开为y~x+h，h是一个由z和w的平方和创建的新变量
function	可以在表达式中用的数学函数，例如log(y)~x+z+w表示通过x、z和w来预测log(y)

对拟合线性模型非常有用的其他函数
函数	用途
Summary（）	展示拟合的详细结果
Coefficients（）	列出拟合模型的模型参数（截距项和斜率）
Cofint（）	提供模型参数的置信区间（默认95%）
Fitted（）	列出拟合模型的预测值
Residuals（）	列出拟合模型的残差值
Anova（）	生成一个拟合模型的方差分析，或者比较两个或更多拟合模型的方差分析表
Vcov（）	列出模型参数的协方差矩阵
AIC（）	输出赤池信息统计量
Plot（）	生成评价拟合模型的诊断图
Predict（）	用拟合模型对新的数据集预测响应变量值

4. 简单线性回归

eg：

[plain]  view plain  copy

1. fit<-lm(weight~height,data=women)  
2. summary(fit)  

在Pr(>|t|)栏，可以看到回归系数（3.45）显著不为0（p<0.001），表明身高每增加1英寸，体重将预期地增加3.45磅

R平方项（0.991）表明模型可以解释体重99.1%的方差，它也是实际和预测值之间的相关系数（R^2=r^2）

残差的标准误（1.53lbs）则可认为模型用身高预测体重的平均误差

F统计量检验所有的预测变量预测响应变量是否都在某个几率水平之上

[plain]  view plain  copy

1. fitted(fit)#拟合模型的预测值  

[plain]  view plain  copy

1. residuals(fit)#拟合模型的残差值  

[plain]  view plain  copy

1. plot(women$height,women$weight,  
2.      xlab="Height （in inches）",  
3.      ylab="Weight（in pounds）")  
4. abline(fit)  

5. 多项式回归

[plain]  view plain  copy

1. fit2<-lm(weight~height+I(height^2),data=women)  
2. summary(fit2)  

[plain]  view plain  copy

1. plot(women$height,women$weight,  
2.      xlab="Height（in inches）",  
3.      ylab="Weight（in lbs）")  
4. lines(women$height,fitted(fit2))  

一般来说，n次多项式生成一个n-1个弯曲的曲线

car包中的scatterplot（）函数，可以很容易、方便地绘制二元关系图

[plain]  view plain  copy

1. scatterplot(weight~height,  
2.             data=women,  
3.             spread=FALSE,  
4.             lty.smooth=2,  
5.             pch=19,  
6.             main="Women Age 30-39",  
7.             xlab="Height (inches)",  
8.             ylab="Weight(lbs.)")  

6.多元线性回归

采用的数据集：state.x77

[plain]  view plain  copy

1. states<-as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])  

检测二变量关系

[plain]  view plain  copy

1. cor(states)  

[plain]  view plain  copy

1. library(car)  
2. scatterplotMatrix(states,spread=FALSE,lty.smooth=2,main="Scatter Plot Matrix")  

scatterplotMatrix（）函数默认在非对角线区域绘制变量间的散点图，并添加平滑（loess）和线性拟合曲线

多元线性回归

[plain]  view plain  copy

1. fit<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)  
2. summary(fit)  

7.有交互项的多元线性回归

[plain]  view plain  copy

1. fit<-lm(mpg~hp+wt+hp:wt,data=mtcars)  
2. summary(fit)  

通过effects包中的effect（）函数，可以用图形展示交互项的结果

[plain]  view plain  copy

1. fit<-lm(mpg~hp+wt+hp:wt,data=mtcars)  
2. summary(fit)  
3.   
4. install.packages("effects")  
5. library(effects)  
6. plot(effect("hp:wt",fit,  
7.             list(wt=c(2.2,3.2,4.2))),multiline=TRUE)  

8.回归诊断

（1）标准方法

[plain]  view plain  copy

1. 简单线性回归  

[plain]  view plain  copy

1. fit<-lm(weight~height,data=women)  
2. par(mfrow=c(2,2))  
3. plot(fit)  

正态性：当预测变量值固定时，因变量成正态颁，则残差图也应是一个均值为0的正态颁。正态Q-Q图是在正态颁对应的值上，标准化残差的概率图，若满足正态假设，则图上的点应该落在吓45度角的直线上，若不是，则违反了正态性假设。

独立性：只能从收集的数据中来验证。

线性：若因变量与自变量线性相关，则残差值与预测（拟合）值就没有任务系统关联，若存在关系，则说明可能城要对回归模型进行调整。

同方差性：若满足不变方差假设，则在位置尺度图（Scale-Location Graph）中，水平线周围的点应随机分布。

二次拟合诊断图

[plain]  view plain  copy

1. fit2<-lm(weight~height+I(height^2),data=women)  
2. par(mfrow=c(2,2))  
3. plot(fit2)  

（2）改进的方法

（car包中的）回归诊断实用函数 

函数	目的
qqPlot（）	分位数比较图
durbinWatsonTest（）	对误差自相关性做Durbin-Watson检验
crPlots（）	成分与残差图
ncvTest（）	对非恒定的误差方差做得分检验
spreadLevelPlot（）	分散水平检验
outlierTest（）	Bonferroni离群点检验
avPlots（）	添加的变量图形
inluencePlot（）	回归影响图
scatterplot（）	增强的散点图
scatterplotMatrix（）	增强的散点图矩阵
vif（）	方差膨胀因子

另gvlma包提供了对所有线性模型进行检验的方法

正态性：

与 plot（）函数相比，qqplot（）函数提供了更为精确的正态假设检验方法，画出了n-p-1个自由度的t分布下的学生化残差图形，n为样本大小，p是回归参数的数目（包括截距项）

eg：

[plain]  view plain  copy

1. library(car)  
2. fit<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)  
3. qqPlot(fit,labels=row.names(states),id.method="identify",simulate=TRUE,main="Q-Q Plot")  

绘制学生残差图的函数

[plain]  view plain  copy

1. residplot<-function(fit,nbreaks=10){  
2.   z<-rstudent(fit)  
3.   hist(z,breaks=nbreaks,freq=FALSE,  
4.        xlab="Studnetized Residual",  
5.        main="Distribution of Errors")  
6.   rug(jitter(z),col="brown")  
7.   curve(dnorm(x,mean=mean(z),sd=sd(z)),  
8.         add=TRUE,col="blue",lwd=2)  
9.   lines(density(z)$x,density(z)$y,  
10.         col="red",lwd=2,lty=2)  
11.   legend("topright",  
12.          legend=c("Normal Curve","Kernel Density Curve"),  
13.          lty=1:2,col=c("blue","red"),cex=0.7)}  
14. residplot(fit)  

误差的独立性：

之前提到可依据收集数据判断因变量是否独立

car包中提供了一个可做Durbin-Watson检验的函数，可检测误差的序列相关性

[plain]  view plain  copy

1. durbinWatsonTest(fit)  

线性：

可通过成分残差图即偏残差图，判断因变量与自变量之间是否呈非线性关系，也可以看是否不同于已设定线性模型的系统偏差，图形可用car包中crPlots（）函数绘制

[plain]  view plain  copy

1. library(car)  
2. crPlots(fit)  

若图形存在非线性，则说明可能对预测变量的函数形式建模不够充分

car包提供了两个有用的函数，可判断误差方差是否恒定

ncvTest（）函数生成一个计分检验，零假设为误差方差不变

spreadLevelPlot（）函数创建一个添加了最佳拟合曲线的散点图，展示标准化残差绝对值与拟合值的关系

检验同方差性：

[plain]  view plain  copy

1. library(car)  
2. ncvTest(fit)  
3. spreadLevelPlot(fit)  

（3）线性模型假设的综合验证

gvlma包中的gvlma（）函数

[plain]  view plain  copy

1. install.packages("gvlma")  
2. library(gvlma)  
3. gvmodel<-gvlma(fit)  
4. summary(gvmodel)  

（4）多重共线性

VIF（Variance Inflation Factor，方差膨胀因子）进行检测

一般原则下，(VIF)^1/2 >2表明存在多重共线性问题

[plain]  view plain  copy

1. library(car)  
2. vif(fit)  
3. sqrt(vif(fit))>2  

9.异常观测值

（1）离群点

离群点指那些模型预测效果不佳的观测点，通常有很大的、或正或负的残差，正残差说明模型低估了响应值，负残差说明高佑了响应值

[plain]  view plain  copy

1. library(car)  
2. outlierTest(fit)  

outlierTest（）函数是根据单个最大（或正或负）残差值的显著性来判断是否有离群点，若不显著，则说明数据集中没有离群点，若显著，则必须删除该离群点，然后再检验是否还有其他离群点存在。

（2）高杠杆值点

高杠杆值观测点，即是与其他预测变量有关的离群点，即它们是由许多异常的预测变量组合起来的，与响应变量值没有关系。

高杠杆值的观测点可通过帽子统计量（hat statistic）判断。对于一个给定的数据集，帽子均值为p/n，其中p是模型估计的参数数目（包含截距项），n是样本量。一般来说，若观测点的帽子值大于帽子均值的2或3倍，则可认定为高杠杆值点。

[plain]  view plain  copy

1. hat.plot<-function(fit){  
2.   p<-length(coefficients(fit))  
3.   n<-length(fitted(fit))  
4.   plot(hatvalues(fit),main="Index Plot of Hat Values")  
5.   abline(h=c(2,3)*p/n,col="red",lty=2)  
6.   identify(1:n,hatvalues(fit),names(hatvalues(fit)))  
7. }  
8. hat.plot(fit)  

（3）强影响点

强影响点，即对模型参数估计值影响有些比例失衡的点。例如，当移除 模型的一个观测点时模型会发生巨大的改变，那么需要检测一下数据中是否存在强影响点。

检测方法

Cook距离，或称为D统计量     Cook's D值大于4/(n-k-1)，则表明它是强影响点，其中n为样本量大小，k是预测变量数目（有助于鉴别强影响点，但并不提供关于这些点如何影响模型的信息）

变量添加图（added variable plot）（弥补了该缺陷）（对于每个预测变量Xk，绘制Xk在其他k-1个预测变量上回归的残差值相对于响应变量在其他k-1个预测变量上回归的残差值的关系图）

[plain]  view plain  copy

1. cutoff<-4/(nrow(states)-length(fit$coefficients)-2)  
2. plot(fit,which=4,cook.levels=cutoff)  
3. abline(h=cutoff,lty=2,col="red")  

[plain]  view plain  copy

1. library(car)  
2. avPlots(fit,ask=FALSE,onepage=TRUE,id.method="identify")  

car包中的influencePlot（）函数，可将离群点、杠杆点和强影响点的信息整合到一幅图形中

[plain]  view plain  copy

1. library(car)  
2. influencePlot(fit,id.method="identify",main="Influence Plot",  
3.               sub="Circle size if proportional to Cook's distance")  

影响图。纵坐标超过2或小于-2的州可被认为是离群点，水平轴超过0.2或0.3的州有高杠杆值（通常为预测值的组合）。圆圈大小与影响成比例，圆圈很大的点可能是对模型估计造成的不成比例影响的强影响点。

10.改进的措施

（1）删除观测点

删除观测点可提高数据集对于 正态假设的拟合度，而强影响点会干扰结果，通常也会被删除。删除最大的离群点或强影响点，模型需要重新拟合，若离群点或强影响点仍然存在，重复以上过程直到获得比较满意的拟合。

对删除观测点应持谨慎态度。

（2）变量变换 

当模型不符合正态性、线性或同方差性假设时，一个或多个变量的变换通常可以改善或调整模型效果。

当模型违反了正态假设时，通常可以对响应变量尝试某种变换。

car包中的powerTransform（）函数

Box-Cox正态变换

[plain]  view plain  copy

1. library(car)  
2. summary(powerTransform(states$Murder))  

（3）增删变量

改变模型的变量会影响模型的拟合度，增加或删除变量

多重共线问题：岭回归

11.选择“最佳”的回归模型

（1）模型比较

anova（）函数可比较两个嵌套模型的拟合优度

嵌套模型即指它的一个些项完全饮食在另一个模型中

用anova（）函数比较

[plain]  view plain  copy

1. fit1<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)  
2. fit2<-lm(Murder~Population+Illiteracy,data=states)  
3. anova(fit2,fit1)  

模型1嵌套在模型2中，检验不显著，基础知识 不需要将Income和Frost添加到线性模型中，可将它们从模型中删除 

AIC（Akaike Information Criterion，赤池信息准则）可用来比较模型，考虑了模型的统计拟合度及用来拟合的参数数目 

AIC值越小的模型要优行选择，说明模型用较少的参数获得了足够的拟合度

[plain]  view plain  copy

1. fit1<-lm(Murder~Population+Illiteracy+Income+Frost,  
2.          data=states)  
3. fit2<-lm(Murder~Population+Illiteracy,data=states)  
4. AIC(fit1,fit2)  

（2）变量选择

逐步回归法（stepwise method）：

向前逐步回归（forward stepwise）每次添加一个预测变量到模型中，直到添加变量不会使模型有所改进为止。

向后逐步回归（backward stepwise）从模型包含所有预测变量开始，一次删除一个变量直到会降低模型质量为止。

向前向后逐步回归（stepwise stepwise 逐步回归）

MASS包中的steAIC（）函数可实现逐步回归模型，依据的是精确AIC准则

[plain]  view plain  copy

1. 后向回归  
2. library(MASS)  
3. fit1<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)  
4. stepAIC(fit,direction="backward")  

全子集回归（all-subsets regression）：

全子集回归，即所有可能的酣篮队支被检验，可选择展示所有可能的结果，也可展示n个不同子集大小（一个、两个或多个预测变量）的最佳模型

可用leaps包中的regsubsets（）函数实现

可通过R平方、调整R平方或Mallows Cp统计量等准则来选择“最佳”模型

R平方是预测变量解释响应变量的程度

调整R平方与之类似，但考虑了模型的参数数目 

Mallows Cp统计量也用来作为逐步回归的判停规则，对于一个好的模型，它的Cp统计量非常迫近于模型的参数数目（包括截距项）

[plain]  view plain  copy

1. install.packages("leaps")  
2. library(leaps)  
3. leaps<-regsubsets(Murder~Population+Illiteracy+Income+Frost,data=states,nbest=4)  
4. plot(leaps,scale="adjr2")  

[plain]  view plain  copy

1. library(car)  
2. subsets(leaps,statistic="cp",main="Cp Plot for All Subsets Regression")  
3. abline(1,1,lty=2,col="red")  

12. 深层次分析

（1）交叉验证

交叉验证即将一定比例的数据挑选出来作为训练样本，另外的样本作为保留样本，先在训练样本上获取回归方程，然后在保留样本上做预测。 由于保留样本不涉及模型及参数的选择，该样本可获得比新数据更为精确的估计。

k重交叉难中，样本被分为k个子样本，轮流将k-1个子样本组合作为训练集，另外1个子样本作为保留集，这样会获得k个预测方程，记录k个保留样本的预测表现结果，然后求其平均值。【当n是观测总数目，k为n时，该方法又称作刀切法（jackknifing）】

bootstrap包中的crossval（）函数可实现k重交叉验证

[plain]  view plain  copy

1. install.packages("bootstrap")  
2. library(bootstrap)  
3. shrinkage<-function(fit,k=10){  
4.   require(bootstrap)  
5.   theta.fit<-function(x,y){lsfit(x,y)}  
6.   theta.predict<-function(fit,x){cbind(1,x)%*%fit$coef}  
7.   x<-fit$model[,2:ncol(fit$model)]  
8.   y<-fit$model[,1]  
9.   results<-crossval(x,y,theta.fit,theta.predict,ngroup=k)  
10.   r2<-cor(y,fit$fitted.values)^2  
11.   r2cv<-cor(y,results$cv.fit)^2  
12.   cat("Original R-square=",r2,"\n")  
13.   cat(k,"Fold Cross-Validated R-square=",r2cv,"\n")  
14.   cat("Change=",r2-r2cv,"\n")  
15. }  
16. fit<-lm(Murder~Population+Income+Illiteracy+Frost,data=states)  
17. shrinkage(fit)  
18. fit2<-lm(Murder~Population+Illiteracy,data=states)  
19. shrinkage(fit2)  

（2）相对重要性

[plain]  view plain  copy

1. zstates<-as.data.frame(scale(states))  
2. zfit<-lm(Murder~Population+Income+Illiteracy+Frost,data=zstates)  
3. coef(zfit)  

相对权重：是对所有可能子模型添加一个预测变量引起的R平方均增加量的一个近似值。

[plain]  view plain  copy

1. relweights<-function(fit,……){  
2.   R<-cor(fit$model)  
3.   nvar<-ncol(R)  
4.   rxx<-R[2:nvar,2:nvar]  
5.   rxy<-R[2:nvar,1]  
6.   svd<-eigen(rxx)  
7.   evec<-svd$vectors  
8.   ev<-svd$values  
9.   delta<-diag(sqrt(ev))  
10.   lambda<-evec%*%delta%*%t(evec)  
11.   lambdasq<-lambda^2  
12.   beta<-solve(lambda)%*%rxy  
13.   rsqrare<-colSums(beta^2)  
14.   rawwgt<-lambdasq%*%beta^2  
15.   import<-(rawwgt/rsquare)*100  
16.   lbls<-names(fit$model[2:nvar])  
17.   rownames(import)<-lbls  
18.   colnames(import)<-"Weight"  
19.   barplot(t(import),names.arg=lbls,  
20.           ylab="% of R-Square",  
21.           xlab="Predictor Variables",  
22.           main="Relative Importance of Predictor Variables",  
23.           sub=paste("R-Square=",round(rsquare,digits=3)),……)  
24.   return(import)  
25. }  
26. fit<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)  
27. relweights(fit,col="lightgrey")