This function implements an association test for multiple alpha diversity measures. The test is based on a mixed-effects model fitted by the `lmer` function from the `lmerTest` package. The function accepts a data object as input and returns a list of tests, one for each alpha diversity index.
Usage
generate_alpha_test_pair(
data.obj,
alpha.obj = NULL,
alpha.name = NULL,
depth = NULL,
time.var,
subject.var,
group.var,
adj.vars = NULL,
change.base = NULL
)Arguments
- data.obj
A list object in a format specific to MicrobiomeStat, which can include components such as feature.tab (matrix), feature.ann (matrix), meta.dat (data.frame), tree, and feature.agg.list (list). The data.obj can be converted from other formats using several functions from the MicrobiomeStat package, including: 'mStat_convert_DGEList_to_data_obj', 'mStat_convert_DESeqDataSet_to_data_obj', 'mStat_convert_phyloseq_to_data_obj', 'mStat_convert_SummarizedExperiment_to_data_obj', 'mStat_import_qiime2_as_data_obj', 'mStat_import_mothur_as_data_obj', 'mStat_import_dada2_as_data_obj', and 'mStat_import_biom_as_data_obj'. Alternatively, users can construct their own data.obj. Note that not all components of data.obj may be required for all functions in the MicrobiomeStat package.
- alpha.obj
An optional list containing pre-calculated alpha diversity indices. If NULL (default), alpha diversity indices will be calculated using mStat_calculate_alpha_diversity function from MicrobiomeStat package.
- alpha.name
A character vector with the names of alpha diversity indices to compute. Options include: "shannon", "simpson", "observed_species", "chao1", "ace", and "pielou".
- depth
An integer specifying the sequencing depth for the "Rarefy" and "Rarefy-TSS" methods. If NULL, no rarefaction is performed.
- time.var
A string representing the time variable's name in the metadata. The default is NULL.
- subject.var
A string specifying the subject variable column in the metadata.
- group.var
A string representing the group variable's name in the metadata.
- adj.vars
A character vector with the names of adjustment variables in the metadata.
- change.base
A value indicating the base level for the time variable. If provided, the specified level will be used as the reference category in the model. Default is NULL, which means the first level of the factor will be used.
Details
The mixed-effects model includes the time variable, group variable, and any additional adjustment variables as fixed effects, and the subject variable as a random effect.
The output is a list of coefficient tables, one for each alpha diversity index. Each table includes the term, estimate, standard error, t value, and p-value for each fixed effect in the model.
Examples
data(peerj32.obj)
generate_alpha_test_pair(
data.obj = peerj32.obj,
alpha.obj = NULL,
time.var = "time",
alpha.name = c("shannon", "simpson", "ace"),
subject.var = "subject",
group.var = NULL
)
#> Warning: It appears the data may not have been rarefied. Please verify.
#> Calculating shannon diversity...
#> Calculating simpson diversity...
#> Calculating ace diversity...
#> Diversity calculations complete.
#> $shannon
#> # A tibble: 2 × 5
#> Term Estimate Std.Error Statistic P.Value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 3.58 0.0226 158. 8.10e-55
#> 2 time2 0.0256 0.0261 0.982 3.37e- 1
#>
#> $simpson
#> # A tibble: 2 × 5
#> Term Estimate Std.Error Statistic P.Value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.951 0.00238 400. 4.96e-74
#> 2 time2 0.00303 0.00299 1.01 3.23e- 1
#>
#> $ace
#> # A tibble: 2 × 5
#> Term Estimate Std.Error Statistic P.Value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 101. 1.86 54.0 4.45e-30
#> 2 time2 3.29 1.45 2.27 3.39e- 2
#>
generate_alpha_test_pair(
data.obj = peerj32.obj,
alpha.obj = NULL,
time.var = "time",
alpha.name = c("shannon", "simpson", "ace"),
subject.var = "subject",
group.var = NULL,
change.base = "2"
)
#> Warning: It appears the data may not have been rarefied. Please verify.
#> Calculating shannon diversity...
#> Calculating simpson diversity...
#> Calculating ace diversity...
#> Diversity calculations complete.
#> $shannon
#> # A tibble: 2 × 5
#> Term Estimate Std.Error Statistic P.Value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 3.60 0.0226 159. 6.19e-55
#> 2 time1 -0.0256 0.0261 -0.982 3.37e- 1
#>
#> $simpson
#> # A tibble: 2 × 5
#> Term Estimate Std.Error Statistic P.Value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.954 0.00238 401. 4.37e-74
#> 2 time1 -0.00303 0.00299 -1.01 3.23e- 1
#>
#> $ace
#> # A tibble: 2 × 5
#> Term Estimate Std.Error Statistic P.Value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 104. 1.86 55.8 1.81e-30
#> 2 time1 -3.29 1.45 -2.27 3.39e- 2
#>
generate_alpha_test_pair(
data.obj = peerj32.obj,
alpha.obj = NULL,
time.var = "time",
alpha.name = c("shannon", "simpson", "ace"),
subject.var = "subject",
group.var = "group"
)
#> Warning: It appears the data may not have been rarefied. Please verify.
#> Calculating shannon diversity...
#> Calculating simpson diversity...
#> Calculating ace diversity...
#> Diversity calculations complete.
#> Complex model failed. Trying a simpler model...
#> Complex model failed. Trying a simpler model...
#> Complex model failed. Trying a simpler model...
#> $shannon
#> # A tibble: 4 × 5
#> Term Estimate Std.Error Statistic P.Value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 3.51 0.0356 98.6 8.55e-47
#> 2 groupPlacebo 0.104 0.0446 2.34 2.47e- 2
#> 3 time2 0.0606 0.0432 1.40 1.76e- 1
#> 4 groupPlacebo:time2 -0.0551 0.0541 -1.02 3.21e- 1
#>
#> $simpson
#> # A tibble: 4 × 5
#> Term Estimate Std.Error Statistic P.Value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.945 0.00384 247. 4.08e-64
#> 2 groupPlacebo 0.00920 0.00481 1.91 6.29e- 2
#> 3 time2 0.00610 0.00501 1.22 2.38e- 1
#> 4 groupPlacebo:time2 -0.00482 0.00628 -0.768 4.52e- 1
#>
#> $ace
#> # A tibble: 4 × 5
#> Term Estimate Std.Error Statistic P.Value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 97.9 3.11 31.5 9.18e-23
#> 2 groupPlacebo 4.27 3.90 1.09 2.83e- 1
#> 3 time2 4.94 2.42 2.04 5.47e- 2
#> 4 groupPlacebo:time2 -2.59 3.03 -0.853 4.04e- 1
#>
generate_alpha_test_pair(
data.obj = peerj32.obj,
alpha.obj = NULL,
time.var = "time",
alpha.name = c("shannon", "simpson", "ace"),
subject.var = "subject",
group.var = "group",
adj.vars = "sex"
)
#> Warning: It appears the data may not have been rarefied. Please verify.
#> Calculating shannon diversity...
#> Calculating simpson diversity...
#> Calculating ace diversity...
#> Diversity calculations complete.
#> Complex model failed. Trying a simpler model...
#> Complex model failed. Trying a simpler model...
#> Complex model failed. Trying a simpler model...
#> $shannon
#> # A tibble: 5 × 5
#> Term Estimate Std.Error Statistic P.Value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 3.53 0.0372 95.0 1.54e-43
#> 2 groupPlacebo 0.0993 0.0437 2.27 2.89e- 2
#> 3 time2 0.0606 0.0432 1.40 1.76e- 1
#> 4 sexmale -0.0575 0.0354 -1.62 1.21e- 1
#> 5 groupPlacebo:time2 -0.0551 0.0541 -1.02 3.21e- 1
#>
#> $simpson
#> # A tibble: 5 × 5
#> Term Estimate Std.Error Statistic P.Value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.947 0.00404 234. 1.02e-59
#> 2 groupPlacebo 0.00877 0.00477 1.84 7.39e- 2
#> 3 time2 0.00610 0.00501 1.22 2.38e- 1
#> 4 sexmale -0.00488 0.00371 -1.32 2.04e- 1
#> 5 groupPlacebo:time2 -0.00482 0.00628 -0.768 4.52e- 1
#>
#> $ace
#> # A tibble: 5 × 5
#> Term Estimate Std.Error Statistic P.Value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 97.0 3.45 28.1 5.35e-20
#> 2 groupPlacebo 4.49 3.96 1.13 2.68e- 1
#> 3 time2 4.94 2.42 2.04 5.47e- 2
#> 4 sexmale 2.45 3.78 0.648 5.25e- 1
#> 5 groupPlacebo:time2 -2.59 3.03 -0.853 4.04e- 1
#>
data("subset_pairs.obj")
generate_alpha_test_pair(
data.obj = subset_pairs.obj,
alpha.obj = NULL,
time.var = "Antibiotic",
alpha.name = c("shannon", "simpson", "ace"),
subject.var = "MouseID",
group.var = "Sex"
)
#> Warning: It appears the data may not have been rarefied. Please verify.
#> Calculating shannon diversity...
#> Calculating simpson diversity...
#> Calculating ace diversity...
#> Diversity calculations complete.
#> Complex model failed. Trying a simpler model...
#> Complex model failed. Trying a simpler model...
#> Complex model failed. Trying a simpler model...
#> $shannon
#> # A tibble: 4 × 5
#> Term Estimate Std.Error Statistic P.Value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 1.55 0.0947 16.3 1.11e-26
#> 2 SexM -0.390 0.135 -2.88 5.20e- 3
#> 3 AntibioticWeek 2 -0.0691 0.115 -0.601 5.51e- 1
#> 4 SexM:AntibioticWeek 2 0.474 0.165 2.88 6.33e- 3
#>
#> $simpson
#> # A tibble: 4 × 5
#> Term Estimate Std.Error Statistic P.Value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.631 0.0331 19.1 4.46e-31
#> 2 SexM -0.130 0.0474 -2.75 7.49e- 3
#> 3 AntibioticWeek 2 -0.0372 0.0412 -0.904 3.71e- 1
#> 4 SexM:AntibioticWeek 2 0.177 0.0589 3.01 4.43e- 3
#>
#> $ace
#> # A tibble: 4 × 5
#> Term Estimate Std.Error Statistic P.Value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 72.8 3.59 20.3 1.01e-33
#> 2 SexM -4.24 5.13 -0.827 4.11e- 1
#> 3 AntibioticWeek 2 -6.19 5.06 -1.23 2.27e- 1
#> 4 SexM:AntibioticWeek 2 8.15 7.23 1.13 2.66e- 1
#>
generate_alpha_test_pair(
data.obj = subset_pairs.obj,
alpha.obj = NULL,
time.var = "Antibiotic",
alpha.name = c("shannon", "simpson", "ace"),
subject.var = "MouseID",
group.var = "Sex",
change.base = "Week 2"
)
#> Warning: It appears the data may not have been rarefied. Please verify.
#> Calculating shannon diversity...
#> Calculating simpson diversity...
#> Calculating ace diversity...
#> Diversity calculations complete.
#> Complex model failed. Trying a simpler model...
#> Complex model failed. Trying a simpler model...
#> Complex model failed. Trying a simpler model...
#> $shannon
#> # A tibble: 4 × 5
#> Term Estimate Std.Error Statistic P.Value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 1.48 0.0947 15.6 1.68e-25
#> 2 SexM 0.0841 0.135 0.621 5.37e- 1
#> 3 AntibioticBaseline 0.0691 0.115 0.601 5.51e- 1
#> 4 SexM:AntibioticBaseline -0.474 0.165 -2.88 6.33e- 3
#>
#> $simpson
#> # A tibble: 4 × 5
#> Term Estimate Std.Error Statistic P.Value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.593 0.0331 17.9 2.15e-29
#> 2 SexM 0.0475 0.0474 1.00 3.19e- 1
#> 3 AntibioticBaseline 0.0372 0.0412 0.904 3.71e- 1
#> 4 SexM:AntibioticBaseline -0.177 0.0589 -3.01 4.43e- 3
#>
#> $ace
#> # A tibble: 4 × 5
#> Term Estimate Std.Error Statistic P.Value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 66.6 3.59 18.6 4.09e-31
#> 2 SexM 3.91 5.13 0.762 4.48e- 1
#> 3 AntibioticBaseline 6.19 5.06 1.23 2.27e- 1
#> 4 SexM:AntibioticBaseline -8.15 7.23 -1.13 2.66e- 1
#>