R/utils.R

In funStatTest: Statistical Testing for Functional Data

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

#' Euclidean norm
#' 
#' @keywords internal
#' @param x Vector of Numeric values
#' 
#'
#'
#' @return Euclidean norm of vector x
#' 
#' @inherit funStatTest_authors author
#' @noRd
norm2 <- function(x) return(sqrt(sum(x^2)))

#' Function to prepare simulated data for plotting
#' 
#' @keywords internal
#' @param mat numeric matrix of dimension `n_point x n_obs` containing 
#' `n_obs` trajectories (in columns) of size `n_point` (in rows).
#' @param sample_id integer value, id of the sample, among `1` or `2`.
#' 
#'
#'
#' @return a `tibble` data frame containing 
#' 
#' @importFrom dplyr mutate row_number bind_rows n
#' @importFrom ggplot2 ggplot aes geom_point geom_line 
#' theme_bw scale_colour_brewer
#' @importFrom stats setNames
#' @importFrom stringr str_c str_extract
#' @importFrom tibble as_tibble
#' @importFrom tidyr pivot_longer
#' @importFrom tidyselect starts_with
#' 
#' @inherit funStatTest_authors author
#' @noRd
prepare_data <- function(mat, sample_id) {
    
    assert_matrix(mat, any.missing = FALSE, min.rows = 1, min.cols = 1)
    qassert(sample_id, "X1[1,2]")
    mat_sample <- c("mat_sample1", "mat_sample2")[sample_id]
    
    return(
        as_tibble(mat, .name_repair = "minimal") %>%
        setNames(str_c("traj", 1:length(names(.)))) %>%
        mutate(sample = sample_id, point = 1L:n()) %>%
        pivot_longer(
            cols = starts_with("traj"), 
            names_to = "trajectory", values_to = "value") %>%
        mutate(
            sample = factor(sample),
            trajectory = as.integer(str_extract(trajectory, "[0-9]+"))
        )
    )
}

#' Compute permutation-based stat values
#' 
#' @keywords internal
#' @inheritParams permut_pval
#' @param verbose boolean, if TRUE, enable verbosity.
#' @param return_perm boolean, if TRUE, the matrix of random permutations used 
#' for the computation is returned.
#' @param n_core integer, number of cores to be used for parallel 
#' computations. Default is 1. **Note:** parallel computing is not enabled 
#' for the moment.
#' 
#'
#'
#' @return a list containing
#' - the `1 x n_perm` matrix of statistic values on permuted samples for each 
#' statistic listed in `stat` input parameter (a `2 x n_perm` matrix for 
#' the `"hkr"` statistic, see [funStatTest::comp_stat()]).
#' 
#' @inherit funStatTest_authors author
#' 
#' @importFrom pbapply pblapply
#' 
#' @seealso [funStatTest::permut_pval()], [funStatTest::comp_stat()]
#' @noRd
permut_stats <- function(
        MatX, MatY, n_perm, stat, verbose = FALSE, return_perm = FALSE, 
        n_core = 1L) {
    
    # global data matrix
    MatW <- cbind(MatX, MatY)
    M <- ncol(MatX)
    N <- ncol(MatY)
    
    # generate permutations
    perm <- replicate(n_perm, sample(N+M))
    
    # temp function to compute stats
    tmp_fun <- function(rep) {
        # permutation
        tmp_perm <- perm[,rep]
        MatXsim <- MatW[,tmp_perm[1:M]]
        MatYsim <- MatW[,tmp_perm[(M+1):(M+N)]]
        
        # compute stat values on permuted samples
        tmp_stat <- comp_stat(MatXsim, MatYsim, stat)
        return(tmp_stat)
    }
    
    # compute permutation-based stat values
    tmp_perm_stat_val <- NULL
    if(verbose) {
        tmp_perm_stat_val <- pblapply(1:n_perm, tmp_fun)
    } else {
        tmp_perm_stat_val <- lapply(1:n_perm, tmp_fun)
    }
    
    
    # reformat permutation-based stat values by stat
    perm_stat_val <- lapply(
        stat,
        function(stat_name) {
            tmp_stat <- Reduce("cbind", lapply(
                1:n_perm,
                function(rep) {
                    return(
                        as.matrix(tmp_perm_stat_val[[rep]][[stat_name]]))
                }
            ))
            return(tmp_stat)
        }
    )
    names(perm_stat_val) <- stat
    
    # return permutation?
    if(return_perm) {
        perm_stat_val[["perm_matrix"]] <- perm
    }
    
    # output
    return(perm_stat_val)
}

#' Compute p-value for a given statistic computed on the original data and on 
#' the permuted data.
#' 
#' @keywords internal
#' 
#' @param original_value a `1 x 1` matrix of statistic values on original data 
#' for each statistic (a `2 x 1` matrix for the `"hkr"` statistic, 
#' see [funStatTest::comp_stat()]).
#' @param perm_value a `1 x n_perm` matrix of statistic values on permuted 
#' samples for each statistic (a `2 x n_perm` matrix for the `"hkr"` 
#' statistic, see [funStatTest::comp_stat()]).
#' 
#'
#'
#' @return a list containing
#' - the `1 x n_perm` matrix of statistic values on permuted samples for each 
#' statistic listed in `stat` input parameter (a `2 x n_perm` for the `"hkr"`
#' statistic, see [funStatTest::comp_stat()]).
#' 
#' @inherit funStatTest_authors author
#' 
#' @seealso [funStatTest::permut_pval()], [funStatTest::comp_stat()]
#' @noRd
pval_stat <- function(original_value, perm_values) {
    
    original_value <- as.matrix(original_value)
    
    assert_matrix(
        original_value, mode = "numeric", 
        min.rows = 1, max.rows = 2, ncols = 1, any.missing = FALSE
    )
    assert_matrix(
        perm_values, mode = "numeric", 
        min.rows = 1, max.rows = 2, min.cols = 1, any.missing = FALSE
    )
    
    assert_set_equal(nrow(original_value), nrow(perm_values))
    
    # number of permuted sample-based stats higher 
    # than original stat values
    # print(str(perm_values))
    o_value_mat <- matrix(
        rep(original_value, ncol(perm_values)), ncol = ncol(perm_values))
    # print(str(o_value_mat))
    tmp_count <- apply(perm_values - o_value_mat > 0, 1, sum)
    tmp_pval <- (tmp_count + 2)/(ncol(perm_values) + 1)
    
    return(tmp_pval)
}