corls_test.R 4.37 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
#' @seealso \code{\link[stats]{p.adjust}}
#'
72 73 74
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
75
corls_test <- function(...,
76
                       permutations = permute::how(nperm = 999),
77
                       p_adjust_method="holm") {
78 79 80 81 82
  eps <- sqrt(sqrt(.Machine$double.eps))
  xs <- list(...)

  if (length(xs) == 1) {
    x <- xs[[1]]
83
    if (is_procmod_frame(x)) {
84 85
      xs <- x
    } else {
86
      xs <- procmod_frame(x)
87 88 89
    }
  }
  else {
90
    xs <- as_procmod_frame(xs)
91 92 93 94 95 96 97 98 99 100 101 102 103
  }

  x_names <- names(xs)

  xs <- ortho(xs)


  cov <- varls(xs, nperm = 0)
  lcov <- cov - eps
  ngreater <- array(0,dim = dim(cov))

  n <- nrow(xs)
  nx <- length(xs)
Eric Coissac committed
104
  pmatrix <- .getPermuteMatrix(permutations, n)
105 106 107 108 109 110 111 112 113 114 115 116 117 118 119

  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
    )

120
    rcov = varls(as_procmod_frame(
121 122 123 124 125 126 127 128 129 130 131 132 133 134 135
        lapply(
          1:nx,
          function(j) xs[[j]][pmatrix[ps[j], ], ]
        )
      ),
      nperm = 0
      )
    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)],
136
           method = p_adjust_method,
137 138 139 140 141 142 143 144 145 146
           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

}