如何执行引导程序并在 R 中找到置信区间答案

【问题标题】：How to perform a bootstrap and find confidence interval in R如何执行引导程序并在 R 中找到置信区间
【发布时间】：2021-02-19 01:33:43
【问题描述】：

我想创建一个自定义引导函数，因为我想更好地了解引导程序在做什么，而且似乎其他引导库无法解决我的问题。

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

现在，我遇到了这种类型的错误：

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

目标：我的目标是让 bootresults 数据集返回 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 的新手。

【问题讨论】：

标签： r confidence-interval statistics-bootstrap

【解决方案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)

【讨论】：