nlme fit : vcov versus summary

Question

I have fitted a model using nlme() from the package nlme.

Now I wish to simulate some prediction intervals, taking into account parameter uncertainty.

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

As far as I am aware, there are two ways of doing this:

vcov(fit)

and

summary(fit)$varFix

These two give the same matrix.

However, if I inspect

diag(vcov(fit))^.5

it is NOT the same as the reported Std Error in summary(fit)

Am I wrong to expect these two to be the same?

Edit: Here is a code example

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 due to a bias correction term; it's documented in ?summary.lme.

adjustSigma: an optional logical value. If ‘TRUE’ and the estimation method used to obtain ‘object’ was maximum likelihood, the residual standard error is multiplied by sqrt(nobs/(nobs - npar)), converting it to a REML-like estimate. This argument is only used when a single fitted object is passed to the function. Default is ‘TRUE’.

If you look inside nlme:::summary.lme (which is the method used to generate the summary of an nlme object as well, since it has class c("nlme", "lme")), you 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 is being scaled by sqrt(n/(n-p)) where n is the number of observations and p the number of fixed-effect parameters. Let's 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

I have to admit that I found the answer in the code and only then looked back to see that it was perfectly clearly stated in the documentation ...