R语言glm模型公式参数非标准求值问题求解

在循环或迭代拟合多个logistic回归模型时，使用glm(mpg_20 ~ poly(horsepower, i), data = Auto)生成的模型对象，其call参数始终保留符号i，而非替换为迭代时的具体数值。以下是三种解决该非标准求值问题的实用方法：

方法1：用bquote()构造已求值的公式 #

bquote()可以将表达式中用.(变量)标记的部分替换为当前环境中的实际值，适合动态生成带参数的公式。

修改后的可复现代码：

# 加载数据
data("Auto", package = "ISLR2")

# 创建二分类响应变量
Auto$mpg_20 <- as.numeric(Auto$mpg < 20)

# 批量拟合模型
mlist <- lapply(
  1:3, 
  \(i) {
    # 构造替换了i值的公式
    eval_formula <- bquote(mpg_20 ~ poly(horsepower, .(i)))
    glm(eval_formula, data = Auto, family = binomial)
  }
)

# 查看模型call参数
lapply(mlist, \(x) x$call)

方法2：用reformulate()从字符串生成公式 #

reformulate()通过字符串拼接的方式构造公式，逻辑直观，尤其适合处理动态变化的项。

代码示例：

mlist <- lapply(
  1:3, 
  \(i) {
    # 拼接右侧项的字符串
    rhs_term <- sprintf("poly(horsepower, %d)", i)
    # 生成完整公式
    eval_formula <- reformulate(rhs_term, response = "mpg_20")
    glm(eval_formula, data = Auto, family = binomial)
  }
)

# 验证call参数
lapply(mlist, \(x) x$call)

方法3：事后修改模型的call参数 #

如果已经完成模型拟合，也可以直接修改模型对象的call部分，将符号i替换为具体数值：

mlist <- lapply(
  1:3, 
  \(i) {
    model <- glm(mpg_20 ~ poly(horsepower, i), data = Auto, family = binomial)
    # 替换call中的i为当前迭代值
    model$call$formula <- do.call(substitute, list(model$call$formula, list(i = i)))
    model
  }
)

# 查看修改后的call
lapply(mlist, \(x) x$call)

预期输出 #

以上三种方法运行后，都会得到符合预期的结果：

[[1]]
glm(formula = mpg_20 ~ poly(horsepower, 1), data = Auto, family = binomial)

[[2]]
glm(formula = mpg_20 ~ poly(horsepower, 2), data = Auto, family = binomial)

[[3]]
glm(formula = mpg_20 ~ poly(horsepower, 3), data = Auto, family = binomial)

注意：由于是logistic回归，建议显式指定family = binomial，避免glm默认使用高斯线性模型族。

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