NLME拟合：vcov与summary对比

Question

我使用nlme()函数从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)

- Ketil B T

---

如果您提供数据集或至少代表性样本，并展示用于拟合的代码，那么您更有可能获得帮助。 - jlhoward

我同意。但数据集不是我的，我推断任何有能力回答的人都曾经使用过nlme，并且已经准备好了nlme拟合。既然我指出的问题在一般情况下应该是与数据无关的，我希望它不会成为问题。也就是说，如果人们不能够确认在他们自己的示例中这两个矩阵的非等价性，那么这将是我做错了什么的很大提示。但如果你认为这会有所帮助，我可以离开并模拟一个数据集。 - Ketil B T

是的，展示你可以发布的数据问题非常重要。 - jlhoward

已编辑帖子以模拟数据。每次运行代码块时，您会得到两个标准差/标准误差估计的不同答案。 - Ketil B T

这个更容易复现：library(nlme); example(nlme); all.equal(sqrt(diag(vcov(fm2))),summary(fm2)$tTable[,"Std.Error"]) - Ben Bolker

Answer 1

这是由于偏差校正项引起的；在?summary.lme中记录了此项。

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

如果您查看nlme ::： summary.lme （它是生成 nlme 对象摘要的方法，因为它具有类 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))，其中n是观测次数，p是固定效应参数的数量。让我们来看一下：

 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

我必须承认，我在代码中找到了答案，然后回头看文档才发现它已经非常明确地说明了......

- Ben Bolker

---

太好了！非常感谢 - 我从没想过要查看摘要文档。也许我有点懒，没有查看代码。我相信我已经接受并赞同了你的答案。 - Ketil B T