## 6.4节例题


################## 例6.4.1 ##################
data <- read.table("/Users/macbook/Desktop/线性模型例题/Chapter6/data/例6.2.3.txt")
data1 <- as.data.frame(data)[,-1]
colnames(data1) <- c("X1","X2","X3","X4","Y")

lm.reg <- lm(Y ~ X1 + X2, data = data1) ## 拟合BIC选模型 
par(ann = FALSE, mfrow = c(1, 2), mar = c(5, 6, 4, 2))
plot(lm.reg,pch = 16, 1)
title(
  xlab = expression(italic(hat(y)[i])),  # 斜体显示 \hat{y}_i
  ylab = expression(italic(hat(e)[i])),  # 斜体显示 \hat{e}_i
  cex.lab = 1.1         # 标签字体大小
)
plot(lm.reg,pch = 16, 2)
title(
  xlab = expression(italic(q[(i)])),  
  ylab = expression(italic(r[(i)])),  
  cex.lab = 1.1         # 标签字体大小
)
y.res <- residuals(lm.reg)
shapiro.test(y.res)  






################## 例6.4.2 ##################
data(mtcars)
model<-lm(mpg~disp+hp, data=mtcars)
summary(model) 
par(mfrow=c(1,1))
plot(model, 1, pch=19, lwd=3)  ## Figure6.4.4(a)
title(
  xlab = expression(italic(hat(y)[i])),  # 斜体显示 \hat{y}_i
  ylab = expression(italic(hat(e)[i])),  # 斜体显示 \hat{e}_i
  cex.lab = 1.1         # 标签字体大小
)
plot(mtcars$disp, model$residual, pch=19, lwd=3)  ## Figure6.4.4(b)
title(
  xlab = expression(italic(x[i1])),  # 斜体显示 \hat{y}_i
  ylab = expression(italic(hat(e)[i])),  # 斜体显示 \hat{e}_i
  cex.lab = 1.1         # 标签字体大小
)
plot(mtcars$hp, model$residual, pch=19, lwd=3)  ## Figure6.4.4(c)
title(
  xlab = expression(italic(x[i2])),  # 斜体显示 \hat{y}_i
  ylab = expression(italic(hat(e)[i])),  # 斜体显示 \hat{e}_i
  cex.lab = 1.1         # 标签字体大小
)

## Spearman秩相关检验
e <- resid(model) #获取残差
abse <- abs(e)
cor.test(mtcars$disp, abse, method="spearman" )
cor.test(mtcars$hp, abse, method="spearman" )

library(lmtest) 
gqtest(model, order.by = ~disp+hp, data = mtcars, fraction = 7) # Goldfeld Quandt test
bptest(model)  # Breusch-Pagan test
bptest(model, ~ disp*hp + I(disp^2) + I(hp^2), data = mtcars) # White检验





################## 例6.4.3 ##################
library(MASS)
data<-read.table("/Users/macbook/Desktop/线性模型例题/Chapter6/data/例6.4.3.txt",header=T)
model1<-lm(Y~X,data=data)  #模型1
summary(model1) 

par(ann=FALSE, mfrow=c(2,2)) ##Fig-645
plot(data$X,data$Y, pch=19 )
title(
  xlab = expression(italic(X)),  
  ylab = expression(italic(Y)),  
  cex.lab = 1.1,         # 标签字体大小
  line = 2.5
)
plot(model1,1, pch=19 )
title(
  xlab = expression(italic(hat(y)[i])),  
  ylab = expression(italic(hat(e)[i])),  
  cex.lab = 1.1,        # 标签字体大小
  line = 2.5
)
plot(model1,2, pch=19 )
title(
  xlab = "理论分位数",  
  ylab = "标准化残差",  
  cex.lab = 1.1,         # 标签字体大小
  line = 2.5
)
#plot(predict(model1), model1$residual,   xlab="fitted value",ylab="Residuals", pch=19, lwd=3)
plot(data$X, model1$residual, xlab="X",ylab="Residuals", pch=19 )
title(
  xlab = expression(italic(X)),  
  ylab = expression(italic(hat(e)[i])),  
  cex.lab = 1.1,         # 标签字体大小
  line = 2.5
)

X2<-(data$X)^2
model3<-lm(log(Y)~X+X2,data=data)  # 模型3
summary(model3) 
par(mfrow=c(1,2))
plot(model3, 1, pch = 19, lwd = 3)  
title(
  xlab = expression(paste(ln(italic(Y)), "的拟合值")),  # 正确组合表达式与中文
  ylab = "普通残差",  
  cex.lab = 1.1,        # 标签字体大小
  line = 2.5)
plot(model3, 2, pch = 19, lwd = 3)
title(
  xlab = expression(paste('理论分位数', (italic(log(Y)~X+X^2)))),  # 修正逗号分隔和括号闭合问题
  ylab = '标准化残差',  
  cex.lab = 1.1,        # 标签字体大小
  line = 2.5)

shapiro.test(model3$residual)






