dopasowanie nlme: vcov versus podsumowanie

Question

Zainstalowałem model używając nlme() z package nlme.

Teraz chcę zasymulować niektóre przedziały predykcyjne, biorąc pod uwagę niepewność parametru.

W tym celu muszę wyodrębnić macierz zależności dla ustalonych efektów.

O ile mi wiadomo, istnieją dwa sposoby, aby to zrobić:

vcov(fit) 

i

summary(fit)$varFix 

Te dwa dają taką samą matrycę.

Jednak gdybym skontrolować

diag(vcov(fit))^.5 

wcale nie jest taka sama, jak donosi Std Błąd w summary(fit)

Mylę się spodziewać tych dwóch będzie tak samo?

Edit: Tutaj jest przykład kodu

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) 

Answer 1

Wynika to okres korekcji bias; jest to udokumentowane w ?summary.lme.

adjustSigma: opcjonalna wartość logiczna. Jeśli "TRUE" i oszacowanie metoda zastosowana do uzyskania "obiektu" było maksymalnym prawdopodobieństwem, standardowy błąd rezydualny jest mnożony przez sqrt (nobs/(nobs - npar)), konwertując go na oszacowanie zbliżone do REML. Ten argument jest używany tylko wtedy, gdy pojedynczy dopasowany obiekt jest przekazywany do funkcji . Domyślna wartość to "PRAWDA".

Jeśli spojrzeć wewnątrz nlme:::summary.lme (która to metoda wykorzystywana do generowania podsumowanie obiektu nlme jak dobrze, ponieważ ma klasę c("nlme", "lme")), widzisz:

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

Oznacza to, że standard błąd jest skalowany przez sqrt(n/(n-p)), gdzie n jest liczbą obserwacji i p liczbą parametrów o stałym efekcie. Powiedzmy to sprawdzić:

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 

Muszę przyznać, że znalazłem odpowiedź w kodzie, a dopiero potem spojrzał, aby zobaczyć, że to było doskonale jasno określone w dokumentacji ...