## 6.3节例题

################## 例6.3.1 ##################
## 使用例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)

## 最优子集选择
install.packages("olsrr")
library(olsrr)
# 使用ols_step_all_possible函数生成所有可能的回归模型
ols_step_all_possible(lm.reg)
best<-ols_step_all_possible(lm.reg)  #（不显示BIC值）
str(best)
best <- best$result # 输出结果

### 表 6.3.1
p <- best$n         #c(rep(1,4),rep(2,6),rep(3,4),4) 
varibles <- best$predictors
R2 <- best$rsquare
adjR2 <- best$adjr
Cp <- best$cp 
PRESS <- best$predrsq 
AIC <- best$aic
BIC <- best$sbc
A <- round(cbind(R2,adjR2, Cp,PRESS, AIC, BIC ),3)
results <- cbind(varibles, A)   

### 图 6.3.1 
# 设置3行2列的图形布局
par(mfrow = c(3, 2), pch = 19, mar = c(5, 5, 4, 2))

# (a) R²图
plot(R2 ~ p,    pch = 19,  
     xlab = "(a)", ylab = expression(italic(R)^2), col = "gray50",
     xaxp = c(1, 4, 3))  # xaxp = c(起点, 终点, 间隔数)，1到4间隔3次即显示1,2,3,4
lines(1:4, tapply(R2, p, max), lwd = 2)

# (b) 调整后R²图
plot(adjR2 ~ p,   pch = 19,   
     xlab = "(b)", ylab = expression(italic(R)[italic(A)]^2), col = "gray50",
     xaxp = c(1, 4, 3))  # 强制显示1-4的整数刻度
lines(1:4, tapply(adjR2, p, max), lwd = 2)

# (c) PRESS图
plot(PRESS ~ p, pch = 19,  
     xlab = "(c)", ylab = expression(PRESS), col = "gray50",
     xaxp = c(1, 4, 3))
lines(1:4, tapply(PRESS, p, min), lwd = 2)

# (d) Cp图
plot(Cp ~ p,   pch = 19, 
     xlab = "(d)", ylab = expression(italic(C)[italic(p)]), col = "gray50",
     xaxp = c(1, 4, 3))
lines(1:4, tapply(Cp,p, min), lwd = 2)

# (e) AIC图
plot(AIC ~ p,  pch = 19, 
     xlab = "(e)", ylab = "AIC", col = "gray50",
     xaxp = c(1, 4, 3))
lines(1:4, tapply(AIC, p, min), lwd = 2)

# (f) BIC图
plot(BIC ~ p, pch = 19,  
     xlab = "(f)", ylab = "BIC", col = "gray50",
     xaxp = c(1, 4, 3))
lines(1:4, tapply(BIC, p, min), lwd = 2)



## 逐步回归分析
lm.aic = step(lm.reg,  direction ="backward") ## AIC准则
lm.aic = step(lm.reg,  direction ="forward")
lm.aic = step(lm.reg,  direction ="both")
summary(lm.aic)                  

lm.bic <- step(lm.reg,k = log(length(data1[,1])),trace = 0) ## BIC准则
summary(lm.bic)






###########（唐年胜，李会琼，2014）变量选择典型案例 ###########
##逐步回归分析
data <- read.table("/Users/macbook/Desktop/线性模型例题/Chapter6/data/tang.txt",header = T)
data1 <- as.data.frame(data)[,-1]

lm.reg <- lm(y ~ x1+x2+x3+x4+x5+x6, data = data1) ##拟合方程
summary(lm.reg)

lm.bic = step(lm.reg, k = log(length(data1[,1])),direction ="both")  ##BIC准则
summary(lm.bic)

### 基于函数regsubsets()变量选择和逐步回归分析 
install.packages("leaps")
library(leaps)

## 定义best()函数：提取最优子集的关键统计量
best <- function(model, ...) {
  # 对输入的模型（lm.reg）执行最优子集选择
  subsets <- regsubsets(formula(model), model.frame(model), ...)
  # 从summary结果中提取统计量：变量数量、入选变量、RSS、R²、调整后R²、Cp、BIC
  subsets <- with(summary(subsets),
                  cbind(p = as.numeric(rownames(which)),  # 模型包含的变量数量
                        which,  # 逻辑矩阵：标记哪些变量被选中（TRUE/FALSE）
                        rss,    # 残差平方和
                        rsq,    # 决定系数R²
                        adjr2,  # 调整后R²
                        cp,     # Mallows' Cp统计量
                        bic))   # BIC信息准则
  return(subsets)  # 返回整理后的结果
}  

