R/resample.R

In LMMstar: Repeated Measurement Models for Discrete Times

### resample.R --- 
##----------------------------------------------------------------------
## Author: Brice Ozenne
## Created: okt 31 2022 (10:09) 
## Version: 
## Last-Updated: nov  8 2023 (18:12) 
##           By: Brice Ozenne
##     Update #: 566
##----------------------------------------------------------------------
## 
### Commentary: 
## 
### Change Log:
##----------------------------------------------------------------------
## 
### Code:

## * resample (documentation)
##' @title Inference via Resampling for Linear Mixed Model
##' @description Non-parametric bootstrap or permutation test for Linear Mixed Models.
##' @name resample
##'
##' @param object a \code{lmm} object.
##' @param type [character] should permutation test (\code{"perm-var"} or \code{"perm-res"}) or non-parametric bootstrap (\code{"boot"}) be used?
##' @param effects [character vector] the variable(s) to be permuted or the effect(s) to be tested via non-parametric bootstrap.
##' @param n.sample [integer] the number of samples used.
##' @param studentized [logical] should a studentized boostrap or permutation test be used?
##' @param level [numeric,0-1] the confidence level of the confidence intervals.
##' @param trace [logical] should the execution of the resampling be traced?
##' @param seed [integer, >0] Random number generator (RNG) state used when starting resampling.
##' @param correction [logical] correction to ensure non-0 p-values, e.g. with permutations the p.value is evaluated as (#more extreme + 1)/(n.sample + 1) instead of (#more extreme)/(n.sample).
##' @param cpus [integer, >0] number of child-processes for parallel evaluation.
##' If \code{NULL} no state is set.
##' @param export.cpus [character vector] name of the variables to export to each cluster.
##' @param ... Not used. For compatibility with the generic method.
##'
##' @details All approach are carried at the cluster level: \itemize{
##' \item Bootstrap: sampling with replacement clusters. If a cluster is picked twice then different cluster id is used for each pick.
##' \item Permutation: permuting covariate values between clusters (this only lead to the null hypothesis when the covariate values are constant within clusters)
##' or assigning new outcome values by, under the null, permuting the normalized residuals, rescaling to residuals, and adding back the permuted fixed effect (any mean effect under H1 would be 0 because of the permutation if the variance-covariance structure is correct). The later procedure originates from Oliver et al (2012).
##' }
##'
##' The studentized bootstrap keep the original estimate and standard error but uses the samples to evaluates the quantiles of the distribution used to form the confidence intervals.
##' The studentized permutation test approximate the distribution of the test statistic under the null (instead of the distribution of the estimate under the null.).
##'
##' P-values for the bootstrap are computed by centering the bootstrap distribution of the estimate or test statistic around 0 and evaluating the frequency at which it takes values more extreme than the observed estimate or test statistics.
##' 
##' @references Oliver E. Lee and Thomas M. Braun (2012), \bold{Permutation Tests for Random Effects in Linear Mixed Models}. \emph{Biometrics}, Journal 68(2).
##'
##' @keywords htest
##' 
##' @examples
##' \dontrun{
##'
##' #### univariate linear regression ####
##' data(gastricbypassW, package = "LMMstar")
##' ## rescale to ease optimization
##' gastricbypassW$weight1 <- scale(gastricbypassW$weight1)
##' gastricbypassW$weight2 <- scale(gastricbypassW$weight2)
##' gastricbypassW$glucagonAUC1 <- scale(gastricbypassW$glucagonAUC1)
##'
##' e.lm <- lmm(weight2~weight1+glucagonAUC1, data = gastricbypassW)
##' model.tables(e.lm)
##'
##' ## non-parametric bootstrap
##' resample(e.lm, type = "boot", effects = c("weight1","glucagonAUC1"), seed = 10)
##' ## permutation test
##' resample(e.lm, type = "perm-var", effects = "weight1", seed = 10) 
##' resample(e.lm, type = "perm-var", effects = "glucagonAUC1", seed = 10)
##' ## using multiple cores
##' resample(e.lm, type = "boot", effects = c("weight1","glucagonAUC1"), cpus = 4)
##' 
##' #### random intercept model ####
##' data(gastricbypassL, package = "LMMstar")
##' gastricbypassL$weight <- scale(gastricbypassL$weight)
##' gastricbypassL$glucagonAUC <- scale(gastricbypassL$glucagonAUC)
##' gastricbypassL$gender <- as.numeric(gastricbypassL$id) %% 2
##' gastricbypassLR <- na.omit(gastricbypassL)
##' 
##' eCS.lmm <- lmm(weight~glucagonAUC+gender, data = gastricbypassLR,
##'                repetition = ~visit|id, structure = "CS")
##' model.tables(eCS.lmm)
##'
##' ## non-parametric bootstrap
##' resample(eCS.lmm, type = "boot", effects = c("glucagonAUC","gender"), seed = 10, trace = FALSE)
##' ## permutation test
##' resample(eCS.lmm, type = "perm-var", effects = "gender", seed = 10)
##' resample(eCS.lmm, type = "perm-res", effects = "glucagonAUC", seed = 10) 
##' }
##' 
#' @export
`resample` <-
  function(object, type, ...) UseMethod("resample")


