偏差校正百分位自举答案

【问题标题】：bias-corrected percentile bootstrap偏差校正百分位自举
【发布时间】：2019-05-21 15:55:48
【问题描述】：

我正在尝试针对 3 种不同的效果大小运行经过偏差校正的百分位引导程序，并运行了前任的一些代码。虽然此代码有效，但运行时间太长了。完成 3 种不同的效果大小中的每一种大约需要 3 天。

它还显示了我不需要的 1000 次迭代中的每一次。理想情况下，我想最大限度地优化运行时间。我还希望控制台上的所有 3 个输出或 word/excel 文档中的所有输出，而不必一次运行它们。

我知道代码有点笨拙，但我们将不胜感激。

我对 R 相当陌生，正在寻找一些关于如何正确进行的建议。代码如下：

library(pscl)
library(boot)

# RNG MODULE FOR TWO_PART HURDLE MODEL

gen.hurdle = function(n, a, b1, b2, c1, c2, i0, i1, i2){

  x = round(rnorm(n),3)
  e = rnorm(n)
  m = round(i0 + a*x + e, 3)

  lambda = exp(i1 + b1*m + c1*x)                       # PUT REGRESSION TERMS FOR THE CONTINUUM PART HERE; KEEP exp()
  ystar = qpois(runif(n, dpois(0, lambda), 1), lambda) # Zero-TRUNCATED POISSON DIST.; THE CONTINUUM PART

  z = i2 + b2*m  + c2*x                                # PUT REGRESSION TERMS FOR THE BINARY    PART HERE
  z_p = exp(z) / (1+exp(z))                            # p(1) = 1-p(0)
  tstar = rbinom(n, 1, z_p)                            # BINOMIAL DIST.         ; THE BINARY    PART

  y= ystar*tstar                                       # TWO-PART COUNT OUTCOME

  return(cbind(x,m,y,z,z_p,tstar))
}


##################################################################################################
# MODEL ##########################################################################################
##################################################################################################

# i=1  ###################################

#small effect
seed=51; n=50  ;a=.18; b=.16; c=.25; i=1
seed=53; n=100 ;a=.18; b=.16; c=.25; i=1
seed=55; n=200 ;a=.18; b=.16; c=.25; i=1
seed=57; n=300 ;a=.18; b=.16; c=.25; i=1
seed=58; n=500 ;a=.18; b=.16; c=.25; i=1
seed=59; n=1000;a=.18; b=.16; c=.25; i=1

#medium effect
seed=61; n=50  ;a=.31; b=.35; c=.25; i=1
seed=63; n=100 ;a=.31; b=.35; c=.25; i=1
seed=65; n=200 ;a=.31; b=.35; c=.25; i=1
seed=67; n=300 ;a=.31; b=.35; c=.25; i=1
seed=68; n=500 ;a=.31; b=.35; c=.25; i=1
seed=69; n=1000;a=.31; b=.35; c=.25; i=1

#large effect
seed=81; n=50  ;a=.52; b=.49; c=.25; i=1
seed=73; n=100 ;a=.52; b=.49; c=.25; i=1
seed=75; n=200 ;a=.52; b=.49; c=.25; i=1
seed=77; n=300 ;a=.52; b=.49; c=.25; i=1
seed=78; n=500 ;a=.52; b=.49; c=.25; i=1
seed=79; n=1000;a=.52; b=.49; c=.25; i=1


#model
set.seed(seed)
iterations = 1000
r = 1000

