1

ここの初心者。Y がイベント数、D が治療、X が対数オフセットであるカウント データに負の二項モデルを当てはめています。

out <- glm.nb(y ~ d + offset(log(x)),data=d1)

D=1 と D=0 の最初の差の信頼区間をブートストラップしたいと思います。私はここまでやってきましたが、それが正しいアプローチであるかどうかはわかりません:

holder <- matrix(NA,1200,1)
out <- out <- glm.nb(y ~ d + offset(log(x)),data=d1)

for (i in 1:1200){
q <- sample(1:nrow(d1), 1)
d2 <- d1[q,]
d1_1 <- d1_2 <- d2
d1_1$d <- 1
d1_2$d <- 0
d1pred <- predict(out,d1_1,type="response")
d2pred <- predict(out,d1_2,type="response")
holder[i,1] <- (d1pred[1] - d2pred[1])
}
mean(holder)

これは最初の違いをブートストラップする正しい方法ですか?

4

1 に答える 1

0

一般的に、あなたのアプローチは問題ありませんが、よりR っぽい方法で行うことができます。まず、ブートストラップに真剣に取り組んでいる場合は、bootライブラリを使用して、よりコンパクトなコード、ループなし、その他多くの高度なオプションの恩恵を受けることができます。

あなたの場合、次のようになります。

## Data generation
N <- 100
set.seed(1)
d1 <- data.frame(y=rbinom(N, N, 0.5),
                 d=rbinom(N, 1, 0.5),
                 x=rnorm(N, 10, 3))
## Model
out <- glm.nb(y ~ d + offset(log(x)), data=d1)

## Statistic function (what we are bootstrapping)
## Returns difference between D=1 and D=0
diff <- function(x,i,model){
  v1 <- v2 <- x[i,]
  v1$d <- 1
  v2$d <- 0
  predict(model,v1,type="response") - predict(model,v2,type="response")
}

## Bootstrapping itself
b <- boot(d1, diff, R=5e3, model=out)
mean(b$t)

b$tブートストラップされた値を保持するようになりました。詳細についてはnames(b)、および/または?bootを参照してください。

ブートストラップは時間のかかる操作であり、bootライブラリの明白な利点の 1 つは並列操作のサポートです。次のように簡単です。

 b <- boot(d1, diff, R=5e3, model=out, parallel="multicore", ncpus=2)

Windows を使用している場合は、parallel="snow"代わりに使用してください。

于 2013-02-06T21:14:56.340 に答える