R/calc_power.R

In li-xinran/RIQITE: Randomization Inference for Quantiles of Individual Treatment Effects

# This code supports doing simulation-based power calculations to
# understand which $s$ values are optimal given specific data contexts.



#library( tidyverse )


# Given a Y0 vector and vector of treatment effects tau, arrange them
# so they have a rough correlation of rho with Y0.  In particular, if
# both datasets are transformed to a normal based on their relative
# quantiles, the correlation will be about rho.
#
# If rho = 0, do nothing.
correlate_Y0_and_tx <- function( Y0, tau, rho = 0 ) {

    if ( rho == 0 ) {
        return( tau )
    }

    ords = MASS::mvrnorm(n = length(Y0),
                         mu = c(0,0),
                         Sigma = matrix( c(1,rho,rho,1), nrow=2 ) )
    ords = ords[ order( ords[,1] ), ]
    ords = ords[ rank(Y0), ]

    tau = sort(tau)
    tau = tau[ rank( ords[,2] ) ]

    tau
}


#' Generate a function that makes treatment effects.
#'
#' The function returned is a function that returns a vector of
#' treatment impacts to go with passed vector of Y0 values.  I.e.,
#' final function is of format: `tx_function( Y0 )`.
#'
#' For the linear option, the impacts are a linear function of Y0 of
#'
#' tau = ATE + rho Y0 + epsilon,
#'
#' with epsilon normally distributed with standard deviation of
#' `tx_scale`.
#'
#' For rexp and rnorm, impacts are generated with the passed
#' distribution, and then ordered so if they were transformed to a
#' normal distribution they would have about a rho correlation with
#' the passed outcomes.
#'
#' For scale, the impact is simply a scaling of the Y0 by "rho". This
#' is equivalent to setting ATE = 0 and tx_scale = 0 with "linear", above).
#'
#' For constant, all parameters other than the ATE are ignored.
#'
#' @param distribution A string name of a distribution.  "constant" is
#'   a constant treatment effect. "scale" is a scaling of Y0.  "rexp"
#'   and "rnorm" are those random distributions, generated independent
#'   of the Y0.  "linear" is a linear function of Y0; see above.
#' @param ATE the average treatment effect.  Ignored by "scale"
#'   option.
#' @param tx_scale Scaling factor for impacts.  For rexp and rnorm,
#'   final treatment effects will have this standard deviation, in
#'   expectation.  to achieve desired standard deviation of impacts
#'   generated by tx_function by (default is 1). Ignored by "constant"
#'   and "scale" options.
#' @param rho Strength of relationship between the control-side
#'   distribution and the individual treatment impacts (assuming
#'   variation in both).  NOTE: Depending on the details of how
#'   individual impacts are generated, rho may not be a strict
#'   correlation.  See details of the impact generation functions for
#'   more.
#' @param ... Extra parameters to pass to distribution function (e.g.,
#'   sd for rnorm).
#' @return A function of the form tx_function( Y0 ).
#'
#' @export
tx_function_factory <- function( distribution, ATE = 0, rho = 0, tx_scale = 1, ... ) {
    dots = list( ... )
    if ( distribution == "constant" ) {
        return( function( Y0 ) {
            rep( ATE, length(Y0) )
        } )
    } else if ( distribution == "scale" ) {
        stopifnot( ATE == 0 )
        return( function( Y0 ) {
            Y0 * rho
        } )
    } else if ( distribution == "rexp") {
        return( function( Y0 ) {
            tau = tx_scale * (rexp( length( Y0 ) ) - 1 ) + ATE
            correlate_Y0_and_tx( Y0, tau, rho )
        } )
    } else if ( distribution == "rnorm") {
        return( function( Y0 ) {
            tau = ATE + tx_scale * (rnorm( length( Y0 ) ) )
            correlate_Y0_and_tx( Y0, tau, rho )
        } )
    } else if ( distribution == "linear" ) {
        return( function( Y0 ) {
            ATE + rho * Y0 + tx_scale * rnorm( length( Y0 ) )
        })
    } else {
        stop( "Unrecognized distribution")
    }
}


