Prepare environment and data

data("WRS_OTU_model")
data("WRS_abiotic_model")
# Obtain prior groupings based on environmental data
abiotic_data <- data.frame(unname(WRS_abiotic_model))
WRS_abiotic_gower <- daisy(abiotic_data, metric = "gower")
abiotic_clust <- hclust(WRS_abiotic_gower, method = "average")
abiotic_cut <- cutree(abiotic_clust, k = 4)
abiotic_groups <- letters[abiotic_cut]
names(abiotic_groups) <- names(abiotic_cut)
# Create community distances & turn into PCoA axis values & Euclidean distances
WRS_bray <- lapply(WRS_OTU_model, vegdist)
WRS_PCoA <- lapply(WRS_bray, wcmdscale, eig = TRUE, x.ret = TRUE)
WRS_points <- lapply(WRS_PCoA, "[[", "points")
WRS_euclid <- lapply(WRS_points, dist)

Linear discriminant analysis (LDA)

# Test for multivariate homogeneity of groups
disper <- lapply(WRS_euclid, betadisper, group = abiotic_groups)
lapply(disper, anova)
lapply(disper, permutest)
# Shrug

WRS_LDA <- lapply(seq_along(WRS_points), function(i) {
  lda(abiotic_groups ~ WRS_points[[i]][, 1:4], CV = TRUE)
})
names(WRS_LDA) <- names(WRS_points)

LDA_class <- structure(.Data = lapply(WRS_LDA, "[[", "class"), class = "list", names = names(WRS_points))
# Convert NA to new class
LDA_class <- lapply(LDA_class, addNA)

WRS_table <- lapply(seq_along(LDA_class), function(i) {
  table(abiotic_groups, LDA_class[[i]])
})
names(WRS_table) <- names(WRS_points)

WRS_LDA_predict <- sapply(seq_along(WRS_table), FUN = function(i) {
  diag(prop.table(WRS_table[[i]], 1))
  })
colnames(WRS_LDA_predict) <- names(WRS_table)
WRS_LDA_predict

LDA_posterior <- structure(.Data = lapply(WRS_LDA, "[[", "posterior"), class = "list", names = names(WRS_points))

Try a multinomial logistic regression instead

archaea_logit <- multinom(abiotic_groups ~ WRS_points$archaea[, 1:4])


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