R/glm_combine.R

In BATSS: Bayesian Adaptive Trial Simulator Software (BATSS) for Generalised Linear Models

#' @name batss.combine
#' @title Combines outputs generated by [batss.glm]
#' @description Combines different evaluations of [batss.glm] considering the same trial design but different sets of seeds. This function is useful when the evaluation of Monte Carlo samples generated by different seeds was split in sets computed by different nodes/cpus. The output of this function is of class 'batss' meaning that the usual generic functions (print, summary, plot) can be used.
#' @param paths Vector indicating the paths to the rdata files containing the outputs of the function [batss.glm] considering the same trial design but different set of seeds. This requires the argument '`extended`' of the function [batss.glm] to be > 0. 
#' @param force a \link[base]{logical} with default `force=FALSE`. Among other checks, [batss.glm] controls that the \link[base]{call}s of the Monte Carlo trials to be combined are \link[base]{identical} and stops if they are not (Note that this check is not bullet proof: such a check, for example, would be able to note that two sets of Monte Carlo trials used a `eff.arm` function named the same way and considered the same optional parameters but would be blind to the fact that they could correspond to two *different* functions). `force=TRUE` forces [batss.glm] to ignore this check. This could be useful if the calls differ due to the `batss` objects to be combined being generated using different versions of [batss.glm].
#' @returns an object of class 'batss'. Refer to the section 'Value' in [batss.glm] for details about this object structure. 
#' @seealso [batss.glm()], the function allowing to simulate Bayesian adaptive trials with GLM endpoint for different seeds. 
#' @export
batss.combine = function(paths, force=FALSE){# paths
    ##
    ## load results 
    ##
    objlist = as.list(rep(NA,length(paths)))
    for(ow in 1:length(paths)){# ow=2
        loaded = load(paths[ow])
        if(length(loaded)==1){
            if(!inherits(get(loaded),"batss")){stop("loaded object not of class 'batss'")}
        }else{
            if(length(loaded)>1){
                class.obj = rep(NA,length(loaded))
                for(oww in 1:length(loaded)){
                    class.obj[oww] = class(get(loaded[oww]))=="batss"
                }
                loaded = loaded[class.obj][1]
            }
        }
        objlist[[ow]]  = get(loaded)
    }
    
    ##
    ## check calls
    ##
    tmp   = lapply(objlist,function(x)as.list(x$call)[names(as.list(x$call))!="R"])
    check = all(sapply(tmp,function(x,ref){identical(x,ref)},ref=tmp[[1]]))
    if(!check){
        if(!force){
            stop("the calls that generated the objects to be combined are different")
        }else{
            warning("the calls that generated the objects to be combined are different")
        }
    }

    ##
    ## combine results
    ##
    out     = objlist[[1]]
    if(is.null(out$call$extended)|(out$call$extended<1)){stop("argument 'extended' of the loaded object is NULL or equal to 0")}
    n.look  = nrow(out$look)
    # seeds
    if(length(objlist)>1){
        for(ow in 2:length(objlist)){
            out$par$seed = c(out$par$seed,objlist[[ow]]$par$seed)
        }
    }
    if(any(table(out$par$seed)>1)){stop("duplicated seeds")}
    # H0 and H1 trials
    Hw      = c("H0","H1")[c(objlist[[1]]$par$H0,objlist[[1]]$par$H1)]
    for(hw in 1:length(Hw)){# hw=1
        n.obj   = length(objlist)
        id.target = objlist[[1]][[Hw[hw]]]$trial[[1]]$target
        estimate  = objlist[[1]][[Hw[hw]]]$estimate
        trial_r   = objlist[[1]][[Hw[hw]]]$trial
        if(length(objlist)>1){
            for(ow in 2:n.obj){
                if(is.null(objlist[[ow]]$call$extended)|(objlist[[ow]]$call$extended<1)){
                        stop("argument 'extended' of the loaded object is NULL or equal to 0")
                    }            
                estimate = abind::abind(estimate,objlist[[ow]][[Hw[hw]]]$estimate,along=3)
                trial_r  = c(trial_r,objlist[[ow]][[Hw[hw]]]$trial)
            }
            }
        tar.p  = batss.res.tp(estimate,id.target)
        tar.g  = batss.res.tg(estimate,id.target)
        eff.p  = batss.res.ep(estimate,id.target,n.look)
        eff.g  = batss.res.eg(estimate,id.target,n.look)
        fut.p  = batss.res.fp(estimate,id.target,n.look)
        fut.g  = batss.res.fg(estimate,id.target,n.look)
        sample = batss.res.s1(trial_r,group=out$par$group$id,
                              type=c(apply(estimate[,"type",,drop=FALSE],2:3,paste0,collapse="")),
                              early=apply(estimate[,"look",,drop=TRUE]<n.look,2,all))
        scenario = batss.res.s2(sample,target=id.target$id)
        out[[Hw[hw]]] = list(estimate = estimate,
                        target   = list(par=tar.p,global=tar.g),
                        efficacy = list(par=eff.p,global=eff.g),
                        futility = list(par=fut.p,global=fut.g),
                        sample=sample,scenario=scenario)
        if(!is.null(out$call$extended)){
            if(out$call$extended>0){
                out[[Hw[hw]]]$trial = trial_r
            }       
        }
    }
    out
}