## * resample.lmm (code)
##' @export
##' @rdname resample
resample.lmm <- function(object, type, effects, n.sample = 1e3, studentized = TRUE,
                         level = 0.95, correction = TRUE,
                         trace = TRUE, seed = NULL, cpus = 1, export.cpus = NULL,
                         ...){

    ## ** check user input
    alpha <- 1-level
    type <- match.arg(type, c("perm-var","perm-res","boot"))
    ## boot: non-parametric bootstrap
    ## perm-var: exchange covariate value between clusters (need to be constant within cluster)
    ## perm-res: outcome becomes permuted normalized residuals under the null to which the (permuted) fixed effect under the null are added.
    if(!is.logical(correction) || correction %in% 0:1 == FALSE){
        stop("Argument \'correction\' must be binary or logical. \n")
    }

    name.meanvar <- attr(object$design$mean,"variable")
    value.meancoef <- coef(object, effects = "mean")
    sd.meancoef <- sqrt(diag(vcov(object, effects = "mean")))

    name.meancoef <- names(value.meancoef)
    if(type == "perm-var"){
        if(any(is.character(effects)==FALSE)){
            stop("Argument \'effects\' should contain character strings refering to a variable for the mean structure. \n")
        }
        if(length(effects)==0){
            stop("Argument \'effects\' should have length at least 1. \n")
        }
        if(any(effects %in% name.meanvar == FALSE)){
            stop("Argument \'effects\' should be one of \"",paste(name.meanvar, collapse = "\", \""),"\". \n")
        }

        M.factor <- attr(stats::terms(object$formula$mean.design),"factor")
        colnames(M.factor) <- setdiff(name.meancoef,"(Intercept)")

        index.keepcoef <- which(colSums(M.factor[effects,,drop=FALSE]!=0)>0)
        name.keepcoef <- names(index.keepcoef)
        effects.vcov <- any(effects %in% manifest(object, effects = c("variance","correlation")))
        
        ## test whether the covariate is constant within cluster
        test.Wvar <- unlist(by(object$data[effects], object$data[[object$cluster$var]], function(iData){sum(!duplicated(iData))>1}, simplify = FALSE))
        if(any(test.Wvar)){
            stop("The covariate(s) indicated by the argument \"effects\" vary within cluster. \n",
                 "This is not support when using permutations. Consider using type=\"boot\" instead of type=\"perm\". \n")
        }
        Uvar <- by(object$data[effects], object$data[[object$cluster$var]], function(iData){iData[1,,drop=FALSE]}, simplify = FALSE)

    }else if(type == "perm-res"){
        if(any(is.character(effects)==FALSE)){
            stop("Argument \'effects\' should contain character strings refering to a parameter of the mean structure. \n")
        }
        if(any(effects %in% name.meancoef == FALSE)){
            stop("Argument \'effects\' should be one of \"",paste(name.meancoef, collapse = "\", \""),"\". \n")
        }
        effects.vcov <- FALSE
        name.keepcoef <- effects
    }else if(type == "boot"){
        if(any(is.character(effects)==FALSE)){
            stop("Argument \'effects\' should contain character strings refering to a parameter of the mean structure. \n")
        }
        if(any(effects %in% name.meancoef == FALSE)){
            stop("Argument \'effects\' should be one of \"",paste(name.meancoef, collapse = "\", \""),"\". \n")
        }
        effects.vcov <- TRUE
        name.keepcoef <- effects
    }

    if(length(cpus)!=1){
        stop("Argument \'cpus\' should have length 1.\n ")
    }else if(identical(cpus,"all")){
        cpus <- parallel::detectCores()
    }else if(!is.numeric(cpus) || cpus <=0 || cpus %% 1 != 0){
        stop("Argument \'cpus\' should be a positive integer or \'all\'.\n ")
    }else if(cpus>1 && cpus>parallel::detectCores()){
        stop("Only ",parallel::detectCores()," CPU cores are available. \n")
    }

    ## dots
    dots <- list(...)
    if(length(dots)>0){
        stop("Unknown argument(s) \'",paste(names(dots),collapse="\' \'"),"\'. \n")
    }

    ## ** initialize
    ## *** data
    data <- object$data
    if(length(object$index.na)>0){
        data <- data[-object$index.na,]
    }
    var.cluster <- object$cluster$var
    vec.Uid <- object$cluster$level
    n.cluster <- object$cluster$n
    n.obs <- NROW(data)
    index.cluster <- object$design$index.cluster
    Vindex.cluster <- unlist(index.cluster)
    for(iCluster in names(index.cluster)){
        names(index.cluster[[iCluster]]) <- rep(iCluster, length(index.cluster[[iCluster]]))
    }
    precompute.moments <- !is.null(object$design$precompute.XX)
  
    var.all <- manifest(object, effects = "all")
    var.mean <- manifest(object, effects = "mean")
    XX.all <- c("XXindexXX", "XXclusterXX", "XXcluster.indexXX", "XXtimeXX", "XXtime.indexXX", "XXstrataXX", "XXstrata.indexXX")

    var.weights <- object$weights$var
    var.outcome <- object$outcome$var

    ## *** parameters
    if(inherits(object$design$vcov,"ID") || inherits(object$design$vcov,"IND")){
        param.init <- NULL ## OLS solver + empirical residual standard deviation can directly find the solution
    }else{
        param.init <- object$param
        if(type == "perm-var"){
            param.init[name.keepcoef] <- 0
        }else if(type == "perm-res"){
            param.init[effects] <- 0
        }
    }
    ## ## *** missing patterns
    ## if(type == "perm-var"){
    ##     Upattern <- object$design$vcov$Upattern
    ##     ## find patterns sharing the same times 
    ##     ls.time.pattern <- unique(Upattern$time)  
    ##     time.pattern <- lapply(ls.time.pattern, paste, collapse="|")
    ##     n.patterns <- length(time.pattern)

    ##     ## gather clusters with patterns sharing the same times 
    ##     match.tp <- match(lapply(Upattern$time,paste,collapse="|"), time.pattern)
    ##     ls.indexcluster.pattern <- as.list(by(Upattern, match.tp, function(iData){unlist(iData$index.cluster)}, simplify = FALSE))
    ##     n.idpatterns <- lapply(ls.indexcluster.pattern,length)
    ## }

    ## *** residuals
    if(type == "perm-res"){
        ## re-estimate the model under the null hypothesis
        call0 <- object$call
        call0$formula <- stats::update(eval(call0$formula), paste0("~.-",effects))
        if(!is.null(object$index.na)){
            call0$data <- object$data.original[-object$index.na,,drop=FALSE]
        }
        object0 <- eval(call0)
        Xbeta0 <- stats::predict(object0, newdata = data[var.mean], se = FALSE, df = FALSE, simplify = TRUE)
        epsilon0.norm <- stats::residuals(object0, type = "normalized")
        OmegaChol0 <- lapply(stats::sigma(object0, chol = TRUE, cluster = as.character(vec.Uid)), FUN = base::t)
    }

    ## ** function
    warperResample <- function(iSample){

        ## *** resample
        if(type == "perm-var"){
            ## permute X-values (constant within cluster)
            iPerm <- sample(n.cluster, replace = FALSE)
            iData <- data
            iData[unlist(index.cluster),effects] <- unlist(mapply(x = index.cluster, y = Uvar[iPerm], FUN = function(x,y){rep(y,length(x))}))
                        
            ## permute X-values between individuals with same missing data pattern
            ## iData <- data
            ## for(iPattern in 1:n.patterns){ ## iPattern <- 1
            ##     iPerm <- sample(n.idpatterns[[iPattern]], replace = FALSE)
            ##     iPatternindex.cluster <-  index.cluster[ls.indexcluster.pattern[[iPattern]]]
            ##     iData[[effects]][unlist(iPatternindex.cluster)] <- data[[effects]][unlist(iPatternindex.cluster[iPerm])]
            ## }
        }else if(type == "perm-res"){            
            iData <- data
            ## permute residuals (including the fixed effect to be tested)
            iPerm <- sample(1:n.obs, replace = FALSE)
            iEpsilon0.norm <- epsilon0.norm[iPerm]
            ## rescale residuals
            iData[[var.outcome]][unlist(index.cluster)] <- unlist(lapply(1:n.cluster, function(iC){OmegaChol0[[iC]] %*% iEpsilon0.norm[index.cluster[[iC]]]}))
            ## add fixed effects
            iData[[var.outcome]] <- iData[[var.outcome]] + Xbeta0[iPerm]
            ## range(iData[[var.outcome]] - data[,var.outcome])
        }else if(type == "boot"){
            ils.Boot <- index.cluster[sample(n.cluster, replace = TRUE)]
            iBoot <- unlist(ils.Boot)
            
            iData <- data[iBoot,,drop=FALSE]
            iData[[var.cluster]] <- unlist(lapply(1:n.cluster, function(iC){rep(iC,length(ils.Boot[[iC]]))}))
            rownames(iData) <- NULL
            ## lmm(Y~X1*X2+X5,data = iData[,c("Y","X1","X2","X5")])
        }

        ## *** update design
        if(effects.vcov){
            ## update design according to the permutation (mean and variance)
            iStructure <- object$design$vcov
            iStructure$var <- NULL
            iStructure$cor <- NULL
            iStructure$param <- NULL
            iStructure$pattern <- NULL
            iStructure$Upattern <- NULL
            iStructure$pair.vcov <- NULL
            iStructure$pair.meanvcov <- NULL
            
            iData2 <- .lmmNormalizeData(iData[,var.all,drop=FALSE],
                                        var.outcome = var.outcome,
                                        var.cluster = attr(var.cluster, "original"),
                                        var.time = attr(object$time$var, "original"),
                                        var.strata = attr(object$strata$var, "original"),
                                        droplevels = TRUE,
                                        initialize.cluster = iStructure$ranef$crossed,
                                        initialize.time = setdiff(iStructure$ranef$vars, iStructure$var.cluster))$data

            iDesign <- .model.matrix.lmm(formula.mean = object$formula$mean.design,
                                         structure = iStructure,
                                         data = iData2,
                                         var.outcome = object$outcome$var,
                                         var.weights = object$weights$var,
                                         precompute.moments = precompute.moments,
                                         drop.X = object$design$drop.X)

        }else if(type == "perm-res"){ ## change in the Y values

            iDesign <- object$design
            iDesign$Y <- iData[[var.outcome]]

            ## update pre-computation              
            if(precompute.moments && NCOL(iDesign$mean)>0){
                if(is.na(var.weights[1])){
                    iwY <- cbind(iData[[var.outcome]])
                }else{
                    iwY <- cbind(iData[[var.outcome]]*sqrt(iData[[var.weights[1]]]))
                }
                iDesign$precompute.XY <- .precomputeXR(X = iDesign$precompute.XX$Xpattern, residuals = iwY, pattern = iDesign$vcov$Upattern$name,
                                                       pattern.ntime = stats::setNames(iDesign$vcov$Upattern$n.time, iDesign$vcov$Upattern$name),
                                                       pattern.cluster = attr(iDesign$vcov$pattern,"list"), index.cluster = iDesign$index.cluster)
            }

        }else{ ## change in the X values
            ## update design matrix according to the permutation (only mean)
            iDesign <- object$design
            iDesign$mean <- model.matrix(object, data = iData[,var.mean,drop=FALSE])
            attr(iDesign$mean, "assign") <- attr(object$design$mean, "assign")
            attr(iDesign$mean, "contrasts") <- attr(object$design$mean, "contrasts")
            attr(iDesign$mean, "variable") <- attr(object$design$mean, "variable")
            attr(iDesign$mean, "terms") <- attr(object$design$mean, "terms")

            ## update pre-computation                
            if(precompute.moments && NCOL(iDesign$mean)>0){
                if(is.na(var.weights[1])){
                    iwX.mean <- iDesign$mean
                    iwY <- cbind(iData[[var.outcome]])
                }else{
                    iwX.mean <- sweep(iDesign$mean, FUN = "*", MARGIN = 1, STATS = sqrt(iData[[var.weights[1]]]))
                    iwY <- cbind(iData[[var.outcome]]*sqrt(iData[[var.weights[1]]]))
                }
                iIndex.cluster <- .extractIndexData(data = iData, structure = iDesign$vcov)$index.cluster
                iDesign$precompute.XX <- .precomputeXX(X = iwX.mean, pattern = iDesign$vcov$Upattern$name, 
                                                       pattern.ntime = stats::setNames(iDesign$vcov$Upattern$n.time, iDesign$vcov$Upattern$name),
                                                       pattern.cluster = attr(iDesign$vcov$pattern,"list"), index.cluster = iIndex.cluster)
                iDesign$precompute.XY <- .precomputeXR(X = iDesign$precompute.XX$Xpattern, residuals = iwY, pattern = iDesign$vcov$Upattern$name,
                                                       pattern.ntime = stats::setNames(iDesign$vcov$Upattern$n.time, iDesign$vcov$Upattern$name),
                                                       pattern.cluster = attr(iDesign$vcov$pattern,"list"), index.cluster = iIndex.cluster)
            }
        }

        ## *** re-estimate
        iEstimate <- try(.estimate(design = iDesign,
                                   time = object$time,
                                   method.fit = object$args$method.fit,
                                   type.information = object$args$type.information,
                                   transform.sigma = object$reparametrize$transform.sigma,
                                   transform.k = object$reparametrize$transform.k,
                                   transform.rho = object$reparametrize$transform.rho,
                                   precompute.moments = precompute.moments, 
                                   optimizer = object$args$control$optimizer,
                                   init = param.init,
                                   n.iter = object$opt$control["n.iter"],
                                   tol.score = object$opt$control["tol.score"],
                                   tol.param = object$opt$control["tol.param"],
                                   n.backtracking = object$opt$control["n.backtracking"],
                                   trace = FALSE))

        if(!inherits(iEstimate,"try-error") && studentized){
            iVcov <- .moments.lmm(value = iEstimate$estimate,
                                  design = iDesign,
                                  time = object$time,
                                  method.fit = object$args$method.fit,
                                  type.information = object$args$type.information,
                                  transform.sigma = object$reparametrize$transform.sigma,
                                  transform.k = object$reparametrize$transform.k,
                                  transform.rho = object$reparametrize$transform.rho,
                                  logLik = FALSE, score = FALSE, information = FALSE, vcov = TRUE, df = FALSE, indiv = FALSE, effects = c("mean"), robust = FALSE,
                                  trace = FALSE, precompute.moments = precompute.moments, method.numDeriv = "simple", transform.names = FALSE)$vcov
            
        }

        ## *** export
        if(inherits(iEstimate,"try-error")){
            return(iEstimate)
        }else if(!studentized){
            return(c(iEstimate$cv, iEstimate$estimate[name.keepcoef]))
        }else if(studentized){
            c(iEstimate$cv, iEstimate$estimate[name.keepcoef], se = sqrt(diag(iVcov)[name.keepcoef]))
        }

    }
    ## res <- warperResample(1)
    ## ** run
    if(trace){
        if(type == "perm-var"){
            if(studentized){
                cat("\tStudentized permutation test with ",n.sample," replicates \n",
                    "\t(permutation of the regressor values between clusters)\n\n", sep = "")
            }else{
                cat("\tPermutation test with ",n.sample," replicates \n",
                    "\t(permutation of the regressor values between clusters)\n\n", sep = "")
            }
        }else if(type == "perm-res"){
            if(studentized){
                cat("\tStudentized permutation test with ",n.sample," replicates \n",
                    "\t(permutation of the normalized residuals under the null)\n\n", sep = "")
            }else{
                cat("\tPermutation test with ",n.sample," replicates \n",
                    "\t(permutation of the normalized residuals under the null)\n\n", sep = "")
            }
        }else if(type == "boot"){
            if(studentized){
                cat("\tNon-parametric studentized bootstrap with ",n.sample," replicates \n\n", sep = "")
            }else{
                cat("\tNon-parametric bootstrap with ",n.sample," replicates \n\n", sep = "")
            }
        }
    }
    
    if(!is.null(seed)){
        tol.seed <- 10^(floor(log10(.Machine$integer.max))-1)
        if(n.sample>tol.seed){
            stop("Cannot set a seed per simulation when considering more than ",tol.seed," similations. \n")
        }
        if(!is.null(get0(".Random.seed"))){ ## avoid error when .Random.seed do not exists, e.g. fresh R session with no call to RNG
            old <- .Random.seed # to save the current seed
            on.exit(.Random.seed <<- old) # restore the current seed (before the call to the function)
        }else{
            on.exit(rm(.Random.seed, envir=.GlobalEnv))
        }
        set.seed(seed)
        test.seed <- TRUE
        seqSeed <- sample.int(tol.seed, n.sample,  replace = FALSE)        
    }else{
        test.seed <- FALSE
        seqSeed <- NULL
    }

    if(cpus==1){
        if (trace > 0) {
            requireNamespace("pbapply")
            method.loop <- pbapply::pblapply
        }else{
            method.loop <- lapply
        }
    }else if(cpus>1){
        cl <- parallel::makeCluster(cpus)
        ## link to foreach
        doParallel::registerDoParallel(cl)
        ## export from user
        if(!is.null(export.cpus)){
            parallel::clusterExport(cl, export.cpus)
        } 
        ## export seed 
        if (test.seed) {
            parallel::clusterExport(cl, varlist = "seqSeed", envir = environment())
        }
        ## export BuyseTest 
        fct2export <- c(".estimate",".precomputeXR",".precomputeXX",".extractIndexData",".model.matrix.lmm",".lmmNormalizeData",".moments.lmm")
        parallel::clusterExport(cl = cl, varlist = fct2export, envir = as.environment(asNamespace("LMMstar")))
        
    }

    if (cpus == 1) {
        ls.sample <- do.call(method.loop,
                             args = list(X = 1:n.sample,
                                         FUN = function(iX){
                                             if(test.seed){set.seed(seqSeed[iX])}
                                             iOut <- warperResample(iX)
                                             return(iOut)
                                         })
                             )
    }else if(cpus > 1){
        if(trace>0){
            pb <- utils::txtProgressBar(max = n.sample, style = 3)          
            progress <- function(n){utils::setTxtProgressBar(pb, n)}
            opts <- list(progress = progress)
        }else{
            opts <- list()
        }

        iX <- NULL ## [:forCRANcheck:] foreach        
        ls.sample <- foreach::`%dopar%`(
                                  foreach::foreach(iX=1:n.sample, .options.snow = opts, .packages = c("LMMstar","nlme")), {
                                      if(test.seed){set.seed(seqSeed[iX])}
                                      iOut <- warperResample(iX)
                                      return(iOut)
                                  })
        if(trace>0){close(pb)}
        parallel::stopCluster(cl)
    }


    ## ** post-process
    M.sample <- do.call(rbind,ls.sample[sapply(ls.sample, inherits, "try-error")==FALSE])
    Mcv.sample <- M.sample[M.sample[,1]==1,-1,drop=FALSE]
    ncv.sample <- NROW(Mcv.sample)

    if(studentized){
        Mcv.sample.se <- Mcv.sample[,(length(name.keepcoef)+1):(2*length(name.keepcoef)),drop=FALSE]
        Mcv.sample <- Mcv.sample[,1:length(name.keepcoef),drop=FALSE]
    }
    
    out <- model.tables(object)[name.keepcoef,,drop=FALSE]*NA
    out$estimate <- as.double(value.meancoef[name.keepcoef])

    if(type %in% c("perm-var","perm-res") ){

        if(studentized){
            Mcv.estimate <- matrix(value.meancoef[name.keepcoef]/sd.meancoef[name.keepcoef], nrow = ncv.sample, ncol = length(name.keepcoef), byrow = TRUE,
                                   dimnames = list(NULL,name.keepcoef))
            Mcv.sample.Wald <- Mcv.sample/Mcv.sample.se

            out$se <- as.double(sd.meancoef[name.keepcoef])
            out$p.value <- (colSums(abs(Mcv.sample.Wald) > abs(Mcv.estimate), na.rm = TRUE)+correction)/(colSums(!is.na(Mcv.sample.Wald))+correction)
        }else{
            Mcv.estimate <- matrix(value.meancoef[name.keepcoef], nrow = ncv.sample, ncol = length(name.keepcoef), byrow = TRUE,
                                   dimnames = list(NULL,name.keepcoef))

            out$p.value <- (colSums(abs(Mcv.sample) > abs(Mcv.estimate), na.rm = TRUE)+correction)/(colSums(!is.na(Mcv.sample))+correction)
        }
    }else if(type == "boot"){
        if(studentized){
            Mcv.sample.Wald0 <- apply(Mcv.sample, MARGIN = 2, FUN = scale, scale = FALSE, center = TRUE)/Mcv.sample.se  ## center around the null
            
            out$se <- as.double(sd.meancoef[name.keepcoef])
            out$lower <- out$estimate + out$se * apply(Mcv.sample.Wald0, MARGIN = 2, FUN = stats::quantile, probs = alpha/2, na.rm = TRUE)
            out$upper <- out$estimate + out$se * apply(Mcv.sample.Wald0, MARGIN = 2, FUN = stats::quantile, probs = 1-alpha/2, na.rm = TRUE)
            out$p.value <- sapply(name.keepcoef, function(iName){
                boot2pvalue(stats::na.omit(out[iName,"estimate"] + out[iName,"se"] * Mcv.sample.Wald0[,iName]),
                            null = 0,
                            estimate = out[iName,"estimate"],
                            alternative = "two.sided",
                            add.1 = correction)
            })
            
        }else{
            Mcv.estimate <- matrix(value.meancoef[name.keepcoef], nrow = ncv.sample, ncol = length(name.keepcoef), byrow = TRUE,
                                   dimnames = list(NULL,name.keepcoef))

            out$se <- apply(Mcv.sample, MARGIN = 2, FUN = stats::sd, na.rm = TRUE)
            out$lower <- apply(Mcv.sample, MARGIN = 2, FUN = stats::quantile, probs = alpha/2, na.rm = TRUE)
            out$upper <- apply(Mcv.sample, MARGIN = 2, FUN = stats::quantile, probs = 1-alpha/2, na.rm = TRUE)
            out$p.value <- sapply(name.keepcoef, function(iName){
                boot2pvalue(stats::na.omit(Mcv.sample[,iName]),
                            null = 0,
                            estimate = out[iName,"estimate"],
                            alternative = "two.sided",
                            add.1 = correction)
            })
        }
    }

    ## ** export
    attr(out,"call") <- match.call()
    attr(out,"args") <- list(type = type, effects = effects, n.sample = n.sample, studentized = studentized, seed = seqSeed)
    attr(out,"M.sample") <- M.sample
    attr(out,"n.sample") <- n.sample
    class(out) <- append("resample",class(out))
    return(out)

}

