是否有一个函数或包可以模拟从 lm() 返回的对象的预测？答案

【问题标题】：Is there a function or package which will simulate predictions for an object returned from lm()?是否有一个函数或包可以模拟从 lm() 返回的对象的预测？
【发布时间】：2013-02-04 17:50:54
【问题描述】：

是否有一个类似于“runif”、“rnorm”等的函数可以为线性模型生成模拟预测？我可以自己编写代码，但代码很难看，我认为这是以前有人做过的事情。

slope = 1.5
intercept = 0
x = as.numeric(1:10)
e = rnorm(10, mean=0, sd = 1)
y = slope * x + intercept + e
fit = lm(y ~ x, data = df)
newX = data.frame(x = as.numeric(11:15))

我感兴趣的是一个类似于下面一行的函数：

sims = rlm(1000, fit, newX)

该函数将根据新的 x 变量返回 1000 次 y 值模拟。

【问题讨论】：

你 Q 中的最后一行让我很困惑。 x 已修复；你的意思是为新的x 数据模拟y（响应）？
对不起，加文，你是对的。我的意思是说响应将被模拟。这已被编辑。
好的，所以您可以查看?simulate，但这仅适用于当前的x。但是您可以更改它 (simulate.lm()) 以使用 newdata = newX 在模型对象上调用 predict() 而不是当前对 fitted() 的调用，然后允许它像使用普通代码一样继续。假设 weights 没有被使用，因为这会使事情复杂化......
我看到这只是来自 stats.stackexchange。写的这个问题属于这里。但我认为如果在那儿以更一般的方式询问它会更有趣，而不是编程细节。只是我的两分钱。

标签： r regression lm

【解决方案1】：

表明 Gavin Simpson 修改stats:::simulate.lm 的建议是可行的。

## Modify stats:::simulate.lm by inserting some tracing code immediately
## following the line that reads "ftd <- fitted(object)" 
trace(what = stats:::simulate.lm,
      tracer = quote(ftd <- list(...)[["XX"]]),
      at = list(6))

## Prepare the data and 'fit' object 
df <- data.frame(x =x<-1:10, y = 1.5*x + rnorm(length(x)))
fit <- lm(y ~ x, data = df)

## Define new covariate values and compute their predicted/fitted values
newX <- 8:1
newFitted <- predict(fit, newdata = data.frame(x = newX))

## Pass in fitted via the argument 'XX'
simulate(fit, nsim = 4, XX = newFitted)
#        sim_1     sim_2       sim_3     sim_4
# 1 11.0910257 11.018211 10.95988582 13.398902
# 2 12.3802903 10.589807 10.54324607 11.728212
# 3  8.0546746  9.925670  8.14115433  9.039556
# 4  6.4511230  8.136040  7.59675948  7.892622
# 5  6.2333459  3.131931  5.63671024  7.645412
# 6  3.7449859  4.686575  3.45079655  5.324567
# 7  2.9204519  3.417646  2.05988078  4.453807
# 8 -0.5781599 -1.799643 -0.06848592  0.926204

这行得通，但这是一种更清洁（并且可能更好）的方法：

## A function for simulating at new x-values
simulateX <- function(object, nsim = 1, seed = NULL, X, ...) {
    object$fitted.values <- predict(object, X)
    simulate(object = object, nsim = nsim, seed = seed, ...)
}
    
## Prepare example data and a fit object
df <- data.frame(x =x<-1:10, y = 1.5*x + rnorm(length(x)))
fit <- lm(y ~ x, data = df)

## Supply new x-values in a data.frame of the form expected by
## the newdata= argument of predict.lm()
newX <- data.frame(x = 8:1)
   
## Try it out
simulateX(fit,  nsim = 4, X = newX)
#       sim_1     sim_2     sim_3     sim_4
# 1 11.485024 11.901787 10.483908 10.818793
# 2 10.990132 11.053870  9.181760 10.599413
# 3  7.899568  9.495389 10.097445  8.544523
# 4  8.259909  7.195572  6.882878  7.580064
# 5  5.542428  6.574177  4.986223  6.289376
# 6  5.622131  6.341748  4.929637  4.545572
# 7  3.277023  2.868446  4.119017  2.609147
# 8  1.296182  1.607852  1.999305  2.598428

【讨论】：

太棒了。第二种解决方案很好很干净。第一个使用了一些我还没有学过的调试技术。优秀。谢谢。
@PirateGrunt -- 很高兴您欣赏 trace() 示例 -- 它是 R 开发和代码探索的真正强大工具。如果您经常使用它，您可能会喜欢the answers to this quesion，它可以更容易地找到与所需代码插入点对应的at= 的值。（我经常使用迈克尔霍夫曼的回答中的print.func()。试试吧，例如，使用print.func(stats:::simulate.lm)）。享受吧！
如果有多个预测变量，这个函数会是什么样子？
我知道这已经七岁了（StackOverflow 很疯狂），但这是错误的。您正在为 simulationX 函数提供新的协变量数据，但您需要提供的是新的响应。如果在模拟X 函数调用中将X = newX 替换为X = predict(fit, new_data = newX)，这将起作用。
@DaltonHance 是的。已编辑我的答案的两个部分以纠正错误。谢谢！