skbio.stats.composition.dirmult_lme#

skbio.stats.composition.dirmult_lme(table, metadata, formula, grouping, pseudocount=0.5, draws=128, p_adjust='holm', seed=None, re_formula=None, vc_formula=None, model_kwargs={}, fit_method=None, fit_converge=False, fit_warnings=False, fit_kwargs={})[source]#

Fit a Dirichlet-multinomial linear mixed effects model.

Added in version 0.7.0.

The Dirichlet-multinomial distribution is a compound distribution that combines a Dirichlet distribution over the probabilities of a multinomial distribution. This distribution is used to model the distribution of species abundances in a community.

To fit the linear mixed effects model we first fit a Dirichlet-multinomial distribution for each sample, and then we compute the fold change and p-value for each feature. The fold change is computed as the slopes from the resulting model. Statistical tests are then performed on the posterior samples, drawn from each Dirichlet-multinomial distribution. The log-fold changes as well as their credible intervals, the p-values and the multiple comparison corrected p-values are reported.

Changed in version 0.7.4: Significantly improved computational efficiency under the default configuration by introducing an analytical kernel fitting all features at once. See Notes for details.

Parameters:
tabletable-like of shape (n_samples, n_features)

A matrix containing count or proportional abundance data of the samples. Counts are recommended over proportions for lower statistical uncertainty. See Notes of dirmult_ttest for details.

metadatapd.DataFrame or 2-D array-like

The metadata for the model. Rows correspond to samples and columns correspond to covariates in the model. Must be a pandas DataFrame or convertible to a pandas DataFrame.

formulastr or generic Formula object

The formula defining the model. Refer to Patsy’s documentation on how to specify a formula.

groupingstr, pd.Series or 1-D array-like

A vector or a metadata column name indicating the assignment of samples to groups. Samples are independent between groups during model fitting.

pseudocountfloat, optional

A non-zero value added to the input counts to ensure that all of the estimated abundances are strictly greater than zero. Default is 0.5.

drawsint, optional

Number of draws from the Dirichlet-multinomial posterior distribution. Default is 128.

p_adjuststr, optional

Method to correct p-values for multiple comparisons. Options are: Bonferroni (“bonf”/”bonferroni”), Holm-Bonferroni (“holm”/”holm-bonferroni”, default), Benjamini-Hochberg (“bh”/”benjamini-hochberg”), and Benjamini-Yekutieli (“by”/”benjamini-yekutieli”), or any method supported by statsmodels’ multipletests function. Case-insensitive. If None, no correction will be performed.

seedint, Generator or RandomState, optional

A user-provided random seed or random generator instance for drawing from the Dirichlet distribution. See details.

re_formulastr, optional

Random coefficient formula. See MixedLM.from_formula for details.

vc_formulastr, optional

Variance component formula. See MixedLM.from_formula for details.

model_kwargsdict, optional

Additional keyword arguments to pass to MixedLM.

fit_methodstr or list of str, optional

Optimization method for model fitting. Can be a single method name, or a list of method names to be tried sequentially. See statsmodels.optimization for available methods. If None (default), the default model is solved in closed form and no iterative optimizer is used. See Notes.

fit_convergebool, optional

If True, model fittings that were completed but did not converge will be excluded from the calculation of final statistics. Default is False. Setting this selects the iterative optimizer, which is what reports convergence.

fit_warningsbool, optional

Issue warnings if any during the model fitting process. Default is False. Warnings are usually issued when the optimization methods do not converge, which is common in the analysis. Default is False.

fit_kwargsdict, optional

Additional keyword arguments to pass to MixedLM.fit.

Returns:
pd.DataFrame

