Parallel implementation of bootCI for the R DCchoice package

bootCI.R

## Computing a bootsrap confidence interval
bootCI.parallel <- function (obj, nboot = 1000, CI = 0.95, individual = NULL, threads = 1) {

    ## Manage any potential incorrect inputs, error out as approrpiate
    if(!inherits(obj, c("sbchoice", "dbchoice", "oohbchoice")))
        stop("the object must be sbchoice, dbchoice, or oohbchoice class")
    if(CI > 1 | CI < 0)
        stop("CI must be between 0 and 1")
    if(!is.null(threads) && !is.numeric(threads))
        stop("threads must be an integer or NULL for autodetection")

    ## If threads is not 1, prepare for parallel processing.  Also, if
    ## threads == NULL, allow registerDoMC to do its own thing.  If
    ## 
    if(threads != 1 || is.null(threads))
        registerDoMC(cores = threads)
    
    tmp.dat <- eval(obj$data, parent.frame())   # retrieving the data from the object
    
    nobs <- obj$nobs
    ind <- 1:nobs
    
    dist <- obj$distribution
    fr <- obj$formula

    if (!is.null(individual)) {
        formula <- formula(obj$formula, lhs = 0, rhs = -2)
        mf.nexX <- model.frame(formula, individual, xlev = obj$xlevels)
        mm.newX <- model.matrix(formula, mf.nexX, contrasts.arg = obj$contrasts)
    }
    
    if(inherits(obj, c("dbchoice", "oohbchoice"))) {
        results <- foreach(i = 1:nboot, .combine = rbind) %dopar% {
            ## determinw the number of rows for bootstrap sample
            ## with replacement and resample
            ind.boot <- sample(ind, nobs, replace = TRUE)  
            boot.dat <- tmp.dat[ind.boot, ]
            
            if(!inherits(obj, "oohbchoice")) {
                suppressWarnings(tmp.re <- dbchoice(fr, data = boot.dat, dist = dist, par = obj$coefficients))
            } else if(inherits(obj, "oohbchoice")) {
                suppressWarnings(tmp.re <- oohbchoice(fr, data = boot.dat, dist = dist, par = obj$coefficients))
            }

            if(tmp.re$convergence){
                if (is.null(individual)) {
                    boot <- DCchoice:::wtp(object = tmp.re$covariates, b = tmp.re$coefficients, bid = tmp.re$bid, dist = tmp.re$dist)
                } else {
                    boot <- DCchoice:::wtp(object = mm.newX, b = tmp.re$coefficients, bid = tmp.re$bid, dist = tmp.re$dist)
                }
            } else {
                i < i - 1        # discard an unconverged trial
            }
        }
    } else if(class(obj) == "sbchoice"){
        results <- foreach(i = 1:nboot, .combine = rbind) %dopar% {
            ind.boot <- sample(ind, nobs, replace = TRUE)
            boot.dat <- tmp.dat[ind.boot, ]
            suppressWarnings(tmp.re <- sbchoice(fr, data = boot.dat, dist = dist, par = obj$coefficients))
            if(tmp.re$glm.out$converged){
                if (is.null(individual)) {
                    boot <- DCchoice:::wtp(object = tmp.re$covariates, b = tmp.re$coefficients, bid = tmp.re$bid, dist = tmp.re$dist)
                } else {
                    boot <- DCchoice:::wtp(object = mm.newX, b = tmp.re$coefficients, bid = tmp.re$bid, dist = tmp.re$dist)
                }
            } else {
                i < i -1        # discard an unconverged trial
            }
        }
    }
    
    results <- as.data.frame(results)
    boot.mean <- unlist(results$meanWTP)
    boot.median <- unlist(results$medianWTP)
    boot.trmean <- unlist(results$trunc.meanWTP)
    boot.adj.trmean <- unlist(results$adj.trunc.meanWTP)
    
    output <- list(mWTP = boot.mean, tr.mWTP = boot.trmean, adj.tr.mWTP = boot.adj.trmean, medWTP = boot.median)
    
    ## sorting the simulation outcomes
    boot.mean <- sort(boot.mean)
    boot.median <- sort(boot.median)
    boot.trmean <- sort(boot.trmean)
    boot.adj.trmean <- sort(boot.adj.trmean)
    
    lb <- 0.5 * (1 - CI)   # lower bound of the empirical distribution
    ub <- CI + lb        # upper bound of the empirical distribution
    
    int <- c(ceiling(nboot*lb), floor(nboot*ub))  # subscripts corresponding to lb and ub
    
    ## 100*CI% bootstrap CI
    CImat <- rbind(boot.mean[int],    # for mean
                   boot.trmean[int],  # for truncated mean
                   boot.adj.trmean[int],  # for truncated mean with adjustment
                   boot.median[int])  # for median
    
    ## the mean estimates in the original outcome
    if (is.null(individual)) {
        tmp.sum <- DCchoice:::wtp(object = obj$covariates, b = obj$coefficients, bid = obj$bid, dist = obj$dist)
    } else {
        tmp.sum <- DCchoice:::wtp(object = mm.newX, b = obj$coefficients, bid = obj$bid, dist = obj$dist)
    }
    
    out <- cbind(c(tmp.sum$meanWTP, tmp.sum$trunc.meanWTP, tmp.sum$adj.trunc.meanWTP, tmp.sum$medianWTP), CImat)
    rownames(out) <- c("Mean", "truncated Mean", "adjusted truncated Mean", "Median")
    colnames(out) <- c("Estimate", "LB", "UB")
    
    if(!is.finite(tmp.sum$meanWTP)){
        out[1, 2:3] <- -999   # return -999 if the original outcome does not meet the condition for a finite mean WTP estimate
        output$mWTP <- -999
    }
    
    output$out <- out
    
    class(output) <- "bootCI"
    return(output)   
}