6

我已经nlme()使用package nlme.

现在我想模拟一些预测区间，考虑到参数的不确定性。

为此，我需要提取固定效应的方差矩阵。

据我所知，有两种方法可以做到这一点：

vcov(fit)

和

summary(fit)$varFix

这两个给出相同的矩阵。

但是，如果我检查

diag(vcov(fit))^.5

它与报告的标准错误不同summary(fit)

我期望这两个相同是错误的吗？

编辑：这是一个代码示例

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)
4

1 回答 1

4

这是由于偏差校正项造成的；它记录在?summary.lme.

adjustSigma：一个可选的逻辑值。如果“真”且用于获得“对象”的估计方法是最大似然，则残差标准误差乘以 sqrt(nobs/(nobs - npar))，将其转换为类似 REML 的估计。此参数仅在将单个拟合对象传递给函数时使用。默认为“真”。

如果您查看内部nlme:::summary.lme（这也是用于生成nlme对象摘要的方法，因为它具有 class c("nlme", "lme")），您会看到：

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

也就是说，标准误差是由观察次数和固定效应参数的sqrt(n/(n-p))数量来衡量的。让我们检查一下：np

 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

我不得不承认我在代码中找到了答案，然后才回头看看它在文档中已经完全清楚地说明了......

于 2014-05-08T20:25:21.967 回答