R/mopet.r

Defines functions summary.fc print.fc fourth.corner.0 fourth.corner summary.mopet print.mopet mopet

Documented in fourth.corner mopet print.fc print.mopet summary.fc summary.mopet

#' Modified permutation test with weighted mean of species attributes
#' 
#' Function performing modified permutation test to calculate significance of relationship between weighted mean of species attributes and other variables.
#' 
#' @name mopet
#' @param M An object of the class \code{wm} 
#' @param env Vector or matrix with variables. See details.
#' @param method Statistical method used to analyse the relationship between M (of class \code{wm} and env), partial match to \code{'lm'}, \code{'aov'}, \code{'cor'}, \code{'kruskal'}, \code{'slope'} and \code{'fourthcorner'}.
#' @param cor.coef Correlation coefficient in case of \code{method = 'cor'}. Partial match to 'pearson', 'spearman' and 'kendal'.
#' @param dependence Should M be dependent variable and env independent (\code{'M ~ env'}), or opposite? Applicable only for \code{method = 'lm'}. Partial match to \code{'M ~ env'} and \code{'env ~ M'}, so to write \code{dep = 'M'} is enough.
#' @param permutations Number of permutations.
#' @param test Vector of character values. Which test should be conducted? Partial match of vector with \code{'standard'}, \code{'modified'}, \code{'sequential'} and \code{'all'}.
#' @param parallel NULL (default) or integer number. Number of cores for parallel calculation of modified permutation test. Maximum number of cores should correspond to number of available cores on the processor.
#' @param object,x object of the class \code{"mopet"}, generated by function \code{mopet}.
#' @param digits number of digits reported by \code{print} method on object of \code{"mopet"} class (default is 3).
#' @param ... Other arguments for \code{print} or \code{summary} functions (not implemented).

#' @export
#' @examples 
#' data (vltava)
#' mean.eiv <- wm (vltava$spe, vltava$ell)
#' mopet (mean.eiv, vltava$env$pH, perm = 49, test = 'all')
#' mopet (mean.eiv, vltava$env$pH, perm = 49, test = 'stand', method = 'cor', cor.coef = 'spearm')
#' @details
#' Currently implemented statistical methods are correlation (\code{'cor'}), linear regression (\code{method = 'lm'}), ANOVA (\code{'aov'}) and Kruskal-Wallis test (\code{'kruskal'}).
#'
#' Argument \code{env} can be vector or matrix with one column. Only in case of linear regression (\code{method = 'lm'}) is possible to use matrix with several variables, which will be all used as independent variables in the model. For ANOVA and Kruskal-Wallis test, make sure that 'env' is \code{factor} (warning will be returned if this is not the case, but the calculation will be conducted). 
#' 
#' Difference between \code{method = 'lm'} and \code{'aov'} is in the format of summary tables, returned by \code{summary.mopet} function. In case of 'aov', this summary is expressed in the traditional language of ANOVA rather than linear models.
#' 
#' Both \code{method = 'lm'} and \code{'slope'} are based on linear regression and calculated by function \code{\link{lm}}, but differ by test statistic: while 'lm' is using F value and is testing the strength of the regression (measured by r2), 'slope' is using the slope of the regression line (b). This statistic is added here for comparison with the fourth corner method. While r2 (and r) is influenced by the issue of compositional autocorrelation, slope of regression is not.
#' 
#' Specific issue related to weighted mean is the case of missing species attributes. In current implementation, species with missing species attributes are removed from sample x species matrix prior to permutation of species attributes among species. 
#' @return  Function \code{mopet} returns list of the class \code{"mopet"} (with \code{print} and \code{summary} methods), which contains the following items:
#' \itemize{
#'  \item \code{real.summaries} summary of the method
#'  \item \code{coefs} model coefficients
#'  \item \code{stat} test statistic
#'  \item \code{orig.P} P-values from the original (parametric) test
#'  \item \code{perm.P} P-values from the not-modified permutation test (permuted are the whole rows of matrix M)
#'  \item \code{modif.P} P-values from the modified permutation test (permuted are species attributes in object M)
#'  \item \code{seq.P} P-values from sequential test (\code{max (perm.P, modif.P)})
#'  \item \code{permutations} number of permutations
#'  \item \code{method} method of calculation
#'  \item \code{dependence} dependence (in case of \code{method = 'lm'})
#'  }
#' @seealso \code{\link{wm}}

