R语言之违背基本假设的几种情况xt4.9

第4章 违背基本假设的几种情况

4.9 表4.11（课本P126）是用电高峰每小时用电量y与每月总用电量x的数据。
（1）用普通最小二乘法建立y与x的回归方程，并画出残差散点图。
（2）诊断该问题是否存在异方差。
（3）如果存在异方差，用幂指数型的权函数建立加权最小二乘回归方程。
（4）用方差稳定变换y’=sqrt(y)消除异方差。
（5）用BOX-COX变换消除异方差。

rm(list=ls())

用户编号=c(1:53)
x=c(679,292,1012,493,582,1156,997,2189,1097,2078,1818,1700,747,2030,1643,414,354,1276,745,435,540,874,1543,1029,710,1434,837,1748,1381,1428,1255,1777,370,2316,1130,463,770,724,808,790,783,406,1242,658,1746,468,1114,413,1787,3560,1495,2221,1526)
y=c(0.79,0.44,0.56,0.79,2.70,3.64,4.73,9.50,5.34,6.85,5.84,5.21,3.25,4.43,3.16,0.50,0.17,1.88,0.77,1.39,0.56,1.56,5.28,0.64,4.00,0.31,4.20,4.88,3.48,7.58,2.63,4.99,0.59,8.19,4.79,0.51,1.74,4.10,3.94,0.96,3.29,0.44,3.24,2.14,5.71,0.64,1.90,0.51,8.33,14.94,5.11,3.85,3.93)
data4.9<-data.frame(用户编号,x,y)
data4.9



# ---- 4.9用电高峰每小时用电量y与每月总用电量x ----
#（1）用普通最小二乘法建立y与x的回归方程，并画出残差散点图。----
data4.9 <- read.csv('D:/rwork/应用回归/习题数据/表4-11.csv',head=TRUE)
attach(data4.9) #将该数据框添加到R的搜索路径中，以便于下面直接调用x和y
lm4.9 <- lm(y~x,data=data4.9) #以y为因变量，x为自变量建立回归方程，并将结果赋给lm4.9
summary(lm4.9) #回归分析
# 求得y^=-0.831+0.0038x
e <- resid(lm4.9) #计算残差
plot(x,e,main='残差散点图')
abline(h=c(0),lty=5) #添加虚直线e=0
# 从残差图看出，误差项具有明显的异方差性，误差随着x的增加呈现出增加的态势。



#（2）诊断该问题是否存在异方差。 ----
abse <- abs(e) #计算残差e的绝对值
cor.test(data4.9$x,abse,alternative='two.sided',method='spearman',conf.level=0.95)
# 等级相关系数rs=0.3175294 ，P值=0.02091，认为残差绝对值|ei|与自变量xi显著相关，存在异方差


rs <- 0.3175294
n <- length(x)
t <- sqrt(n-2*rs)/sqrt(1-rs^2)
t
alpha <- 0.05
qt(1-alpha/2,n-2)
# t值=7.631294＞t0.025(51)=2.007584，说明xi与|ei|之间存在系统关系，异方差性问题存在。



#（3）如果存在异方差，用幂指数型的权函数建立加权最小二乘回归方程。 ----
s <- seq(-2,2,0.5) #生成序列-2.0,-1.5,-1.0,-0.5,0.0,0.5,1.0,1.5,2.0
result1 <- vector(length=9,mode='list')
#生成一个列表向量，以存储下面循环过程中的回归方程估计的对数似然统计量结果
result2 <- vector(length=9,mode='list')
#生成一个列表向量，以存储下面循环过程中所建立回归方程的估计系数及显著性检验等结果
for (j in 1:9){
  w <- data4.9$x^(-s[j]) #计算权向量
  lm4.9_3 <- lm(y~x,weight=w,data4.9) #使用加权最小二乘估计建立回归方程
  result1[[j]] <- logLik(lm4.9_3) #将第j次计算的对数似然统计量保存在result1的第j个元素中
  result2[[j]] <- summary(lm4.9_3) #将第j次建立的回归方程的结果保存在result2的第j个元素中
}
result1
# 由输出结果可知m取第8个值即m=1.5时对数似然函数达到极大，因而幂指数m的最优取值为1.5
result2[8] #输出m=1.5时使用加权最小二乘法得到的回归方程
# 此时R^2=0.6591，F值=98.6，加权后的回归方程为y2^=-0.683+0.004x


# 计算加权后的残差，并对残差绝对值和自变量做等级相关系数分析，
# 得到rs=0.321，P值为0.019<0.05，即加权最小二乘法没有消除异方差，只是消除异方差的不良影响，从而对模型进行一点改进。
# ???这时候怎么计算加权后的残差？？？



#（4）用方差稳定变换消除异方差。 ----
y1_hat <- sqrt(y)
lm4.9_4 <- lm(y1_hat~x,data=data4.9)
summary(lm4.9_4)
# 得到y'^=0.5822e+0.00095ex，R^2=0.6485，F值=94.08

e4 <- resid(lm4.9_4)
abse4 <- abs(e4) #计算残差e的绝对值
cor.test(data4.9$x,abse4,alternative='two.sided',method='spearman',conf.level=0.95)
# 得到rs=-0.1740042，P值=0.2121>0.05，说明异方差已经消除。




# ----（5）用BOX-COX变换消除异方差。 ----
#install.packages('MASS') #安装MASS包
library(MASS) #加载MASS包
bc4.9 <- boxcox(y~x,data=data4.9,lambda=seq(-2,2,0.01))
#λ的取值区间为[-2,2]上步长为0.01的值，bc4.9中保存了λ的值及其对应的对数似然函数值
lambda <- bc4.9$x[which.max(bc4.9$y)] #将使对数似然函数值达到最大的λ赋给lambda
lambda
y_bc <- (data4.9$y^lambda-1)/lambda #计算变换后的y值
lm4.9_bc <- lm(y_bc~x,data=data4.9) #使用变换后的y值建立回归方程
summary(lm4.9_bc)
abse_bc <- abs(resid(lm4.9_bc)) #计算残差的绝对值
cor.test(data4.9$x,abse_bc,method='spearman') #计算残差与x的相关系数

# 似然函数取值最大的λ=0.55
# 将因变量进行BOX-COX变换后的回归结果得，y^(0.55)=-0.89+0.002x
# 将y^(0.55)=(y^(0.55)-1)/0.55带入上式得，y^=(0.5105+0.0011x)^(1/0.55)
# 并求得残差绝对值与x的等级相关系数t检验的P值=0.2935，在显著性水平为0.05时都不显著，
#  故可认为异方差被消除。
# 另经BOX-COX变换后的R^2=0.6578，F值=98.04；加权最小二乘的R^2=0.6591，F值=98.6。
#  这说明用BOX-COX变换和加权最小二乘估计都能消除异方差，
#  但对于本题的数据使用加权最小二乘的拟合效果要略好。


detach(data4.9) #与attach()相对应，将数据框从搜索路径中移除

参考课本：应用回归分析（R语言版），何晓群编著