R中的更新函数未更新模型。

Question

我在 R 中生成了一些数据。

n <- 1000; p <- 30 

X <- matrix(rnorm(n*p), nrow = n, ncol = p)
beta <- c(rep(1, 10), rep(0, 10), rep(-2, 10))
y <- X %*% beta + rnorm(1000) 

接下来，我想对X的列从1到30运行一次逐步回归，将y作为因变量。首先只包括截距，然后只包括截距和第一列，接着添加第二列、第三列等等。我编写了以下代码。

model <- lm(y~1)
for(i in 1:30){ 
    model <- update(model, ~.+X[, i])
    print(model) 
}

现在输出结果显示，每次迭代时，回归的是y关于截距和X[, i]（即X的第i列），而不是之前的列，尽管我每一步都在更新。例如，当i = 4时，模型是一个关于截距和X[, 4]（即X的第4列）的回归，而不是所有的1、2、3、4列。这是为什么呢？

- gtoques

Answer 1

试试这个

model <- lm(y~1)
for(i in 1:30){ 
    model <- update(model, ~.+X[, 1:i])
    print(model) 
}

- Leonardo

---

但是使用update的意义在哪里呢？我是否可以在循环中总是说model <- lm(y~X[, 1:i])而不使用它？ - gtoques

是的，我认为没有必要使用“update”和“model <- lm(y~X[, 1:i])”，因为它们会得到相同的结果。 - Leonardo

Answer 2

你提出的代码无法运行是因为R看到公式的方式以及在评估i之前更新公式的事实。

相关更新方法的源代码可以通过在命令行中运行update.default来查看。你会看到一些错误检查后，它运行call$formula <- update(formula(object), formula.)，这将调用update.formula()函数。update.formula()看到你想要将术语X[, i]添加到公式中，并执行此操作。但是update.formula()此时不会评估i的值，而是依赖于“惰性求值”。如果我们展开循环，就可以更清楚地看到这一点。

form <- y ~ 1
form
#> y ~ 1
i <- 1
form <- update.formula(form, ~. +X[, i])
form
#> y ~ X[, i]
i <- 2
form <- update.formula(form, ~. +X[, i])
form
#> y ~ X[, i]

该公式正在使用符号X[, i]进行更新，然后简化以删除重复的符号。这种惰性计算很有用，因为它意味着我不需要实际定义上述代码中X和y的内容即可运行。R相信我会在尝试使用它们之前创建适当的对象。

在update()更新公式后，它会eval()更新后的调用。此时会评估i并使用其当前值。因此，实际上，即使它不尝试更改公式，下面的循环也会给出与您的循环完全相同的输出。每次lm()运行时，它都会查找要使用的i的当前值。

for(i in 1:30){ 
  model <- lm(y ~ X[, i])
  print(model)
}

为了达到您想要的效果，您可以在 lm() 函数之外以编程方式创建公式，而不是使用 update() 函数。像这样，

n <- 1000; p <- 30 
X <- matrix(rnorm(n*p), nrow = n, ncol = p)
beta <- c(rep(1, 10), rep(0, 10), rep(-2, 10))
y <- X %*% beta + rnorm(1000) 

xnames <- sapply(list(1:ncol(X)), function(x) paste0("X",x))
colnames(X) <- xnames
dat <- data.frame(y,X)

for(i in 1:30){ 
  form <- as.formula(paste0("y ~ ", paste(xnames[1:i], collapse = "+")))
  model <- lm(form, data = dat)
  print(model)
}

---

编辑：

阅读完这篇文章https://notstatschat.rbind.io/2022/06/23/getting-strings-into-code-in-base-r/后，执行公式操作的另一种方法是使用bquote()。这样做的好处是模型摘要包含了正确的公式。

for(i in 1:30){ 
  model <- eval(bquote(update(model, ~. + .(as.name(xnames[[i]])))))
  print(model) 
}

- DavidLukeThiessen