results = matrix(,iterations,4)
for (iiii in 1:iterations){

  data  = data.frame(gen.hurdle(n, a, b, b, c, c, i, i, i))
  data0 = data.frame(gen.hurdle(n, a, 0, 0, c, c, i, i, i))

  p_0     =1-mean(data$z_p)
  p_0_hat =1-mean(data$tstar)
  p_0_b0     =1-mean(data0$z_p)
  p_0_hat_b0 =1-mean(data0$tstar)

  # power

  boot= matrix(,r,8)
  for(jjjj in 1:r){

    #power
    fit1      = lm(m ~ x, data=data)
    fit2      = hurdle(formula = y ~ m + x, data=data, dist = "poisson", zero.dist = "binomial") 

    a_hat     = summary(fit1)$coef[2,1]
    b1_hat    = summary(fit2)[[1]]$count[2,1]
    b2_hat    = summary(fit2)[[1]]$zero[2,1]
    ab1_hat   = prod(a_hat,b1_hat)
    ab2_hat   = prod(a_hat,b2_hat)

    boot.data = data[sample(nrow(data), replace = TRUE), ]
    boot.data$y[1] = if(prod(boot.data$y) > 0) 0 else boot.data$y[1]

    boot.fit1 = lm(m ~ x, data=boot.data)
    boot.fit2 = hurdle(formula = y ~ m + x, data=boot.data, dist = "poisson", zero.dist = "binomial") 

    boot.a    = summary(boot.fit1)$coef[2,1]
    boot.b1   = summary(boot.fit2)[[1]]$count[2,1]
    boot.b2   = summary(boot.fit2)[[1]]$zero[2,1]
    boot.ab1  = prod(boot.a,boot.b1)
    boot.ab2  = prod(boot.a,boot.b2)

    #type I error
    fit3       = lm(m ~ x, data=data0)
    fit4       = hurdle(formula = y ~ m + x, data=data0, dist = "poisson", zero.dist = "binomial")  

    a_hat_b0   = summary(fit3)$coef[2,1]
    b1_hat_b0  = summary(fit4)[[1]]$count[2,1]
    b2_hat_b0  = summary(fit4)[[1]]$zero[2,1]
    ab1_hat_b0 = prod(a_hat_b0,b1_hat_b0)
    ab2_hat_b0 = prod(a_hat_b0,b2_hat_b0)

    boot.data0 = data0[sample(nrow(data0), replace = TRUE), ]
    boot.data0$y[1] = if(prod(boot.data0$y) > 0) 0 else boot.data0$y[1]

    boot.fit3  = lm(m ~ x, data=boot.data0)
    boot.fit4  = hurdle(formula = y ~ m + x, data=boot.data0, dist = "poisson", zero.dist = "binomial")  

    boot.a_b0   = summary(boot.fit3)$coef[2,1]
    boot.b1_b0  = summary(boot.fit4)[[1]]$count[2,1]
    boot.b2_b0  = summary(boot.fit4)[[1]]$zero[2,1]
    boot.ab1_b0 = prod(boot.a_b0,boot.b1_b0)
    boot.ab2_b0 = prod(boot.a_b0,boot.b2_b0)

    boot[jjjj,] = c(ab1_hat,    ab2_hat,    boot.ab1,    boot.ab2,
                    ab1_hat_b0, ab2_hat_b0, boot.ab1_b0, boot.ab2_b0)

  }

  z0.1 = qnorm((sum(boot[,3] > boot[,1])+sum(boot[,3]==boot[,1])/2)/r)
  z0.2 = qnorm((sum(boot[,4] > boot[,2])+sum(boot[,4]==boot[,2])/2)/r)
  z0.1_b0 = qnorm((sum(boot[,7] > boot[,5])+sum(boot[,7]==boot[,5])/2)/r)
  z0.2_b0 = qnorm((sum(boot[,8] > boot[,6])+sum(boot[,8]==boot[,6])/2)/r)

  alpha=0.05 # 95% limits
  z=qnorm(c(alpha/2,1-alpha/2)) # Std. norm. limits

  p1    = pnorm(z-2*z0.1) # bias-correct & convert to proportions
  p2    = pnorm(z-2*z0.2) 
  p1_b0 = pnorm(z-2*z0.1_b0)
  p2_b0 = pnorm(z-2*z0.2_b0) 

  ci1    = quantile(boot[,3],p=p1) # Bias-corrected percentile lims
  ci2    = quantile(boot[,4],p=p2)
  ci1_b0 = quantile(boot[,7],p=p1_b0)
  ci2_b0 = quantile(boot[,8],p=p2_b0)

  sig.ab1 = if(prod(ci1) > 0) 1 else 0
  sig.ab2 = if(prod(ci2) > 0) 1 else 0
  sig.ab1_b0 = if(prod(ci1_b0) > 0) 1 else 0
  sig.ab2_b0 = if(prod(ci2_b0) > 0) 1 else 0



  #results
  results[iiii,] = c(sig.ab1, sig.ab2, sig.ab1_b0, sig.ab2_b0)

  message(paste0(iiii, " / iterations"))
  flush.console()
}

i
n
a
b
iterations
#bootstrap how many
r

#power of ab1
mean(results[,1])
#power of ab2
mean(results[,2])
#type I error of ab1
mean(results[,3])
#type I error of ab2
mean(results[,4])

在我看来，每个效果单独运行的问题在于将每个循环的结果命名为“结果”。我不知道在不干扰循环的情况下打印所有结果（针对每个效果大小）的最佳方法。

