2013-02-19 81 views
9

是否有一个单一的函数,类似于“runif”,“rnorm”等会产生线性模型的模拟预测?我可以自己编写代码,但代码很难看,而且我认为这是以前做过的事情。是否有函数或包将模拟从lm()返回的对象的预测?

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) 

这个函数将返回1000个模拟y值的基础上,新的x变量。

+1

最后一行在你的Q具有我迷惑。 'x'是固定的;你的意思是模拟新'x'数据的'y'(响应)吗? – 2013-02-19 21:37:05

+0

对不起,加文,你是对的。我的意思是说,答复将被模拟。这已被编辑。 – PirateGrunt 2013-02-19 21:39:16

+2

好的,所以你可以看看'?simulate',但只适用于当前的'x'。但是你可以改变它('simulate.lm()')以'newdata = newX'在模型对象上调用'predict()',而不是当前调用'fits()',然后让它继续正常的代码。假设'重量'没有使用,因为这会使问题复杂化... – 2013-02-19 21:41:37

回答

8

显示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(5)) 

## 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) 
newX <- 8:1 

## Pass in new x-values via the argument 'XX' 
simulate(fit, nsim = 4, XX = newX) 
#  sim_1 sim_2  sim_3 sim_4 
# 1 8.047710 8.647585 7.9798728 8.400672 
# 2 6.398029 7.714972 7.9713929 7.813381 
# 3 5.469346 5.626544 4.8691962 5.282176 
# 4 4.689371 4.310656 4.2029540 5.257732 
# 5 4.628518 4.467887 3.6893648 4.018744 
# 6 2.724857 4.280262 2.8902676 4.347371 
# 7 1.532617 2.400321 2.4991168 3.357327 
# 8 1.300993 1.379705 0.1740421 1.549881 

这样的作品,但是这是一个更清洁(可能更好)的方法:

## A function for simulating at new x-values 
simulateX <- function(object, nsim=1, seed=NULL, X, ...) { 
    object$fitted.values <- X 
    simulate(object=object, nsim=nsim, seed=seed, ...) 
} 

## Prepare a fit object and some new x-values 
df <- data.frame(x=x<-1:10, y=1.5*x + rnorm(length(x))) 
fit <- lm(y ~ x, data = df)  
newX <- 8:1 

## Try it out 
simulateX(fit, nsim = 4, X = newX) 
#  sim_1 sim_2  sim_3  sim_4 
# 1 8.828988 6.890874 7.397280 8.1605794 
# 2 6.162839 8.174032 3.612395 7.7999466 
# 3 5.861858 6.351116 3.448205 4.3721326 
# 4 5.298132 4.448778 2.006416 5.7637724 
# 5 7.260219 4.015543 3.063622 4.2845775 
# 6 3.107047 4.859839 6.202650 -1.0956775 
# 7 1.501132 1.086691 -1.273628 0.4926548 
# 8 1.197866 1.573567 2.137449 0.9694006 
+0

太棒了。第二种解决方案非常干净。第一个使用一些我还没学过的调试技术。优秀。谢谢。 – PirateGrunt 2013-02-20 13:51:10

+0

@PirateGrunt - 非常感谢'trace()'示例 - 它是R开发和代码探索的一个真正的powertool。如果你使用它很多,你可能会欣赏[这个问题的答案](http://stackoverflow.com/questions/11319161/what-is-a-fast-way-to-set-debugging-code-at -a-given-line-in-a-function),这使得更容易找到对应于所需代码插入点的at =的值。 (我经常使用Michael Hoffman的回答中的'print.func()',例如,用'print.func(stats ::: simulate.lm)''试试。请享用! – 2013-02-20 14:11:54