## ## * resample.mlmm
## ##' @export
## ##' @rdname resample
## resample.mlmm <- function(object, type, method = NULL, cluster = NULL, n.sample = 1e3,
##                           level = 0.95, correction = TRUE,
##                           trace = TRUE, seed = NULL, cpus = 1,
##                           ...){

##     ## ** normalize arguments
##     ## alternative
##     if(object$glht[[1]][[1]]$alternative!="two.sided"){
##         stop("Can only perform two sided tests. \n")
##     }
    
##     ## n.sample
##     if(length(n.sample)!=1){
##         stop("Argument \'n.sample\' should have lenght 1. \n")
##     }
##     if(!is.numeric(n.sample) || n.sample<0 || n.sample %% 1 >0 ){
##         stop("Argument \'n.sample\' should be a positive integer. \n")
##     }
    
##     ## cluster
##     if(is.null(cluster)){
##         ls.nameCluster <- lapply(object$model, function(iModel){attr(manifest(iModel),"cluster")})
##         if(length(unique(unlist(ls.nameCluster)))==1 && all(lengths(ls.nameCluster)==1)){
##             cluster <- unname(ls.nameCluster[[1]])
##         }else{
##             stop("Argument \'cluster\' cannot be guessed from the object. \n")
##         }
        
##     }
##     if(length(cluster)!=1){
##         stop("Argument \'cluster\' should have lenght 1. \n")
##     }
##     if(!is.character(cluster)){
##         stop("Argument \'cluster\' should be a character. \n")
##     }

