Skip to contents

Test beta diversity for associations with metadata.

Source: R/stats_table.r
bdiv_stats.Rd

A convenience wrapper for bdiv_table() + stats_table().

Usage

bdiv_stats(
  biom,
  regr = NULL,
  stat.by = NULL,
  bdiv = "Bray-Curtis",
  weighted = TRUE,
  tree = NULL,
  within = NULL,
  between = NULL,
  split.by = NULL,
  transform = "none",
  test = "emmeans",
  fit = "gam",
  at = NULL,
  level = 0.95,
  alt = "!=",
  mu = 0,
  p.adj = "fdr"
)

Arguments

biom

An rbiom object, such as from as_rbiom(). Any value accepted by as_rbiom() can also be given here.

regr

Dataset field with the x-axis (independent; predictive) values. Must be numeric. Default: NULL

stat.by

Dataset field with the statistical groups. Must be categorical. Default: NULL

bdiv

Beta diversity distance algorithm(s) to use. Options are: "Bray-Curtis", "Manhattan", "Euclidean", "Jaccard", and "UniFrac". For "UniFrac", a phylogenetic tree must be present in biom or explicitly provided via tree=. Default: "Bray-Curtis"

Multiple/abbreviated values allowed.

weighted

Take relative abundances into account. When weighted=FALSE, only presence/absence is considered. Default: TRUE

Multiple values allowed.

tree

A phylo object representing the phylogenetic relationships of the taxa in biom. Only required when computing UniFrac distances. Default: biom$tree

within, between

Dataset field(s) for intra- or inter- sample comparisons. Alternatively, dataset field names given elsewhere can be prefixed with '==' or '!=' to assign them to within or between, respectively. Default: NULL

split.by

Dataset field(s) that the data should be split by prior to any calculations. Must be categorical. Default: NULL

transform

Transformation to apply. Options are: c("none", "rank", "log", "log1p", "sqrt", "percent"). "rank" is useful for correcting for non-normally distributions before applying regression statistics. Default: "none"

test

Method for computing p-values: 'wilcox', 'kruskal', 'emmeans', or 'emtrends'. Default: 'emmeans'

fit

How to fit the trendline. 'lm', 'log', or 'gam'. Default: 'gam'

at

Position(s) along the x-axis where the means or slopes should be evaluated. Default: NULL, which samples 100 evenly spaced positions and selects the position where the p-value is most significant.

level

The confidence level for calculating a confidence interval. Default: 0.95

alt

Alternative hypothesis direction. Options are '!=' (two-sided; not equal to mu), '<' (less than mu), or '>' (greater than mu). Default: '!='

mu

Reference value to test against. Default: 0

p.adj

Method to use for multiple comparisons adjustment of p-values. Run p.adjust.methods for a list of available options. Default: "fdr"

Value

A tibble data.frame with fields from the table below. This tibble object provides the $code operator to print the R code used to generate the statistics.

FieldDescription
.meanEstimated marginal mean. See emmeans::emmeans().
.mean.diffDifference in means.
.slopeTrendline slope. See emmeans::emtrends().
.slope.diffDifference in slopes.
.h1Alternate hypothesis.
.p.valProbability that null hypothesis is correct.
.adj.p.p.val after adjusting for multiple comparisons.
.effect.sizeEffect size. See emmeans::eff_size().
.lowerConfidence interval lower bound.
.upperConfidence interval upper bound.
.seStandard error.
.nNumber of samples.
.dfDegrees of freedom.
.statWilcoxon or Kruskal-Wallis rank sum statistic.
.t.ratio.mean / .se
.r.sqrPercent of variation explained by the model.
.adj.r.r.sqr, taking degrees of freedom into account.
.aicAkaike Information Criterion (predictive models).
.bicBayesian Information Criterion (descriptive models).
.loglikLog-likelihood goodness-of-fit score.
.fit.pP-value for observing this fit by chance.

See also

Other beta_diversity: bdiv_boxplot(), bdiv_clusters(), bdiv_corrplot(), bdiv_heatmap(), bdiv_ord_plot(), bdiv_ord_table(), bdiv_table(), distmat_stats()

Other stats_tables: adiv_stats(), distmat_stats(), stats_table(), taxa_stats()

Examples

    library(rbiom)
    
    biom <- rarefy(hmp50)
      
    bdiv_stats(biom, stat.by = "Sex", bdiv = c("bray", "unifrac"))[,1:7]
#> # Model:    gam(.distance ~ Sex, method = "REML")
#> # A tibble: 6 × 7
#>   .bdiv       Sex                   .mean.diff .h1    .p.val .adj.p .effect.size
#>   <fct>       <chr>                      <dbl> <fct>   <dbl>  <dbl>        <dbl>
#> 1 Bray-Curtis Female - Male            0.0562  != 0  0.00693 0.0337       0.245 
#> 2 Bray-Curtis Male - Female vs Male   -0.0462  != 0  0.0112  0.0337      -0.222 
#> 3 UniFrac     Female - Female vs M…   -0.0247  != 0  0.0784  0.157       -0.112 
#> 4 UniFrac     Male - Female vs Male   -0.0205  != 0  0.249   0.373       -0.101 
#> 5 Bray-Curtis Female - Female vs M…    0.00998 != 0  0.473   0.568        0.0457
#> 6 UniFrac     Female - Male           -0.00421 != 0  0.839   0.839       -0.0183
    
    bdiv_stats(biom, stat.by = "Body Site", split.by = "==Sex")[,1:6]
#> # Model:    gam(.distance ~ `Body Site`, method = "REML")
#> # A tibble: 150 × 6
#>    Sex    `Body Site`                         .mean.diff .h1     .p.val   .adj.p
#>    <fct>  <chr>                                    <dbl> <fct>    <dbl>    <dbl>
#>  1 Female Buccal mucosa - Mid vagina vs Stool     -0.793 != 0  1.48e-62 6.68e-61
#>  2 Female Buccal mucosa - Buccal mucosa vs M…     -0.793 != 0  1.62e-62 6.68e-61
#>  3 Female Buccal mucosa - Anterior nares vs …     -0.792 != 0  1.65e-62 6.68e-61
#>  4 Female Buccal mucosa - Mid vagina vs Sali…     -0.792 != 0  1.78e-62 6.68e-61
#>  5 Female Saliva - Mid vagina vs Stool            -0.517 != 0  4.81e-38 1.04e-36
#>  6 Female Saliva - Buccal mucosa vs Mid vagi…     -0.517 != 0  5.15e-38 1.04e-36
#>  7 Female Saliva - Anterior nares vs Mid vag…     -0.516 != 0  5.44e-38 1.04e-36
#>  8 Female Saliva - Mid vagina vs Saliva           -0.516 != 0  5.54e-38 1.04e-36
#>  9 Female Buccal mucosa vs Saliva - Mid vagi…     -0.198 != 0  4.09e-36 6.82e-35
#> 10 Female Buccal mucosa vs Mid vagina - Bucc…      0.197 != 0  4.94e-36 7.41e-35
#> # ℹ 140 more rows