#' @export
mopet <- function (M, env, method = c('lm'), cor.coef = c('pearson'), dependence = "M ~ env", permutations = 499, test = "modified", parallel = NULL)
{
  METHOD <- c('lm', 'aov', 'cor', 'kruskal', 'slope', 'fourthcorner')
  COR.COEF <- c('pearson', 'spearman', 'kendall')
  TEST <- c('none', 'standard', 'modified', 'sequential', 'all')
  DEPENDENCE <- c("M ~ env", "env ~ M")
  method <- match.arg (method, METHOD)
  cor.coef <- match.arg (cor.coef, COR.COEF)
  test <- match.arg (test, TEST, several.ok = T)
  dependence <- match.arg (dependence, DEPENDENCE)
  if (!is.wm (M) & ("modified" %in% test || "sequential" %in% test || "all" %in% test)) stop ("Object M must be of 'wm' class")
  env <- as.matrix (env)
  sweeping.sign <- if (method == 'cor' & cor.coef == 'spearman') "<=" else ">="
  orig.P <- NULL
  perm.P <- NULL
  modif.P <- NULL
  seq.P <- NULL
  sitspe <- attr (M, 'sitspe')
  speatt <- attr (M, 'speatt')
  fun <-  switch (method, 
                  lm = if (dependence == 'M ~ env') expression (lm (M ~ ., as.data.frame (as.matrix (env)))) else expression (lm (env ~ ., as.data.frame (as.matrix (M)))),
                  aov = expression (apply (as.matrix (M), 2, FUN = function (m) aov (m ~ env))),
                  cor = expression (apply (as.matrix (M), 2, FUN = function (m) cor.test (m, env, method = cor.coef))),
                  kruskal = expression (apply (as.matrix (M), 2, FUN = function (m) kruskal.test (m, env))),
                  #slope = expression (apply (as.matrix (M), 2, FUN = function (m) lm (scale (m, center = weighted.mean (m, rowSums (sitspe))) ~ ., as.data.frame (scale (env), center = weighted.mean (env, rowSums (sitspe))), weights = rowSums (sitspe)))),
                  slope = expression ({speatt <- attr (M, 'speatt'); sitspe <- attr (M, 'sitspe'); apply (speatt, 2, FUN = function (q) lm (wm (sitspe, scale (q, center = weighted.mean (q, colSums (sitspe)))) ~ scale (env, center = weighted.mean (env, rowSums (sitspe)))))}),
                  fourthcorner = expression (weimea::fourth.corner (R = R, L = L, Q = Q)))
  summ <- switch (method,
                  lm = expression (summary (obj)),
                  aov = expression (lapply (obj, FUN = function (x) summary (x))),
                  cor = expression (lapply (obj, FUN = function (x) x)), 
                  kruskal = expression (lapply (obj, FUN = function (x) x)),
                  slope = expression (lapply (obj, FUN = function (x) summary (x))),
                  fourthcorner = expression (obj))
  coefs <- switch (method, 
                   lm = expression (coef (obj)),
                   aov = expression (lapply (obj, FUN = function (x) coef (x))),
                   cor = expression (lapply (obj, FUN = function (x) x$estimate)),
                   kruskal = expression (lapply (obj, FUN = function (x) x$parameter)),
                   slope = expression (lapply (obj, FUN = function (x) coef (x))),
                   fourthcorner = expression (obj))
  stat <- switch (method,
                  lm = expression ({temp <- anova (obj)$"F value"[1]; names (temp) <- "F value"; temp}),
                  aov = expression (lapply (obj, FUN = function (x) {temp <- summary (x)[[1]]$"F value"[1]; names (temp) <- "F value"; temp})),
                  cor = expression (lapply (obj, FUN = function (x) x$statistic)),
                  kruskal = expression (lapply (obj, FUN = function (x) x$statistic)),
                  slope = expression (lapply (obj, FUN = function (x) {temp <- coef (summary (x))[2,1]; names (temp) <- "b"; temp})),
                  fourthcorner = expression (obj))
  signi <- switch (method,
                   lm = expression (anova (obj)$"Pr(>F)"[1]),
                   aov = expression (lapply (obj, FUN = function (x) summary (x)[[1]]$"Pr(>F)"[1])),
                   cor = expression (lapply (obj, FUN = function (x) x$p.value)),
                   kruskal = expression (lapply (obj, FUN = function (x) x$p.value)),
                   slope = expression (lapply (obj, FUN = function (x) coef (summary (x))[2,4])),
                   fourthcorner = NULL)
  tail <- switch (method,
                  lm = 'one',
                  aov = 'one',
                  cor = 'two',
                  kruskal = 'one',
                  slope = 'two',
                  fourthcorner = 'two')
  res.real.fun <- if (method == 'fourthcorner') list (with (list (R = env, L = sitspe, Q = speatt), eval (fun))) else 
                  if (method %in% c('aov', 'kruskal', 'cor', 'slope') || (method == 'lm' & dependence == 'env ~ M')) apply (as.matrix (env), 2, FUN = function (i) with (list (M = M, env = i), eval (fun))) else
                  if (method == 'lm' & dependence == 'M ~ env') apply (as.matrix (M), 2, FUN = function (i) with (list (M = i, env = env), eval (fun))) 
  real.summaries <- lapply (res.real.fun, FUN = function (obj) with (list (obj = obj), eval (summ)))
  coefs <- lapply (res.real.fun, FUN = function (obj) with (list (obj = obj), eval (coefs)))
  if (method != 'fourthcorner') orig.P <- lapply (res.real.fun, FUN = function (obj) with (list (obj = obj), eval (signi)))
  res.real.stat <- lapply (res.real.fun, FUN = function (obj) with (list (obj = obj), eval (stat)))
  
  if ('modified' %in% test || 'all' %in% test || 'sequential' %in% test) 
  {
    res.temp.stat.modif <-  
      if (method == 'fourthcorner') randomize (M, permutations = permutations, parallel = parallel, FUN = function (mat) {lapply (list (with (list (R = env, L = attr (mat, 'sitspe'), Q = attr (mat, 'speatt')), eval (fun))), FUN = function (obj) with (list (obj = obj), eval (stat)))}) else
      if (method %in% c('aov', 'kruskal', 'cor', 'slope') || (method == 'lm' & dependence == 'env ~ M')) randomize (M, permutations = permutations, parallel = parallel, FUN = function (mat) lapply (apply (as.matrix (env), 2, FUN = function (i) with (list (M = mat, env = i), eval (fun))), FUN = function (obj) with (obj, eval (stat)))) else
      if (method == 'lm' & dependence == 'M ~ env')  randomize (M, permutations = permutations, parallel = parallel, FUN = function (mat) lapply (apply (as.matrix (mat), 2, FUN = function (i) with (list (M = i, env = env), eval (fun))), FUN = function (obj) with (obj, eval (stat)))) 
    res.temp.stat.modif <- rbind (matrix (unlist (res.temp.stat.modif), nrow = length (res.temp.stat.modif), byrow = T), unlist (res.real.stat))      
    modif.P <- (colSums (sweep (
      if (tail == 'one') res.temp.stat.modif else abs (res.temp.stat.modif), 2,
      if (tail == 'one') unlist (res.real.stat) else abs (unlist (res.real.stat)), sweeping.sign)))/(permutations+1)
  }
  
  if ('standard' %in% test || 'all' %in% test || 'sequential' %in% test)
  {
    res.temp.stat.stand <- 
      if (method == 'fourthcorner') lapply (lapply (1:permutations, FUN = function (i) as.matrix (env)[sample (nrow (env)),]), FUN = function (env) lapply (with (list (R = env, L = sitspe, Q = speatt), eval (fun)), FUN = function (obj) with (list (obj = obj), eval (stat)))) else
      if (method %in% c('aov', 'kruskal', 'cor', 'slope') || (method == 'lm' & dependence == 'env ~ M'))      lapply (lapply (1:permutations, FUN = function (i) as.matrix (env)[sample (nrow (env)),]), FUN = function (env.r) lapply (apply (as.matrix (env.r), 2, FUN = function (i) with (list (M = M, env = i), eval (fun))), FUN = function (obj) with (obj, eval (stat)))) else
      if (method == 'lm' & dependence == 'M ~ env') lapply (lapply (1:permutations, FUN = function (i) as.matrix (M)[sample (nrow (M)),]), FUN = function (mat) lapply (apply (as.matrix (mat), 2, FUN = function (i) with (list (M = i, env = env), eval (fun))), FUN = function (obj) with (obj, eval (stat))))
    res.temp.stat.stand <- rbind (matrix (unlist (res.temp.stat.stand), nrow = length (res.temp.stat.stand), byrow = T), unlist (res.real.stat))      
    perm.P <- (colSums (sweep (
      if (tail == 'one') res.temp.stat.stand else abs (res.temp.stat.stand), 2,
      if (tail == 'one') unlist (res.real.stat) else abs (unlist (res.real.stat)), sweeping.sign)))/(permutations+1)
  }
  if ('sequential' %in% test || 'all' %in% test) seq.P <- apply (cbind (perm.P, modif.P), 1, max)    
  if (!('modified' %in% test) & !('standard' %in% test) & !('all' %in% test)) {modif.P <- NULL; perm.P <- NULL}
  result <- list (real.summaries = real.summaries, coefs = coefs, stat = res.real.stat, orig.P = orig.P, perm.P = perm.P, modif.P = modif.P, seq.P = seq.P, permutations = permutations, method = method, dependence = dependence)
  class (result) <- 'mopet'
  return (result)
}