##     ## data    
##     object.call <- attr(object, "call")
##     by <- object.call$by
##     level.by <- names(object$model)
##     data <- try(eval(object.call$data), silent = TRUE)

##     if(inherits(data, "try-error")){ ## reconstruct dataset based on what was stored in the object
##         object.model <- object$model
##         n.model <- length(object.model)
##         object.manifest <- lapply(object.model, manifest)
##         object.data <- lapply(object.model, "[[", "data")

##         if(any(sapply(object.data, function(iData){cluster %in% names(iData)})==FALSE)){
##             stop("Argument \'cluster\' could not be found in the dataset(s) stored in the object. \n")
##         }
##         object.Udata <- mapply(FUN = function(iData,iName,iBy){
##             iOut <- cbind(iBy,iData[c(cluster,iName)])
##             names(iOut)[1] <- by
##             return(iOut)
##         }, iData = object.data, iName = object.manifest, iBy = level.by, SIMPLIFY = FALSE)
##         data <- object.Udata[[1]]

##         if(n.model>1){
##             for(iModel in 2:n.model){ ## iModel <- 3
##                 data <- merge(data, object.Udata[[iModel]], by = intersect(names(data),c(object.manifest[[iModel]],by,cluster)), all=TRUE)
##             }
##         }
    
##     }else{
##         if(cluster %in% names(data)==FALSE){
##             stop("Argument \'cluster\' could not be found in the dataset. \n")
##         }
##     }
##     data <- as.data.frame(data)