#' Generate data for power simulation
#'
#'
#' Given a desired sample size and a data generating mechanism,
#' generate a dataset of Y0s and individual treatment effects.
#'
#'
#' @param n Sample size to test
#' @param Y0_distribution Either a function to generate a set of Y0
#'   values, or a list of discrete values to bootstrap sample from
#'   (with replacement).  Function takes single parameter of sample
#'   size. Default is normal.
#' @param tx_function Either a string name of an implemented
#'   distribution (rexp, rnorm, linear, etc.) or a function that takes
#'   Y0 and returns a set of treatment effects.  Default is a constant tx effect.
#'   tx_constant.
#' @param ... Extra parameters to generate treatment impacts with.
#' @return Tibble (dataframe) with two columns, one of Y0s and one
#'   (tau) of treatment effects.
#'
#' @import tibble
#' @export
generate_finite_data <- function( n,
                                  Y0_distribution = rnorm,
                                  tx_function = "constant",
                                  ... ) {
    Y0s = NA

    stopifnot( !is.null( Y0_distribution ) )

    if ( is.function(Y0_distribution) ) {
        Y0s = Y0_distribution( n )
    } else {
        stopifnot( is.numeric(Y0_distribution) )
        Y0s = sample( Y0_distribution, n, replace=TRUE )
    }

    if ( is.character( tx_function ) ) {
        tx_function = tx_function_factory( tx_function, ... )
    }
    tau = tx_function( Y0s )

    tibble( Y0 = Y0s, tau = tau )
}






#' Summarize set of power simulation results
#'
#' Summarize the CI_precision simulation to compare the performance of
#' various Stephenson rank statistics.
#'
#' Calculates three types of precision.  First, we summarize the
#' distribution of the confidence limit for the number of units with
#' effects passing the threshold cutoff 'c' (default of c=0). We also,
#' for each rank, calculate the median confidence limits for the
#' individual treatment effect for that rank. Finally, we calculate
#' the power for testing whether the individual treatment effect
#' passes the threshold cutoff for each rank.
#'
#' @param result A data frame with one row for each quantile in the
#'   result object passed (e.g., the result from `calc_power_finite()`
#'   with `targeted_power=FALSE`).
#'
#' @param alternative A character takes value "greater" and "less",
#'   indicating the direction of alternative hypotheses.
#' @param c A number specifying the cutoff.
#' @param quantile_n_CI What quantile of the distribution of lower
#'   confidence limits for the number of significant units to
#'   calculate.  NA (the default) means calculate mean instead of
#'   median. 0.5 would report the median size of the confidence
#'   interval for number of significant units.
#'
#' @export
summary_power_results <- function(result,
                                  alternative = "greater",
                                  c = 0,
                                  quantile_n_CI = NA ) {



    # if ( !is.matrix( result ) && !is.null( result$s ) ) {
    #     # Recursively call this method on each group defined by s
    #     alt_full <- result %>% group_by( s ) %>%
    #         nest()
    #
    #     alt <- alt_full$data %>% set_names( alt_full$s ) %>%
    #         map_df( summary_power_results, .id = "s" ) %>%
    #         mutate( s = as.numeric(s) )
    #     return( alt )
    # }

    if ( alternative == "greater") {
        n <- rowSums( result > c )
        pows <- apply( result > c, 2, mean )
        cis <- apply( result, 2, median )
    } else {
        n <- rowSums( result < c )
        pows <- apply( result < c, 2, mean )
        cis <- apply( result, 2, median )
    }

    n_val = NA
    if ( is.na( quantile_n_CI ) ) {
        n_val = mean( n )
    } else {
        n_val = quantile( n, quantile_n_CI )
    }

    res <- tibble( k = as.numeric( colnames(result) ),
                   power = pows,
                   q_ci = cis,
                   n = n_val )

    return( res )

}






