Hi! I’m trying to obtain the value of the intercept for each of the features in my MaAsLin2 results. Is it possible to generate a column for this in the output? I have fitted a few models with interaction terms where this would be especially helpful for following up with additional analyses and plots. Thanks!
If you set Maaslin2(..., save_models = TRUE)
, among the outputs you’ll get an rds file containing the model fits. You can extract the intercepts from there. Using the example in the package:
fit_data <- Maaslin2(input_data,
input_metadata,
output = '~/tmp',
transform = "AST",
fixed_effects = c('diagnosis', 'dysbiosisnonIBD','dysbiosisUC','dysbiosisCD', 'antibiotics', 'age'),
random_effects = c('site', 'subject'),
normalization = 'NONE',
reference = 'diagnosis,nonIBD',
standardize = FALSE,
save_models = TRUE)
models <- readRDS("~/tmp/fits/models.rds")
models[[1]] |> summary() |> coefficients()
# Estimate Std. Error df t value Pr(>|t|)
#(Intercept) 2.655611e-02 0.0092456815 121.0273 2.8722722 0.004813519
#diagnosisCD -1.812594e-02 0.0088601782 125.6027 -2.0457758 0.042866199
#diagnosisUC 4.091809e-03 0.0096240661 123.9921 0.4251643 0.671453859
#dysbiosisnonIBDYes -6.719608e-03 0.0071908428 1551.8470 -0.9344674 0.350208208
#dysbiosisUCYes -1.056876e-02 0.0054566951 1498.1532 -1.9368420 0.052952474
#dysbiosisCDYes 6.020841e-04 0.0041147376 1567.5584 0.1463238 0.883684589
#antibioticsYes -6.128646e-03 0.0033373554 1526.8368 -1.8363779 0.066496111
#age 4.813489e-05 0.0001984042 122.9060 0.2426102 0.808711741
Thanks for your quick response! I just tried that and unfortunately received this error:
Error in Maaslin2(input_data = asv.relabun, input_metadata = input_metadata, :
unused argument (save_models = TRUE)
In case it’s helpful, this is what I am running. It normally works fine without the “save_models = TRUE” parameter:
mod3 <- Maaslin2(
input_data = asv.relabun, input_metadata = input_metadata, output = "output_file",
fixed_effects = c("timepoint_type","GroupHL_01","GroupHL_int"),
random_effects = "participant_id",
reference = "timepoint_type,Pre",
normalization = "NONE",
save_models = TRUE)
Oops, I didn’t realize that option hadn’t made it to the BioC version yet. You can get it by installing the development version from github: devtools::install_github("biobakery/Maaslin2")
That worked, and this is exactly what I was looking for. Thanks so much for your help!
1 Like