# 应用函数到lm.reg模型，保留3位小数
round(best(lm.reg), 3)


## 最优子集回归分析（regsubsets()）
# 对data1数据集中的y与所有自变量（y~.）执行最优子集选择
regit.full <- regsubsets(y ~ ., data = data1)
reg.summary <- summary(regit.full)  # 提取结果摘要
reg.summary  # 打印摘要（包含各模型的统计量）

# 查看摘要中包含的统计量名称
names(reg.summary)  # 输出："which" "rsq" "adjr2" "rss" "cp" "bic"等

# 提取具体统计量
reg.summary$rsq    # 各模型的R²
reg.summary$adjr2  # 调整后R²（修正自由度对R²的影响）
reg.summary$bic    # BIC值（贝叶斯信息准则）


## 可视化统计量，确定最优变量数量
# 设置图形布局为2行2列
par(mfrow = c(2, 2))

# （1）残差平方和（RSS）：值越小，模型拟合越好
plot(reg.summary$rss, xlab = "Number of Variables", ylab = "RSS", type = "l")

# （2）调整后R²：值越大越好，标记最大值点
plot(reg.summary$adjr2, xlab = "Number of Variables", ylab = "Adjusted RSq", type = "l")
id.adj <- which.max(reg.summary$adjr2)  # 找到调整后R²最大的模型索引
points(id.adj, reg.summary$adjr2[id.adj], col = "red", cex = 2, pch = 20)  # 标记红色点

# （3）Mallows' Cp：接近变量数量p+1时模型较优，标记最小值点
plot(reg.summary$cp, xlab = "Number of Variables", ylab = "Cp", type = "l")
id.cp <- which.min(reg.summary$cp)  # 找到Cp最小的模型索引
points(id.cp, reg.summary$cp[id.cp], col = "red", cex = 2, pch = 20)

# （4）BIC：值越小越好，标记最小值点
id.bic <- which.min(reg.summary$bic)  # 找到BIC最小的模型索引
plot(reg.summary$bic, xlab = "Number of Variables", ylab = "BIC", type = "l")
points(id.bic, reg.summary$bic[id.bic], col = "red", cex = 2, pch = 20)

## ##由图可知，推荐最优变量为2

par(mfrow = c(2, 2))  # 2行2列布局
# 按不同准则（R²、调整后R²、Cp、BIC）展示变量入选情况
plot(regit.full, scale = "r2")     # 按R²排序，深色表示入选该模型
plot(regit.full, scale = "adjr2")  # 按调整后R²排序
plot(regit.full, scale = "Cp")     # 按Cp排序
plot(regit.full, scale = "bic")    # 按BIC排序

# 提取包含2个变量的最优模型的系数
coef(regit.full, 2)  # 输出：2个自变量的回归系数（含截距）

# 向前逐步回归：从空模型开始，每次添加一个使模型最优的变量（最多19个）
regfit.fwd <- regsubsets(y ~ ., data = data1, nvmax = 19, method = "forward")
summary(regfit.fwd)  # 查看各步骤的模型结果
coef(regfit.fwd, 2)  # 提取包含2个变量的模型系数

# 向后逐步回归：从全模型开始，每次剔除一个对模型影响最小的变量（最多19个）
regfit.bwd <- regsubsets(y ~ ., data = data1, nvmax = 19, method = "backward")
summary(regfit.bwd)  # 查看各步骤的模型结果
coef(regfit.bwd, 2)  # 提取包含2个变量的模型系数





################## 例6.3.2 ##################
## 使用例6.2.2数据
install.packages('latex2exp')
install.packages('glmnet')
library(glmnet)
library(latex2exp)
data <- read.table("/Users/macbook/Desktop/线性模型例题/Chapter6/data/例6.2.2.txt", header = T)
x<-model.matrix(Y~.,data)[,-1]; y=data$Y