#' Calculate power or properties of quantile confidence intervals for
#' a finite-sample dataset.
#'
#' Given a fixed schedule of potential outcomes and a proportion to
#' assign to treatment, calculate (via simulation) either the power to
#' reject a null hypothesis or the distribution of all the confidence
#' intervals.
#'
#' If targeted_power = FALSE, this method will estimate (via
#' simulation) the distributions all simultaneous confidence intervals
#' that can then be post processed to determine their variability, and
#' also the CI for the number of positive units, etc.
#'
#' Otherwise, it will, by simulation, calculate the chance of
#' rejecting a null by repeatedly rerandomizing the data and
#' conducting a hypothesis test on the results, via using
#' `pval_quantile()` (or `pval_bound()` if `use_pval_bound` is TRUE)
#'
#' @inheritParams pval_quantile
#'
#' @param alpha Alpha for testing (default of 0.05)
#' @param R Number of permutation replicates for power calculation.
#' @param percentile Which quantile is of interest (default is 1, the
#'   max).
#' @param c Threshold for null (see pval_quantile)
#' @param If TRUE, calculates the overall bound on impacts using
#'   `pval_bound()` (which allows for the Diff in Means test statistic
#'   and other stats that do not need to be distribution free).
#'   Otherwise uses `pval_quantile()`.
#' @param quantile_n_CI What quantile of the distribution of number of
#'   significant units should be reported?
#' @return If targeted_power=TRUE, return a small tibble with
#'   statistics of the power calculation (true effect being tested,
#'   power to detect, sample size, etc.).  sd_Y0, sd_Y1, sd_tau are
#'   the standard deviations of the potential outcomes and individual
#'   treatment impacts, and r is the correlation of Y0 and the
#'   treatment effects.  Otherwise, if summarise = FALSE, return a
#'   matrix with each column being a simulation run and each row being
#'   a quantile, with the listed confidence interval limit for each
#'   entry.  Otherwise a tibble with a row for each quantile specified
#'   in k.vec.  The "n" column is a repeated value of the median
#'   number of quantiles found significant across simulations.#'
#' @import tibble
#' @export
calc_power_finite <- function( Y0, tau, p_tx, R = 100,
                               percentile = 1,
                               alpha = 0.05,
                               c = 0,
                               alternative = "greater",
                               method.list = list(name = "Stephenson", s = 10),
                               score = NULL,
                               stat.null = NULL,
                               use_pval_bound = FALSE,
                               nperm = 1000,
                               k.vec = NULL,
                               quantile_n_CI = NA,
                               targeted_power = TRUE,
                               summarise = TRUE ) {


    n = length(Y0)
    stopifnot( length(tau) == n )

    k = NA
    if ( targeted_power ) {
        stopifnot( percentile >= 0 && percentile <= 1 )
        k = round( percentile * n )
    }

    if ( is.null( stat.null ) && !(targeted_power && use_pval_bound) ) {
        stat.null = null_dist(n, floor( n * p_tx ), method.list = method.list, nperm = nperm )
    }

    rps = map( 1:R, ~ {
        Z = as.numeric( sample(n) <= p_tx * n )
        Yobs = Y0 + ifelse( Z, tau, 0 )

        if ( !targeted_power ) {
            ci_quantile( Z, Yobs, k.vec = k.vec,
                         alternative = alternative,
                         method.list = method.list,
                         score = score,
                         stat.null = stat.null,
                         nperm = nperm,
                         alpha = alpha,
                         switch = FALSE )
        } else if ( use_pval_bound ) {
            pval = pval_bound( Z, Yobs, c = c,
                               alternative = alternative,
                               method.list = method.list,
                               nperm = nperm )
            data.frame( pval = pval )

        } else {
            pval = pval_quantile( Z, Yobs, k = k, c = c,
                                  alternative = alternative,
                                  method.list = method.list,
                                  score = score,
                                  stat.null = stat.null,
                                  nperm = nperm,
                                  switch = FALSE )
            data.frame( pval = pval )
        }
    } )

    res = bind_rows( rps )

    if ( !targeted_power ) {
        # We are looking at confidence intervals, not pvalues
        if ( alternative == "greater" ) {
            res <- matrix( res$lower, ncol = R )
        } else {
            res <- matrix( res$upper, ncol = R )
        }
        res = t(res)

        if ( is.null(k.vec) ) {
            colnames(res) = 1:ncol(res)
        } else {
            colnames(res) = k.vec
        }

        if ( summarise ) {
            return( summary_power_results( res,
                                           alternative = alternative,
                                           quantile_n_CI = quantile_n_CI,
                                           c = c ) )
        } else {
            return( res )
        }
    } else {
        true_tau = sort(tau)[k]
        tibble( n = n, k = k, tau_k = true_tau,
                sd_Y0 = sd( Y0 ),
                sd_Y1 = sd( Y0 + tau ),
                sd_tau = sd( tau ),
                r = ifelse( sd_Y0 > 0 & sd_tau > 0, cor( tau, Y0 ), NA ),
                power = mean(res$pval <= alpha) )
    }
}