##     ## type
##     if(tolower(type) == "perm"){
##         type <- "permutation"
##     }else if(tolower(type) == "boot"){
##         type <- "bootstrap"
##     }
##     type <- match.arg(type, c("permutation","bootstrap"))
##     if(type=="permutation"){

##         effect.param <- object$univariate$parameter
##         if(any(object$univariate$type!="mu")){
##             stop("Can only test mean parameters when using permutations. \n")
##         }
##         ls.variable.perm <- mapply(FUN = function(iParam,iObject){ ## iParam <- "X2b" ## iObject <- object$model[[1]]
##             iX <- iObject$design$mean 
##             attr(iX,"variable")[attr(iX,"assign")[colnames(iX)==iParam]]
##         }, iParam = effect.param, iObject = object$model, SIMPLIFY = FALSE)
##         variable.perm <- unname(unique(unlist(ls.variable.perm)))
##         if(length(variable.perm)>1){
##             stop("Permutation test not available for effects depending on more than one variable. \n",
##                  "Can only permute a single variable. \n")
##         }

##         ## test whether the covariate is constant within cluster
##         test.Wvar <- unlist(by(data[[variable.perm]], data[[cluster]], function(iData){sum(!duplicated(iData))>1}, simplify = FALSE))
##         if(any(test.Wvar)){
##             stop("The covariate(s) indicated by the argument \"effects\" vary within cluster. \n",
##                  "This is not support when using permutations. Consider using type=\"boot\" instead of type=\"perm\". \n")
##         }
##         Uvar <- by(data[[variable.perm]], data[[cluster]], function(iData){iData[1]}, simplify = FALSE)
##     }else if(type=="bootstrap"){
##         calc.pvalue <- try(requireNamespace("BuyseTest")) ## to get the p-values
##     }

