## 6.2节例题

################## 例6.2.1 ##################

## 数据准备
y <- c(243, 261, 244, 285, 202, 180, 183, 207, 216, 160, 104, 110)  # 因变量Y
X1 <- c(1.5, 1.5, 1.5, 1.5, 2.0, 2.0, 2.0, 2.0, 2.5, 2.5, 2.5, 2.5 )  # 自变量X1
X2 <- c(6.0, 6.0, 9.0, 9.0, 7.5, 7.5, 7.5, 7.5, 9.0, 9.0, 6.0, 6.0)  # 自变量X2
X3 <- c(1315, 1315, 1890, 1890, 1575, 1575, 1575, 1575, 1315, 1315, 1890, 1890)  # 自变量X3
X0 <- rep(1, 12)  # 常数项（截距）列，全为1
X <- cbind(X0, X1, X2, X3)  # 构造设计矩阵X（包含截距项）

## 最小二乘估计
beta <- solve(t(X) %*% X) %*% t(X) %*% y  # 计算回归系数的LS估计
beta  # 输出估计的回归系数

## 计算协方差矩阵
D <- solve(t(X) %*% X)  # 计算(X^TX)^{-1}（系数的协方差矩阵的倍数）
D * 10^4  

# 计算系数的标准误
s <- sqrt(diag(D) * 435.86)  # 计算每个系数的标准误
# 总平方和、残差平方和、回归平方和
SST <- sum((y - mean(y))^2 )  # 总平方和（总变差）
SSE <- sum((y - X %*% beta )^2)  # 残差平方和（未被解释的变差）
RSS <- SST - SSE  # 回归平方和（被解释的变差）
SS <- c(RSS, SSE, SST)  # 组合成向量
# 自由度
f1 <- ncol(X) - 1  # 回归自由度（自变量个数，p-1，p为参数总数）
f2 <- nrow(X) - ncol(X)  # 残差自由度（n-p）
freed <- c(f1, f2, nrow(X)-1 )  # 总自由度（n-1）

## 均方和与F统计量
MSS <- SS / freed  # 均方和（平方和/自由度）
F <- c(MSS[1]/MSS[2], NA, NA)  # F统计量（用于检验整体回归显著性）
F <- round(F, 2)
fpv <- 1 - pf(F, f1, f2)  # F统计量的P值

## 回归系数的 t 检验
t <- beta / s  # 计算每个系数的t统计量（系数/标准误）
tpvalue <- 2 * (1 - pt(abs(t), f2))  # 计算t检验的双侧P值
cbind(beta, s, t, tpvalue)  # 组合展示系数、标准误、t值、P值

## 输出方差分析表
cbind(SS, freed, MSS, F, fpv)  # 展示平方和、自由度、均方和、F值、P值


## 构建多元线性回归模型
data <- read.table("/Users/macbook/Desktop/线性模型例题/Chapter6/data/例6.2.1.txt",header = T) 
reg <- lm(Y ~ X1+X2+X3,data = data)
summary(reg)     
result <- summary(reg)

## 构建简化模型：移除不显著的自变量X3
summary(lm(Y ~ X1+X2, data = data))   

## OLS 回归扩展功能
install.packages("olsrr")
library(olsrr)
ols_regress(Y ~ X1+X2, data = data)





################## 例6.2.2 ##################
data <- read.table("/Users/macbook/Desktop/线性模型例题/Chapter6/data/例6.2.2.txt",header = T) 

reg <- lm(Y ~ X1+X2+X3+X4+X5+X6, data = data)
summary(reg)

round(cor(model.frame(reg)),3)  ##相关系数






################## 例6.2.3 ##################
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 ~ .,data = data1) ##拟合方程
summary(lm.reg)

##  相关系数系数矩阵和两两变量散点图
data1.names <- c(expression(italic(Y)), 
                 expression(italic(X)[1]), 
                 expression(italic(X)[2]), 
                 expression(italic(X)[3]), 
                 expression(italic(X)[4]))
round(cor(model.frame(lm.reg)),3)  ##相关系数
## 图6.2.1
pairs(model.frame(lm.reg), labels = data1.names, pch = 19, col = "grey50")





################## 例6.2.4 ##################
install.packages("cowplot")
library(ggplot2)
library(cowplot)
data <- read.table("/Users/macbook/Desktop/线性模型例题/Chapter6/data/例6.2.4.txt",header = T) 
result1 <- summary(lm(Y1 ~ X1,data = data))
result2 <- summary(lm(Y2 ~ X2,data = data))
result3 <- summary(lm(Y3 ~ X3,data = data))
result4 <- summary(lm(Y4 ~ X4,data = data))

rbind(result1$coef,result2$coef,result3$coef,result4$coef) ###四个模型回归系数

F<-c(result1$fstatistic[1],result2$fstatistic[1],
     result3$fstatistic[1],result4$fstatistic[1])   ###四个模型回归方程的F检验统计量  
R2<-c(result1$r.squared,result2$r.squared,
      result3$r.squared,result4$r.squared)  ###四个模型的R^2
p <-1- pf(F,1,9)   ##F的p值
t(cbind(F,p,R2) )




library(ggplot2)
library(cowplot)

p1 <- ggplot(data = data, aes(x = X1, y = Y1)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE, color = "gray50") +
  theme_classic() +
  xlim(0, 20) + 
  ylim(0, 20) +
  xlab(expression(italic(X[1]))) +  # x轴标签设为斜体X₁
  ylab(expression(italic(Y[1]))) +  # y轴标签设为斜体Y₁
  theme(
    axis.title.x = element_text(face = "italic"),  # 确保x轴标题斜体
    axis.title.y = element_text(face = "italic")   # 确保y轴标题斜体
  )


p2 <- ggplot(data = data, aes(x = X2, y = Y2)) + 
  geom_point() +
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE, color = "gray50") +
  theme_classic() +
  xlim(0, 20) + 
  ylim(0, 20) +
  xlab(expression(italic(X[2]))) +  # 斜体X₂
  ylab(expression(italic(Y[2]))) +  # 斜体Y₂
  theme(
    axis.title.x = element_text(face = "italic"),
    axis.title.y = element_text(face = "italic")
  )


p3 <- ggplot(data = data, aes(x = X3, y = Y3)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE, color = "gray50") +
  theme_classic() +
  xlim(0, 20) + 
  ylim(0, 20) +
  xlab(expression(italic(X[3]))) +  # 斜体X₃
  ylab(expression(italic(Y[3]))) +  # 斜体Y₃
  theme(
    axis.title.x = element_text(face = "italic"),
    axis.title.y = element_text(face = "italic")
  )


p4 <- ggplot(data = data, aes(x = X4, y = Y4)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE, color = "gray50") +
  theme_classic() +
  xlim(0, 20) + 
  ylim(0, 20) +
  xlab(expression(italic(X[4]))) +  # 斜体X₄
  ylab(expression(italic(Y[4]))) +  # 斜体Y₄
  theme(
    axis.title.x = element_text(face = "italic"),
    axis.title.y = element_text(face = "italic")
  )

# 组合图形
plot_grid(p1, p2, p3, p4, nrow = 2)




