6

nlme()を使用してモデルをフィッティングしましたpackage nlme。

ここで、パラメーターの不確実性を考慮して、いくつかの予測区間をシミュレートしたいと考えています。

この目的のために、固定効果の分散行列を抽出する必要があります。

私が知る限り、これを行うには2つの方法があります。

vcov(fit)

と

summary(fit)$varFix

これら 2 つは同じ行列を与えます。

しかし、調べてみると

diag(vcov(fit))^.5

これは、報告された Std Error と同じではありません。summary(fit)

これら2つが同じであると期待するのは間違っていますか?

編集：これはコード例です

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: オプションの論理値。'TRUE' で、'object' を取得するために使用された推定方法が最尤法であった場合、残差標準誤差に sqrt(nobs/(nobs - npar)) を掛けて、REML のような推定値に変換します。この引数は、単一の適合オブジェクトが関数に渡される場合にのみ使用されます。デフォルトは「TRUE」です。

内部を見ると(オブジェクトには class があるためnlme:::summary.lme、オブジェクトの概要を生成するためにも使用されるメソッドです)、次のことがわかります。nlmec("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 に答える