##     ## cpus
##     if(length(cpus)!=1){
##         stop("Argument \'cpus\' should have length 1.\n ")
##     }else if(identical(cpus,"all")){
##         cpus <- parallel::detectCores()
##     }else if(!is.numeric(cpus) || cpus <=0 || cpus %% 1 != 0){
##         stop("Argument \'cpus\' should be a positive integer or \'all\'.\n ")
##     }else if(cpus>1 && cpus>parallel::detectCores()){
##         stop("Only ",parallel::detectCores()," CPU cores are available. \n")
##     }

##     ## dots
##     dots <- list(...)
##     if(length(dots)>0){
##         stop("Unknown argument(s) \'",paste(names(dots),collapse="\' \'"),"\'. \n")
##     }

##     ## ** recover dataset
##     method.merge <- c("average", "pool.fixse", "pool.se", "pool.gls", "pool.gls1", "pool.rubin")
##     keep.cols <- c("by","parameter","estimate","se","df","statistic","lower","upper","null","p.value")
        
##     Ucluster <- sort(unique(unlist(data[[cluster]])))
##     n.cluster <- length(Ucluster)
##     index.cluster <- split(1:NROW(data), data[[cluster]], identity)
##     template <- list(multivariate = object$multivariate,
##                      univariate = object$univariate)
##     template$multivariate[,c("statistic","df.num","df.denom","p.value")] <- NULL
##     template$univariate[,c("estimate","se","df","statistic","lower","upper","p.value")] <- NULL
                     
