R/ipr_boot.R

In ipr: Iterative Proportional Repartition Algorithm

ipr_boot <- function(X, y, method = "normal", conf.level = 0.95, B = 100, ...){

  # Point estimate
  est <- ipr(X,y,...)$coef

  # Check bootstrap method
  if(!(method %in% c("normal","basic","perc","all"))){
    stop("Unknown method for IPR bootstrap. Should be 'normal', 'basic', 'perc', or 'all'.")
  }

  # Initialize vectors with lower bounds and upper bounds
  lb_norm <- rep(NA,length(est))
  ub_norm <- rep(NA,length(est))
  lb_basic <- rep(NA,length(est))
  ub_basic <- rep(NA,length(est))
  lb_perc <- rep(NA,length(est))
  ub_perc <- rep(NA,length(est))

  # build a matrix with IPR results at each simulation
  n <- dim(X)[1]
  p <- dim(X)[2]
  res <- matrix(NA,nrow=B,ncol=p)
  alpha <- 1-conf.level

  for(j in 1:B){
    ind <- sample(1:n,n,replace=T)
    Xind <- X[ind,]
    yind <- y[ind]
    res[j,] <- ipr(Xind, yind, print.it = FALSE, ...)$coef
  }

  # Standard error for each disease cost, used for normal method
  sd <- apply(res,2,sd, na.rm=T)
  lb_norm <- est - qnorm(1-alpha/2) * sd
  ub_norm <- est + qnorm(1-alpha/2) * sd
  
  # Percentile method
  lb_perc <- apply(res,2,quantile,probs=alpha/2, na.rm=T)
  ub_perc <- apply(res,2,quantile,probs=1-alpha/2, na.rm=T)

  # Basic method
  lb_basic <- 2*est - ub_perc
  ub_basic <- 2*est - lb_perc

  # Collect the results accordingly to the method
  if(method == "normal"){
    finallist <- list(estimate=est,lower=lb_norm, upper=ub_norm)
  }

  if(method == "perc"){
    finallist <- list(estimate=est,lower=lb_perc, upper=ub_perc)
  }

  if(method == "basic"){
    finallist <- list(estimate=est,lower=lb_basic, upper=ub_basic)
  }

  if(method == "all"){
    finallist <- list(estimate=est,lower_normal = lb_norm, upper_normal = ub_norm,
                      lower_perc = lb_perc, upper_perc = ub_perc,
                      lower_basic = lb_basic, upper_basic = ub_basic)
  }

  return(finallist)
}