2013-02-19 11 views
9

Gibt es eine einzelne Funktion, ähnlich wie "runf", "rnorm" und dergleichen, die simulierte Vorhersagen für ein lineares Modell erzeugen? Ich kann es selbst programmieren, aber der Code ist hässlich und ich nehme an, dass dies etwas ist, was jemand vorher getan hat.Gibt es eine Funktion oder ein Paket, das Vorhersagen für ein von lm() zurückgegebenes Objekt simuliert?

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

Was ich bin interessiert, ist eine Funktion, die unten wie die Linie aussieht:

sims = rlm(1000, fit, newX) 

Diese Funktion 1000 Simulationen von y-Werte, basierend auf den neuen x-Variablen zurückkehren würde.

+1

Die letzte Zeile in Ihrem Q hat mich verwirrt. 'x' ist festgelegt; meinst du simulieren 'y' (die Antwort) für die neuen' x' Daten? –

+0

Entschuldigung, Gavin, du hast Recht. Ich wollte sagen, dass die Antworten simuliert werden. Dies wurde bearbeitet. – PirateGrunt

+2

OK, Sie könnten also '' simulate' 'sehen, aber das funktioniert nur für das aktuelle 'x'. Aber Sie könnten es ändern ('simulate.lm()'), um 'propective()' auf dem Modellobjekt mit 'newdata = newX' anstelle des aktuellen Aufrufs von' fitted() 'aufzurufen und es dann wie mit fortfahren zu lassen der normale Code. Angenommen, "Gewichte" wurde nicht verwendet, da dies die Angelegenheit komplizieren würde ... –

Antwort

8

Zeigen, dass Gavin Simpsons Vorschlag, stats:::simulate.lm zu modifizieren, ein praktikabler ist.

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

das funktioniert, aber dies ist ein sauberer (und besser wahrscheinlich) Ansatz:

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

Das ist fantastisch. Die zweite Lösung ist schön und sauber. Der erste verwendet einige Debugging-Techniken, die ich noch nicht gelernt habe. Ausgezeichnet. Vielen Dank. – PirateGrunt

+0

@PirateGrunt - Froh, dass Sie das 'trace()' Beispiel schätzen - es ist ein echtes Powertool für R-Entwicklung und Code-Exploration. Wenn Sie es viel verwenden, können Sie [die Antworten auf diese Frage] schätzen (http://stackoverflow.com/questions/11319161/what-isa-fast-way-to-set-debugging-code-at -a-given-line-in-a-function), die es wesentlich einfacher machen, den Wert von at = 'zu finden, der dem gewünschten Codeeinfügepunkt entspricht. (Ich verwende oft 'print.func()' aus Michael Hoffmans Antwort. Probieren Sie es aus, zum Beispiel mit 'print.func (stats ::: simulate.lm)'). Genießen! –