目录
一、数据
首先到网址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,比之前具有更好的结果。