#' @export
#' @rdname mopet
print.mopet <- function (x, digits = max(3, getOption("digits") - 3), ...)
{
  symnum.pval <- function (pval) symnum( pval, corr = FALSE, na = FALSE, cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1), symbols = c("***", "**", "*", ".", " "))
  mat.coef <- matrix (format (unlist (x$coefs), digits = digits), nrow = if (x$method == 'lm' || x$method == 'slope') length (x$coefs) else length (unlist (x$coefs)))
  rownames (mat.coef) <- if (x$method == 'lm') names (x$coefs) else 
    if (x$method == 'fourthcorner') apply (expand.grid (rownames (x$coefs[[1]]), colnames (x$coefs[[1]])), 1, FUN = function (na) paste(rev (na), collapse = ' / ')) else
    apply (expand.grid (names (x$coefs[[1]]), names (x$coefs)), 1, FUN = function (na) paste(rev (na), collapse = ' / '))
  colnames (mat.coef) <- if (x$method == 'lm') names (x$coef[[1]]) else if (x$method == 'fourthcorner') '4thcorner' else names (x$coef[[1]][[1]])
  mat.stat <- matrix (format (unlist (x$stat), digits = digits))
  colnames (mat.stat) <- if (x$method == 'lm') names (x$stat[[1]]) else names (x$stat[[1]][[1]])
  if (!is.null (x$orig.P)) mat.orig.P <- cbind (orig.P = format.pval (unlist (x$orig.P), digits = digits), symnum.pval (unlist (x$orig.P)))
  if (!is.null (x$perm.P)) mat.perm.P <- cbind (perm.P = format.pval (unlist (x$perm.P), digits = digits), symnum.pval (unlist (x$perm.P)))
  if (!is.null (x$modif.P)) mat.modif.P <- cbind (modif.P = format.pval (unlist (x$modif.P), digits = digits), symnum.pval (unlist (x$modif.P)))
  if (!is.null (x$seq.P)) mat.seq.P <- cbind (seq.P = format.pval (x$seq.P, digits = digits), symnum.pval (unlist (x$seq.P)))
  print.default (cbind (mat.coef, if (x$method != 'fourthcorner') mat.stat, if (!is.null (x$orig.P)) mat.orig.P, if (!is.null (x$perm.P)) mat.perm.P, if (!is.null (x$modif.P)) mat.modif.P, if (!is.null (x$seq.P)) mat.seq.P), quote = F, right = T, ...)
}

