Nothing
### ate-bootstrap.R ---
##----------------------------------------------------------------------
## Author: Brice Ozenne
## Created: apr 11 2018 (17:05)
## Version:
## Last-Updated: Mar 7 2022 (08:36)
## By: Thomas Alexander Gerds
## Update #: 343
##----------------------------------------------------------------------
##
### Commentary:
##
### Change Log:
##----------------------------------------------------------------------
##
### Code:
## * calcBootATE
## generate a boot object for the ate function that will be used to compute CI and p.values
calcBootATE <- function(args, n.obs, fct.pointEstimate, name.estimate,
handler, B, seed, mc.cores, cl,
verbose){
# {{{ prepare arguments
## hard copy of the dataset before bootstrap
ls.data <- list(object.event = NULL,
object.treatment = NULL,
object.censor = NULL)
ls.package <- list(object.event = NULL,
object.treatment = NULL,
object.censor = NULL)
for(iModel in c("object.event","object.treatment","object.censor")){ ## iModel <- "object.event"
if(inherits(args[[iModel]]$call$data,"name")){
data.tempo <- eval(args[[iModel]]$call$data)
if(inherits(data.tempo,"function")){
stop("The dataset in argument \'",iModel,"\' has the same name has an existing R function.\n",
"This creates confusion when we try to extract it from the call. Please rename the dataset.\n")
}else{
ls.data[[iModel]] <- data.table::as.data.table(data.tempo)
}
}else{
ls.data[[iModel]] <- args[[iModel]]$call$data
}
if(inherits(args[[iModel]]$call$formula,"name")){
formula.tempo <- eval(args[[iModel]]$call$formula)
if(inherits(formula.tempo,"function")){
stop("The formula in argument \'",iModel,"\' has the same name has an existing R function.\n",
"This creates confusion when we try to extract it from the call. Please rename the formula.\n")
}else{
args[[iModel]]$call$formula <- formula.tempo
}
}
## package to be exported to cluster
if(inherits(args[[iModel]]$call[[1]],"name")){
iSource <- utils::find(deparse(args[[iModel]]$call[[1]]))
iFound <- grepl("package:",iSource)
if(any(iFound)){
ls.package[[iModel]] <- gsub("package:","",iSource[iFound])
}else{
NULL
}
}else{
if(!is.null(environment(args[[iModel]]$call[[1]]))){
ls.package[[iModel]] <- utils::packageName(environment(args[[iModel]]$call[[1]]))
}else{
ls.package[[iModel]] <- NULL
}
}
}
## packages and functions to be exported to the cluster
add.Package <- unique(c("riskRegression","prodlim","data.table","parallel","survival",unlist(ls.package)))
add.Fct <- c("SurvResponseVar","predictRisk.coxphTD","predictRisk.CSCTD","ATE_COMPARISONS")
## if cluster already defined by the user
no.cl <- is.null(cl)
if( (no.cl[[1]] == FALSE) && (mc.cores[[1]] == 1) ){ ## i.e. the user has not initialized the number of cores
mc.cores <- length(cl)
}
## seed
if (!missing(seed)){
set.seed(seed)
}
bootseeds <- sample(1:1000000,size=B,replace=FALSE)
## allArgs <- c("warperBootATE","data","n.obs","fct.pointEstimate",
## "object","treatment","contrasts","times","cause","landmark",
## "n.contrasts","levels","TD","name.estimate","formula","dots")
# }}}
# {{{ warper
warperBootATE <- function(index, args, fct.pointEstimate, name.estimate){
## models for the conditional mean
for(iModel in c("object.event","object.treatment","object.censor")){
if(!is.null(args[[iModel]])){
args[[iModel]]$call$data <- ls.data[[iModel]][index] ## resample dataset
args[[iModel]] <- try(eval(args[[iModel]]$call),silent=TRUE) ## refit model
if (inherits(x=args[[iModel]],what="try-error")){
iBoot <- c(paste0("Failed to fit model ",iModel," on the bootstrap sample", sep = ""),
args[[iModel]]
)
class(iBoot) <- "try-error"
return(iBoot)
}
}
}
## compute ate
if(is.data.frame(args$mydata)){
args$mydata <- args$mydata[index,]
}else{
args$mydata <- args$mydata[index]
}
if(any(args$data.index != 1:NROW(args$mydata))){
args$data.index <- which(index %in% args$data.index)
}
iBoot <- try(do.call(fct.pointEstimate, args), silent = TRUE)
## export
if(inherits(x=iBoot,what="try-error")){ ## error handling
out <- rep(NA, length(name.estimate), name.estimate)
attr(out,"error") <- iBoot
return(out)
}else{
return(c(iBoot$meanRisk$estimate,iBoot$diffRisk$estimate,iBoot$ratioRisk$estimate))
}
}
# }}}
## bootstrap
if(no.cl && mc.cores == 1){
if(verbose){pb <- txtProgressBar(max = B, style = 3,width=30)}
boots <- lapply(1:B, FUN = function(b){
if(verbose>0){setTxtProgressBar(pb, b)}
set.seed(bootseeds[[b]])
warperBootATE(index = sample(1:n.obs, size = n.obs, replace = TRUE),
args = args,
fct.pointEstimate = fct.pointEstimate,
name.estimate = name.estimate)
})
if(verbose>0){close(pb)}
}else if (handler=="foreach"){ #
# {{{ foreach
if(no.cl){
if(verbose>0){
cl <- parallel::makeCluster(mc.cores, outfile = "")
}else{
cl <- parallel::makeCluster(mc.cores)
}
}
doParallel::registerDoParallel(cl)
## progress bar
if(verbose){pb <- txtProgressBar(max = B, style = 3,width=30)}
b <- NULL ## [:forCRANcheck:] foreach
boots <- foreach::`%dopar%`(foreach::foreach(b = 1:B, .packages = add.Package, .export = add.Fct), { ## b <- 1
if(verbose>0){setTxtProgressBar(pb, b)}
set.seed(bootseeds[[b]])
warperBootATE(index = sample(1:n.obs, size = n.obs, replace = TRUE),
args = args,
fct.pointEstimate = fct.pointEstimate,
name.estimate = name.estimate)
})
if(verbose>0){close(pb)}
if(no.cl){parallel::stopCluster(cl)}
# }}}
}else if(handler=="mclapply"){
# {{{ mclapply
boots <- parallel::mclapply(1:B, mc.cores = mc.cores, FUN = function(b){
set.seed(bootseeds[[b]])
warperBootATE(index = sample(1:n.obs, size = n.obs, replace = TRUE),
args = args,
fct.pointEstimate = fct.pointEstimate,
name.estimate = name.estimate)
})
# }}}
}
## output
M.bootEstimate <- do.call(rbind,boots)
if(all(is.na(M.bootEstimate))){
stop(paste0("Error in all bootstrap samples: ", attr(boots[[1]],"error")[1]))
}
colnames(M.bootEstimate) <- name.estimate
return(list(boot = M.bootEstimate,
bootseeds = bootseeds))
}
######################################################################
### ate-bootstrap.R ends here
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.