##     warperResample <- function(iSample){

##         ## *** resample
##         if(type == "permutation"){
##             ## permute X-values (constant within cluster)
##             iPerm <- sample(n.cluster, replace = FALSE)
##             iData <- data
##             for(iCluster in 1:n.cluster){ ## iCluster <- 1
##                 iData[index.cluster[[iCluster]],variable.perm] <- Uvar[[iPerm[iCluster]]] ## seems to properly expand Uvar over multiple timepoints
##             }
##         }else if(type == "bootstrap"){
##             browser()
##             ils.Boot <- index.cluster[sample(n.cluster, replace = TRUE)]
##             iBoot <- unlist(ils.Boot)
##             iBoot.cl <- unlist(mapply(x = 1:n.cluster, y = lengths(ils.Boot), function(x,y){rep(x,y)}, SIMPLIFY = FALSE))

##             iData <- data[iBoot,,drop=FALSE]
##             iData[[cluster]] <- iBoot.cl
##             rownames(iData) <- NULL            
##         }

##         ## *** re-estimate
##         iCall <- object.call
##         iCall$trace <- 0
##         iCall$data <- iData
##         iEstimate <- eval(iCall)

##         ## *** export
##         ## use merge in case effects for some categories are not estimated
##         iEstimate$univariate2 <- model.tables(iEstimate, columns = keep.cols, method = method)

##         iOut <- template
##         if(is.null(method) || all(method %in% method.merge == FALSE)){
##             iOut$multivariate <- merge(iOut$multivariate, iEstimate$multivariate, by = c("type","test","null"), all = TRUE)
##             iOut$univariate <- merge(iOut$univariate, iEstimate$univariate2, by = c("by","parameter"), all = TRUE)
##         }else{
##             iOut$univariate <- iEstimate$univariate2
##         }
##         return(iOut)
##     }

##     ## ** run
##     if(trace){
##         if(type == "permutation"){
##             cat("\tPermutation test with ",n.sample," replicates \n",
##                 "\t(permutation of ",variable.perm," values between clusters)\n\n", sep = "")
##         }else if(type == "boot"){
##             cat("\tNon-parametric bootstrap with ",n.sample," replicates \n\n", sep = "")
##         }
##     }
    
##     ## seed
##     if(!is.null(seed)){
##         tol.seed <- 10^(floor(log10(.Machine$integer.max))-1)
##         if(n.sample>tol.seed){
##             stop("Cannot set a seed per simulation when considering more than ",tol.seed," similations. \n")
##         }
##         if(!is.null(get0(".Random.seed"))){ ## avoid error when .Random.seed do not exists, e.g. fresh R session with no call to RNG
##             old <- .Random.seed # to save the current seed
##             on.exit(.Random.seed <<- old) # restore the current seed (before the call to the function)
##         }else{
##             on.exit(rm(.Random.seed, envir=.GlobalEnv))
##         }
##         set.seed(seed)
##         test.seed <- TRUE
##         seqSeed <- sample.int(tol.seed, n.sample,  replace = FALSE)        
##     }else{
##         seqSeed <- NULL
##         test.seed <- FALSE
##     }

