R建模之回归（一）

3种常见的回归模型：

线性回归（预测连续型变量比如婴儿出生体重），逻辑回归（预测二元变量比如过低出生体重与正常出生体重），泊松分布（计数比如每年或每个国家过低出生体重婴儿人数）

我们以gamlss.data包提供的usair数据集进行研究，US空气污染数据集。我们希望预测根据城市面积（以人口规模/千人为统计依据）估计的空气污染程度（这里也就是数据集中的x3），空气污染以每立方米空气中二氧化硫的含量（毫克）为指标（也就是数据集中的y）

首先查看一下数据集的数据结构：

library(gamlss.data)           #导入数据集所在的包
data(usair)
mode1.0<-lm(y~x3,data=usair)   #线性拟合函数
summary(mode1.0)               #获取描述性统计量

上述代码运行结果如下：

拟合的方程x3系数为0.02，常数项系数为：17.87，将回归结果可视化：

plot(y~x3,data=usair,cex.lab=1.5)  #y和x3的点图，并设置坐标轴标签的缩放倍数
abline(mode1.0,col="red",lwd=2.5)  #添加拟合结果直线，并设置颜色和线条宽度
legend('bottomright',legend='y~x3',lty=1,col='red',lwd=2.5,title='Regression line')  #添加图例，第一个参数表式图例的位置，第二个表示内容，第三个线的类型，第四个颜色第五个宽度，第六个图例的标题

图中直线很好地拟合了数据。我们使用最小二乘法来实现最佳拟合，所以该模型也被称为普通最小二乘回归（OLS）

那么上面的直线是如何得到的呢？是通过使残差和最小来找到最佳的拟合曲线，残差指的是观测值（散点图中的原始数据点）与预测理论值（回归直线上对应的值）的差值：

usair$prediction<-predict(mode1.0)    #预测结果保存在prediction变量里面
usair$residual<-resid(mode1.0)        #残差
plot(y~x3,data=usair,cex.lab=1.5)     #画原始数据图像
abline(mode1.0,col='red',lwd=2.5)       #拟合曲线
segments(usair$x3,usair$y,usair$x3,usair$prediction,col='blue',lty=2)       #对线段作图，
legend('bottomright',legend=c('y~x3','residuals'),lty=c(1,2),col=c('red','blue'),lwd=2.5,title='Regression line')  #添加图例

以上就是对于单变量的拟合。另一方面，如果我们希望将人口数与工业化的影响区分考虑，建立更复杂的模型，需要引入变量x2，该变量为工人数大于20的工厂个数。可以有两种方法对模型进行更新：

mode1.1<-lm(y~x3+x2,data=usair)
summary(mode1.1)

mode1.1<-update(mode1.0,.~.,+x2)
summary(mode1.1)

从图中可以看到，x3的系数是-0.05661，而第一个模型中x3的系数为0.02，这说明，当引入x2这一变量后，y和x3的关系变成了负相关，这表示人口数每上升1000，空气中so2浓度下降0.05661单位，该结论具有统计显著效应。当我们固定工业化水平，人口规模和空气污染之间的关联就变为负相关，这意味着城市工业化水平不变，城市规模增加，则污染物扩散的面积也增加了。由此，可以得到结论，变量x2可能是一个混杂因子，它歪曲了y和x3之间的关联关系。

基于上述模型，我们可以预测变量的响应值：比如：预测一个拥有300 000居民和160个工厂，平均每个工厂人数都超过20的城市二氧化硫浓度：

可以看出效果还是不错的。

如果有两个预测变量，那么回归直线就是一个平面了，使用scatterplot3d展示：

library(scatterplot3d)
plot3d<-scatterplot3d(usair$x3,usair$x2,usair$y,pch=19,type='h',       #绘制3d图形
     highlight.3d=TRUE,main='3-D Scatterplot')
 plot3d$plane3d(mode1.1,lty='solid',col='red')

一般是看二位投影比较好：

