R语言（一）

2023-10-19 21:22| 来源: 网络整理| 查看: 265

一、数据

二、简单拟合

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，比之前具有更好的结果。

【本文地址】

公司简介

联系我们