如何在R中利用leaps.AIC选择最优预测变量子集构建薪资线性模型

任务背景 #

现有一份包含16个预测变量的数据集，需以salary为响应变量（不使用对数形式），利用leaps包中的leaps命令筛选最优变量子集，构建线性模型。同时使用leaps.AIC函数计算不同规模最优子集模型的AIC和BIC值，辅助模型选择。

完整实现代码 #

# 安装leaps包（已安装则注释此行）
install.packages("leaps")
# 加载leaps包
library(leaps)    

# 定义计算最优子集模型AIC和BIC的函数
leaps.AIC=function(X,y){
  # 参数说明：X为n×k的预测变量矩阵，y为长度n的响应变量向量
  # 函数功能：针对每个变量数量j（1到k），筛选出R²最优的子集模型，计算对应AIC和BIC
  
  # 格式转换，确保y为向量、X为矩阵
  y=as.vector(y)
  X=as.matrix(X)
  num=ncol(X)
  
  # 生成leaps结果，每个变量数量仅保留1个最优模型（按R²排序）
  out=leaps(X,y,method='r2',nbest=1)
  
  # 初始化AIC和BIC存储向量
  aic=1:num
  bic=1:num
  
  # 循环计算每个变量数量下最优模型的AIC和BIC
  for(j in 1:num){
    # 提取当前最优模型包含的变量列索引
    cols=(1:num)[out$which[j,]==TRUE]
    # 拟合线性模型
    fit=lm(y~X[,cols])
    # 计算AIC
    aic[j]=AIC(fit)
    # 计算BIC（BIC的AIC公式中k取log(n)，n为样本量）
    bic[j]=AIC(fit,k=log(length(y)))
  }
  
  # 打印结果
  print("AIC值")
  print(aic)
  print("BIC值")
  print(bic)
  
  # 返回结果列表
  list(aic=aic,bic=bic)
}

# 读取本地CSV数据（请将文件路径替换为你的实际路径）
data <- read.csv("你的数据集文件名.csv")

# 分离响应变量和预测变量：提取salary作为响应变量，其余16个变量作为预测变量
y <- data$salary
X <- data[, !names(data) %in% "salary"]

# 调用函数计算各模型的AIC和BIC
model_metrics <- leaps.AIC(X, y)

# 可选：可视化AIC和BIC随变量数量的变化，辅助选择最优模型
plot(1:ncol(X), model_metrics$aic, type="b", col="blue", xlab="变量数量", ylab="准则值", main="AIC与BIC随变量数量变化")
lines(1:ncol(X), model_metrics$bic, type="b", col="red")
legend("topright", legend=c("AIC", "BIC"), col=c("blue", "red"), lty=1)

模型选择说明 #

• AIC准则：选择AIC值最小的模型，倾向于兼顾拟合优度和模型复杂度
• BIC准则：选择BIC值最小的模型，相比AIC更倾向于简洁模型（对变量数量惩罚更重）
• 结合可视化结果，找到曲线的最低点对应的变量数量，即可确定最优的变量子集规模，再通过out$which提取对应的变量，最终构建线性模型

内容的提问来源于stack exchange，提问作者MKHN15