Takes the output of dist_calc function. Or use with the result of the permanova function to ensure the results correspond to exactly the same input data. Runs betadisper for all categorical variables in variables argument. See help('betadisper', package = 'vegan').

dist_bdisp(
  data,
  variables,
  method = c("centroid", "median")[[1]],
  complete_cases = TRUE,
  verbose = TRUE
)

Arguments

data

psExtra output from dist_calc

variables

list of variables to use as group

method

centroid or median

complete_cases

drop samples with NAs in any of the variables listed

verbose

sends messages about progress if true

Value

psExtra containing betadisper results

Examples

library(phyloseq)
library(vegan)
#> Loading required package: permute
#> Loading required package: lattice
#> This is vegan 2.6-4
data("dietswap", package = "microbiome")

# add some missings to demonstrate automated removal
sample_data(dietswap)$sex[3:6] <- NA
# create a numeric variable to show it will be skipped with a warning
dietswap <- ps_mutate(dietswap, timepoint = as.numeric(timepoint))

# straight to the betadisp
bd1 <- dietswap %>%
  tax_agg("Genus") %>%
  dist_calc("aitchison") %>%
  dist_bdisp(variables = c("sex", "bmi_group", "timepoint")) %>%
  bdisp_get()
#> Dropping samples with missings: 4
#> Warning: Variable 'timepoint' is skipped as it cannot be used for grouping (class = 'numeric')
bd1$sex
#> $model
#> 
#> 	Homogeneity of multivariate dispersions
#> 
#> Call: vegan::betadisper(d = distMat, group = meta[[V]], type = method)
#> 
#> No. of Positive Eigenvalues: 123
#> No. of Negative Eigenvalues: 0
#> 
#> Average distance to centroid:
#> female   male 
#> 10.022  9.059 
#> 
#> Eigenvalues for PCoA axes:
#> (Showing 8 of 123 eigenvalues)
#>  PCoA1  PCoA2  PCoA3  PCoA4  PCoA5  PCoA6  PCoA7  PCoA8 
#> 5106.7 2075.2 1360.1 1281.4 1061.7  890.1  665.4  567.5 
#> 
#> $anova
#> Analysis of Variance Table
#> 
#> Response: Distances
#>            Df Sum Sq Mean Sq F value    Pr(>F)    
#> Groups      1  50.14  50.144  20.976 7.853e-06 ***
#> Residuals 216 516.35   2.391                      
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> $tukeyHSD
#>   Tukey multiple comparisons of means
#>     95% family-wise confidence level
#> 
#> Fit: aov(formula = distances ~ group, data = df)
#> 
#> $group
#>                   diff       lwr        upr   p adj
#> male-female -0.9632687 -1.377813 -0.5487242 7.9e-06
#> 
#> 
# quick vegan plotting methods
plot(bd1$sex$model, label.cex = 0.5)

boxplot(bd1$sex$model)


# compute distance and use for both permanova and dist_bdisp
testDist <- dietswap %>%
  tax_agg("Genus") %>%
  dist_calc("bray")

PERM <- testDist %>%
  dist_permanova(
    variables = c("sex", "bmi_group"),
    n_processes = 1, n_perms = 99
  )
#> Dropping samples with missings: 4
#> 2024-04-03 12:20:16.206291 - Starting PERMANOVA with 99 perms with 1 processes
#> 2024-04-03 12:20:16.358622 - Finished PERMANOVA
str(PERM, max.level = 1)
#> Formal class 'psExtra' [package "microViz"] with 15 slots

bd <- PERM %>% dist_bdisp(variables = c("sex", "bmi_group"))
bd
#> psExtra object - a phyloseq object with extra slots:
#> 
#> phyloseq-class experiment-level object
#> otu_table()   OTU Table:         [ 130 taxa and 218 samples ]
#> sample_data() Sample Data:       [ 218 samples by 8 sample variables ]
#> tax_table()   Taxonomy Table:    [ 130 taxa by 3 taxonomic ranks ]
#> 
#> psExtra info:
#> tax_agg = "Genus" 
#> 
#> bray distance matrix of size 218 
#> 0.7639533 0.731024 0.7283254 0.6637252 0.7437293 ...
#> 
#> permanova:
#> Permutation test for adonis under reduced model
#> Marginal effects of terms
#> Permutation: free
#> Number of permutations: 99
#> 
#> vegan::adonis2(formula = formula, data = metadata, permutations = n_perms, by = by, parallel = parall)
#>            Df SumOfSqs      R2      F Pr(>F)   
#> sex         1    0.361 0.00933 2.1539   0.17   
#> bmi_group   2    2.377 0.06143 7.0888   0.01 **
#> Residual  214   35.874 0.92720                 
#> Total     217   38.691 1.00000                 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> betadisper:
#> sex bmi_group