# Fit a MLM to collection of finite power calculations see how much
# different datasets impacts the power.
calc_power_ICC <- function( rps, iter_per_set ) {
    require( lme4 )
    rps$gid = 1:nrow(rps)
    rps$k = round( iter_per_set * rps$power )

    if ( sd( rps$k ) == 0 ) {
        tibble( ICC = NA,
                pow_l = NA,
                pow_h = NA )
    } else {

        M1 = glmer( cbind( k, iter_per_set-k ) ~ 1 + (1|gid), data=rps,
                    family=binomial )
        #display( M1 )
        V = VarCorr( M1 )$gid[1,1]
        V = round( V, digits = 4 )
        ICC =  V / ( V + pi^(2/3) )
        rng = boot::inv.logit( fixef(M1)[[1]] + c( -sqrt(V), sqrt(V) ) )
        ICC = round( ICC, digits = 4 )
        tibble( ICC = ICC,
                pow_l = rng[[1]],
                pow_h = rng[[2]] )
    }
}



#' Calculate power for the noninferiority null
#'
#' Given a specified DGP for control outcomes and treatment impacts,
#' calculate power to detect effects.  Calculates superpopulation
#' power via simulation, where we simulate series of datasets of size
#' n, randomized into treatment and control, and then test for
#' detection of effects. Function returns proportion of times effects
#' are detected.
#'
#' Each individual dataset is tested iter_per_set times, giving very
#' rough power calculations for individual sets that are then
#' aggregated; this provides the ability to do variance decomposition
#' to see how sensitive power is to individual dataset characteristics
#' (for datasets of the same family as defined by the DGP) vs.
#' randomization imbalance.
#'
#' @inheritParams generate_finite_data
#' @inheritParams calc_power_finite
#'
#' @param iter_per_set Number of permutations per finite dataset.
#' @param calc_ICC Calculate how much power varies across finite
#'   datasets.
#' @import tidyverse
#' @export
calc_power <- function( n,
                        Y0_distribution = rnorm,
                        tx_function = "constant",
                        p_tx = 0.5,
                        R = 100, iter_per_set = 10,
                        percentile = 1,
                        alpha = 0.05,
                        c = 0,
                        k.vec = NULL,
                        alternative = "greater",
                        method.list = list(name = "Stephenson", s = 10),
                        targeted_power = TRUE,
                        use_pval_bound = FALSE,
                        score = NULL,
                        stat.null = NULL,
                        nperm = 1000,
                        calc_ICC = FALSE,
                        quantile_n_CI = NA,
                        summarise = TRUE,
                        ... ) {

    stopifnot( !is.null( Y0_distribution ) )

    one_run <- function( Rint ) {
        dat = generate_finite_data( n = n,
                                    Y0_distribution,
                                    tx_function,
                                    ... )

        calc_power_finite( Y0 = dat$Y0, tau = dat$tau, p_tx = p_tx,
                           R = Rint,
                           percentile = percentile, alpha = alpha,
                           c = c,
                           alternative = alternative, method.list = method.list,
                           score = score, stat.null = stat.null, nperm = nperm,
                           targeted_power = targeted_power,
                           use_pval_bound = use_pval_bound,
                           k.vec = k.vec,
                           quantile_n_CI = quantile_n_CI,
                           summarise = summarise )
    }

    n_block = round( R / iter_per_set )
    rps = map( 1:n_block, ~ one_run( iter_per_set ) )

    if ( !summarise ) {
        # Return raw results and stop.
        if ( !targeted_power ) {
            rps = do.call( rbind, rps )
            return( rps )
        } else {
            return( bind_rows( rps, .id = "runID" ) )
        }
    }

    rps <- bind_rows( rps )

    if ( targeted_power ) {
        aggrps <- rps %>%
            summarize( n = mean(n),
                       k = mean(k),
                       sd_tau = sd( .data$tau_k ),
                       tau_k = mean(tau_k),
                       SEpower = sd( .data$power ) / sqrt( n() ),
                       power = mean(power) ) %>%
            relocate( tau_k, sd_tau, power, SEpower, .after = k )

        if ( calc_ICC ) {
            aggrps = cbind( aggrps,
                            calc_power_ICC(rps, iter_per_set) )
        }
    } else {
        aggrps <- rps %>% group_by( k ) %>%
            summarise( SEpower = sd( power ) / sqrt( n() ),
                       power = mean( power ),
                       q_ci = mean( q_ci ),
                       n = mean( n ) ) %>%
            relocate( SEpower, .after = power )
        aggrps$k = as.numeric( aggrps$k )
    }
    aggrps$method = paste0( method.list, collapse = "-" )
    aggrps
}







