2013-02-19 83 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