nlme fit: vcov versus summary

Ho installato un modello utilizzando nlme() da package nlme.nlme fit: vcov versus summary

Ora desidero simulare alcuni intervalli di predizione, tenendo conto dell'incertezza dei parametri.

A tal fine, ho bisogno di estrarre la matrice di varianza per gli effetti fissi.

Per quanto ne so, ci sono due modi per farlo:

vcov(fit)

summary(fit)$varFix

Questi due danno la stessa matrice.

Tuttavia, se ispezionare

diag(vcov(fit))^.5

non è lo stesso del Errore Std riportato in summary(fit)

Sono sbagliato aspettarsi questi due essere lo stesso?

Edit: Ecco un esempio di codice

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)

fonte

2014-05-07 Ketil B T

Si hanno maggiori probabilità per ottenere aiuto se fornisci il set di dati, o almeno un campione rappresentativo, e mostra il codice che hai usato per ottenere l'adattamento. – jlhoward

Sono d'accordo. Ma il set di dati non è mio e ho pensato che chiunque fosse in grado di rispondere avrebbe usato nlme in passato e quindi avrebbe potuto facilmente adattarsi. Poiché sto segnalando un problema che dovrebbe essere indipendente dai dati in generale, avevo sperato che non sarebbe diventato un problema. Cioè, se le persone non riescono a confermare la non uguaglianza delle due matrici nei loro stessi esempi, sarebbe un grande indizio che io stia facendo qualcosa di sbagliato. Ma posso andare via e simulare un set di dati se pensi che possa essere d'aiuto. –

Sì, è importante dimostrare il problema con i dati che è possibile pubblicare. – jlhoward

Ciò è dovuto ad un termine di correzione di polarizzazione; è documentato in ?summary.lme.

adjustSigma: un valore logico opzionale. Se 'TRUE' e il metodo di stima utilizzato per ottenere 'oggetto' era la massima verosimiglianza, l'errore standard residuo viene moltiplicato per sqrt (nobs/(nobs - npar)), convertendolo in una stima simile a REML. Questo argomento viene utilizzato solo quando un singolo oggetto montato viene passato alla funzione . L'impostazione predefinita è "VERO".

Se si guarda dentro nlme:::summary.lme (che è il metodo utilizzato per generare il sommario di un oggetto nlme pure, dal momento che ha classe c("nlme", "lme")), si vede:

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

Cioè, lo standard l'errore viene ridimensionato da sqrt(n/(n-p)) dove n è il numero di osservazioni e p il numero di parametri a effetto fisso. Diamo Check this out:

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

devo ammettere che ho trovato la risposta nel codice e solo allora guardato indietro per vedere che era perfettamente chiaramente indicato nella documentazione ...

fonte

2014-05-08 20:25:21

Grande! Grazie mille - non ho mai pensato di cercare nella documentazione di riepilogo. E forse sono colpevole di pigrizia nel non guardare nel codice. Credo di aver accettato e messo in discussione la tua risposta. –

risposta

Problemi correlati