Skip to main content
`Ctrl`+`K`

# 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={})

  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.

  This function uses the `MixedLM` class from statsmodels.

  Note

  Because the analysis iteratively runs many numeric optimizations, it can take longer than usual to finish. Please allow extra time for completion.

  Parameters:
    table table_like of shape (n_samples, n_features)

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

    metadata pd.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.

    formula str or generic Formula object

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

    grouping str, 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.

    pseudocount float, 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.

    draws int, optional

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

    p_adjust str, optional

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

    seed int, Generator or RandomState, optional

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

    re_formula str, optional

      Random coefficient formula. See `MixedLM.from_formula` for details.

    vc_formula str, optional

      Variance component formula. See `MixedLM.from_formula` for details.

    model_kwargs dict, optional

      Additional keyword arguments to pass to `MixedLM`.

    fit_method str 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, a default list of methods will be tried.

    fit_converge bool, optional

      If True, model fittings that were completed but did not converge will be excluded from the calculation of final statistics. Default is False.

    fit_warnings bool, 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_kwargs dict, 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`

  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.122148  0.403737
      1        Y1  treatment   128 -0.744061 -3.401978  1.581917  0.252057
      2        Y2       time   128  0.210769 -1.122148  1.532255  0.403737
      3        Y2  treatment   128  0.744061 -1.581917  3.401978  0.252057

           qvalue  Signif
      0  0.644470   False
      1  0.440581   False
      2  0.644470   False
      3  0.440581   False