permuRDAv1.6.R

#documentation start
#=============================================================================
# File data
# creator: Christiane Hassenrück
# acknowledgements: Alban Ramette, Marianne Jacob
# primary authority: Christiane Hassenrück
# other authorities: 
#=============================================================================
# File contents
# wrapper for rda function to calculate R², df, F, p for each explanatory variable
# can also handle restricted permutations
#
# input: 
# Y - community matrix
# env - environmental data (type: dataframe, one column per variable, data type per column either numeric or factor, if factor no empty levels)
# cov - is there a conditioning variable that should always be considered? (logical) 
# cond - conditioning variable to always be applied
# CTRL - permutation control (generated by how(), strata must have equal number of observations, no emptly levels in factor)
# pure - should pure parts of the environemntal factors be calculated while conditioning for all other factors (default),
#        or: if (pure == "FALSE") only the total parts (no conditioning)
#
# output:
# table with columns: source of variation (parameter), adj.R², F, P-value, covariance, residual variation
# first rows: pure parts, last row: complete model
# or: if (pure=="FALSE") the source of variation (parameter), adj.R², F, P-value of each environemtnal factor individually 
# 
# dependencies:
# requires R packages vegan and permute
#=============================================================================
#documentation end

permuRDA <- function(Y, env, cond, 
                     cov = FALSE, CTRL = how(within = Within(type = "free"), plots = Plots(type = "free"), nperm=999),
                     pure = TRUE) {
  if (pure == TRUE) {
    result <- data.frame(matrix(NA, ncol(env) + 1, 8))
    colnames(result) <- c("parameter", "Radj", "F", "df-N", "df-D", "p", "cov", "Res")
    envmm <- model.matrix(~., data = as.data.frame(env))[, -1]
    
    if (cov == FALSE) {
      R <- rda(Y, envmm)
      result[nrow(result), 2] <- round(RsquareAdj(R)$adj.r.squared, 3)
      result[nrow(result), 8] <- round(1 - result[nrow(result), 2], 3)
    } else {
      condmm <- model.matrix(~., data = as.data.frame(cond))[, -1]
      R <- rda(Y, envmm, condmm)
      result[nrow(result), 2] <- round(RsquareAdj(R)$adj.r.squared, 3)
      V <- varpart(Y, envmm, condmm)
      result[nrow(result), 7:8] <- round(V$part$indfract$Adj.R.squared[c(2, 4)], 3)
    }
    
    A <- anova(R, permutations = CTRL)
    result[nrow(result), 3:6] <- round(c(A$F[1],
                                         A$Df[1:2],
                                         A$"Pr(>F)"[1]), 3) 
    result[nrow(result), 1] <- c("all")
        
    for (j in 1:ncol(env)) {
      X <- env[, j]
      result[j, 1] <- colnames(env)[j]
      
      if (cov == FALSE) {
        Z <- env[, -j]
      } else {
        Z <- data.frame(env[, -j], cond)
      }
      
      X1 <- model.matrix(~., data = as.data.frame(X))[, -1]
      Z1 <- model.matrix(~., data = as.data.frame(Z))[, -1]
      R <- rda(Y, X1, Z1)
      
      A <- anova(R, permutations = CTRL)
      V <- varpart(Y, X1, Z1)
      
      result[j, 2:8] <- round(c(V$part$indfract$Adj.R.squared[1],
                                A$F[1],
                                A$Df[1:2],
                                A$"Pr(>F)"[1],
                                V$part$indfract$Adj.R.squared[2],
                                V$part$indfract$Adj.R.squared[4]), 3)
    }
    return(result)
  }
  
  if (pure == "FALSE"){
    result <- data.frame(matrix(NA, ncol(env) + 1, 6))
    colnames(result) <- c("parameter", "Radj", "F", "df-N", "df-D", "p")
    envmm <- model.matrix(~., data = as.data.frame(env))[, -1]
    
    R <- rda(Y, envmm)
    result[nrow(result), 2] <- round(RsquareAdj(R)$adj.r.squared, 3)
        
    A <- anova(R, permutations = CTRL)
    result[nrow(result), 3:6] <- round(c(A$F[1],
                                         A$Df[1:2],
                                         A$"Pr(>F)"[1]), 3) 
    result[nrow(result), 1] <- c("all")
    
    for (j in 1:ncol(env)){
      X <- env[, j]
      result[j, 1] <- colnames(env)[j]
           
      X1 <- model.matrix(~., data = as.data.frame(X))[, -1]
      R <- rda(Y, X1)
      
      A <- anova(R, permutations = CTRL)      
      result[j, 2:6] <- round(c(RsquareAdj(R)$adj.r.squared,
                                A$F[1],
                                A$Df[1:2],
                                A$"Pr(>F)"[1]),3)
    }
    return(result)
  }
}