Nlme fit: vcov vs resume

Question

I installed the model using nlme() from package nlme .

Now I want to model some prediction intervals, taking into account the uncertainty of the parameter.

To this end, I need to extract the variance matrix for fixed effects.

As far as I know, there are two ways to do this:

 vcov(fit) 

and

 summary(fit)$varFix 

These two give the same matrix.

However, if I check

 diag(vcov(fit))^.5 

this does NOT match the reported std error in summary(fit)

Am I really mistaken that these two are the same?

Edit: Here is sample code

 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

This is related to the term bias correction; It is documented in ?summary.lme .

adjustSigma: optional boolean. If the "TRUE" and evaluation method used to obtain the object was the maximum probability, the residual standard error is multiplied by sqrt (nobs / (nobs - npar)), converting it to a REML-like evaluation. This argument is used only when one installed object is passed a function. The default value is TRUE.

If you look inside nlme:::summary.lme (this is the method used to generate a summary of the nlme object, since it has a class c("nlme", "lme") ), you will see:

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

That is, the standard error scales sqrt(n/(np)) , where n is the number of observations and p number of parameters of the fixed effect. Check this:

  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 

I must admit that I found the answer in the code and only then looked again to make sure that it is completely clearly stated in the documentation ...