anova.pm.R 2.31 KB
Newer Older
1 2 3 4 5 6 7
#' @include mprocuste.R

#' Anova function for `pm` objects.
#'
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
8
anova.pm = function(m, permutations = how(nperm = 999),...) {
9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

  mf   = m$model
  coef = m$coefficients

  sdxy = sqrt(diag(m$cov))

  mt = terms(mf)

  reponse = attr(mt,"response")
  explicatives = attr(mt,'term.labels')
  nXs = length(explicatives)

  coef.norm = coef * sdxy[explicatives] / sdxy[reponse]


24
  cors = m$cov / outer(sdxy,sdxy)
25

26 27
  #print(reponse)
  #print(explicatives)
28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

  cors.Y = cors[reponse,explicatives]

  SCT  = m$SCT
  SCR = sum(m$residuals^2)

  SCE = SCT * coef.norm * cors.Y

  names(SCE) = explicatives

  SCET = sum(SCE)

  SCIT = SCT - SCET - SCR

  SumSq  = c(SCE,Residuals=SCR)
  Df     = c(rep(1,nXs),m$df.residual)
  MeanSq = SumSq/Df
  Fvalue = c(MeanSq[-length(MeanSq)]/MeanSq[length(MeanSq)],NA)
  Pvalue = 1-pf(Fvalue,
                Df[-length(MeanSq)],
                Df[length(MeanSq)])

  result = data.frame(Df,
                      SumSq,
                      MeanSq,
                      Fvalue,
                      Pvalue)

  colnames(result)=c("Df","Sum Sq","Mean Sq","F value","Pr(>F)")


  attr(result, "heading")= paste("Analysis of Variance Table\n\nResponse:",
                                 all.vars(m$terms)[1]
  )
  class(result)=c("anova",class(result))

  return(result)
}
66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94

#' Generate permutation matrix according to a schema.
#' 
#' The permutation schema is defined using the `how` function.
#' The implementation of this function is currectly extraced 
#' from the VEGAN package and just copied pasted here to avoid
#' dependency on an hidden vegan function.
#'
getPermuteMatrix = function (perm, N, strata = NULL) 
{
  if (length(perm) == 1) {
    perm <- how(nperm = perm)
  }
  if (!missing(strata) && !is.null(strata)) {
    if (inherits(perm, "how") && is.null(getBlocks(perm))) 
      setBlocks(perm) <- strata
  }
  if (inherits(perm, "how")) 
    perm <- shuffleSet(N, control = perm)
  else {
    if (!is.integer(perm) && !all(perm == round(perm))) 
      stop("permutation matrix must be strictly integers: use round()")
  }
  if (is.null(attr(perm, "control"))) 
    attr(perm, "control") <- structure(list(within = list(type = "supplied matrix"), 
                                            nperm = nrow(perm)), class = "how")
  perm
}