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-ish方式来做。首先，如果您对引导很认真，您可以使用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库的明显优势之一是支持并行操作。这很简单：

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

如果您使用的是 Windows，请parallel="snow"改用。

于 2013-02-06T21:14:56.340 回答