| Title: | Linear Models for Sequence Count Data |
|---|---|
| Description: | Provides scalable generalized linear and mixed effects models tailored for sequence count data analysis (e.g., analysis of 16S or RNA-seq data). Uses Dirichlet-multinomial sampling to quantify uncertainty in relative abundance or relative expression conditioned on observed count data. Implements scale models as a generalization of normalizations which account for uncertainty in scale (e.g., total abundances) as described in Nixon et al. (2025) <doi:10.1186/s13059-025-03609-3> and McGovern et al. (2025) <doi:10.1101/2025.08.05.668734>. |
| Authors: | Justin Silverman [aut, cre], Greg Gloor [aut], Kyle McGovern [aut, ctb] |
| Maintainer: | Justin Silverman <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 1.3.0 |
| Built: | 2026-07-24 09:11:06 UTC |
| Source: | https://github.com/jsilve24/aldex3 |
ALDEx3 Linear Models
aldex( Y, X, data = NULL, method = "lm", nsample = 2000, scale = NULL, streamsize = 8000, n.cores = detectCores() - 1, return.pars = c("X", "estimate", "std.error", "p.val", "p.val.adj", "logComp", "logScale"), p.adjust.method = "BH", test = "t.HC3", onesided = FALSE, ... )aldex( Y, X, data = NULL, method = "lm", nsample = 2000, scale = NULL, streamsize = 8000, n.cores = detectCores() - 1, return.pars = c("X", "estimate", "std.error", "p.val", "p.val.adj", "logComp", "logScale"), p.adjust.method = "BH", test = "t.HC3", onesided = FALSE, ... )
Y |
an (D x N) matrix of sequence count, N is number of samples, D is number of taxa or genes |
X |
either a formula (in which case DATA must be non-null) or a model matrix of dimension P x N (P is number of linear model covariates). If a formula is passed it should not include the target Y, e.g., should simply be "~condition-1" (note the lack of the left hand side). If using lme4, this should be the formula including random effects. |
data |
a data frame for use with formula, must have N rows |
method |
(default lm) The regression method; "lm": linear regression,
"lme4": linear mixed effects regression with lme4; "nlme": linear mixed
effects models with nlme, REQUIRES "random" argument representing random
effects be passed into |
nsample |
number of monte carlo replicates |
scale |
A scale-model function, such as |
streamsize |
(default 8000) memory footprint (approximate) at which to use streaming. This should be thought of as the number of Mb for each streaming chunk. If DNnsample*8/1000000 is less than streamsize then no streaming will be performed. Note, to conserve memory, samples from the Dirichlet and scale models will not be returned if streaming is used. Streaming can be turned off by setting streamsize=Inf. |
n.cores |
(default detectCores()-1) If method is 'lme4' or 'blmm',
use this many cores for feature-level parallelism across taxa/features.
For |
return.pars |
what results should be returned, see return section below. |
p.adjust.method |
(default BH) The method for multiple hypothesis test
correction. See |
test |
(default t.HC3), "t", t test is performed for each covariate (fast); "t.HC0" Heteroskedasticsity-Robust Standard Errors used (HC0; White's; slower); "t.HC3" (default) Heteroskedasticsity-Robust Standard Errors used (HC3; unlike HC0, this includes a leverage adjustment and is better for small sample sizes or when there are data with high leverage; slowest). To learn more about these, loko at Long and Ervin (2000) Using Heteroscedasticity Consistent Standard Errors in the Linear Regression Model, The American Statistician. |
onesided |
(default: FALSE) if sided return p-values for two-sided test. Otherwise if "lower" or "upper" return one-sided test corresponding to test that estimate is negative or positive respectively. |
... |
Additional scale-model arguments, such as |
Sequence counts determine the relative composition of each sample, but not
its total abundance or scale. The scale argument tells ALDEx3 what is
known about that missing quantity and how uncertain that knowledge is.
Choose a scale model based on the information available:
Use sample.sm for external measurements such as flow
cytometry, qPCR, or spike-ins for individual samples.
Use coefficient.sm when prior knowledge concerns a
group difference, treatment effect, time trend, or another model term.
Use clr.sm when CLR normalization is a reasonable
center but its implied scale differences are uncertain.
Use tss.sm when the starting assumption is no
systematic scale difference between groups.
Write a custom scale-model function when the built-in models do not
match the information available for the study. The "Writing a custom
scale model" example below describes the required function interface
and N x nsample output.
a list with elements controled by parameter resturn.pars. Options
include: - X: P x N covariate matrix - estimate: (P x D x nsample) array
of linear model estimates - std.error: (P x D x nsample) array of standard
error for the estimates - p.val: (P x D) matrix, unadjusted p-value for
two-sided t-test - p.val.adj: (P x D) matrix, p-value for two-sided t-test
adjusted according to p.adj.method - logScale: (N x S) matrix, samples
of the log scale from the scale model - logComp: (D x N x S) array,
samples of the log composition from the multinomial-Dirichlet - streaming:
boolean, detnote if streaming was used. - random.effects (Pr x N x S): if
using mixed effects models, return all Pr random-effects covariance terms
and the residual variance. For method = "blmm", correlated
random-effect covariance terms are returned when present, matching the
lme4 naming convention as closely as possible. Note, logScale and
logComp are not returned if streaming is active.
Justin Silverman, Kyle McGovern
Y <- matrix(1:110, 10, 11) condition <- c(rep(0, 5), rep(1, 6)) data <- data.frame(condition=condition) ## Use CLR as the center and allow coefficient-scale uncertainty of SD 0.5. res <- aldex(Y, ~condition, data, nsample=2000, scale=clr.sm, gamma=0.5) ## Writing a custom scale model ## A custom model must return one log2-scale draw per sample and Monte ## Carlo replicate. Here, external_scale supplies a known value per sample. known_scale <- function(X, logComp, external_scale) { nsample <- dim(logComp)[3] matrix(rep(external_scale, nsample), nrow=ncol(X), ncol=nsample) } res <- aldex(Y, ~condition, data, nsample=2000, scale=known_scale, external_scale=rep(0, ncol(Y)))Y <- matrix(1:110, 10, 11) condition <- c(rep(0, 5), rep(1, 6)) data <- data.frame(condition=condition) ## Use CLR as the center and allow coefficient-scale uncertainty of SD 0.5. res <- aldex(Y, ~condition, data, nsample=2000, scale=clr.sm, gamma=0.5) ## Writing a custom scale model ## A custom model must return one log2-scale draw per sample and Monte ## Carlo replicate. Here, external_scale supplies a known value per sample. known_scale <- function(X, logComp, external_scale) { nsample <- dim(logComp)[3] matrix(rep(external_scale, nsample), nrow=ncol(X), ncol=nsample) } res <- aldex(Y, ~condition, data, nsample=2000, scale=known_scale, external_scale=rep(0, ncol(Y)))
ALDEx2-inspired effect diagnostics for an ALDEx3 binary contrast
aldex.effect(object, contrast)aldex.effect(object, contrast)
object |
An object returned by |
contrast |
The exact model coefficient or unambiguous binary data column to summarize. |
aldex.effect is a helper for effect-size diagnostics that Greg Gloor and
ALDEx2 users commonly inspect when evaluating pairwise effects. It is not a
replacement for summary.aldex, which reports model estimates, standard
errors and adjusted p-values. This function summarizes a single binary
ALDEx3 model contrast using Cohen's-d-based Monte Carlo diagnostics.
For a selected contrast row k, feature d, Monte Carlo draw s, and
binary group indicators x = object$X[k, ], let B[d, s] be the fitted
ALDEx3 coefficient stored in object$estimate[k, d, s]. Let
logW[d, i, s] = object$logComp[d, i, s] + object$logScale[i, s] be the
reconstructed log abundance for sample i. aldex.effect requires
logComp and logScale, so it is unavailable when aldex streams large
jobs and omits those arrays. For each feature and Monte Carlo draw, the
pooled within-group variance is
((n0 - 1) * var(logW[d, x == 0, s]) + (n1 - 1) * var(logW[d, x == 1, s])) / (n0 + n1 - 2), where n0 and n1 are the
group sizes.
The returned columns are calculated as follows:
estimate: the mean of B[d, s] over Monte Carlo draws.
pooled.SD: the mean of the pooled within-group standard deviation,
sqrt(pooled variance), over Monte Carlo draws.
cohens.d: the mean of B[d, s] / pooled.SD[d, s] over Monte Carlo
draws. This is Cohen's-d-based and is not the same as the historical
robust ALDEx2 effect statistic based on between/within dispersion.
overlap: an ALDEx2-style directional uncertainty diagnostic. It counts
the Monte Carlo Cohen's d draws below and above zero, smooths those two
counts with a 0.5 pseudocount, converts them with the ALDEx/Aitchison
mean, and returns the smaller smoothed sign probability. Values near 0
indicate that nearly all Monte Carlo effect draws have the same sign.
Values near 0.5 indicate uncertainty about direction. overlap is not a
p-value, a confidence interval probability, or a literal area of overlap
between two density curves.
A data.frame with columns parameter, entity, estimate,
pooled.SD, cohens.d, and overlap.
Greg Gloor, Justin Silverman
Not designed to create realistic data. Does not add any noise to linear regression! True W is in CLR Corrdinates
aldex.lm.sim.clr(D = 10, N = 11, P = 2, depth = 10000)aldex.lm.sim.clr(D = 10, N = 11, P = 2, depth = 10000)
D |
number of taxa/genes |
N |
number of samples |
P |
number of covariates |
depth |
sum of counts for each multinomial draw |
a list with elements Y, X, W, and Lambda
Justin Silverman
Includes random intercpt, option to include random slope, second random intercept, and time-correlation.
aldex.mem.sim( D, days, subjects, depth = 10000, location = 1, random_slope = FALSE, corr = 0, rho_ar1 = 0, sd_resid = 0.1 )aldex.mem.sim( D, days, subjects, depth = 10000, location = 1, random_slope = FALSE, corr = 0, rho_ar1 = 0, sd_resid = 0.1 )
D |
number of taxa/genes |
days |
num days (i.e., repeated measurements) for each subject |
subjects |
num of subjects to simulate |
depth |
sum of counts for each multinomial draw |
location |
second random intercept, if 0 ignore, else simulate num of locations |
random_slope |
If true, simulate a random slope for each subject. |
corr |
The correlation between slope/random intercept |
rho_ar1 |
The ar1 time-correlation, if 0, don't simulate |
sd_resid |
The residual error time-correlation. |
Kyle McGovern
Simple plots for an ALDEx3 object
aldex.plot( object, plot = c("volcano", "effect", "MA", "water"), contrast = NULL, threshold = 0.05, min.diff = 0.5, cohen = 0.5, sig.col = rgb(1, 0, 0, 0.5), water.show = 5, water.col = c("red", "blue"), water.names = TRUE, ... )aldex.plot( object, plot = c("volcano", "effect", "MA", "water"), contrast = NULL, threshold = 0.05, min.diff = 0.5, cohen = 0.5, sig.col = rgb(1, 0, 0, 0.5), water.show = 5, water.col = c("red", "blue"), water.names = TRUE, ... )
object |
An object of class |
plot |
type of plot (default='volcano') |
contrast |
the name of the comparison to plot, must be provided |
threshold |
FDR significance threshold (default=0.05) |
min.diff |
(default=0.5) used for MA (display), and waterfall (cutoff) plots |
cohen |
Cohen's d threshold for effect plot (slope) |
sig.col |
color for significant features |
water.show |
number features to show in waterfall plot |
water.col |
vector of colors for waterfall plot |
water.names |
logical, show names of features for waterfall plot |
... |
additional plot parameters (cex) |
Provides volcano, effect, MA and waterfall plots from an ALDEx3 object.
This method plots combinations of adjusted p-values from object$p.val.adj,
posterior estimates, standard errors and log abundance values
averaged across Monte Carlo samples.
The result is returned as a single plot of the desired type. Effect plots
use the ALDEx2-inspired diagnostics returned by aldex.effect. The contrast
must identify exactly one binary model coefficient. Effect diagnostics and
MA plots require logComp and logScale, which are not returned when
aldex streams large jobs.
the desired plot
Greg Gloor
Calculate p-values adjusting for changes in sign as described by Nixon et al. (2024) in Beyond Normalization: Incorporating Scale Uncertainty in Microbiome and Gene Expression Analysis (internal only)
aldex.pvals(p.lower, p.upper, p.adjust.method, onesided = FALSE)aldex.pvals(p.lower, p.upper, p.adjust.method, onesided = FALSE)
p.lower |
A P x D x S matrix for P covariates, D taxa/genes, and S monte carlo samples representing the lower tail p.values |
p.upper |
A P x D x S matrix for P covariates, D taxa/genes, and S monte carlo samples representing the upper tail p.values |
p.adjust.method |
An adjutment method for p.adjust |
onesided |
(default: FALSE) if sided return p-values for two-sided test. Otherwise if "lower" or "upper" return one-sided test corresponding to test that estimate is negative or positive respectively. |
A list with P x D matrices with the non-adjusted and adjusted p-values.
Justin Silverman, Kyle McGovern
Nixon G, Gloor GB, Silverman JD (2025). "Incorporating scale uncertainty in microbiome and gene expression analysis as an extension of normalization". Genome Biology. doi:10.1186/s13059-025-03609-3
Fast approximate mixed-effects engine for ALDEx3.
blmm(logW, formula, data, n.cores = 1L)blmm(logW, formula, data, n.cores = 1L)
logW |
numeric array (N x D x S) |
formula |
lme4 mixed-effects formula |
data |
data.frame for formula evaluation |
n.cores |
number of worker processes for feature-level parallelism.
If |
Internal engine called by aldex when method = "blmm".
Uses an lme4-style profiled Gaussian LMM formulation. One batched
anchor REML problem is fit per feature across all S MC draws; draw-specific
local covariance updates are obtained via a single Newton step; exact
conditional GLS fixed-effect solves are then run per draw. The approximation
is only in the variance-component step — fixed effects and standard errors
are exact conditional on the updated covariance parameters.
Features for which the approximate path fails are silently re-fit with the
exact lme4 engine and a consolidated warning is issued naming the
affected features and their error messages.
named list of arrays (P x D x S): estimate, std.error, df,
p.lower, p.upper, and (Pr x D x S) random.eff — same layout as
sr.mem.
Use this model when centered log-ratio (CLR) normalization is a reasonable
starting point, but you do not want to assume that it recovers the sample
scales exactly. For example, in a case-control study, CLR may provide a
plausible center while gamma allows the true difference in total
abundance between cases and controls to depart from that center.
clr.sm(X, logComp, gamma = 0.5)clr.sm(X, logComp, gamma = 0.5)
X |
Model matrix passed automatically by |
logComp |
Monte Carlo log-compositions passed automatically by
|
gamma |
Non-negative scalar. Standard deviation of the Gaussian random
shift that relaxes the CLR assumption about scale. |
Set gamma = 0 to use the CLR-implied scales without additional
uncertainty. A value such as gamma = 0.5 allows each scale-model
coefficient to vary with SD 0.5 on the log2 scale. Larger values express
less confidence in the CLR assumption.
The uncertainty applies to model coefficients, not independently to every
sample. Samples with the same covariate values therefore receive the same
random shift in a Monte Carlo draw. This differs from sample.sm,
which can draw a separate scale for each sample.
A numeric matrix of dimension N x nsample giving Monte Carlo samples
of the log-scale for each sample (rows) and each Monte Carlo draw (columns).
Justin Silverman
Nixon G, Gloor GB, Silverman JD (2025). "Incorporating scale uncertainty in microbiome and gene expression analysis as an extension of normalization". Genome Biology. doi:10.1186/s13059-025-03609-3
# Allow moderate uncertainty around CLR-implied scale differences. Y <- matrix(seq_len(60), nrow = 10) condition <- factor(c("control", "control", "control", "case", "case", "case")) fit <- aldex(Y, ~ condition, data.frame(condition), nsample = 20, scale = clr.sm, gamma = 0.5)# Allow moderate uncertainty around CLR-implied scale differences. Y <- matrix(seq_len(60), nrow = 10) condition <- factor(c("control", "control", "control", "case", "case", "case")) fit <- aldex(Y, ~ condition, data.frame(condition), nsample = 20, scale = clr.sm, gamma = 0.5)
Use this model when prior information describes a treatment, group, time, or other model coefficient rather than an individual sample. Samples with the same covariate values share the same scale shift in each Monte Carlo draw.
coefficient.sm( X, logComp, c.mu = NULL, c.cor = NULL, c.sd = NULL, c.cov = NULL )coefficient.sm( X, logComp, c.mu = NULL, c.cor = NULL, c.sd = NULL, c.cov = NULL )
X |
Model matrix passed automatically by |
logComp |
Monte Carlo log-compositions passed automatically by
|
c.mu |
Expected log2 scale change for each model coefficient, in the
same order as the rows of |
c.cor |
Deprecated |
c.sd |
SD of the log2 scale change for each model coefficient. Must
have length |
c.cov |
A |
Suppose an effect plot from a control-versus-treatment experiment is
centered near -8, even though prior knowledge suggests that most features
should not change. One way to represent that knowledge is to shift the
treatment samples upward by 8 log2 units relative to the controls. With
control as the reference level, use c.mu = c(0, 8): the first value
leaves the common intercept unchanged, and the second adds 8 to treatment.
Setting c.sd = c(0, 0.5) expresses an SD of 0.5 around that treatment
shift. If the plot were centered near +8, use -8 instead. Adding 8 to both
groups would not recenter the plot because it would leave their difference
unchanged.
More generally, c.mu gives the expected log2 scale change for each
model coefficient and c.sd gives its SD. Use c.cov when the
coefficient uncertainties are correlated. For each Monte Carlo draw,
ALDEx3 samples a coefficient vector and applies it to the model matrix
X; mathematically, and the sample
scales are .
Exactly one of c.sd or c.cov must be provided. The deprecated
c.cor argument may be used in place of c.cov during the
compatibility period.
A numeric matrix of dimension N x nsample giving Monte Carlo
draws of the log2 scale for each sample (rows) across nsample draws
(columns).
Kyle McGovern
# Three control samples followed by three treatment samples. condition <- factor( c("control", "control", "control", "treatment", "treatment", "treatment"), levels = c("control", "treatment") ) X <- t(model.matrix(~ condition)) logComp <- array(0, c(2, 6, 10)) logScale <- coefficient.sm( # Shift treatment upward by 8, with SD 0.5 around that shift. X, logComp, c.mu = c(0, 8), c.sd = c(0, 0.5) )# Three control samples followed by three treatment samples. condition <- factor( c("control", "control", "control", "treatment", "treatment", "treatment"), levels = c("control", "treatment") ) X <- t(model.matrix(~ condition)) logComp <- array(0, c(2, 6, 10)) logScale <- coefficient.sm( # Shift treatment upward by 8, with SD 0.5 around that shift. X, logComp, c.mu = c(0, 8), c.sd = c(0, 0.5) )
Function to compute cohensd on the results provided by the aldex function
cohensd(m, var)cohensd(m, var)
m |
the output of a call to |
var |
if aldex was called with X being a pre-computed model matrix,
this var should be an integer corresponding to a binary covariate
indicating the desired effect size to calculate (an effect size between
two groups indicated by the binary covariate). For example, if the third
covariate in the model is an indicator denoting health (0) and disease (1)
then set |
WARNING: this function is experimental and requires users read the documentation fully.
A (D x nsample)-matrix of Cohen's d statistics for the variable of interest.
Justin Silverman
Takes a list l, where every element is itself a list of 3D arrays. Each array should have matching first and second dimentions. This function combines those arrays along the third dimension.
combine.streams(l)combine.streams(l)
l |
a list, where every element is itself a list of 3D arrays. |
a list of 3D arrays
Justin Silverman
Tailored for ALDEx3 where covariates are shared between massive numbers of linear regressions where only Y is changing.
fflm(Y, X, test = "t.HC3")fflm(Y, X, test = "t.HC3")
Y |
a numeric array (N x D X S) where D is the number of taxa/genes, N is the number of samples, and S is the number of posterior samples |
X |
a numeric matrix (N x P) where P is number of covariates |
test |
(default t.HC3), "t", t test is performed for each covariate (fast); "t.HC0" Heteroskedasticsity-Robust Standard Errors used (HC0; White's; slower); "t.HC3" (default) Heteroskedasticsity-Robust Standard Errors used (HC3; unlike HC0, this includes a leverage adjustment and is better for small sample sizes or when there are data with high leverage; slowest). To learn more about these, loko at Long and Ervin (2000) Using Heteroscedasticity Consistent Standard Errors in the Linear Regression Model, The American Statistician. |
A list of (P x D x S)-arrays with the OLS point estimates, the standard errors, and the two-sided p-values for each coefficient (P), of each model fit to each taxa (D) and each posterior sample (S)
Justin Silverman
This dataset contains a matrix and a data frame: genus-level microbiome profiles and corresponding sample metadata from a Crohn's disease case–control cohort. The dataset is used in examples and vignettes throughout the package.
gut_crohns_datagut_crohns_data
A named list of one matrix and a dataframe:
countsA matrix with read counts with 195 rows and 95 columns)
metadataA data frame with 95 rows and 7 columns containing subject-level covariates.
The counts matrix has one row per sample and one column per genus.
The metadata data frame has one row per sample with metadata, critically Health.status either CD or Control, and Average cell count per gram frozen feces.
Vandeputte D, Falony G, Vieira-Silva S, Wang J, Sailer M, Theis S, Raes J (2017). "Quantitative microbiome profiling links gut community variation to microbial load". Nature, 551, 507–511. doi:10.1038/nature24460
Closure operation
miniclo(X)miniclo(X)
X |
should be a D x N matrix, closes along the D dimension |
D x N matrix with columns summing to 1
Justin Silverman
This dataset contains a matrix and a data frame: genus-level microbiome profiles and corresponding sample metadata from an oral microbiome perturbation study. 28 participants' oral microbiomes were measured before, 15 minutes after, and 2 hours after perturbation with either a water control, antiseptic mouthwash, alchohol-free mouthwash, or soda. The dataset is used in examples and vignettes throughout the package.
oral_mouthwash_dataoral_mouthwash_data
A named list of one matrix and a dataframe:
countsA matrix with read counts with 116 rows and 81 columns)
metadataA data frame with 81 rows and 38 columns containing sample-level covariates.
The counts matrix has one row per sample and one column per genus.
The metadata data frame has one row per sample with metadata, critically the participant ID, flow cytometry average (and replicate) cells, time_c (time points), and treat (the perturbations )
Marotz C, Morton JT, Navarro P, Coker J, Belda-Ferre P, Knight R (2021). "Quantifying live microbial load in human saliva samples over time reveals stable composition and dynamic load". mSystems, 6(3), e01182-21. doi:10.1128/mSystems.01182-20
Test object contains elements
req(obj, names)req(obj, names)
obj |
object (list type) |
names |
names of elements required to be in the list |
Justin Silverman
Function for sampling Dirichlet random variables (base-2 normalized via log-sum-exp for stability)
rLogDirichlet(n, alpha)rLogDirichlet(n, alpha)
n |
number of samples |
alpha |
D-vector of Dirichlet parameters |
a D x n matrix of samples
Justin Silverman
Use this model when each sample has its own external scale measurement, such
as a flow-cytometry cell count, qPCR concentration, or spike-in estimate.
Put the log2 measurement for each sample in s.mu. Use s.sd to
describe the measurement uncertainty for each sample.
sample.sm( X, logComp, s.mu = NULL, s.var = NULL, s.cor = NULL, s.sd = NULL, s.cov = NULL )sample.sm( X, logComp, s.mu = NULL, s.var = NULL, s.cor = NULL, s.sd = NULL, s.cov = NULL )
X |
Model matrix passed automatically by |
logComp |
Monte Carlo log-compositions passed automatically by
|
s.mu |
Expected log2 scale for each sample, in the same order as the
columns of the count matrix. Must have length |
s.var |
Deprecated numeric vector of per-sample log2-scale variances.
Use |
s.cor |
Deprecated |
s.sd |
SD of the log2 scale for each sample. Must have length |
s.cov |
An |
For example, if four samples have estimated total cell counts of
1e8, 2e8, 8e7, and 1.5e8, set
s.mu = log2(c(1e8, 2e8, 8e7, 1.5e8)). If each measurement has an SD
of 0.5 on the log2 scale, set s.sd = rep(0.5, 4).
Use s.sd when measurement errors are independent. Use s.cov
instead when errors are correlated, for
example because several samples share a calibration standard. The diagonal
of s.cov contains variances, and its off-diagonal entries contain
covariances. Thus independent SDs of 0.5 can also be written as
s.cov = diag(0.5^2, 4).
Supply exactly one of s.sd or s.cov. With s.sd, each
sample receives its own independent scale draw. Samples in the same group
do not automatically share uncertainty as they do in coefficient-based
models.
A numeric matrix of dimension N x nsample giving Monte Carlo
draws of the log2 scale for each sample (rows) across nsample draws
(columns).
Kyle McGovern
# External total-cell estimates for four samples. cell_count <- c(1e8, 2e8, 8e7, 1.5e8) X <- matrix(1, 1, 4) logComp <- array(0, c(2, 4, 10)) logScale <- sample.sm( X, logComp, s.mu = log2(cell_count), s.sd = rep(0.5, 4) )# External total-cell estimates for four samples. cell_count <- c(1e8, 2e8, 8e7, 1.5e8) X <- matrix(1, 1, 4) logComp <- array(0, c(2, 4, 10)) logScale <- sample.sm( X, logComp, s.mu = log2(cell_count), s.sd = rep(0.5, 4) )
Implementation of SR-MEM: scale-reliant mixed effects models.
sr.mem(logW, formula, data, n.cores, method, mem.args)sr.mem(logW, formula, data, n.cores, method, mem.args)
logW |
a numeric array (N x D X S) where D is the number of taxa, N is the number of samples, and S is the number of posterior samples |
formula |
an lme4 formula with fixed and random effects for lmer |
data |
A data.frame with the random and fixed effects items in formula |
n.cores |
The number of cores for parallelization. If n.cores=1, no parallelization is used |
method |
The mem method to use: either nlme or lme4 |
mem.args |
Additional arguments including random or correlation |
A list of (P x D x S)-arrays with the fixed effect point estimates, the standard errors, degrees of freedom and the lower and upper p-values for each coefficient (P), of each model fit to each taxa (D) and each posterior sample (S) and a (Pr x D x S)-array with (Pr) random effects.
Kyle McGovern
Summarize an ALDEx3 result object
## S3 method for class 'aldex' summary(object, ignore.intercept = TRUE, ...)## S3 method for class 'aldex' summary(object, ignore.intercept = TRUE, ...)
object |
An object of class |
ignore.intercept |
(default=TRUE), ignore intercept when creating summary table |
... |
Additional arguments (currently ignored). |
Provides a summary of the adjusted p-values, estimates, and standard errors from an ALDEx3 result object.
This method extracts adjusted p-values from object$p.val.adj, along with
posterior estimates and standard errors averaged across Monte Carlo samples.
The result is returned as a long-format data.frame suitable for downstream
analysis or visualization.
A data.frame with columns parameter, entity,
p.val.adjusted, estimate, and std.error.
Justin Silverman
Total-sum scaling (TSS) starts from the assumption that total scale does not change systematically between experimental groups. Use this model when that is a reasonable starting point but you want to express uncertainty about it. For example, if a treatment is expected to redistribute features without changing total biomass, TSS centers the treatment-versus-control scale difference at zero.
tss.sm(X, logComp, gamma = 0.5)tss.sm(X, logComp, gamma = 0.5)
X |
Model matrix passed automatically by |
logComp |
Monte Carlo log-compositions passed automatically by
|
gamma |
Non-negative scalar. Standard deviation of the Gaussian random
shift applied to the scale-model coefficients (in log2 space).
|
Set gamma = 0 to enforce no scale difference. Setting
gamma = 0.5 keeps the expected difference at zero but allows each
scale-model coefficient to vary with SD 0.5 on the log2 scale. Larger values
express less confidence that total scale is unchanged.
ALDEx3 applies this uncertainty to model coefficients: for each draw,
, and the sample scales are
. Samples with the same covariate values therefore share a
random shift. This differs from sample.sm, which can draw
independent uncertainty for every sample.
A numeric matrix of dimension N x nsample giving Monte Carlo
draws of the log2 scale for each sample (rows) across nsample draws
(columns).
Justin Silverman
Nixon G, Gloor GB, Silverman JD (2025). "Incorporating scale uncertainty in microbiome and gene expression analysis as an extension of normalization". Genome Biology. doi:10.1186/s13059-025-03609-3
# Center the group-scale difference at zero, with moderate uncertainty. Y <- matrix(seq_len(60), nrow = 10) condition <- factor(c("control", "control", "control", "treatment", "treatment", "treatment")) fit <- aldex(Y, ~ condition, data.frame(condition), nsample = 20, scale = tss.sm, gamma = 0.5)# Center the group-scale difference at zero, with moderate uncertainty. Y <- matrix(seq_len(60), nrow = 10) condition <- factor(c("control", "control", "control", "treatment", "treatment", "treatment")) fit <- aldex(Y, ~ condition, data.frame(condition), nsample = 20, scale = tss.sm, gamma = 0.5)