#' @export
#' @rdname mopet
summary.mopet <- function (object, ...)
{
  len <- length (object$real.summaries)
  cat ('\nSummary of mopet function:\n\n')
  for (i in seq (1, len))
  {
    cat ('\nOriginal result for variable', names (object$real.summaries)[i], ':\n')
    print (object$real.summaries[[i]])
    cat ('\n------------------------------------------------')
    if (!is.null (object$perm.P[i])) cat ('\nStandard permutation test: P = ', object$perm.P[i])
    if (!is.null (object$modif.P[i]))cat ('\nModified permutation test: P = ', object$modif.P[i])
    if (!is.null (object$seq.P[i]))  cat ('\nSequential permutation test: P = ', object$seq.P[i])
    cat ('\nPermutation results based on', object$permutations, 'permutations')
    cat ('\n************************************************\n')
  }

}

#' @export
fourth.corner <- function (R, L, Q)
{
  R <- as.matrix (R)
  L <- as.matrix (L)
  Q <- as.matrix (Q)
  
  missing.Q <- ifelse (any (is.na (Q)), TRUE, FALSE)
  missing.R <- ifelse (any (is.na (R)), TRUE, FALSE)
    
  if (!missing.Q & !missing.R) r.h <- fourth.corner.0 (R, L, Q)
  if (missing.Q & !missing.R)
    r.h <- apply (Q, 2, FUN = function (q) {L <- L[,!is.na (q)]; q <- q[!is.na (q)]; fourth.corner.0 (R, L, q)})
  if (!missing.Q & missing.R)
    r.h <- t(apply (R, 2, FUN = function (r) {L <- L[!is.na (r),]; r <- r[!is.na (r)]; fourth.corner.0 (r, L, Q)}))
  if (missing.Q & missing.R)
    r.h <- apply (Q, 2, FUN = function (q) apply (R, 2, FUN = function (r) {L <- L[!is.na (r), !is.na (q)]; q <- q[!is.na (q)]; r <- r[!is.na (r)]; fourth.corner.0 (r, L, q)}))
  if (is.null (colnames (r.h))) colnames (r.h) <- colnames (Q)
  if (is.null (rownames (r.h))) rownames (r.h) <- colnames (R)
  class (r.h) <- 'fc'
  return (r.h)
}  
  
  
fourth.corner.0 <- function (R, L, Q)  # this can be used only in case of no missing values in Q
{
  W.k <- diag (colSums (L))  
  Q <- as.matrix (Q)
  R <- as.matrix (R)
  L <- as.matrix (L)
  W.n <- diag (rowSums (L))
  I.k <- matrix (rep (1, ncol (L)))
  I.n <- matrix (rep (1, nrow (L)))
  I.g <- matrix (rep (1, ncol (Q)))
  trace.m <- function (M) sum (diag (M))
  
  Q.std.1 <- (Q - I.k %*% t(I.k) %*% W.k %*% Q * (1/trace.m (W.k)))
  Q.std.2 <- (I.k %*% ((t(I.k) %*% W.k %*% (Q * Q) - ((t(I.k) %*% W.k %*% Q) * (t(I.k) %*% W.k %*% Q))/trace.m (W.k)) / trace.m (W.k)) ^ 0.5)
  Q.std <- Q.std.1 / Q.std.2
  
  R.std.1 <- (R - I.n %*% t(I.n) %*% W.n %*% R * (1/trace.m (W.n)))
  R.std.2 <- (I.n %*% ((t(I.n) %*% W.n %*% (R * R) - ((t(I.n) %*% W.n %*% R) * (t(I.n) %*% W.n %*% R))/trace.m (W.n)) / trace.m (W.n)) ^ 0.5)
  R.std <- R.std.1 / R.std.2
  
  #Q.L <- apply (Q.std, 2, FUN = function (Q.std.1) I.n %*% t (Q.std.1) * L)
  X.P <- L %*% Q.std / W.n %*% I.n %*% t(I.g)
  
  r.h <- t(apply (R.std, 2, FUN = function (R.h.std) (t(R.h.std) %*% W.n %*% R.h.std) ^-1 %*% t(R.h.std) %*% W.n %*% X.P))
  colnames (r.h) <- colnames (Q)
  return (r.h)
}

#' @export
print.fc <- function (x, ...) print.default (as.matrix (unclass (x)))

#' @export
summary.fc <- function (object, ...) print (object)
zdealveindy/weimea documentation built on Dec. 5, 2017, 11:25 p.m.