R/power.R

In funStatTest: Statistical Testing for Functional Data

# WARNING - Generated by {fusen} from dev/flat_package.Rmd: do not edit by hand

#' Simulation-based experiment for power analysis
#' 
#'
#' 
#' @description Computation of the statistical power (i.e. risk to reject 
#' the null hypothesis when it is false) associated to any statistics
#' described in the package based on simulation permutation-based p-values
#' computations. See Smida et al (2022) for more details.
#' 
#' @details
#' The `c_val` input argument should be strictly positive so that the null
#' hypothesis is not verified when simulating the data (i.e. so that the two
#' sample are not generated from the same probability distribution).
#' 
#' @param n_simu integer value, number of simulations to compute the 
#' statistical power.
#' @param alpha numerical value, between 0 and 1, type I risk level to reject 
#' the null hypothesis in the simulation. Default value is `5%`.
#' @inheritParams permut_pval
#' @inheritParams simul_data
#' 
#' @return a list with the following elements:
#' - `power_res`: a list of named numeric value corresponding to the 
#' power values for each statistic listed in `stat` input.
#' - `pval_res`: a list of named numeric values corresponding to the p-values 
#' for each simulation and each statistic listed in the `stat` input.
#' - `simu_config`: information about input parameters used for simulation, 
#' including `n_simu`, ` c_val`, ` distrib`, ` delta_shape`, ` n_point`, 
#' ` n_obs1`, ` n_obs2`.
#' 
#' @inherit funStatTest_authors author
#' 
#' @references
#' Zaineb Smida, Lionel Cucala, Ali Gannoun & Ghislain Durif (2022) 
#' A median test for functional data, Journal of Nonparametric Statistics, 
#' 34:2, 520-553,  
#' \doi{10.1080/10485252.2022.2064997}, 
#' [hal-03658578](https://hal.science/hal-03658578)
#' 
#' @seealso [funStatTest::stat_mo()], [funStatTest::stat_med()], 
#' [funStatTest::stat_wmw()], [funStatTest::stat_hkr()], 
#' [funStatTest::stat_cff()], [funStatTest::comp_stat()]
#' 
#' @export
#' 
#' @examples
#' # simulate a few small data for the example
#' res <- power_exp(
#'     n_simu = 20, alpha = 0.05, n_perm = 100, 
#'     stat = c("mo", "med", "wmw", "hkr", "cff"), 
#'     n_point = 25, n_obs1 = 4, n_obs2 = 5, c_val = 10, delta_shape = "constant", 
#'     distrib = "normal", max_iter = 10000, verbose = FALSE
#' )
#' res$power_res
power_exp <- function(
        n_simu = 100, alpha = 0.05, n_perm = 100, stat = c("mo", "med"), 
        n_point = 100, n_obs1 = 50, n_obs2 = 50, 
        c_val = 1, delta_shape = "constant", distrib = "normal", 
        max_iter = 10000, verbose = FALSE) {
    
    assert_count(n_simu, positive = TRUE)
    qassert(alpha, "N1[0,1]")
    assert_count(n_perm, positive = TRUE)
    assert_subset(stat, c("mo", "med", "wmw", "hkr", "cff"), empty.ok = FALSE)
    assert_logical(verbose, len = 1)
    assert_count(n_point, positive = TRUE)
    assert_count(n_obs1, positive = TRUE)
    assert_count(n_obs2, positive = TRUE)
    qassert(c_val, c("N1(0,)", "N1(,0)"))
    assert_choice(
        delta_shape, c("constant", "linear", "quadratic"), null.ok = FALSE)
    assert_count(max_iter, positive = TRUE)
    
    # simulate data and compute p-values
    tmp_pval_res <- lapply(
        1:n_simu, 
        function(sim_id) {
            simu_data <- simul_data(
                n_point, n_obs1, n_obs2, c_val, delta_shape, distrib
            )
            
            MatX <- simu_data$mat_sample1
            MatY <- simu_data$mat_sample2
            
            tmp_res <- permut_pval(MatX, MatY, n_perm, stat, verbose)
            return(tmp_res)
        }
    )
    
    # reformat p-values results by stat
    pval_res <- lapply(
        stat,
        function(stat_name) {
            tmp_stat <- Reduce("cbind", lapply(
                1:n_simu,
                function(sim_id) {
                    return(
                        as.matrix(tmp_pval_res[[sim_id]][[stat_name]]))
                }
            ))
            return(tmp_stat)
        }
    )
    names(pval_res) <- stat
    
    # compute power
    power_res <- lapply(
        stat,
        function(stat_name) {
            tmp_count <- apply(pval_res[[stat_name]] < alpha, 1, sum)
            tmp_power <- tmp_count / n_simu
            return(tmp_power)
        }
    )
    names(power_res) <- stat
    
    return(
        lst(
            power_res, 
            pval_res,
            simu_config = lst(
                n_simu, c_val, distrib, delta_shape, n_point, n_obs1, n_obs2
            )
        )
    )
}