2013-02-19 24 views
9

Có một hàm duy nhất, tương tự như "runif", "rnorm" và các hàm tương tự sẽ tạo ra các dự đoán mô phỏng cho một mô hình tuyến tính không? Tôi có thể tự viết mã, nhưng mã là xấu và tôi cho rằng đây là điều mà ai đó đã làm trước đây.Có chức năng hoặc gói nào sẽ mô phỏng các dự đoán cho một đối tượng được trả về từ lm() không?

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)) 

Những gì tôi quan tâm là một chức năng giống như các dòng dưới đây:

sims = rlm(1000, fit, newX) 

chức năng đó sẽ quay trở lại năm 1000 mô phỏng các giá trị y, dựa trên x biến mới.

+1

Dòng cuối cùng trong Q của bạn đã cho tôi bối rối. 'x' là cố định; bạn có nghĩa là mô phỏng 'y' (phản hồi) cho dữ liệu' x' mới? –

+0

Xin lỗi, Gavin, bạn đã đúng. Tôi có nghĩa là để nói rằng các phản ứng sẽ được mô phỏng. Điều này đã được chỉnh sửa. – PirateGrunt

+2

OK, Vì vậy, bạn có thể nhìn vào '? Mô phỏng' nhưng chỉ hoạt động với' x' hiện tại. Nhưng bạn có thể thay đổi nó ('simulate.lm()') để gọi 'predict()' trên đối tượng mô hình với 'newdata = newX' thay vì lệnh gọi hiện tại thành' fitted() 'và sau đó cho phép nó tiếp tục mã bình thường. Giả sử 'trọng số' không được sử dụng vì điều đó sẽ làm phức tạp vấn đề ... –

Trả lời

8

Cho thấy đề xuất sửa đổi của Gavin Simpson stats:::simulate.lm là một khả thi.

## 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 

đó làm việc, nhưng đây là một trình dọn dẹp (và có khả năng tốt hơn) cách tiếp cận:

## 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

Thật tuyệt vời. Giải pháp thứ hai là tốt đẹp và sạch sẽ. Người đầu tiên sử dụng một số kỹ thuật gỡ lỗi mà tôi chưa học được. Xuất sắc. Cảm ơn bạn. – PirateGrunt

+0

@PirateGrunt - Vui mừng khi bạn đánh giá cao ví dụ 'trace()' - đó là một powertool thực sự cho phát triển R và khám phá mã. Nếu bạn sử dụng nó nhiều, bạn có thể đánh giá cao [câu trả lời cho câu hỏi này] (http://stackoverflow.com/questions/11319161/what-is-a-fast-way-to-set-debugging-code-at -a-given-line-in-a-function), giúp tìm giá trị của 'at =' tương ứng với điểm chèn mã mong muốn dễ dàng hơn nhiều. (Tôi thường sử dụng 'print.func()' từ câu trả lời của Michael Hoffman. Hãy thử nó, ví dụ, với 'print.func (stats ::: simulate.lm)'). Thưởng thức! –