A table of features and covariates, their log-fold changes and other relevant statistics.

  • FeatureID: Feature identifier, i.e., dependent variable.

  • Covariate: Covariate name, i.e., independent variable.

  • Reps: Number of Dirichlet-multinomial posterior draws that supported the reported statistics, i.e., the number of successful model fittings on this feature. Max: draws (if none failed). Min: 0 (in which case all statistics are NaN).

  • Log2(FC): Expected log2-fold change of abundance from the reference category to the covariate category defined in the formula. The value is expressed in the center log ratio (see clr) transformed coordinates. The reported value is the average of all of the log2-fold changes computed from each of the posterior draws.

  • CI(2.5): 2.5% quantile of the log2-fold change. The reported value is the minimum of all of the 2.5% quantiles computed from each of the posterior draws.

  • CI(97.5): 97.5% quantile of the log2-fold change. The reported value is the maximum of all of the 97.5% quantiles computed from each of the posterior draws.

  • pvalue: p-value of the linear mixed effects model. The reported value is the average of all of the p-values computed from each of the posterior draws.

  • qvalue: Corrected p-value of the linear mixed effects model for multiple comparisons. The reported value is the average of all of the q-values computed from each of the posterior draws.

  • Signif: Whether the covariate category is significantly differentially abundant from the reference category. A feature-covariate pair marked as “True” suffice: 1) The q-value must be less than or equal to the significance level (0.05). 2) The confidence interval (CI(2.5)..CI(97.5)) must not overlap with zero.

See also

dirmult_ttest
statsmodels.formula.api.mixedlm
statsmodels.regression.mixed_linear_model.MixedLM

Notes

The default model has one random intercept per group, for which the REML fit has a closed form. In that case all features of a replicate are fitted at once by maximizing the one-dimensional profile likelihood over the ratio of the group variance to the residual variance, instead of running one numeric optimization per feature through MixedLM.fit. The estimates, standard errors, p-values and confidence intervals are those of MixedLM, computed from the same REML profile likelihood and the same observed information matrix.

Supplying re_formula, vc_formula, model_kwargs, fit_method, fit_kwargs or fit_converge selects a model, an optimizer, or a convergence filter that has no closed form, and falls back to fitting each feature separately with MixedLM.

Changed in version 0.7.4: The default random intercept model is fitted for all features at once. Results may differ slightly from earlier versions on replicates where the iterative optimizer of MixedLM stopped short of the optimum. Any other configuration falls back to running one numeric optimization per feature per draw using MixedLM, which can take much longer to finish.

Examples

>>> import pandas as pd
>>> from skbio.stats.composition import dirmult_lme
>>> table = pd.DataFrame(
...     [[1.00000053, 6.09924644],
...      [0.99999843, 7.0000045],
...      [1.09999884, 8.08474053],
...      [1.09999758, 1.10000349],
...      [0.99999902, 2.00000027],
...      [1.09999862, 2.99998318],
...      [1.00000084, 2.10001257],
...      [0.9999991, 3.09998418],
...      [0.99999899, 3.9999742],
...      [1.10000124, 5.0001796],
...      [1.00000053, 6.09924644],
...      [1.10000173, 6.99693644]],
...     index=['u1', 'u2', 'u3', 'x1', 'x2', 'x3', 'y1', 'y2', 'y3', 'z1',
...            'z2', 'z3'],
...     columns=['Y1', 'Y2'])
>>> metadata = pd.DataFrame(
...     {'patient': [1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4],
...      'treatment': [1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2],
...      'time': [1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3],},
...     index=['x1', 'x2', 'x3', 'y1', 'y2', 'y3', 'z1', 'z2', 'z3', 'u1',
...            'u2', 'u3'])
>>> result = dirmult_lme(table, metadata, formula='time + treatment',
...                      grouping='patient', seed=0, p_adjust='sidak')
>>> result
  FeatureID  Covariate  Reps  Log2(FC)   CI(2.5)  CI(97.5)    pvalue  \
0        Y1       time   128 -0.210769 -1.532255  1.123699  0.404007
1        Y1  treatment   128 -0.744061 -3.401965  1.500992  0.250282
2        Y2       time   128  0.210769 -1.123699  1.532255  0.404007
3        Y2  treatment   128  0.744061 -1.500992  3.401965  0.250282

     qvalue  Signif
0  0.644793   False
1  0.437922   False
2  0.644793   False
3  0.437922   False
On this page