复杂统计方法R语言——简单回归
简单回归
数据来源：http://www.statsci.org/data/general/cofreewy.html
1.读入数据
setwd("D:/数学建模/寒假美赛集训/R统计")
w=read.table("COfreewy.txt",header=T,encoding = "utf-8")
2.线性回归
a=lm(CO~.,w) #a=lm(CO~Traffic+Wind,w)#线性回归 #print(w) summary(a)
第一列为回归系数，第二列为标准差，一般都<0.05较好，第四列为p值。原假设为：系数为0。p<0.05，则拒绝原假设。若p>0.05，可以考虑去除该变量再做一次回归，以使每一个p值都<0.05。
3.AIC逐步回归
b=step(a,direction = "backward")#逐步回归 summary(b)
通过使用逐步回归，我们选择了变量，舍弃了Hour变量，每一个p值都<0.05。
4.残差检验
shapiro.test(b$res)#做残差的正态性检验
$y-\hat{y}$为残差，应符合正态分布，均值为0。如果均值不为0，则存在系统误差。
5.画出qq图
qqnorm(b$res);qqline(b$res)#查看回归残差的qq图
qq图有两个作用：1、检验一组数据是否服从某一分布。2、检验两个分布是否服从同一分布。qq图全称是quantile-quantile plot，从名称中可以了解到是和分位数相关的图。qq图没有过原点，说明仍存在一定系统误差。
6.各变量散点图
plot(w) #pairs(w)
7.傅里叶技术变换
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/24)*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)
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)#查看回归残差的qq图
所有代码:
setwd("D:/数学建模/寒假美赛集训/R统计")
w=read.table("COfreewy.txt",header=T,encoding = "utf-8")
attach(w)#把数据变成全局变量
print(CO)
head(w)#查看前6条数据
a=lm(CO~.,w)
#a=lm(CO~Traffic+Wind,w)#线性回归
#print(w)
summary(a)
b=step(a,direction = "backward")#逐步回归
summary(b)
shapiro.test(b$res)#做残差的正态性检验
qqnorm(b$res);qqline(b$res)#查看回归残差的qq图
par(mfrow=c(2,3))#同时显示6张图
plot(Traffic,Wind)
plot(CO,Wind)
plot(Traffic,CO)
plot(Hour,Wind)
plot(CO,Hour)
plot(Hour,Traffic)
plot(w)
#pairs(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/24)*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)
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)#对残差的正态性检验
par(mfrow=c(1,1))
qqnorm(b$res);qqline(b$res)#查看回归残差的qq图