长度显然来自于大量的迭代和缺乏足够快的 RAM 来处理它们。这是可以补救的吗？与在程序运行时无法获得所有效果大小的结果相比，时间问题几乎没有那么重要。

【问题讨论】：

看起来还需要很多改进。如果代码有效，则可能是Code Review 的问题无论如何，您可以将引导程序包装到function 并使用replicate 而不是for 循环进行迭代。将data$variable替换为data["variable"]，不要多次使用summary(boot.fit2)，一次进行并将结果存储在一个矩阵中，而不是使用.lm.fit而不是lm，尽量使用整个colMeans数据框而不是每个变量的mean，使用数字矩阵而不是数据框等
@jay.sf 感谢您抽出宝贵时间查看此内容。我不知道 Code Review 甚至是一个论坛！我将在那里交叉发布。我承认它需要做很多工作。如果你有时间，你能解释一下如何使引导程序成为一个函数吗？
我投票结束这个问题，因为这个网站太宽泛了，不能在这里on-topic，它是cross posted on CR
@SᴀᴍOnᴇᴌᴀ 关闭它还是版主更好？
@JBacon 你应该会在标签下看到一个 delete 链接，你可以用它来删除这篇文章

标签： r statistics-bootstrap pscl

【解决方案1】：

为了部分地给你答案，这里是如何引导一个函数。

考虑这些数据，

summary(beaver2[c("time", "activ")]); nrow(beaver2)
#      time          activ     
# Min.   :   0   Min.   :0.00  
# 1st Qu.:1128   1st Qu.:0.00  
# Median :1535   Median :1.00  
# Mean   :1446   Mean   :0.62  
# 3rd Qu.:1942   3rd Qu.:1.00  
# Max.   :2350   Max.   :1.00  
# [1] 100

还有这个回归：

summary(with(beaver2, lm(temp ~ activ)))$coe
#               Estimate Std. Error    t value      Pr(>|t|)
# (Intercept) 37.0968421 0.03456240 1073.32955 2.895092e-201
# activ        0.8062224 0.04389429   18.36736  1.682051e-33

使用for 循环，我们会像这样引导：

set.seed(528)

for.res <- rep(NA, 5e2)
for (i in 1:5e2) {
  X <- beaver2[sample(seq(nrow(beaver2)), replace=TRUE), ]
  y0 <- with(X, mean(temp[activ == 0]))
  y1 <- with(X, mean(temp[activ == 1]))
  for.res[i] <- y1 - y0
}

cat("beta:", b <- mean(for.res), "\tse:", s <-  sd(for.res), "\tt:", t <- b/s,
    "\tp:", 2 * pt(-abs(t), df = nrow(beaver2) - 1))
# beta: 0.8041428   se: 0.04089321  t: 19.66446     p: 5.761707e-36

而对于一个函数，我们会像这样引导：

bootFun <- function(X) {
  y0 <- with(X, mean(temp[activ == 0]))
  y1 <- with(X, mean(temp[activ == 1]))
  return(yh1 - yh0)
}

set.seed(528)
boot.res <- replicate(5e2, bootFun(X=beaver2[sample(seq(nrow(beaver2)), 
                                              replace=TRUE), ]))


cat("beta:", b <- mean(boot.res), "\tse:", s <-  sd(boot.res), "\tt:", t <- b/s,
    "\tp:", 2 * pt(-abs(t), df = nrow(beaver2) - 1))
# beta: 0.8049935   s: 0.04255446   t: 18.91678     p: 1.203037e-34

在上面的例子中，replicate 相对于for 的速度提高了大约 10%：

microbenchmark::microbenchmark(
  replicate=replicate(5e2, bootFun(beaver2[sample(seq(nrow(beaver2)), 
                                              replace=TRUE), ])),
  for.loop={
    for (i in 1:5e2) {
      X <- beaver2[sample(seq(nrow(beaver2)), replace=TRUE), ]
      y0 <- with(X, mean(temp[activ == 0]))
      y1 <- with(X, mean(temp[activ == 1]))
      for.res[i] <- y1 - y0
    }
  })
# Unit: milliseconds
#      expr      min       lq      mean   median       uq      max neval cld
# replicate 88.37455 89.12988  93.36896 89.66325  99.8587 133.1689   100  a 
#  for.loop 95.88837 96.65481 102.40125 97.07489 107.5985 295.3379   100   b

【讨论】：