Prepare data and environment

data("WRS_abiotic_model")
list2env(WRS_abiotic_model, envir = .GlobalEnv)
data("WRS_OTU_model")
list2env(WRS_OTU_model, envir = .GlobalEnv)

# Choose PCoA axes and extract eigenvectors
# Step 1: make Gower distances
WRS_abiotic_gower <- lapply(WRS_abiotic_model, daisy, metric = "gower")
# Step 2: make ordinations
WRS_PCoA <- lapply(WRS_abiotic_gower, wcmdscale, eig = TRUE, x.ret = TRUE)
# Determine number of useful axes using broken stick
taxa_screeplot <- function(PCoA_list) {
  for(i in seq_along(PCoA_list)){ # assumes all elements have same length as first element. I should abstract this later probably.
    # Create broken-stick model
    eig <- PCoA_list[[i]]$eig
    n <- length(eig)
    bsm <- data.frame(j = seq(1:n), p = 0)
    bsm$p[1] <- 1 / n
    for(x in 2 : n) {
      bsm$p[x] <- bsm$p[x - 1] + (1 / (n + 1 - x))
    }
    bsm$p <- 100 * bsm$p / n
    # Initialise plot
    barplot(t(cbind(100 * eig / sum(eig), bsm$p[n:1])),
            beside = TRUE, col = c("bisque", 2), las = 2)
    # Plot legend
    legend("topright", c("% eigenvalue", "Broken stick model"),
           pch = 15, col = c("bisque", 2), bty = "n", title = names(PCoA_list)[i])
  }
}
taxa_screeplot(WRS_PCoA)
# I'll take the first three axes of each class
PCoA_axes <- lapply(WRS_PCoA, "[[", "points")
PCoA_axes <- lapply(seq_along(PCoA_axes), function(i) {
  PCoA_axes[[i]] <- as.data.frame(PCoA_axes[[i]][, 1:3])
})
names(PCoA_axes) <- names(WRS_PCoA)
PCoA_axes <- data.frame(unlist(PCoA_axes, recursive = FALSE))

Create dbRDAs

dbRDAs <- lapply(seq_along(WRS_OTU_model), function(i) {
  capscale(WRS_OTU_model[[i]] ~ ., data = PCoA_axes, distance = "bray")
})
names(dbRDAs) <- names(WRS_OTU_model)
dbRDA_summaries <- lapply(dbRDAs, summary)
lapply(dbRDAs, anova)

lapply(seq_along(dbRDAs), function(i) {
  plot(dbRDAs[[i]], type = "points", display = c("wa", "cn"), main = names(dbRDAs)[i])
  text(dbRDAs[[i]], cex = 1.5)
})

inertias <- sapply(dbRDA_summaries, "[", c("constr.chi", "tot.chi"))
as.numeric(inertias[1, ]) / as.numeric(inertias[2, ])
# Proportion is approximately 0.32 - 0.36 for all taxa, although this is an over-estimate.


mixtrak/wellington.aquifer.ecol documentation built on Nov. 30, 2017, 4:25 a.m.