################## 例6.4.4 ##################
library(MASS)
data<-read.table("/Users/macbook/Desktop/线性模型例题/Chapter6/data/例6.4.4.txt",header=T)
model<-lm(Y~X,data=data)   
summary(model)
par(ann=FALSE)
plot(model,1, pch=19, lwd=3) 
title(
  xlab = expression(italic(hat(y)[i])),  # 斜体显示 \hat{y}_i
  ylab = expression(italic(hat(e)[i])),  # 斜体显示 \hat{e}_i
  cex.lab = 1.1,         # 标签字体大小
  line = 2.6
)
plot(model,2, pch=19, lwd=3)
title(
  xlab = expression(italic(q[(i)])),  
  ylab = expression(italic(r[(i)])),  
  cex.lab = 1.1,         # 标签字体大小
  line = 2.6
)
shapiro.test(residuals(model))


## Box-Cox变换
library(MASS) 
# 寻找最优lambda值
# 定义函数：计算给定 lambda下Box-Cox变换后模型的SSE（残差平方和）
boxcox.sse <- function(lambda, model) {
  # 从模型对象中提取自变量和因变量
  mf <- model.frame(model)
  y <- model.response(mf)
  x <- model.matrix(model)
  K2 <- prod(y) ^ (1 / length(y))  
  K1 <- 1 / (lambda * K2 ^ (lambda - 1))  
  if (lambda != 0) {
    W <- K1 * (y ^ lambda - 1)  
  } else {
    W <- K2 * log(y)  
  }
  new_model <- lm(W ~ x - 1) # - 1 表示移除截距项
  return(deviance(new_model))  
}

# 设定要尝试的 lambda 取值
lambda <- c(-2, -1, -0.5, 0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 1, 2)  
# 使用 sapply 遍历每个 lambda 值，调用 boxcox.sse 函数计算对应的 SSE
SSE <- sapply(lambda, boxcox.sse, model)  
# 将 lambda 取值和对应的 SSE 整理成矩阵
result_matrix <- t(cbind('lambda' = lambda, 'SSE' = round(SSE, 2)))  
# 输出结果矩阵：表6.4.3
result_matrix  

## 绘制Box-Cox变换后的SSE随lambda的变化图:图6.4.8
par(ann=TRUE, mfrow = c(1, 1))
lambda <- seq(-1, 1.2, by = 0.1)  # 生成 lambda 序列
# 计算每个 lambda 对应的 SSE
SSE <- sapply(lambda, boxcox.sse, model)  # 直接传入完整模型
# 绘制 lambda 与 SSE 的关系图
plot(lambda, SSE, type = "l", lwd = 2, 
     xlab = expression(lambda), ylab ="SSE",
     col = "blue")
# 添加垂直参考线
best_lambda <- lambda[which.min(SSE)]  # 自动找到最优 lambda
abline(v = best_lambda, lty =3 , lwd = 2, col = "red")

## boxcox()函数寻找最优lambda值
b <- boxcox(model)
lambda <- b$x[which.max(b$y)]
# [1] 0.5454545


## 变换后的经验回归方程
model2<-lm(sqrt(Y)~X,data=data)   # 取lambda=0.5
summary(model2)
par(ann=FALSE, mfrow= c(1,2))
plot(model2,1, pch=19, lwd=3) 
title(
  xlab = expression(italic(hat(z)[i])),  
  ylab = expression(italic(hat(e)[i])),  
  cex.lab = 1.1,         # 标签字体大小
  line = 2.6
)
plot(model2,2, pch=19, lwd=3)
title(
  xlab = expression(italic(q[(i)])),  
  ylab = expression(italic(r[i])),  
  cex.lab = 1.1,         # 标签字体大小
  line = 2.6
)
e<-residuals(model2)
shapiro.test(e)







################## 例6.4.5 ##################
install.packages("zoo")
install.packages("lmtest")
install.packages("orcutt")
library(zoo)
library(lmtest)
data = read.table("/Users/macbook/Desktop/线性模型例题/Chapter6/data/例6.4.5.txt", header=TRUE) 
lm.reg=lm(y~x, data)
summary(lm.reg)
e<-lm.reg$residuals
e2<-e[-1]   ##e_t
e1<-e[-31]  ##e_{t-1}
## 图6.4.11
par(ann=TRUE, mfrow=c(2,2))
plot(data$x,data$y, xlab=expression(italic(x[t])),ylab=expression(italic(y[t])), pch=19, lwd=2)
plot(predict(lm.reg), e, xlab=expression(italic(hat(y)[t])),ylab=expression(italic(hat(e)[t])), pch=19, lwd=1)
plot(data$i,e, pch=19,  xlab =expression(italic("t")), ylab = expression(hat(italic("e"))[italic("t")]),  lwd = 1)
plot(e2, e1, pch=19, xlab =expression(hat(italic("e"))[italic("t-1")]), ylab = expression(hat(italic("e"))[italic("t")]), lwd = 1)

## D-W检验
dwtest(lm.reg, alternative ="greater")  # Durbin-Watson test statistic
library(orcutt)
cochrane.orcutt(lm.reg)