if ( FALSE ) {
    pow = calc_power( 500,
                      tx_function = "rexp",
                      p_tx = 0.5 )
    pow
}




#### Exploring stephenson s statistics ####


#' Explore power and precision across a range of s values
#'
#' Conduct power simulation against specific finite-sample dataset to
#' see how power changes as a function of s.
#'
#' @param s Numeric list of s values to test using the Stephenson
#'   rank.
#' @param Y0 Either a name of function to generate Y0s, or an
#'   empirical distribution represented as a list of values.
#' @param p_tx Proportion of units to be treated.
#' @param R Number of simulation iterations to do.
#' @param percentile For targeted power, what percentile to target.
#' @param k.vec List of ranks to calculate power for (if not targeted
#'   power).
#' @param alpha Significance level of the test.
#' @param c Threshold bound for treatment impacts for hypothesis
#'   testing.
#' @param stat.nulls Optional list of the precomputed null
#'   distributions of the statistics.
#' @param k.vec Which quantiles to examine (if targeted_power =
#'   FALSE).
#' @param targeted_power TRUE means focus on passed percentile and
#'   calculate power for that specific quantile.  FALSE means
#'   calculate simultenaous confidence intervals for all quantiles,
#'   leaving power caculation to post-simulation exploration.
#'
#' @inheritParams summary_power_results
#'
#' @import dplyr purrr
#' @export
explore_stephenson_s_finite <- function( s = c(2,5,10,30),
                                         Y0, tau, p_tx,
                                         R = 100,
                                         percentile = 1,
                                         k.vec = NULL,
                                         alpha = 0.05,
                                         c = 0,
                                         quantile_n_CI = NA,
                                         alternative = "greater",
                                         stat.nulls = NULL,
                                         nperm = 1000,
                                         targeted_power = TRUE ) {
    n = length(Y0)
    stopifnot( n == length( tau ) )
    stopifnot( p_tx*n >= 1 && p_tx*n <= n-1 )
    if ( !is.null( stat.nulls ) ) {
        stopifnot( length(stat.nulls) == length(s) )
        names(stat.nulls) = s
    }

    rs = s %>% set_names(s) %>%
        map_df( function( ss ) {
            sn = NULL
            if ( !is.null(stat.nulls) ) {
                sn = stat.nulls[[ as.character(ss) ]]
            }
            calc_power_finite( Y0 = Y0, tau = tau, p_tx = p_tx, R = R,
                               percentile = percentile,
                               alpha = alpha, c = c,
                               method.list = list(name = "Stephenson", s = ss),
                               alternative = alternative,
                               stat.null = sn, nperm = nperm,
                               targeted_power = targeted_power,
                               k.vec = k.vec,
                               quantile_n_CI = quantile_n_CI,
                               summarise = TRUE )
        }, .id="s" )

    rs$s = as.numeric(rs$s)

    return( rs )
}



