WdStar is an R package for
- Are you using PERMANOVA and concerned about how heteroscedasticity and unbalanced sample sizes may affect your analyses?
- Have you ever questioned the reliability of community-wide microbiome analysis results?
- Are you looking for a global test for your multivariate dataset?
$W_d^*$ addresses these concerns and is
- robust to heteroscedasticity;
- handles multi-level factors and stratification;
- allows for multiple post hoc testing scenarios;
- allows for adjustment of covariates; and
- can be used with any distance or dissimilarity matrix.
Preprint: Covariate-adjusted multivariate analysis for omics data
- Hamidi B, Fanning L, Wallace K, & Alekseyenko AV. Submitted to Briefings in Bioinformatics on June 23, 2026. Not currently accepted or under review. DOI forthcoming.
- Code repository
See the goodness-of-fit and covariate-adjustment methods inventory for implementation notes and comparisons with vegan, fast.adonis, DISTLM, MultANOVA, and microbiome-analysis packages.
Source installation of the WdStar R package is available directly from GitHub using remotes for R 3.6 or later:
install.packages("remotes")
remotes::install_github("alekseyenko/WdStar", force = TRUE)
library(WdStar)
packageVersion("WdStar")Until a CRAN release is available, install WdStar from GitHub as shown above.
For detailed and complex examples please refer to our publication repositories, which contain Markdown files with application datasets and code.
The following is a simple example using the mtcars dataset to assess the effect of gear on mpg, cyl, and disp (first three variables of the dataset):
# Load dataset
data(mtcars)
# The outcome could be a single variable or multiple variables (such as multidimensional omics data).
### This is an example of outcome with a single variable (`mpg`):
dm <- dist(mtcars$mpg, method="euclidean")
### This is an example of outcome with multiple variables (`mpg`, `cyl`, and `disp`):
dm <- dist(mtcars[1:3], method="euclidean")
# Grouping/independent variable. You could use multiple variables here too.
f <- factor(mtcars$gear)
# Basic multivariate test example ###########
#############################################
WdS.test(dm=dm, f=f)
## By default, the unadjusted test's goodness.of.fit reports the
## distance-based pseudo-R-squared for the tested factor.
unadjusted_res <- WdS.test(dm=dm, f=f)
unadjusted_res$goodness.of.fit
## Use goodness="none" to skip goodness-of-fit calculations.
WdS.test(dm=dm, f=f, goodness="none")
# Stratified example ########################
#############################################
strata <- factor(mtcars$vs)
WdS.test(dm=dm, f=f, strata=strata)
# Covariate adjustment/elimination examples #
#############################################
## Right-hand side adjustment formula to specify adjustment covariates.
formula <- ~ wt + as.factor(am)
## Adjustment example 1: pass unadjusted `dm` and formula to WdS.test()
WdS.test(dm=dm, f=f, formula=formula, formula_data=mtcars) ## Perform adjusted test
## By default, goodness.of.fit reports the distance-based semi-partial
## pseudo-R-squared for the tested factor after adjustment.
## Additional components computed along the way are stored but not printed.
res <- WdS.test(dm=dm, f=f, formula=formula, formula_data=mtcars)
res$goodness.components
## Eigenvalue and tolerance diagnostics for residual distance matrices are
## stored separately from goodness-of-fit values.
res$distance.diagnostics
## Request all available goodness-of-fit components, including factor-only,
## adjustment-only, full-model, semi-partial, and partial pseudo-R-squared.
WdS.test(dm=dm, f=f, formula=formula, formula_data=mtcars, goodness="all")
## Interpretation of adjusted components:
## - adjustment: variation explained by the adjustment variables alone.
## - full: variation explained jointly by adjustment variables and the tested factor.
## - semi.partial: additional total variation explained by the tested factor after adjustment.
## - partial: adjustment-residual variation explained by the tested factor.
## Negative values can occur if the residual distances contain more variation
## than the original distance matrix.
## Adjustment example 2: Create the adjusted distance matrix `a.dm` outside the function
a.dm <- a.dist(dm=dm, formula=formula, formula_data=mtcars)
WdS.test(dm=a.dm, f=f) ## Perform adjusted test with `a.dm`; input diagnostics are preserved.
attr(a.dm, "distance.diagnostics")
## Store raw eigenvalues only when deeper diagnostics are needed.
a.dm.with.eigenvalues <- a.dist(
dm=dm,
formula=formula,
formula_data=mtcars,
keep.eigenvalues=TRUE
)
length(attr(a.dm.with.eigenvalues, "distance.diagnostics")$eigenvalues[[1]])
## Distance-based pseudo-R-squared can also be computed directly
dist.goodness.of.fit(dm=dm, dm_residual=a.dm)
## Taxa/ASV importance can be ranked by adjusting for one taxon at a time.
## The abundance values are used exactly as supplied.
taxa_table <- data.frame(
ASV1 = mtcars$wt,
ASV2 = mtcars$hp,
ASV3 = mtcars$qsec,
row.names = rownames(mtcars)
)
taxonomy_table <- data.frame(
Genus = c("TaxonA", "TaxonB", "TaxonC"),
Species = c("species1", "species2", "species3"),
row.names = colnames(taxa_table)
)
WdS.taxa.importance(
dm=dm,
f=f,
taxa_table=taxa_table,
taxa_are_rows=FALSE,
taxonomy_table=taxonomy_table,
nrep=9
)
## Add sample-level terms to every taxon-specific adjustment model.
## For example, this evaluates each taxon together with a paired subject term.
sample_data <- data.frame(
Subject_ID = factor(rep(seq_len(16), each = 2)),
row.names = rownames(mtcars)
)
WdS.taxa.importance(
dm=dm,
f=f,
taxa_table=taxa_table,
taxa_are_rows=FALSE,
taxonomy_table=taxonomy_table,
formula=~ Subject_ID,
formula_data=sample_data,
rank.by="adjustment.goodness.of.fit",
nrep=9
)Further examples are provided in the package documentation and may be accessed by running the following commands:
?WdS.test
?a.dist
?dist.goodness.of.fit
?WdS.taxa.importanceWe welcome feature requests and bug reports and kindly ask you to submit them via our GitHub issue tracker.