##     if(cpus==1){
##         if (trace > 0) {
##             requireNamespace("pbapply")
##             method.loop <- pbapply::pblapply
##         }else{
##             method.loop <- lapply
##         }
##     }else if(cpus>1){
##         cl <- parallel::makeCluster(cpus)
##         ## link to foreach
##         doParallel::registerDoParallel(cl)
##         ## export from user
##         if(!is.null(export.cpus)){
##             parallel::clusterExport(cl, export.cpus)
##         } 
##         ## export seed 
##         if (test.seed) {
##             parallel::clusterExport(cl, varlist = "seqSeed", envir = environment())
##         }
##         ## export BuyseTest 
##         fct2export <- NULL
##         parallel::clusterExport(cl = cl, varlist = fct2export, envir = as.environment(asNamespace("LMMstar")))
        
##     }
##     if (cpus == 1) {
##         ls.sample <- do.call(method.loop,
##                              args = list(X = 1:n.sample,
##                                          FUN = function(iX){ ## iX <- 1
##                                              if(test.seed){set.seed(seqSeed[iX])}
##                                              iOut <- warperResample(iX)
##                                              return(iOut)
##                                          })
##                              )
##     }else if(cpus > 1){
##         if(trace>0){
##             pb <- utils::txtProgressBar(max = n.sample, style = 3)          
##             progress <- function(n){utils::setTxtProgressBar(pb, n)}
##             opts <- list(progress = progress)
##         }else{
##             opts <- list()
##         }

##         iX <- NULL ## [:forCRANcheck:] foreach        
##         ls.sample <- foreach::`%dopar%`(
##                                   foreach::foreach(iX=1:n.sample, .options.snow = opts, .packages = c("LMMstar","nlme")), {
##                                       if(test.seed){set.seed(seqSeed[iX])}
##                                       iOut <- warperResample(iX)
##                                       return(iOut)
##                                   })
##         if(trace>0){close(pb)}
##         parallel::stopCluster(cl)
##     }
    
##     ## ** post process
##     keep.cols2 <- c("estimate", "se", "df", "statistic", "lower", "upper", "null", "p.value")
    
##     out <- list()
##     sample.multivariate <- do.call(cbind,lapply(ls.sample, function(iSample){iSample$multivariate$statistic}))
##     sample.univariate <- do.call(cbind,lapply(ls.sample, function(iSample){iSample$univariate$statistic}))

##     if(type == "permutation"){
##         obs.multivariate <- matrix(object$multivariate[,"statistic"], nrow = NROW(object$multivariate), ncol = NCOL(sample.multivariate))
##         obs.univariate <- matrix(object$univariate[,"statistic"], nrow = NROW(object$univariate), ncol = NCOL(sample.univariate), byrow = FALSE)
        
##         if(is.null(method) || all(method %in% method.merge == FALSE)){
##             out$multivariate <- object$multivariate
##             out$multivariate$df.num <- NA
##             out$multivariate$df.denom <- NA
##             out$multivariate$p.value <- (rowSums(abs(sample.multivariate) > abs(obs.multivariate), na.rm = TRUE)+1)/(rowSums(!is.na(sample.multivariate))+1)
##             attr(out$univariate,"sample") <- sample.multivariate
##         }

##         if(is.null(method) || all(method %in% method.merge == FALSE)){
##             out$univariate <- object$univariate
##         }else{
##             out$univariate <- model.tables(object, method = method, columns = keep.cols2)
##         }
##         out$univariate$se <- NA
##         out$univariate$df <- NA
##         out$univariate$p.value <- (rowSums(abs(sample.univariate) > abs(obs.univariate), na.rm = TRUE)+1)/(rowSums(!is.na(sample.univariate))+1)
##         attr(out$univariate,"sample") <- sample.univariate

##     }else if(type == "bootstrap"){
##         alpha <- 1-level

##         if(is.null(method) || all(method %in% method.merge == FALSE)){
##             out$univariate <- object$univariate
##         }else{
##             out$univariate <- model.tables(object, method = method, columns = keep.cols2)
##         }

##         out$univariate$se <- apply(sample.univariate, MARGIN = 1, FUN = stats::sd, na.rm = TRUE)
##         out$univariate$statistic <- out$univariate$estimate/out$univariate$se
##         out$univariate$df <- NA
##         out$univariate$lower <- apply(sample.univariate, MARGIN = 1, FUN = stats::quantile, probs = alpha/2, na.rm = TRUE)
##         out$univariate$upper <- apply(sample.univariate, MARGIN = 1, FUN = stats::quantile, probs = 1-alpha/2, na.rm = TRUE)

##         if(!inherits(calc.pvalue,"try-error")){
##             out$univariate$p.value <- sapply(1:NROW(sample.univariate), function(iE){
##                 boot2pvalue(stats::na.omit(sample.univariate[iE,]),
##                             null = 0,
##                             estimate = out$univariate[iE,"estimate"],
##                             alternative = "two.sided",
##                             add.1 = correction)
##             })
##         }
##         attr(out$univariate,"sample") <- sample.univariate
##     }
    
##     ## ** export
##     attr(out,"seed") <- seqSeed
##     return(out)
## }

##----------------------------------------------------------------------
### resample.R ends here