nlme fit: vcv versus summarum

Ich habe ein Modell mit nlme() aus dem package nlme ausgestattet.nlme fit: vcv versus summarum

Jetzt möchte ich einige Prädiktionsintervalle unter Berücksichtigung der Parameterunsicherheit simulieren.

Zu diesem Zweck muss ich die Varianzmatrix für die festen Effekte extrahieren.

Soweit mir bekannt ist, gibt es zwei Möglichkeiten, dies zu tun:

vcov(fit)

und

summary(fit)$varFix

Diese beiden die gleiche Matrix geben.

Wenn ich jedoch

diag(vcov(fit))^.5

es prüfen ist nicht das gleiche wie die gemeldeten Std Fehler bei summary(fit)

Bin ich falsch diese beide gleich sein zu erwarten?

Edit: Hier ist ein Codebeispiel

require(nlme) 

f=expression(exp(-a*t)) 
a=c(.5,1.5) 
pts=seq(0,4,by=.1) 

sim1=function(t) eval(f,list(a=a[1],t))+rnorm(1)*.1 
y1=sapply(pts,sim1) 

sim2=function(t) eval(f,list(a=a[2],t))+rnorm(1)*.1 
y2=sapply(pts,sim2) 

y=c(y1,y2) 
t=c(pts,pts) 
batch=factor(rep(1:2,82)) 
d=data.frame(t,y,batch) 

nlmeFit=nlme(y~exp(-a*t), 
    fixed=a~1, 
    random=a~1|batch, 
    start=c(a=1), 
    data=d 
) 

vcov(nlmeFit) 
summary(nlmeFit)$varFix 
vcov(nlmeFit)^.5 
summary(nlmeFit)

Quelle

2014-05-07 Ketil B T

Sie sind eher um Hilfe zu erhalten, wenn Sie Ihren Datensatz oder zumindest eine repräsentative Stichprobe bereitstellen und den Code anzeigen, den Sie für die Anpassung verwendet haben. – jlhoward

Ich stimme zu. Aber der Datensatz ist nicht meins und ich dachte, dass jemand, der wahrscheinlich antworten könnte, in der Vergangenheit nlme verwendet hätte und daher schnell verfügbar wäre. Da ich auf ein Thema hinweise, das im Allgemeinen datenunabhängig sein sollte, hatte ich gehofft, dass es kein Problem werden würde. Das heißt, wenn Menschen die Nicht-Gleichheit der beiden Matrizen in ihren eigenen Beispielen nicht bestätigen können, wäre das ein ziemlich großer Hinweis, dass ich etwas falsch mache. Aber ich kann weggehen und einen Datensatz simulieren, wenn Sie denken, dass das helfen wird. –

Ja, es ist wichtig, das Problem mit Daten aufzuzeigen, die Sie posten können. – jlhoward

Dies ist auf eine Vorspannung Korrekturterm; es ist in ?summary.lme dokumentiert.

adjustSigma: ein optionaler logischer Wert. Wenn 'TRUE' und die Schätzung , die zum Erhalten des 'Objekts' verwendet wurde, maximale Wahrscheinlichkeit war, wird der Restfehler mit sqrt (nobs/(nobs - npar)) multipliziert und in eine REML-ähnliche Schätzung umgewandelt. Dieses Argument wird nur verwendet, wenn ein einzelnes angepasstes Objekt an die -Funktion übergeben wird. Standard ist 'TRUE'.

Wenn Sie in nlme:::summary.lme aussehen (was ist das Verfahren verwendet auch die Zusammenfassung eines nlme Objekt zu erzeugen, da es Klasse c("nlme", "lme")), sehen Sie:

... 
stdFixed <- sqrt(diag(as.matrix(object$varFix))) 
... 
if (adjustSigma && object$method == "ML") 
    stdFixed <- stdFixed * sqrt(object$dims$N/(object$dims$N - 
     length(stdFixed)))

Das heißt, der Standard Der Fehler wird durch sqrt(n/(n-p)) skaliert, wobei n die Anzahl der Beobachtungen und p die Anzahl der Parameter mit festem Effekt ist. Lassen Sie uns diese Besuche:

library(nlme) 
fm1 <- nlme(height ~ SSasymp(age, Asym, R0, lrc), 
      data = Loblolly, 
      fixed = Asym + R0 + lrc ~ 1, 
      random = Asym ~ 1, 
      start = c(Asym = 103, R0 = -8.5, lrc = -3.3)) 
summary(fm1)$tTable[,"Std.Error"] 
##  Asym   R0  lrc 
## 2.46169512 0.31795045 0.03427017 

nrow(Loblolly) ## 84 
sqrt(diag(vcov(fm1)))*sqrt(84/(84-3)) 
##  Asym   R0  lrc 
## 2.46169512 0.31795045 0.03427017

Ich muss zugeben, dass ich die Antwort im Code gefunden und sah nur dann zurück, um zu sehen, dass es vollkommen klar in der Dokumentation angegeben wurde ...

Quelle

2014-05-08 20:25:21

Großartig! Vielen Dank - ich hätte nie gedacht, in der zusammenfassenden Dokumentation nachzusehen. Und vielleicht bin ich der Faulheit schuldig, weil ich nicht in den Code geschaut habe. Ich glaube, ich habe Ihre Antwort akzeptiert und aufgewertet. –

nlme fit: vcv versus summarum

Antwort

Verwandte Themen