## path.plot函数
path.plot =function(lam, beta.hat){
  plot(log(lam), beta.hat[1,],type='o',ylim= c(min(beta.hat),
                                               max(beta.hat)), col=3, pch=1,lty=1,lwd=2,
       xlab=expression(log(lambda)),ylab=expression(hat(beta)[i](lambda)), xpd = T)
  lines(log(lam), beta.hat[2,],type='o',col=2,pch=3,lty=3,lwd=2)
  lines(log(lam), beta.hat[3,],type='o',col=1,pch=2,lty=2,lwd=2)
  lines(log(lam), beta.hat[4,],type='o',col=4,pch=4,lty=4,lwd=2)
  lines(log(lam), beta.hat[5,],type='o',col=5,pch=5,lty=5,lwd=2)
  lines(log(lam), beta.hat[6,],type='o',col=6,pch=6,lty=6,lwd=2)
  legend(x = -1.5,          # x 轴位置（log(λ) 坐标，根据你的数据范围调整）
         y = -0.7, 
         legend=c(TeX('$\\textit{X}_{1}$'),TeX('$\\textit{X}_{2}$'),TeX('$\\textit{X}_{3}$'),TeX('$\\textit{X}_{4}$'),TeX('$\\textit{X}_{5}$'),TeX('$\\textit{X}_{6}$')), 
         col=c(3,2,1,4,5,6), pch=c(1,3,2,4,5,6),lty=c(1,3,2,4,5,6),lwd=2,cex = 1.4,
         text.width = NULL,  # 调整文本宽度，避免文字重叠
         y.intersp = 1, # 增加图例内部行间距
         x.intersp = 0.8 # 增加图例内部列间距
         )
}

## Lasso估计
fit_lasso = glmnet(x, y, alpha = 1, nlambda = 20)
lam = fit_lasso$lambda
beta.hat = as.matrix(fit_lasso$beta)

par(mar = c(5, 5, 4, 2) ) 
path.plot(lam, beta.hat)     ## 绘制Lasso估计的路径图
## 用函数cv.glmnet()选择最优的lambda
set.seed(2021)
cv.lasso = cv.glmnet(x, y, alpha = 1)
par(ann = FALSE)  # 关闭所有默认标签（标题、坐标轴标签）
plot(cv.lasso)    ## 绘制交叉验证误差图
title(
  xlab = expression(-log(lambda)),  
  ylab = 'MSE',  
  cex.lab = 1.1         
)

cv.lasso$lambda.min
cv.lasso$lambda.1se

coef(cv.lasso, s="lambda.min")
coef(cv.lasso, s="lambda.1se")


### 用程序包msgps中的函数mspgs()进行自适应lasso分析
install.packages("msgps")
library(msgps)
# 拟合自适应 Lasso 模型
alasso_fit = msgps(x, y, penalty = "alasso", gamma = 1, lambda = 0)
## 用函数 summary() 汇总结果，并输出结果
summary(alasso_fit)

## 绘制自适应 Lasso 估计的路径图
# 将绘图区域设置为 1 行 2 列，后续绘制两幅图
par(mfrow = c(1, 2))
par(ann = FALSE)
# 绘制GCV准则的路径图
plot(alasso_fit, criterion = "gcv")
# 通过title()函数统一设置标题和坐标轴标签
title(
  main = "GCV",          # 添加标题"GCV"
  xlab = expression(italic(t)),  # x轴标签：斜体t
  ylab = "标准化系数",        # y轴标签
  cex.lab = 1.1,         # 坐标轴标签字体大小
  cex.main = 1.2         # 标题字体大小（可选）
)

# 绘制基于 BIC 准则的路径图
plot(alasso_fit, criterion = "bic")
# 通过title()函数统一设置标题和坐标轴标签
title(
  main = "BIC",          # 添加标题"BIC"
  xlab = expression(italic(t)),  # x轴标签：斜体t
  ylab = "标准化系数",        # y轴标签
  cex.lab = 1.1,         # 坐标轴标签字体大小
  cex.main = 1.2         # 标题字体大小
)




### SCAD估计
install.packages('ncvreg')
library(ncvreg); library(latex2exp)
fit_SCAD = ncvreg(x, y, penalty = "SCAD", nlambda = 20)
lam = fit_SCAD$lambda
beta.hat = fit_SCAD$beta[-1,]
par(mfrow = c(1, 1))
par(ann = FALSE)
path.plot(lam, beta.hat)     ## 绘制SCAD估计的路径图
title(
  xlab = expression(log(lambda)),  
  ylab = '系数',  
  cex.lab = 1.1         
)


#### 用函数cv.glmnet()选择最优的lambda
set.seed(2021)
cv.SCAD = cv.ncvreg(x, y, penalty = "SCAD")
plot(cv.SCAD)                ## 绘制交叉验证误差图
title(
  xlab = expression(log(lambda)),  
  ylab = '交叉验证误差',  
  cex.lab = 1.1         
)
summary(cv.SCAD)             ## 汇总SCAD回归分析结果
cv.SCAD$lambda.min
# [1] 0.1119247
fit = cv.SCAD$fit;  plot(fit)
beta = fit$beta[,cv.SCAD$min]

beta  #### 输出回归系数的SCAD估计：

cv.SCAD$lambda[36]

fit$beta[,36]