par(mfrow=c(1,2))                                  #画一行两列的图形
 plot(y~x3,data=usair,main='2D projection for x3')
 abline(mode1.1,col='red',lwd=2.5)                   #添加回归直线
plot(y~x2,data=usair,main='2D projection for x2' )
abline(lm(y~x2+x3,data=usair),col='red',lwd=2.5)      

#添加回归直线

5.3 模型假定 

以下：（完全引用）

  采用标准预测技术的线性回归模型能够对输出变量、预测变量以及两者之间的相关性进行假设：
  1）Y是一连续变量（非二元变量，定类变量或定序变量）
  2）错误（残差）具有统计无关性；
  3）在Y和每个X之间具有随机的线性关联
  4）固定X，Y服从正态分布

  5）不论X的固定值是多少，Y拥有相同的方差

假定2）的特例发生在趋势分析中，如果我们使用时间作为预测变量，由于连续年份是相互依赖的，因此相互之间的错误不会独立无关。

 假定3）的特例是变量之间相互关系不是完全线性相关，而是背离了趋势直线。假定4和假定5要求Y的条件分布服从正态分布，假定5也被称为同方差假定。如果违背了该假定，线性回归模型则存在异方差性。

写了这么多，我自己也不是很理解，但是下面的程序出来就能明白一些了：

library(Hmisc)
library(ggplot2)
library(gridExtra)
set.seed(7) 
x <-sort(rnorm(1000,10,100))[26:975]
y <-x*500+rnorm(950,5000,20000)
df <-data.frame(x=x,y=y,cuts=factor(cut2(x,g=5)),resid=resid(lm(y~x)))
scatterP1<-ggplot(df,aes(x=x,y=y))+
           geom_point(aes(colour=cuts,fill=cuts),shape=1,show.legend=FALSE)+
           geom_smooth(method=lm,level=0.99)
plot_left<-ggplot(df,aes(x=y,fill=cuts))+
           geom_density(alpha=.5)+coord_flip()+scale_y_reverse()
plot_right<- ggplot(data=df,aes(x=resid,fill=cuts))+
           geom_density(alpha=.5)+coord_flip()
grid.arrange(plot_left,scatterP1,plot_right,ncol=3,nrow=1,widths=c(1,3,1))

下面介绍去掉孤立点的方法，去掉孤立点后有可能使得假定成立，可以使用gvlma包来快速验证上面五个模型假定是否都能满足，

library(gvlma)
gvlma(mode1.1)

五个假定只有三个能满足。当我们把第31个观测点去掉，再对R方进行分析：

从结果可以看出去掉31个观测值，拟合效果变得更好。

下面简要介绍评价回归线的拟合效果（AIC）

      尽管我们知道结果所得的趋势直线在所有可能的线性直线中拟合效果最佳，但是我们并不知道其对于实际数据的拟合效果。这个时候可以通过验证零假设来求得回归参数的显著性，零假设假定回归参数等于零。

F检验适用于回归参数确实为零的情况。p检验适用于一般情况的回归显著性检验，p值小于0.05即说明回归直线具有显著性。

除了显著的F值，残差是对拟合误差的度量。R平方值将以上这些值统一为一个度量。

R平方值是指：响应变量的变异中有多少百分比是有回归模型解释的；数学表示是预测Y值得方差除以观测Y值的方差。

mode1.2明显比mod1.0R平方大，同种第一个是R平方值，第二个是R平方调整值，但是R平方调整值不能作为非嵌套模型筛选的标准。如果用户的模型为非嵌套模型，可以考虑使用赤池信息量准则（AIC），优先考虑AIC值小的那一个，如果两个模型的差小于2，这两个模型差别不大：

summary(mode1.3 <-update(mode1.2,.~.-x2+x1))$coefficients
summary(mode1.4 <-update(mode1.2,.~.-x3+x1))$coefficients
AIC(mode1.3,mode1.4)       #AIC准则，赤池信息量准则