corls_test.R 4.69 KB
Newer Older
1
#' @include covls.R
2
#' @import permute
3 4 5 6 7 8
#' @author Christelle Gonindard-Melodelima
#' @author Eric Coissac
NULL

#' Generate permutation matrix according to a schema.
#'
9 10 11 12 13
#' The permutation schema is defined using the `how` function.
#' The implementation of this function is inspired
#' from the VEGAN package and reproduced here to avoid an extra
#' dependency on an hidden vegan function.
#'
Eric Coissac committed
14 15 16 17 18 19 20 21
#' @param permutations a list of control values for the permutations as returned
#'              by the function \code{\link[permute]{how}}, or the number of
#'              permutations required.
#' @param n numeric; the number of observations in the sample set.
#'              May also be any object that nobs knows about;
#'              see \code{\link[permute]{nobs}} methods.
#' @param strata A factor, or an object that can be coerced to a
#'               factor via as.factor, specifying the strata for permutation.
22
#'
Eric Coissac committed
23
#' @note Internal function do not use.
24
#'
Eric Coissac committed
25 26
#' @rdname internal.getPermuteMatrix
.getPermuteMatrix = function(permutations, n, strata = NULL)
27
{
Eric Coissac committed
28 29
  if (length(permutations) == 1) {
    permutations <- permute::how(nperm = permutations)
30 31
  }
  if (!missing(strata) && !is.null(strata)) {
Eric Coissac committed
32
    if (inherits(permutations, "how") && is.null(permute::getBlocks(permutations)))
33
      setBlocks(permutations) <- strata
34
  }
Eric Coissac committed
35 36
  if (inherits(permutations, "how"))
    permutations <- permute::shuffleSet(n, control = permutations)
37
  else {
Eric Coissac committed
38
    if (!is.integer(permutations) && !all(permutations == round(permutations)))
39 40
      stop("permutation matrix must be strictly integers: use round()")
  }
Eric Coissac committed
41 42 43 44
  if (is.null(attr(permutations, "control")))
    attr(permutations, "control") <- structure(list(within = list(type = "supplied matrix"),
                                            nperm = nrow(permutations)), class = "how")
  permutations
45 46 47 48
}



Eric Coissac committed
49 50 51 52 53
#' Monte-Carlo Test on the sum of the singular values of a procustean rotation.
#'
#' performs a Monte-Carlo Test on the sum of the singular values of a
#' procustean rotation (see \code{\link[ade4]{procuste.rtest}}).
#'
54
#' @param ...   the set of matrices or a \code{\link[ProcMod]{procmod_frame}}
Eric Coissac committed
55 56 57 58
#'              object.
#' @param permutations a list of control values for the permutations as returned
#'              by the function \code{\link[permute]{how}}, or the number of
#'              permutations required.
59 60
#' @param p_adjust_method the multiple test correction method used
#'              to adjust p values. \code{p_adjust_method}
61
#'              belongs one of the folowing values: \code{"holm"},
Eric Coissac committed
62 63 64 65 66 67 68
#'              \code{"hochberg"}, \code{"hommel"}, \code{"bonferroni"},
#'              \code{"BH"}, \code{"BY"}, \code{"fdr"},  \code{"none"}.
#'              The default is,set to \code{"holm"}.
#'
#' @references {
#'  \insertRef{Jackson:95:00}{ProcMod}
#' }
69
#'
70 71 72 73 74 75 76 77 78 79 80 81
#' @examples
#' A <- simulate_matrix(10,3)
#' B <- simulate_matrix(10,5)
#' C <- simulate_correlation(B,10,r2=0.6)
#'
#' # Computes the correlation matrix
#' data <- procmod_frame(A = A, B = B, C = C)
#'
#' corls_test(A, B, permutations = 100)
#' corls_test(B, C, permutations = 100)
#' corls_test(data, permutations = 100)
#'
82 83
#' @seealso \code{\link[stats]{p.adjust}}
#'
84 85 86
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
87
corls_test <- function(...,
88
                       permutations = permute::how(nperm = 999),
89
                       p_adjust_method="holm") {
90 91 92 93 94
  eps <- sqrt(sqrt(.Machine$double.eps))
  xs <- list(...)

  if (length(xs) == 1) {
    x <- xs[[1]]
95
    if (is_procmod_frame(x)) {
96 97
      xs <- x
    } else {
98
      xs <- procmod_frame(x)
99 100 101
    }
  }
  else {
102
    xs <- as_procmod_frame(xs)
103 104 105 106 107 108 109
  }

  x_names <- names(xs)

  xs <- ortho(xs)


110
  cov <- varls(xs, nrand = 0)
111 112 113 114 115
  lcov <- cov - eps
  ngreater <- array(0,dim = dim(cov))

  n <- nrow(xs)
  nx <- length(xs)
Eric Coissac committed
116
  pmatrix <- .getPermuteMatrix(permutations, n)
117 118 119 120 121 122 123 124 125 126 127 128 129 130 131

  if (ncol(pmatrix) != n) {
    stop(gettextf(
      "'permutations' have %d columns, but data have %d observations",
      ncol(pmatrix), n
    ))
  }

  npermutation <- nrow(pmatrix)
  for (i in seq_len(npermutation)) {
    ps <- sample(1:npermutation,
      size = nx,
      replace = FALSE
    )

132
    rcov = varls(as_procmod_frame(
133 134 135 136 137
        lapply(
          1:nx,
          function(j) xs[[j]][pmatrix[ps[j], ], ]
        )
      ),
138
      nrand = 0
139 140 141 142 143 144 145 146 147
      )
    ngreater <- ngreater + (
       rcov >= lcov)
  }

  p_values <- ngreater / npermutation
  diag(p_values) <- 0

  c_p_values <- p.adjust(p_values[upper.tri(p_values,diag = FALSE)],
148
           method = p_adjust_method,
149 150 151 152 153 154 155 156 157 158
           n = (nx - 1) * nx / 2)

  p_values[upper.tri(p_values,diag = FALSE)] <- c_p_values
  p_values <- as.matrix(Matrix::forceSymmetric(p_values, uplo = "U"))
  colnames(p_values) <- x_names
  rownames(p_values) <- x_names

  p_values

}