if ( FALSE ) {
    Y0 = rnorm( 500 )
    tau = rexp( 500 ) - 1
    explore_stephenson_s_finite( Y0 = Y0, tau = tau, p_tx = 0.5, R = 10, nperm = 100, targeted_power = FALSE )
    explore_stephenson_s_finite( Y0 = Y0, tau = tau, p_tx = 0.5, R = 10, nperm = 100, targeted_power = TRUE )

    explore_stephenson_s_finite( Y0 = Y0, tau = tau, p_tx = 0.5, nperm = 10, targeted_power = FALSE,
                                 summarise = FALSE )

}




#' Explore power across a range of s values
#'
#' Conduct power simulation against a specific data generating process to
#' see how power changes as a function of s in that scenario.
#'
#' @inheritParams explore_stephenson_s_finite
#'
#' @param s Numeric list of s values to test using the Stephenson rank.
#' @import tidyverse
#' @import future
#' @import parallel
#' @import furrr
#' @export
explore_stephenson_s <- function( s = c(2,5,10,30),
                                  n,
                                  Y0_distribution = rnorm,
                                  tx_function = "constant",
                                  p_tx = 0.5,
                                  R = 100, iter_per_set = 10,
                                  percentile = 1,
                                  alpha = 0.05,
                                  c = 0,
                                  k.vec = NULL,
                                  quantile_n_CI = NA,
                                  targeted_power = TRUE,
                                  alternative = "greater",
                                  nperm = 1000,
                                  parallel = FALSE,
                                  n_workers = NULL,
                                  calc_ICC = FALSE,
                                  verbose = FALSE,
                                  ... ) {


    # generate null distributions
    m = floor( n * p_tx )
    Z.perm = assign_CRE(n, m, nperm)
    stat.null.all = list()
    method.list.all = list()
    for( s.ind in seq_along(s) ) {
        method.list.all[[s.ind]] = list( name = "Stephenson", s = s[s.ind] )
        stat.null.all[[s.ind]] = null_dist(n, m, method.list = method.list.all[[s.ind]], Z.perm = Z.perm )
    }

    if ( verbose ) {
        cat( glue::glue( "Null distributions calculated. {nperm} permutations across {length(s)} s-values." ), "\n" )
    }

    one_run <- function( Rint ) {
        dat = generate_finite_data( n = n,
                                    Y0_distribution = Y0_distribution,
                                    tx_function = tx_function, ... )
        explore_stephenson_s_finite( s = s,
                                     Y0 = dat$Y0, tau = dat$tau,
                                     p_tx = p_tx,
                                     R = Rint,
                                     percentile = percentile, alpha = alpha,
                                     c = c, k.vec = k.vec,
                                     quantile_n_CI = quantile_n_CI,
                                     alternative = alternative,
                                     stat.nulls = stat.null.all,
                                     nperm = nperm,
                                     targeted_power = targeted_power )
    }

    n_blocks = round( R / iter_per_set )

    if ( !parallel ) {
        if ( verbose ) {
            cat( glue::glue( "Doing {iter_per_set * n_blocks} iterations in {n_blocks} blocks of {iter_per_set} iterations each." ), "\n" )
        }
        rps = map( rep( iter_per_set, n_blocks ), one_run )
    } else {
        require( future )
        require( parallel )
        require( furrr )
        if ( is.null( n_workers ) ) {
            n_workers = parallel::detectCores() - 1
        }
        future::plan(multisession, workers = n_workers )
        if ( verbose ) {
            cat( glue::glue( "Spreading {n_blocks} jobs of {iter_per_set} iterations each over {n_workers} cores." ), "\n" )
        }
        rps = future_map( rep( iter_per_set, n_blocks ),
                          one_run,
                          .options = furrr_options(seed = NULL), .progress = verbose )
        future::plan("default")
    }
    rps <- bind_rows( rps, .id = "gid" )

    if ( targeted_power ) {
        aggrps <- rps %>% group_by( s ) %>%
            summarise( EvarY0 = mean( sd_Y0^2 ),
                       EvarY1 = mean( sd_Y1^2 ),
                       Evartau = mean( sd_tau^2 ),
                       Er = mean( r ),
                       SEpower = sd( power ) / sqrt( n() ),
                       power = mean( power ) ) %>%
            relocate( SEpower, .after = power )

        if ( calc_ICC ) {
            rpg = rps %>% group_by( s ) %>% nest()
            mods = rpg$data %>%
                set_names(rpg$s) %>%
                map_df( calc_power_ICC,
                        iter_per_set = iter_per_set,
                        .id = "s" ) %>%
                mutate( s = as.numeric( s ) )

            aggrps <- left_join( aggrps, mods, by="s" )
        }
    } else {
        aggrps <- rps %>% group_by( s, k ) %>%
            summarise( SEpower = sd( power ) / sqrt( n() ),
                       power = mean( power ),
                       q_ci = mean( q_ci ),
                       n = mean( n ) ) %>%
            relocate( SEpower, .after = power )
        aggrps$s = as.numeric( aggrps$s )
        aggrps$k = as.numeric( aggrps$k )

        if ( calc_ICC ) {
            rpg = rps %>% group_by( s, k ) %>% nest()
            mods = rpg$data %>%
                map_df( calc_power_ICC,
                        iter_per_set = iter_per_set )
            mods$s = rpg$s
            mods$k = rpg$k
            aggrps <- left_join( aggrps, mods, by=c( "s", "k" ) )
        }
    }


    aggrps
}







