R语言（一）——横截面数据回归

目录

一、数据

二、简单拟合

1.回归

2.残差分析

三、模型修正

1.数据分析

2.回归

---

一、数据

首先到网址http://www.statsci.org/data/general/cofreewy.html获取cofreewy.txt数据

二、简单拟合

1.回归

rm(list=ls())
#setwd("")  #设置路径
w=read.table("cofreewy.txt",header=T)
a=lm(CO~.,w)  #利用全部自变量
summary(a) 

将CO作为因变量，其余作为自变量，得到回归方程，Estimate为各项系数，其中R²约为0.95左右，但是发现因子的P值过大了，进一步进行逐步回归。

这里运用倒向逐步回归（backward），就是不断剔除一个变量，利用其余变量进行回归；接着剔除两个变量....每次计算AIC值，看剔除哪些变量的情况下能达到更好的效果。

b=step(a,direction = "backward") #逐步回归,与lm配合
summary(b)

此时P值都小于0.05，且R²和调整后R²也稍稍提高了一些，看起来具有更好的效果，此时拟合的回归直线是：

2.残差分析

shapiro.test(b$res)  #残差的正态性检验
qqnorm(b$res);qqline(b$res)

残差分析用文字和图像两种方式呈现，可以看到文字中P值>0.05，因此不够合格。而从图像上，合格的残差结果应该具备两点：（1）直线经过原点（2）直线在第一象限呈45°角，显然也不符合。因此需要对数据进一步修正。

三、模型修正

1.数据分析

可以通过plot(w)或pairs(w)将变量间的关系用散点图显示，可以发现Hour和Wind之间似乎呈现正弦曲线关系，可以将其加入考虑。CO和Wind、Traffic的关系比较复杂，很难用线性关系表示，利用变量的二次项和三次项考察对CO的影响

2.回归

attach(w)#变量存入内存
cor(cbind(CO,Traffic,Tsq=Traffic^2,Tcub=Traffic^3,  #相关系数矩阵
          Hour,Hsq=Hour^2,Hcub=Hour^3,
          Wind,Wsq=Wind^2,Wcub=Wind^3))
a=lm(CO~Traffic+Wind+I(Wind^2)+I(Wind^3)+sin((2*pi/25)*Hour)
     +cos((2*pi/24)*Hour)+sin((4*pi/24)*Hour)+cos((4*pi/24)*Hour))
b=step(a) #逐步回归，按照AIC选择变量
summary(b);  #anova(b)
#shapiro.test(b$res)

改变自变量内容后，进行回归分析。观察各项的P值，其中Wind的三次项和Hour的sin项的P值都过大，去掉之后重新拟合，代码为:

b1=lm(CO~Traffic+Wind+I(Wind^2)+
      +cos((2*pi/24)*Hour)+cos((4*pi/24)*Hour))
summary(b1)
anova(b1)  #方差分析表
shapiro.test(b1$res)  #对残差正态性检验
qqnorm(b$res);qqline(b$res)

最终回归结果： 

方差分析表：

最终选择出的模型：

其中在对残差的正态性检验中，P值为0.6396，不能拒绝其来自正太总体。

系数R²和调整的R²分别为0.9913和0.9889，比之前具有更好的结果。