anova.pm.R 2.3 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

  cors.Y = cors[reponse,explicatives]

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

34
  SCE = SCT * abs(coef.norm * cors.Y)
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

  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

#' Generate permutation matrix according to a schema.
68
#'
69
#' The permutation schema is defined using the `how` function.
70
#' The implementation of this function is currectly extraced
71 72 73
#' from the VEGAN package and just copied pasted here to avoid
#' dependency on an hidden vegan function.
#'
74
getPermuteMatrix = function (perm, N, strata = NULL)
75 76 77 78 79
{
  if (length(perm) == 1) {
    perm <- how(nperm = perm)
  }
  if (!missing(strata) && !is.null(strata)) {
80
    if (inherits(perm, "how") && is.null(getBlocks(perm)))
81 82
      setBlocks(perm) <- strata
  }
83
  if (inherits(perm, "how"))
84 85
    perm <- shuffleSet(N, control = perm)
  else {
86
    if (!is.integer(perm) && !all(perm == round(perm)))
87 88
      stop("permutation matrix must be strictly integers: use round()")
  }
89 90
  if (is.null(attr(perm, "control")))
    attr(perm, "control") <- structure(list(within = list(type = "supplied matrix"),
91 92 93 94
                                            nperm = nrow(perm)), class = "how")
  perm
}