#### Testing code ####

if ( FALSE ) {

    library( tidyverse )
    library( RIQITE )
    dat = generate_finite_data( n = 100,
                                Y0_distribution = rnorm )
    dat

    result = calc_CI_precision_finite( dat$Y0, dat$tau, 0.5, R = 10 )


    Y0 = dat$Y0
    tau = dat$tau


    ess = explore_stephenson_s( s = c( 2, 5, 10 ),
                                n = 500,
                                tx_function = "rexp",
                                p_tx = 0.5,
                                nperm = 100 )
    ess

    ess2 = explore_stephenson_s( s = c( 2, 4, 5, 7, 10, 30 ),
                                 n = 500,
                                 tx_function = "rnorm",
                                 ATE = 0.2,
                                 p_tx = 0.5 )
    ess2

    ess3 = explore_stephenson_s( s = c( 2, 4, 5, 7, 10, 30 ),
                                 n = 500,
                                 tx_function = "constant",
                                 ATE = 0.2,
                                 p_tx = 0.5 )
    ess3

    ess4 = explore_stephenson_s( s = c( 2, 4, 5, 7, 10, 30 ),
                                 n = 500,
                                 Y0_distribution = rexp,
                                 tx_function = "constant",
                                 ATE = 0,
                                 p_tx = 0.5 )
    ess4



}