首先從負二項差異自舉

Question

新手在這裏。我安裝上計數數據負二項式模型，其中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) 

這是引導第一個區別的正確方法嗎？

Answer 1

通常情況下，您的方法沒問題，但您可以在更多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"。