-
Notifications
You must be signed in to change notification settings - Fork 34
Open
Description
I noticed that mia::getRDA is internally using vegan::anova.cca for significance testing (this is also routinely used for significance testing with vegan dbrda, although as a separate step).
This is same than PERMANOVA test as far as I can see.
The getRDA and getPERMANOVA are therefore using the same test. The results also seem to be the same.
Do we need both, or could we just keep getRDA (and remove getPERMANOVA) and then document that users can obtain the PERMANOVA p-values from getRDA?
Having both might be convenient to have both just because both are commonly used but this adds maintenance burden and does not add anything.
# load library
library(mia)
library(vegan)
library(DT)
# load some tse object
data("GlobalPatterns")
tse <- GlobalPatterns
tse <- transformAssay(tse, method = "relabundance")
# perform permanova
perm_mia <- getPERMANOVA(
tse,
assay.type = "relabundance",
method = "bray",
formula = x ~ SampleType,
permutations = 99
)
perm_mia_df <- as.data.frame(perm_mia)
datatable(perm_mia_df,
options = list(pageLength = 10),
caption = paste("Results from mia::getPERMANOVA"))
## Using mia::getRDA
# Calculate dbRDA
rda_res <- getRDA(
tse, assay.type = "relabundance", method = "bray",
formula = x ~ SampleType, permutations = 99)
# Significance results are similar to PERMANOVA
rda_mia <- attr(rda_res, "significance")
rda_mia_df <- as.data.frame(rda_mia)
datatable(rda_mia_df,
options = list(pageLength = 10),
caption = paste("Results from mia::getRDA"))
## Using vegan; dbrda+anova.cca
response_matrix <- t(assay(tse, "relabundance"))
predictor_matrix <- colData(tse)[, "SampleType", drop = FALSE]
#predictor_matrix["SampleType"] <- lapply(predictor_matrix["SampleType"], as.factor)
dbrda_vegan <- dbrda(response_matrix ~ ., data = predictor_matrix,
distance = "bray", na.action = na.omit)
# Perform permutational analysis
perm_vegan <- anova.cca(dbrda_vegan, by = "margin",
permutations = 99)
perm_vegan_df <- as.data.frame(perm_vegan)
datatable(perm_vegan_df,
options = list(pageLength = 10),
caption = paste("Results from vegan; dbrda+anova.cca"))
Reactions are currently unavailable
Metadata
Metadata
Assignees
Labels
No labels