如何执行 bootstrap 并在 R 中找到置信区间

Question

我想创建一个自定义 bootstrap 函数，因为我想更好地了解 bootstrap 正在做什么，而其他 bootstrap 库似乎无法解决我的问题.

问题：我想创建自己的 wald 置信区间函数，它接收 bootstrap 数据，输出置信区间，测试置信区间是否在一个范围内，并获得覆盖率。

现在，我收到此类错误：

Error in bootresults[i,}<-waldCI(y=bootdata[i], n=numTrials):number of 
  items to replace is not a multiple of replacement length

目标：我的目标是让bootresults数据集达到return 4列（p值，一个显示上限，下限，以及p是否在区间内）和得到一个类似于这个的图表：

瓦尔德区间图

代码：

set.seed(42)
samples10 <- list()
i <- 1
while(i < 100) {
  sample10[[i]] <- rbinom(1500, size=10, prob=i*.01)  ## rows=1500 ;columns=10
  i <- i + 1
}
sample10 <- data.frame(samples10)
colnames(sample10) <- c(seq(.01, .99, .01)) ## p-values

waldconfidenceinterval <- function(y, n, alpha=0.05) {
  p <- colSums(y)/(n*200)
  sd <- sqrt(p*((1 - p)/(n*200)))
  z <- qnorm(c(alpha/2, 1 - alpha/2))
  ci <- p + z*sd
  return(ci)
}

B <- 200
numTrials <- 10
bootresults <- matrix(ncol=length(sample10), nrow=B)  ## rows=200, cols=99
                                                      ## empty matrix in the beginning
set.seed(42)

for(i in seq_len(B)) {
  bootdata <- sample10[sample(B, replace=T), ]
  bootresults[i, ] <- waldCI(y=bootdata[i], n=numTrials)
  ## Pseudocode:
  # boot_test_data$in_interval <- 
  #   ifelse(boot_test_data$lower1 < i/100 & i/100 < boot_test_data$upper1, 1, 0)
  # coverage[i] <- sum(boot_test_data$in_interval) / length(boot_test_data$in_interval)
}

非常感谢任何帮助，因为我是 R 的新手。

Answer 1

您似乎想要初始化一个三维数组 bootresults 而不是二维矩阵。在你的 waldCI() 中你可以使用 colMeans.

waldCI <- function(y, alpha=0.05) {
  p <- colMeans(y)
  se <- sqrt(p*(1 - p)/nrow(y))
  z <- qnorm(1 - alpha/2)
  ci <- p + z*se %*% cbind(lower=-1, upper=1)
  return(ci)
}

B <- 200
numTrials <- 10
## initialize array
bootresults1 <- array(dim=c(ncol(samples10), 4, B), 
                     dimnames=list(c(), c("p.values", "lower", "upper", "in.int"), c()))

set.seed(42)
for(i in seq_len(B)) {
  samp <- samples10[sample(nrow(samples10), numTrials, replace=F), ]
  ci <- waldCI(samp)
  bootresults1[,,i] <- cbind(p.values, ci, in.int=ci[, 1] < p.values & p.values < ci[, 2])
}

coverage <- rowMeans(bootresults[,4,])
plot(p.values, coverage, type="l", main="My Plot")

类似的方法，但更像 R-ish：

p.values <- seq(.01, .99, .01)
set.seed(42)
samples10 <- `colnames<-`(sapply(p.values, function(pr) rbinom(1.5e3, 1, pr)), p.values)

BOOT <- function(numTrials, ...) {
  samp <- samples10[sample(nrow(samples10), numTrials, replace=F), ]
  ci <- waldCI(samp, ...)
  cbind(p.values, ci, in.int=ci[, 1] < p.values & p.values < ci[, 2])
}

B <- 200
numTrials <- 10

set.seed(42)
bootresults2 <- replicate(B, BOOT(numTrials=10))

stopifnot(all.equal(bootresults1, bootresults2))

数据：

请注意，我使用 rbinom(..., size=1, ...) 创建了您的样本数据。使用“p”作为对象名称表明数据应该是二项式的。

set.seed(42)
samples10 <- matrix(nrow=1500, ncol=99, dimnames=list(c(), c(seq(.01, .99, .01))))
i <- 1
while (i < 100) {
  samples10[, i] <- rbinom(1500, size=1, prob=i*.01)  ## rows=1500 ;columns=10
  i <- i + 1
}

没有 while 循环，您可以继续向量化：

p.values <- seq(.01, .99, .01)
set.seed(42)
samples10 <- `colnames<-`(sapply(p.values, function(pr) rbinom(1.5e3, 1, pr)), p.values)

如何执行 bootstrap 并在 R 中找到置信区间

How to perform a bootstrap and find confidence interval in R

r

confidence-interval

statistics-bootstrap