R/functions_main_gmtfd.R

In gmtFD: General Multiple Tests for Univariate and Multivariate Functional Data

#' 
#' General multiple tests for univariate and multivariate functional data
#'
#' The function \code{gmtFD()} calculates the statistical tests based on globalizing and supremum 
#' pointwise Hotelling's \eqn{T^2}-test statistics (GPH and SPH) for the global null hypothesis and 
#' multiple local null hypotheses. Respective \eqn{p}-values are obtained by a parametric bootstrap strategy.
#'
#' @param x a list of \eqn{k} elements corresponding to groups. Each element representing a group 
#' is a list of \eqn{p} elements corresponding to functional variables, and each such element
#' (representing a functional variable) is a matrix of size \eqn{n_i\times ntp} of descrete observations 
#' in design time points. \eqn{ntp} denotes a number of design time points. 
#' @param h contrast matrix. For contrast matrices based on Dunnett’s and Tukey’s contrasts, 
#' it can be created by the \code{contr_mat()} function from the package \code{GFDmcv} (see examples).
#' @param blocks_contrasts a vector with blocks of contrasts labels. The integer labels (from 1 to a number 
#' of blocks of contrasts) should be used in non-decreasing order. For univariate case (\eqn{p=1}), 
#' this is the vector representing group labels, e.g., for \eqn{k=3} and Tukey's contrasts, we have \code{c(1, 2, 3)}.
#' For multivariate case (\eqn{p>1}), we have more possibilities. Mainly, we compare samples (groups) rather 
#' than particular functional variables, and then we use the vector grouping all variables in samples
#' taking into account contrasts, e.g., for \eqn{k=4}, \eqn{p=2} and Tukey's contrasts,
#' we have \code{c(1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6)}, where pairs \code{(i,i)} for \code{i=1,...,6} 
#' represent the two variables of observations in the \code{i}th contrast. See examples.
#' @param n_boot number of bootstrap samples.
#' @param alpha significance level.
#' @param parallel a logical indicating whether to use parallelization.
#' @param n_cores if \code{parallel = TRUE}, a number of processes used in parallel computation.
#' Its default value means that it will be equal to the number of cores of a computer used.
#' @param multi_gen a logical indicating of whether to use separate multiple generations of Gaussian processes
#' for the parametric bootstrap tests. The default is FALSE, which means that the processes will be
#' generated once in a big matrix. This method is much faster, but for larger \eqn{n=n_1+\dots+n_k} and \eqn{p}
#' the generated data can be too large for RAM. In such a case, we suggest using separate generation
#' (\code{multi_gen = TRUE}), which is slower, but possible to calculate.
#'
#' @details The function \code{gmtFD()} concerns the tests for the heteroscedastic
#' contrast testing problem for univariate and multivariate functional data. The details are 
#' presented in Munko et al. (2023, 2024), but here we present some summary of the problem 
#' and its solutions implemented in the package.
#'
#' Suppose we have \eqn{k} independent functional samples \eqn{\mathbf{x}_{i1},\dots,\mathbf{x}_{in_i}},
#' which consist of independent and identically distributed \eqn{p}-dimensional stochastic processes defined
#' on interval \eqn{\mathcal{T}} with mean function vector \eqn{\boldsymbol{\eta}_i} and covariance function
#' \eqn{\boldsymbol{\Gamma}_i} for each \eqn{i\in\{1,\dots,k\}}. Note that the covariance functions
#' of the different groups may differ from each other, i.e., heteroscedasticity is explicitly allowed.
#'
#' We consider the null and alternative hypothesis
#' \deqn{\mathcal H_0: \mathbf{H}\boldsymbol{\eta}(t) =\mathbf{0}_r \text{ for all }
#' t\in \mathcal{T} \quad \text{vs.} \quad \mathcal H_1: \mathbf{H}\boldsymbol{\eta}(t) \neq
#' \mathbf{0}_r \text{ for some } t\in \mathcal{T},} where \eqn{\mathbf{H} \in \mathbb{R}^{r \times pk}}
#'  denotes a known contrast matrix, i.e., \eqn{\mathbf{H}\mathbf{1}_{pk} = \mathbf{0}_r}, and 
#'  \eqn{\boldsymbol{\eta} := (\boldsymbol{\eta}_1^{\top},\dots,\boldsymbol{\eta}_k^{\top})^{\top}} 
#'  is the vector of the mean functions. The formulation of this testing framework is very general
#'  and contains many special cases like the analysis of variance for univariate and multivariate functional data
#'  (FANOVA and FMANOVA) problems. 
#'  
#'  For univariate functional data (\eqn{p=1}), we may choose \eqn{\mathbf{H} = \mathbf{P}_k} 
#'  for the one-way FANOVA problem to test the null hypothesis of no main effect, where 
#'  \eqn{\mathbf{P}_k:=\mathbf{I}_k-\mathbf{J}_k/k} with \eqn{\mathbf{I}_k \in\mathbb{R}^{k\times k}} 
#'  denoting the unit matrix and \eqn{\mathbf{J}_k := \mathbf{1}_k\mathbf{1}_k^{\top} \in\mathbb{R}^{k\times k}}
#'   denoting the matrix of ones. However, there are different possible choices of the
#'   contrast matrix \eqn{\mathbf{H}} which lead to this global null hypothesis.
#'   Many-to-one comparisons can be considered by choosing Dunnett's contrast matrix
#'  \eqn{\mathbf{H} = [-\mathbf{1}_{k-1}, \mathbf{I}_{k-1}]}, where the mean
#'  functions \eqn{\eta_2,\dots,\eta_k} are compared to the mean function \eqn{\eta_1}
#'   of the first group regarding the different contrasts. To compare all pairs
#'   of mean functions \eqn{\eta_{i_1},\eta_{i_2}, i_1,i_2 \in\{1,\dots,k\}} with
#'   \eqn{i_1 \neq i_2}, the Tukey's contrast matrix:
#'  \deqn{\mathbf{H} = \begin{bmatrix}
#'  -1 & 1 & 0 & 0 & \cdots & \cdots & 0 \\
#'  -1 & 0 & 1 & 0 &\cdots & \cdots & 0 \\
#'  \vdots  & \vdots &\vdots & \vdots & \ddots & \vdots & \vdots  \\
#'  -1 & 0 & 0 & 0& \cdots & \cdots & 1\\
#'  0 & -1 & 1 & 0& \cdots & \cdots & 0 \\
#'  0 & -1 & 0 & 1& \cdots & \cdots & 0 \\
#'  \vdots  & \vdots & \vdots  & \vdots & \ddots & \vdots & \vdots \\
#'  0 & 0 & 0 & 0 & \cdots & -1 & 1
#'  \end{bmatrix}  \in \mathbb{R}^{k(k-1)/2 \times k}} can be used.
#'  
#'  For multivariate functional data (\eqn{p>1}), we extend the matrices for \eqn{p=1} given above
#'  as follows: \eqn{\mathbf{H}=\mathbf{H}_1\otimes \mathbf{I}_p}, where \eqn{\mathbf{H}_1} is the
#'  contrast matrix for the univariate case. For more examples of contrast matrices in different
#'  settings, see Merle et al. (2023, 2024).
#'
#'  For this testing problem, we consider the pointwise Hotelling's \eqn{T^2}-test statistic
#'  \deqn{\mathrm{PH}_{n,\mathbf{H}}(t):= n(\mathbf{H}\boldsymbol{\widehat\eta}(t))^{\top}
#'  (\mathbf{H}\mathbf{\widehat \Gamma}(t,t) \mathbf{H}^{\top})^+
#'  \mathbf{H} \boldsymbol{\widehat\eta}(t)}
#'  for all \eqn{t \in\mathcal{T}}, where 
#'  \eqn{\boldsymbol{\widehat\eta} := (\boldsymbol{\widehat\eta}_1^{\top},\dots,\boldsymbol{\widehat\eta}_k^{\top})^{\top}}
#'  denotes the vector of all mean function estimators, \eqn{\mathbf{A}^+} denotes the
#'  Moore-Penrose inverse of the matrix \eqn{\mathbf{A}}, and
#'  \deqn{\boldsymbol{\widehat{\Gamma}}(t,s):= \mathrm{diag}\left( \frac{n}{n_1}\widehat{\boldsymbol{\Gamma}}_1(t,s),
#'  \ldots,\frac{n}{n_k}\widehat{\boldsymbol{\Gamma}}_k(t,s)\right),} \eqn{n=n_1+\dots+n_k}, 
#'  \eqn{\widehat{\boldsymbol{\Gamma}}_i(t,s)} is the sample covariance function for the \eqn{i}-th group, 
#'  \eqn{i\in\{1,\dots,k\}}. Based on this pointwise Hotelling's \eqn{T^2}-test statistic, we construct the globalizing
#'  and supremum of pointwise Hotelling's \eqn{T^2}-test (GPH and SPH) statistics by integrating and supremum over 
#'  the pointwise Hotelling's \eqn{T^2}-test statistic, that is
#'  \deqn{I_{n}(\mathbf{H}) := \int_{\mathcal{T}} \mathrm{PH}_{n,\mathbf{H}}(t) \,\mathrm{ d }t,\ \ 
#'  T_{n}(\mathbf{H}) := \mathrm{sup}_{t\in\mathcal{T}} \mathrm{PH}_{n,\mathbf{H}}(t).}
#'  We consider the parametric bootstrap test based on these test statistics. However, for better post hoc
#'  analysis, we also consider the multiple contrast testing procedures. The main idea of multiple contrast
#'  tests is to split up the global null hypothesis with matrix
#'  \eqn{\mathbf{H}= [\mathbf{H}_1^{\top}, \dots, \mathbf{H}_R^{\top}]^{\top}} into \eqn{R} matrices 
#'  \eqn{\mathbf{H}_{\ell}\in\mathbb{R}^{r_{\ell}\times k}} with \eqn{\mathrm{rank}(\mathbf{H}_{\ell})\geq 1}, 
#'  \eqn{\ell\in\{1,\dots,R\}}, where \eqn{R,r_{\ell}\in\{1,\dots,r\}}, and \eqn{\sum_{\ell=1}^Rr_{\ell}=r}.
#'  This leads to the multiple testing problem with null hypotheses
#'  \deqn{\mathcal H_{0,{\ell}} : \; \mathbf{H}_{\ell} \boldsymbol{\eta}(t) =
#'  \mathbf{0}_{r_{\ell}} \;\text{ for all }t\in\mathcal{T},	\text{for }\ell\in \{1,\ldots,R\}.}
#'  To verify this family of null hypotheses, we adopt two approaches. First, we simply apply
#'  the above test to each hypothesis \eqn{\mathcal H_{0,{\ell}}}, and the resulting \eqn{p}-values
#' are then corrected by the Bonferroni's method. However, this approach, denoted in the package as
#' GPH and SPH, may give conservative test and loss of power. Thus, we also consider the test adopting the
#' idea for the construction of simultaneous confidence bands proposed by Buhlmann (1998).
#' This test is denoted by mGPH and mSPH in the package and is a more powerful solution than the GPH and SPH
#' procedures, which was shown in Munko et al. (2023, 2024).
#'
#' Note that the value of the test statistics for the mGPH and mSPH tests for global hypotheses are
#' equal to \deqn{\max_{\ell\in\{1,\ldots,R\}}\frac{I_{n}(\mathbf{H}_{\ell})}{q_{GPH,\ell,\widetilde{\beta}}^{\mathcal{P}}}\text{ and }
#' \max_{\ell\in\{1,\ldots,R\}}\frac{T_{n}(\mathbf{H}_{\ell})}{q_{SPH,\ell,\widetilde{\beta}}^{\mathcal{P}}},}
#' respectively, where \eqn{q_{GPH,\ell,\widetilde{\beta}}^{\mathcal{P}}} and \eqn{q_{SPH,\ell,\widetilde{\beta}}^{\mathcal{P}}} 
#' are the quantiles calculated using the adaptation of the method by Buhlmann (1998). 
#' The critical value for it is always 1.
#'
#' Please have a look at a summary function designed for this package. It can be used
#' to simplify the output of \code{gmtFD()} function.
#'
#' @return A list of class \code{multifmanova} containing the following components:
#' \item{res_global}{a data frame containing the results for testing the global null hypothesis,
#' i.e., test statistics and \eqn{p}-values.}
#' \item{res_multi}{all results of multiple contrasts tests for particular hypothesis in
#' a contrast matrix \code{h}, i.e., test statistics, critical values and \eqn{p}-values.}
#' \item{k}{a number of groups.}
#' \item{p}{a number of variables.}
#' \item{ntp}{a number of design time points.}
#' \item{n}{a vector of sample sizes.}
#' \item{h}{an argument \code{h}.}
#' \item{n_boot}{an argument \code{n_boot}.}
#' \item{alpha}{an argument \code{alpha}.}
#' \item{multi_gen}{an argument \code{multi_gen}.}
#'
#' @examples
#' # Some of the examples may run some time.
#'
#' # Canadian weather data set
#' # There are three samples of mean temperature and precipitations for
#' # fifteen weather stations in Western Canada, another fifteen in Eastern Canada, 
#' # and the remaining five in Northern Canada.
#' 
#' # one functional variable - temperature
#' library(fda)
#' data_set_t <- t(CanadianWeather$dailyAv[,, "Temperature.C"])
#' # number of samples
#' k <- 3
#' # number of variables
#' p <- 1
#' # preparing data set
#' gr_label <- rep(c(1, 2, 3), c(15, 15, 5))
#' data_set <- list(list(data_set_t[gr_label == 1, ]),
#'                  list(data_set_t[gr_label == 2, ]),
#'                  list(data_set_t[gr_label == 3, ]))
#' # trajectories of mean temperature and precipitation
#' oldpar <- par(mar = c(4, 4, 2, 0.1))
#' matplot(t(data_set_t), type = "l", col = gr_label, lty = 1,
#'         xlab = "Day", ylab = "Temperature (C)",
#'         main = "Canadian weather data set")
#' legend("bottom", legend = c("Western Canada", "Eastern Canada", "Northern Canada"),
#'        col = 1:3, lty = 1)
#' par(oldpar)
#' 
#' # Tukey's contrast matrix
#' h_tukey <- GFDmcv::contr_mat(k, type = "Tukey")
#' h_tukey_m <- kronecker(h_tukey, diag(p))
#' # vector of blocks of contrasts labels
#' blocks_contrasts <- rep(1:(nrow(h_tukey_m) / p), each = p)
#' \donttest{
#' # testing without parallel computing
#' res <- gmtFD(data_set, h_tukey_m, blocks_contrasts)
#' summary(res, digits = 3)}
#' # plots for pointwise Hotelling's T^2-test statistics
#' oldpar <- par(mfrow = c(2, 2), mar = c(4, 2, 2, 0.1))
#' plot(ph_test_statistic(data_set, h_tukey_m), type = "l",
#'      ylim = c(0, max(ph_test_statistic(data_set, h_tukey_m))),
#'      main = "Global hypothesis", xlab = "Day")
#' plot(ph_test_statistic(data_set, matrix(h_tukey_m[1, ], 1)), type = "l",
#'      ylim = c(0, max(ph_test_statistic(data_set, h_tukey_m))),
#'      main = "Contrast 1", xlab = "Day")
#' plot(ph_test_statistic(data_set, matrix(h_tukey_m[2, ], 1)), type = "l",
#'      ylim = c(0, max(ph_test_statistic(data_set, h_tukey_m))),
#'      main = "Contrast 2", xlab = "Day")
#' plot(ph_test_statistic(data_set, matrix(h_tukey_m[3, ], 1)), type = "l",
#'      ylim = c(0, max(ph_test_statistic(data_set, h_tukey_m))),
#'      main = "Contrast 3", xlab = "Day")
#' par(oldpar)
#' \dontshow{
#' data_set <- list(list(data_set_t[gr_label == 1, 1:5]),
#'                  list(data_set_t[gr_label == 2, 1:5]),
#'                  list(data_set_t[gr_label == 3, 1:5]))
#' # testing without parallel computing
#' res <- gmtFD(data_set, h_tukey_m, blocks_contrasts)
#' summary(res, digits = 3)
#' # testing with parallel computing
#' library(doParallel)
#' res <- gmtFD(data_set, h_tukey_m, blocks_contrasts, parallel = TRUE, n_cores = 2)
#' summary(res, digits = 3)}
#' # Dunnett's contrast matrix
#' h_dunnett <- GFDmcv::contr_mat(k, type = "Dunnett")
#' h_dunnett_m <- kronecker(h_dunnett, diag(p))
#' # vector of blocks of contrasts labels
#' blocks_contrasts <- rep(1:(nrow(h_dunnett_m) / p), each = p)
#' \donttest{
#' # testing without parallel computing
#' res <- gmtFD(data_set, h_dunnett_m, blocks_contrasts)
#' summary(res, digits = 3)
#' # plots for pointwise Hotelling's T^2-test statistics
#' oldpar <- par(mfrow = c(3, 1), mar = c(4, 2, 2, 0.1))
#' plot(ph_test_statistic(data_set, h_dunnett_m), type = "l",
#'      ylim = c(0, max(ph_test_statistic(data_set, h_dunnett_m))),
#'      main = "Global hypothesis", xlab = "Day")
#' plot(ph_test_statistic(data_set, matrix(h_dunnett_m[1, ], 1)), type = "l",
#'      ylim = c(0, max(ph_test_statistic(data_set, h_dunnett_m))),
#'      main = "Contrast 1", xlab = "Day")
#' plot(ph_test_statistic(data_set, matrix(h_dunnett_m[2, ], 1)), type = "l",
#'      ylim = c(0, max(ph_test_statistic(data_set, h_dunnett_m))),
#'      main = "Contrast 2", xlab = "Day")
#' par(oldpar)}
#' 
#' # two functional variables - temperature and precipitation
#' library(fda)
#' data_set_t <- t(CanadianWeather$dailyAv[,, "Temperature.C"])
#' data_set_p <- t(CanadianWeather$dailyAv[,, "Precipitation.mm"])
#' # number of samples
#' k <- 3
#' # number of variables
#' p <- 2
#' # preparing data set
#' gr_label <- rep(c(1, 2, 3), c(15, 15, 5))
#' data_set <- list(list(data_set_t[gr_label == 1, ], data_set_p[gr_label == 1, ]),
#'                  list(data_set_t[gr_label == 2, ], data_set_p[gr_label == 2, ]),
#'                  list(data_set_t[gr_label == 3, ], data_set_p[gr_label == 3, ]))
#' # trajectories of mean temperature and precipitation
#' oldpar <- par(mfrow = c(1, 2), mar = c(4, 4, 2, 0.1))
#' matplot(t(data_set_t), type = "l", col = gr_label, lty = 1,
#'         xlab = "Day", ylab = "Temperature (C)",
#'         main = "Canadian weather data set")
#' legend("bottom", legend = c("Western Canada", "Eastern Canada", "Northern Canada"),
#'        col = 1:3, lty = 1)
#' matplot(t(data_set_p), type = "l", col = gr_label, lty = 1,
#'         xlab = "Day", ylab = "Precipitation (mm)",
#'         main = "Canadian weather data set")
#' legend("topleft", legend = c("Western Canada", "Eastern Canada", "Northern Canada"),
#'        col = 1:3, lty = 1)
#' par(oldpar)
#' 
#' # Tukey's contrast matrix
#' h_tukey <- GFDmcv::contr_mat(k, type = "Tukey")
#' h_tukey_m <- kronecker(h_tukey, diag(p))
#' # vector of blocks of contrasts labels
#' blocks_contrasts <- rep(1:(nrow(h_tukey_m) / p), each = p)
#' \donttest{
#' # testing without parallel computing
#' res <- gmtFD(data_set, h_tukey_m, blocks_contrasts)
#' summary(res, digits = 3)}
#' # plots for pointwise Hotelling's T^2-test statistics
#' oldpar <- par(mfrow = c(2, 2), mar = c(4, 2, 2, 0.1))
#' plot(ph_test_statistic(data_set, h_tukey_m), type = "l",
#'      ylim = c(0, max(ph_test_statistic(data_set, h_tukey_m))),
#'      main = "Global hypothesis", xlab = "Day")
#' plot(ph_test_statistic(data_set, h_tukey_m[1:2, ]), type = "l",
#'      ylim = c(0, max(ph_test_statistic(data_set, h_tukey_m))),
#'      main = "Contrast 1", xlab = "Day")
#' plot(ph_test_statistic(data_set, h_tukey_m[3:4, ]), type = "l",
#'      ylim = c(0, max(ph_test_statistic(data_set, h_tukey_m))),
#'      main = "Contrast 2", xlab = "Day")
#' plot(ph_test_statistic(data_set, h_tukey_m[5:6, ]), type = "l",
#'      ylim = c(0, max(ph_test_statistic(data_set, h_tukey_m))),
#'      main = "Contrast 3", xlab = "Day")
#' par(oldpar)
#' 
#' # Dunnett's contrast matrix
#' h_dunnett <- GFDmcv::contr_mat(k, type = "Dunnett")
#' h_dunnett_m <- kronecker(h_dunnett, diag(p))
#' # vector of blocks of contrasts labels
#' blocks_contrasts <- rep(1:(nrow(h_dunnett_m) / p), each = p)
#' \donttest{
#' # testing without parallel computing
#' res <- gmtFD(data_set, h_dunnett_m, blocks_contrasts)
#' summary(res, digits = 3)}
#' # plots for pointwise Hotelling's T^2-test statistics
#' oldpar <- par(mfrow = c(3, 1), mar = c(4, 2, 2, 0.1))
#' plot(ph_test_statistic(data_set, h_dunnett_m), type = "l",
#'      ylim = c(0, max(ph_test_statistic(data_set, h_dunnett_m))),
#'      main = "Global hypothesis", xlab = "Day")
#' plot(ph_test_statistic(data_set, matrix(h_dunnett_m[1, ], 1)), type = "l",
#'      ylim = c(0, max(ph_test_statistic(data_set, h_dunnett_m))),
#'      main = "Contrast 1", xlab = "Day")
#' plot(ph_test_statistic(data_set, matrix(h_dunnett_m[2, ], 1)), type = "l",
#'      ylim = c(0, max(ph_test_statistic(data_set, h_dunnett_m))),
#'      main = "Contrast 2", xlab = "Day")
#' par(oldpar)
#' \donttest{
#' # testing with parallel computing
#' library(doParallel)
#' blocks_contrasts <- rep(1:(nrow(h_tukey_m) / p), each = p)
#' res <- gmtFD(data_set, h_tukey_m, blocks_contrasts, parallel = TRUE, n_cores = 2)
#' summary(res, digits = 3)}
#' 
#' @references Buhlmann P. (1998) Sieve bootstrap for smoothing in nonstationary time series.
#' Annals of Statistics 26, 48-83.
#'
#' Dunnett C. (1955) A multiple comparison procedure for comparing several treatments
#' with a control. Journal of the American Statistical Association 50, 1096-1121.
#'
#' Munko M., Ditzhaus M., Pauly M., Smaga L., Zhang J.T. (2023) General multiple tests for functional data. 
#' Preprint https://arxiv.org/abs/2306.15259
#'
#' Munko M., Ditzhaus M., Pauly M., Smaga L. (2024) Multiple comparison procedures for simultaneous inference 
#' in functional MANOVA. Preprint https://arxiv.org/abs/2406.01242
#'
#' Tukey J.W. (1953) The problem of multiple comparisons. Princeton University.
#'
#' @importFrom stats p.adjust quantile
#' @import doParallel
#' @import foreach
#' @import MASS
#' @import Matrix
#' @import GFDmcv
#' @import fda
#'
#' @export
gmtFD <- function(x, h, blocks_contrasts, n_boot = 1000, alpha = 0.05,                                         
                   parallel = FALSE, n_cores = NULL, multi_gen = FALSE) {
  if (parallel) {
    # #' @importFrom utils installed.packages
    # if (!("doParallel" %in% rownames(installed.packages()))) {
    #   stop("Please install package 'doParallel'")
    # }
    # require(foreach, quietly = TRUE)
    requireNamespace("doParallel", quietly = TRUE)
    requireNamespace("foreach", quietly = TRUE)
    if (is.null(n_cores)) {
      n_cores <- parallel::detectCores()
    } else if (!(n_cores > 1)) {
      stop("n_cores should be greater than 1")
    }
    cl <- parallel::makePSOCKcluster(n_cores)
    on.exit(parallel::stopCluster(cl))
    doParallel::registerDoParallel(cl)
  }
  
  # Sprawdzamy, czy wektor jest posortowany niemalejąco
  if (!all(diff(blocks_contrasts) >= 0)) {
    stop("Vector blocks_contrasts is not sorted in non-descending order.")
  }
  
  num_hyp <- length(unique(blocks_contrasts))
  mc <- mean_cov_m_cpp(x)
  gph <- gsph_f_m(h, mc)
  gph_d <- matrix(NA, num_hyp, 2)
  for (i in seq_len(num_hyp)) {
    gph_d[i, ] <- gsph_f_m(matrix(h[blocks_contrasts == i, ], nrow = sum(blocks_contrasts == i)), mc)
  }
  ntp <- mc$ntp
  n_vec <- mc$n_vec
  kk <- mc$kk
  nn <- mc$nn
  pp <- mc$pp
  if (!all.equal(rowSums(h), rep(0, nrow(h)))) {
    stop("matrix h is not a contrast matrix - sums in rows are not equal to zero")
  }
  
  # bootstrap
  if (!parallel) {
    # parametric bootstrap
    gph_pb <- matrix(NA, n_boot, 2)
    gph_pb_d_i <- matrix(NA, num_hyp, n_boot)
    gph_pb_d_s <- matrix(NA, num_hyp, n_boot)
    if (multi_gen == TRUE) {
      for (i_b in seq_len(n_boot)) {
        x_pb_l <- vector("list", kk)
        labels <- rep(1:pp, ntp)
        for (i_k in seq_len(kk)) {
          x_pb <- MASS::mvrnorm(n_vec[i_k], numeric(ntp * pp), mc$gr_cov[[i_k]])
          for (i_x in seq_len(pp)) {
            x_pb_l[[i_k]][[i_x]] <- x_pb[, labels == i_x]
          }
        }
        mc_p <- mean_cov_m_cpp(x_pb_l)
        gph_pb[i_b, ] <- gsph_f_m(h, mc_p)
        for (i in seq_len(num_hyp)) {
          temp <- gsph_f_m(matrix(h[blocks_contrasts == i, ], nrow = sum(blocks_contrasts == i)), mc_p)
          gph_pb_d_i[i, i_b] <- temp[1]
          gph_pb_d_s[i, i_b] <- temp[2]
        }
      }
    } else {
      x_pb_temp <- vector("list", kk)
      for (i_k in seq_len(kk)) {
        x_pb_temp[[i_k]] <- MASS::mvrnorm(n_vec[i_k] * n_boot, numeric(ntp * pp), mc$gr_cov[[i_k]])
      }
      for (i_b in seq_len(n_boot)) {
        x_pb_l <- vector("list", kk)
        labels <- rep(1:pp, ntp)
        for (i_k in seq_len(kk)) {
          x_pb <- x_pb_temp[[i_k]][((i_b - 1) * n_vec[i_k] + 1):(i_b * n_vec[i_k]), ]
          for (i_x in seq_len(pp)) {
            x_pb_l[[i_k]][[i_x]] <- x_pb[, labels == i_x]
          }
        }
        mc_p <- mean_cov_m_cpp(x_pb_l)
        gph_pb[i_b, ] <- gsph_f_m(h, mc_p)
        for (i in seq_len(num_hyp)) {
          temp <- gsph_f_m(matrix(h[blocks_contrasts == i, ], nrow = sum(blocks_contrasts == i)), mc_p)
          gph_pb_d_i[i, i_b] <- temp[1]
          gph_pb_d_s[i, i_b] <- temp[2]
        }
      }
    }
  } else {
    # parametric bootstrap
    if (multi_gen == TRUE) {
      RNGkind(kind = "Mersenne-Twister", normal.kind = "Inversion")
      test_stat_boot <- foreach(i_b = 1:n_boot, .combine = rbind,
                                .packages = c("MASS", "gmtFD", "Rcpp"),
                                .noexport = c("gr_cov_cpp")) %dopar%
        {
          x_pb_l <- vector("list", kk)
          labels <- rep(1:pp, ntp)
          for (i_k in seq_len(kk)) {
            x_pb <- MASS::mvrnorm(n_vec[i_k], numeric(ntp * pp), mc$gr_cov[[i_k]])
            for (i_x in seq_len(pp)) {
              x_pb_l[[i_k]][[i_x]] <- x_pb[, labels == i_x]
            }
          }
          
          mc_p <- mean_cov_m_cpp(x_pb_l)
          gph_pb <- gsph_f_m(h, mc_p)
          gph_pb_d_i <- numeric(num_hyp)
          gph_pb_d_s <- numeric(num_hyp)
          for (i in seq_len(num_hyp)) {
            temp <- gsph_f_m(matrix(h[blocks_contrasts == i, ], nrow = sum(blocks_contrasts == i)), mc_p)
            gph_pb_d_i[i] <- temp[1]
            gph_pb_d_s[i] <- temp[2]
          }
          
          c(gph_pb, gph_pb_d_i, gph_pb_d_s)
        }
      gph_pb <- test_stat_boot[, 1:2]
      gph_pb_d_i <- t(test_stat_boot[, 3:(num_hyp + 2)])
      gph_pb_d_s <- t(test_stat_boot[, (num_hyp + 3):(2 * num_hyp + 2)])
    } else {
      x_pb_temp <- vector("list", kk)
      for (i_k in seq_len(kk)) {
        x_pb_temp[[i_k]] <- MASS::mvrnorm(n_vec[i_k] * n_boot, numeric(ntp * pp), mc$gr_cov[[i_k]])
      }
      RNGkind(kind = "Mersenne-Twister", normal.kind = "Inversion")
      test_stat_boot <- foreach(i_b = 1:n_boot, .combine = rbind,
                                .packages = c("MASS", "gmtFD", "Rcpp"),
                                .noexport = c("gr_cov_cpp")) %dopar%
        {
          x_pb_l <- vector("list", kk)
          labels <- rep(1:pp, ntp)
          for (i_k in seq_len(kk)) {
            x_pb <- x_pb_temp[[i_k]][((i_b - 1) * n_vec[i_k] + 1):(i_b * n_vec[i_k]), ]
            for (i_x in seq_len(pp)) {
              x_pb_l[[i_k]][[i_x]] <- x_pb[, labels == i_x]
            }
          }
          
          mc_p <- mean_cov_m_cpp(x_pb_l)
          gph_pb <- gsph_f_m(h, mc_p)
          gph_pb_d_i <- numeric(num_hyp)
          gph_pb_d_s <- numeric(num_hyp)
          for (i in seq_len(num_hyp)) {
            temp <- gsph_f_m(matrix(h[blocks_contrasts == i, ], nrow = sum(blocks_contrasts == i)), mc_p)
            gph_pb_d_i[i] <- temp[1]
            gph_pb_d_s[i] <- temp[2]
          }
          
          c(gph_pb, gph_pb_d_i, gph_pb_d_s)
        }
      gph_pb <- test_stat_boot[, 1:2]
      gph_pb_d_i <- t(test_stat_boot[, 3:(num_hyp + 2)])
      gph_pb_d_s <- t(test_stat_boot[, (num_hyp + 3):(2 * num_hyp + 2)])
    }
  }
  gph_pb_d_b_i <- numeric(num_hyp)
  gph_pb_d_b_s <- numeric(num_hyp)
  gph_pb_d_cv_i <- numeric(num_hyp)
  gph_pb_d_cv_s <- numeric(num_hyp)
  for (i in seq_len(num_hyp)) {
    gph_pb_d_b_i[i] <- mean(gph_pb_d_i[i, ] > gph_d[i, 1])
    gph_pb_d_b_s[i] <- mean(gph_pb_d_s[i, ] > gph_d[i, 2])
    gph_pb_d_cv_i[i] <- quantile(gph_pb_d_i[i, ], 1 - alpha / num_hyp)
    gph_pb_d_cv_s[i] <- quantile(gph_pb_d_s[i, ], 1 - alpha / num_hyp)
  }
  
  # global hypothesis
  beta_i <- gph_pb_d_b_i[which.min(gph_pb_d_b_i)]
  beta_s <- gph_pb_d_b_s[which.min(gph_pb_d_b_s)]
  data_order_i <- t(apply(gph_pb_d_i, 1, sort))
  data_order_s <- t(apply(gph_pb_d_s, 1, sort))
  j_low_b_i <- ceiling(beta_i * n_boot) # - 1
  j_low_b_s <- ceiling(beta_s * n_boot) # - 1
  if (n_boot == j_low_b_i) {
    A_i <- gph_pb_d_i / data_order_i[, 1]
  } else {
    A_i <- gph_pb_d_i / data_order_i[, n_boot - j_low_b_i]
  }
  if (n_boot == j_low_b_s) {
    A_s <- gph_pb_d_s / data_order_s[, 1]
  } else {
    A_s <- gph_pb_d_s / data_order_s[, n_boot - j_low_b_s]
  }
  mgph_p_val_gl_i <- mean(apply(A_i, 2, max) > 1)
  mgph_p_val_gl_s <- mean(apply(A_s, 2, max) > 1)
  
  res_global <- data.frame(gsph_ts = gph,
                           gph_pval = c(mean(gph_pb[, 1] > gph[1]), mean(gph_pb[, 2] > gph[2])),
                           mgsph_ts = c(max(gph_d[, 1] / crit_values_ls(gph_pb_d_i, alpha)),
                                        max(gph_d[, 2] / crit_values_ls(gph_pb_d_s, alpha))),
                           mgph_pval = c(mgph_p_val_gl_i, mgph_p_val_gl_s))
  
  # multi
  mgph_p_val_i <- numeric(num_hyp)
  mgph_p_val_s <- numeric(num_hyp)
  for (i_beta in seq_len(num_hyp)) {
    beta_i <- gph_pb_d_b_i[i_beta]
    beta_s <- gph_pb_d_b_s[i_beta]
    data_order_i <- t(apply(gph_pb_d_i, 1, sort))
    data_order_s <- t(apply(gph_pb_d_s, 1, sort))
    j_low_b_i <- ceiling(beta_i * n_boot)
    j_low_b_s <- ceiling(beta_s * n_boot)
    if (n_boot == j_low_b_i) {
      A_i <- gph_pb_d_i / data_order_i[, 1]
    } else {
      A_i <- gph_pb_d_i / data_order_i[, n_boot - j_low_b_i]
    }
    if (n_boot == j_low_b_s) {
      A_s <- gph_pb_d_s / data_order_s[, 1]
    } else {
      A_s <- gph_pb_d_s / data_order_s[, n_boot - j_low_b_s]
    }
    mgph_p_val_i[i_beta] <- mean(apply(A_i, 2, max) > 1)
    mgph_p_val_s[i_beta] <- mean(apply(A_s, 2, max) > 1)
  }
  
  res_multi <- data.frame(ts = gph_d,
                          gph_cv = matrix(c(gph_pb_d_cv_i, gph_pb_d_cv_s), ncol = 2),
                          gph_pval = matrix(c(p.adjust(gph_pb_d_b_i, "bonferroni"), 
                                              p.adjust(gph_pb_d_b_s, "bonferroni")), ncol = 2),
                          mgph_cv = matrix(c(crit_values_ls(gph_pb_d_i, alpha), 
                                             crit_values_ls(gph_pb_d_s, alpha)), ncol = 2),
                          mgph_pval = matrix(c(mgph_p_val_i, mgph_p_val_s), ncol = 2))
  colnames(res_multi) <- c("gph_ts", "sph_ts",
                           "gph_cv", "sph_cv",
                           "gph_pval", "sph_pval",
                           "mgph_cv", "msph_cv",
                           "mgph_pval", "msph_pval")
  
  result <- list(res_global = res_global, res_multi = res_multi,
                 k = kk, p = pp, ntp = ntp, n = n_vec,
                 h = h, blocks_contrasts = blocks_contrasts, n_boot = n_boot, 
                 alpha = alpha, multi_gen = multi_gen)
  class(result) <- "multifmanova"
  return(result)
}