Questions about choosing analysis method Masslin3

Downstream analysis and statistics MaAsLin
1

Hi,

There are many choices for analysis methods, normalizations, and transformations in Masslin3. I don’t know which will be the optimal methods for my data analysis. My taxa table includes relative abundance values in percentage (ranged 0-100 not 0-1) at species level. Values of each taxa is summed up to 100. By checking the histogram of some species in my data, over 70% of the relative abundance values of a certain species are 0 (highly skewed, and zero inflation). Several bigger values above 5% in some species (outliers?).

Based on the tutorial, I think the optimal analysis method should be CPLM, while I saw the below in the tutorial of Masslin2: Regarding whether to use LM , CPLM , or other models, intuitively, CPLM or a zero-inflated alternative should perform better in the presence of zeroes, but based on our benchmarking, we do not have evidence that CPLM is significantly better than LM in practice. Which method should I try?

Following the analysis method selected from the above, which normalization is suitable? Is the default method TSS good enough?

For the transformation, is AST the best choice? Do I need to scale my original data from 0-100 to 0-1 before conducting AST in Masslin3()?

I would appreciate it if someone could provide any suggestions. Thank you!

2

Hi,

The defaults in MaAsLin 3 (TSS normalization, no filtering, log transforming, split prevalence and abundance analysis) will probably be best for your data as they are for most scenarios. TSS normalization (relative abundances), log (base 2) transformations, and splitting prevalence (presence/absence) from abundance (how much is there, if it’s there) make for very clean interpretations, and we’ve extensively benchmarked these and found them to perform well. To answer your specific question, the TSS normalization will take care of converting your 0-100s to 0-1s.

You might decide to turn on or off the median_comparison_abundance parameter based on reading this section, but I’d recommend leaving median_comparison_prevalence off.

Let me know if you have additional questions.

Will

3

Hi Will,

Thank you for the reply! I tried the models using TSS and LOG in MaAslin3. The features I want to test are time, treatment and their interactions. I used two different formulas. One is ~ Time + Treatment + Time:Treatment + (1|SubjectID). I have six samples from the same subject. Then I created a new interaction column in my metadata named as Trt_time to manually combine the time and treatment. Then run another model with ~ Trt_time + (1|SubjectID). While the two models have different outputs of species (some overlapped), p and q values in significant results for the interaction. Why? I set the reference level for the model 1 as reference = c(‘Treatment,TrtA’ , ‘Time,baseline’) and model 2 as reference = c(‘Trt_time,TrtA_baseline’).

Another question for the model 2: If I set the reference level: reference = c(‘Trt_time,TrtA_baseline’). Any significant associations in the final output were generated by only comparing others with the reference level, right? Then I created a contrast matrix to compare different treatments at the same time points (e.g., trtA_baseline vs trtB_baseline). In my significant results for the model ~ Trt_time + (1|SubjectID), I has one row/species having qval_joint<0.01 for trtB_baseline, I believe this siginificant result is against my reference ‘trtA_baseline’. While in my contrast test output, all the p and q values in the rows of abundance and prevalence results are NA. Is there any possible error? I created the contrast matrix using ‘name’ column in the model result output.

4

Regarding your two models, assuming both time and treatment are categorical variables, both of those models will fit a coefficient to each category. However, what that coefficient means and how significance is tested is very different.

In the first model, [the time groups are tested against the baseline time] for the baseline treatment group, [the treatment groups are tested against the baseline treatment] for the baseline time group, and the interactions then test whether the treatment effect differs from its effect at the reference time. It’s quite complicated and requires being very careful about your interpretations.

In the second model, all groups are simply tested against the single baseline time-treatment group (which may or may not be what you care about).

Regarding your second question for model 2, that is correct: everything is compared against the overall reference. Regarding the NA contrast test, can you provide a reproducible example that I can test on my end? It’s hard to say whether the issue is in defining the contrast matrix, in trying to get coefficients for treatment-time groups that are never present, or in the MaAsLin 3 code itself.

Will

5

For the answer 1, I would like to double check if I understand correctly. I have two levels in treatment : TrtA (reference) and TrtB. I have three levels in time: day0 (reference), day10, day22. If I build a model as below:

DAmodel1= maaslin3(input_data = otu_df,
input_metadata = metadata,
output = “output_damodel1”,
formula = ‘~ Treatment + Time + Treatment:Time + (1|SubjectID)’,
normalization = “TSS”,
transform = “LOG”,
reference = c(“Treatment,TrtA”, ‘Time,day0’),
standardize = FALSE,
augment = TRUE,
min_prevalence = 0.1,
max_significance = 0.1,
median_comparison_abundance = TRUE,
median_comparison_prevalence = FALSE,
max_pngs = 250)

For the main treatment effects: It compared TrtB with the reference TrtA across all the times? You mentioned [the treatment groups are tested against the baseline treatment] for the baseline time group. It makes me a little confused, the baseline time group here is only day0 group or regardless the different time points?
For the main time effects: My understanding is to compare day10 with day0 and compare day22 with day0 regardless different treatments here.
For the interaction effects: Compare TrtA_day10 vs TrtA_day0, TrtA_day22 vs TrtA_day0, TrtB_day10 vs TrtB_day0, and TrtB_day22 vs TrtB_day0, right?

For the model2, my R codes are below: Based on the levels of treatments and times, I have six different levels in Trt_time.
DAmodel2= maaslin3(input_data = otu_df,
input_metadata = metadata,
output = “output_damodel2”,
formula = ‘~ Trt_time + (1|SubjectID)’,
reference = c(‘Trt_time,TrtA_day0’),
normalization = “TSS”,
transform = “LOG”,
standardize = FALSE,
augment = TRUE,
min_prevalence = 0.1,
max_significance = 0.1,
median_comparison_abundance = TRUE,
median_comparison_prevalence = FALSE,
max_pngs = 250)

contrast_mat ← matrix(c(
1, -1, 0, 0, 0, 0, # TrtA_day0 vs. TrtB_day0
0, 0, 1, -1, 0, 0, # TrtA_day10 vs. TrtB_day10
0, 0, 0, 0, 1, -1, # TrtA_day22 vs. TrtB_day22
0, 1, 0, -1, 0, 0, # TrtB_day0 vs. TrtB_day10
0, 1, 0, 0, 0, -1 # TrtB_day0 vs. TrtB_day22
), ncol = 6, byrow = TRUE)

colnames(contrast_mat) ← c(“Trt_timeTrtA_day0”, “Trt_timeTrtB_day0”,
“Trt_timeTrtA_day10”, “Trt_timeTrtB_day10”,
“Trt_timeTrtA_day22”, “Trt_timeTrtB_day22”)

rownames(contrast_mat) ← c(“Baseline_Comparison”,
“10day_Comparison”,
“22day_Comparison”,
“TrtB_day0_vs_4day”,
“TrtB_day0_vs_22day”)

contrast_out ← maaslin_contrast_test(maaslin3_fit = DAmodel2,
contrast_mat = contrast_mat)

significant_abundance_results ← contrast_out$fit_data_abundance$results %>%
filter(qval_individual <= 0.05 | qval_joint <= 0.05) #####???Got 0 rows. All the q and p values are NAs, same for the prevalence results.

6

For model 1, the coefficient that’s reported as TrtB will be the difference between TrtA and TrtB only in the baseline time group (day0). It’s not comparing them across all time points because of the interaction effect you have. If you want to compare across all time points, you should use ~ Treatment + Time + (1|SubjectID), which assumes that the treatment effect is the same at each time but that each time has a difference from day0 that applies equally to both treatment groups.

The time effect as you have it here works analogously, so if you want to compare day10 with day0 and day22 with day0 regardless of treatment, you should not use the interaction term.

The interaction terms you have right now will be:
TrtB:day10, TrtB:day22 which gives the difference in the day10 coefficient for the TrtB group relative to the TrtA group and likewise for day22. The other terms will appear as main effects in your output because they involve either the reference time or the reference group. You can check this in the name column of your results.

Are the matrix names you’re using for the contrast matrix actually all in the names column of your maaslin3 output table? Trt_timeTrtA_day0 for instance should be the baseline (intercept) and therefore not show up in the output table. Does the error column in the contrast fits say anything like Predictors not in the model had non-zero contrast values? Some error should be reported.

If it’s still unclear why the contrast tests are not working, can you provide a small example of the data that I can run so I can evaluate what the issue is?

Will

7

Thank you so much for the detailed explanation, then now I understand why I got different outputs than what I expected. And you are right, I just noticed that I had ‘Predictors not in the model had non-zero contrast values’ in the error column because I set wrong names in the contrast matrix. I also asked a statistician in my department. She suggested not only using an interaction term or a variable representing interactions (like Trt_time in my model2) to build a model. So I think I will stick with the formular like ~ Treatment + Time + Treatment:Time. While based on your explanation, I don’t think how I set up the model and reference now can compare the pairs I want.

To reclarify my study aims: Two treatment including TrtA(Control) and TrtB(Interested treatment). I used Control and Interest instead of TrtA and B for better description in the below. Three time points: day0(Baseline), day10 and day22. The comparisons I want to get in the output are:
Control vs Interest at day0,
Control vs Interest at day10,
Control vs Interest at day22,
Interest_day10 vs day0,
Interest_day22 vs day0,
Interest_day22 vs day10. I want to build a correct model and may need do level wise comparisons by a contrast test based on the correct model. How do I set up the formula and reference in the model? Thank you!

8

Glad you got that figured out with the contrast matrix!

To make those comparisons, you’ll want to fit ~ Treatment + Time + Treatment:Time with Control and day0 as the references and then analyze as follows (the coefficients here won’t literally be what’s in the name column because I’m not sure what your levels are called, but it should be close):

  1. Control vs. interest at day0: TreatmentInterest
  2. Control vs. Interest at day10: TreatmentInterest + TreatmentInterest:Timeday10
  3. Control vs. Interest at day22: TreatmentInterest + TreatmentInterest:Timeday22
  4. Interest_day10 vs day0: Timeday10 + TreatmentInterest:Timeday10
  5. Interest_day22 vs day0: Timeday22 + TreatmentInterest:Timeday22
  6. Interest_day22 vs day10: Timeday22 + TreatmentInterest:Timeday22 - TreatmentInterest:Timeday10 - Timeday10

Put the things that involve sums or differences in the contrast matrix to get proper standard errors and p-values. Linear models with interactions can get thorny :slight_smile:

You can check that these are as expected by doing the addition/subtraction on the coefficients in this toy example:

df <- data.frame('Time' = factor(c('day0', 'day10', 'day22', 'day0', 'day10', 'day22'), 
                                   levels = c('day0', 'day10', 'day22')), 
                 'Treatment' = factor(c('Control', 'Interest', 'Control', 'Interest', 'Control', 'Interest'), 
                                   levels = c('Control', 'Interest')), 
                 'outcome' = c(1, 3, 7, 11, 16, 15))
lm(outcome ~ Time + Treatment + Time:Treatment, df)

Also, it might be easier to do all of this by repeatedly fitting the MaAsLin 3 model with different references and just taking the comparisons you care about. It’ll take longer to run, but the risk of messing something up is lower, and you’ll be more confident in what you did.

9

Thank you for the info, it’s really helpful! Now I got the correct models and outputs for the contrast tests. I would like to make sure if I interpret the significant results correctly. To check the significant differences associated with prevalence or abundance, I look at the coefficients and q values. For example to compare Control vs. interest at day0: TreatmentInterest, if the model is ‘prevalence’, the coefficient is positive and q_individual is very small, it means the prevalence of the correspondering species in Interest group is significantly more than that in Control, right? Then if the model is ‘abundance’, it means the relative abundance of the species in Interest group is significantly higher than that in Control.

Questions 2 about the q_joint, what I understand is it indicates the significant association with either prevalence or abundance. So for the same species with the same value and ‘name’, the q_joint values are the same no matter it’s prevalence or abundance model? Then e.g., if I get a q_individual>0.05 but q_joint<0.05 in a prevalence model, it suggests that the abundance model for the same species and same ‘name’ has a small q individual value <0.05? I am trying to understand the relationship and difference between individual and joint q values.

Question 3 is I got the some big coefficients above 10 (2 absolute values even more than 15, see as below), which are suspicious based on the tutorial. But I didn’t receive any errors, the metadata are not collinear, and 15 samples in each group. I got q values = 0 (both individual and joint) with the big coefficients. How to explain the 0 value here, very different so the q value is extremely small? Or is it still reliable to trust these results?

image
image995×241 6.28 KB

Appreciate for your kind help!

10

Your first paragraph is correct.

The joint q-value answers the question: “Is at least one of the associations (prevalence or abundance) significant?” Therefore, the joint q-value will be the same for an association between the prevalence and abundance. If the prevalence q-value is above 0.05 and the joint q-value is below 0.05, the abundance q-value will likely be below 0.05. However, this isn’t strictly true because the joint q-values are generated from FDR correction over the joint p-values, and the individual q-values are generated from FDR correction over the individual p-values separately.

I’m not sure if you’re using a MaAsLin version that’s up to date enough to warn about it, but your very large coefficients with tiny p-values on prevalence associations in a random intercept model are strongly suggestive of constrained model fits that can be misleading. I’d run your model with small_random_effects=TRUE as described in the wiki here and see whether things hold up. Even though the wiki recommends at least 4 samples per group, you might need even more with a model containing as many terms as yours.

11

Thank you for the info! I got no big coefficients above 10 and zero q values after adding small_random_effects=TRUE.

12

Sounds right - glad to hear!

13

Hi Will,
Thank you for support with every questions here.
Since my issue is very closely related to this, I think I might post the question here instead of making new post.
I am planning to fit the model with interactions (2 x 4) and run pairwise comparison, which is very similar to the OP example. I am testing the contrast comparison with model fitting with 1 factor first.

maaslin3_1 ← phyloseq%>%
ps_filter(Protein_sources %in% ‘Canola’) %>%
ps_filter(!(Piglet_number %in% c(‘202’,‘212’))) %>%
tax_fix(min_length = 3,unknowns = c(‘NA’,“uncultured”), sep = " unidentified ", anon_unique = TRUE, suffix_rank = “current”) %>%
tax_transform(trans = ‘identity’,rank = ‘Species’)%>%
otu_get() %>%
maaslin3(input_metadata = metadata3_filtered ,output = ‘maaslin3’,formula = ‘~ Enzymes_treatment +(1|Batch) + (1|Room)’,
min_abundance = 0.05,cores = 2,small_random_effects=F, verbosity = ,save_models = T)

contrast_mat ← matrix(rep(0), ncol = 3, nrow = 6, byrow = TRUE)

colnames(contrast_mat) ← maaslin3_1$fit_data_abundance$results$name %>% as.factor() %>% levels()
rownames(contrast_mat) ← c(‘All vs Mul’,‘All vs Protease’,‘Mul vs Pro’,‘All vs Con’,‘Mul vs Con’,‘Pro vs Con’)

contrast_mat[1,] ← c(1,-1,0)
contrast_mat[2,] ← c(1,0,-1)
contrast_mat[3,] ← c(0,1,-1)
contrast_mat[4,] ← c(1,0,0)
contrast_mat[5,] ← c(0,1,0)
contrast_mat[6,] ← c(0,0,1)

contrast_mat
Enzymes_treatmentAll Enzymes_treatmentMultigrains Enzymes_treatmentProtease
All vs Mul 1 -1 0
All vs Protease 1 0 -1
Mul vs Pro 0 1 -1
All vs Con 1 0 0
Mul vs Con 0 1 0
Pro vs Con 0 0 1
contrast_out ← maaslin_contrast_test(maaslin3_fit = maaslin3,
contrast_mat = contrast_mat)

The issue I had is, when I compare the results between the original maaslin3 fit and the contrast fit in case of comparison with reference, the coef, nul.hypo, stdder are the same but the pvalue, qvalue (individual and joint) are different (qvalue in contrast was lower). Below is an examples from the output.

original abundance fit: uncultured Gemmiger sp. /Enzymes_treatment/ All/ Enzymes_treatmentAll/ 6.966441 / 0.13719521 /1.965344 / 0.007439222 / 0.4562723/ 0.0148231/ 0.5039855 / /linear/ 23 / 14

contrast fit: uncultured Gemmiger sp./All vs Con/6.966441127 /0.1371952/1.9653441/ 8.479825e-04/ 2.426535e-02/ 1.695246e-03/ 2.660232e-02/ NA/ abundance/ 23/14

Am I missing something in here? I though the re-run comparison with the reference group in contrast test should yield the same result as the original fit?
Thank you for your help.
Lenny

14

Hi Lenny,

Would you be able to provide a reproducible example (subset of the data, possibly deidentified, that recreates this issue)? The fact that the q-values are different is not necessarily surprising to me since the FDR is over a different set of p-values (all covariates originally vs. just the contrast tests here). However, I’m not sure why the p-values are different if the coefficient, standard error, and null hypothesis are the same. There’s one monte-carlo step that might be giving some variation, but the results show a larger difference than I’d expect. As an initial try, you could test with different seeds and see how much variation in the p-value you’re seeing.

Will

15

Hi Will,
Thank you for your prompt reply. Do you mean different set.seed()? I tried and results are still the same.

My experimental design consisted of 2 different feed (Ca and So) mix with 4 enzymes (Phytase-only -Control, Pro, Mul, All), then feed to animal, which made a classic 2x4 factorial designs. Batch and Room was added as random effect. Before fit the whole data to maaslin3 with interaction model, I tested with one feed first (Ca). I subset the data with samples from Ca and set Enzymes as factor with Batch and Room are random effect. I used formula = ‘~ Enzymes_treatment +(1|Batch) + (1|Room)’.The code is below:

maaslin3_1 ← maaslin3(input_data = subset_count, input_metadata = metadata3_filtered ,output = ‘maaslin3’,formula = ‘~ Enzymes_treatment +(1|Batch) + (1|Room)’,reference = c(‘Enzymes_treatments,Phytase_only’), min_abundance = 0.05,cores = 2,small_random_effects=F, save_models = T)

contrast matrix was established with column names from column “name” of the maaslin3_1 ( Pro, Mul, All), and row is the contrast that I want to investigate: ‘All vs Mul’,‘All vs Protease’,‘Mul vs Pro’,‘All vs Con’,‘Mul vs Con’,‘Pro vs Con’. I put 1 and -1 to the pair that I want to compare and 0 to rest. For the comparison between treatments vs Con, I leave 1 in the column of treatment that I want to re-run, then 0, 0 to other column.

contrast_mat ← matrix(rep(0), ncol = 3, nrow = 6, byrow = TRUE)

colnames(contrast_mat) ← maaslin3_1$fit_data_abundance$results$name %>% as.factor() %>% levels()
rownames(contrast_mat) ← c(‘All vs Mul’,‘All vs Pro’,‘Mul vs Pro’,‘All vs Control’,‘Mul vs Control’,‘Pro vs Control’)

contrast_mat[1,] ← c(1,-1,0)
contrast_mat[2,] ← c(1,0,-1)
contrast_mat[3,] ← c(0,1,-1)
contrast_mat[4,] ← c(1,0,0)
contrast_mat[5,] ← c(0,1,0)
contrast_mat[6,] ← c(0,0,1)

I have subset the results from both maaslin3_1 and constrast as below:
maaslin3_1

feature metadata value name coef null_hypothesis stderr pval_individual qval_individual pval_joint qval_joint error model N N_not_zero
1 uncultured Gemmiger sp. Enzymes_treatment All Enzymes_treatmentAll 6.966441127 0.137195212 1.965344122 0.007365382 0.451743401 0.014676514 0.499001484 NA linear 23 14
3 Escherichia coli Enzymes_treatment Multigrains Enzymes_treatmentMultigrains 7.629424946 0.131611244 0.18771331 0.063535109 0.995774868 0.123033508 0.999705319 NA linear 23 8
4 Floccifex porci Enzymes_treatment All Enzymes_treatmentAll -4.631491308 0.137195212 1.638787302 0.024396598 0.995774868 0.048198001 0.999705319 NA linear 23 12
6 Limosilactobacillus reuteri Enzymes_treatment Multigrains Enzymes_treatmentMultigrains -3.051987987 0.131611244 1.333893365 0.064941839 0.995774868 0.125666236 0.999705319 NA linear 23 14
7 Paratractidigestivibacter faecalis Enzymes_treatment Multigrains Enzymes_treatmentMultigrains 5.479961949 0.131611244 2.038042088 0.05631977 0.995774868 0.109467624 0.999705319 NA linear 23 10
8 Prevotella sp. P3-122 Enzymes_treatment Multigrains Enzymes_treatmentMultigrains 3.535368459 0.131611244 1.371536577 0.054758432 0.995774868 0.106518377 0.999705319 NA linear 23 16

contrast

feature test coef null_hypothesis stderr pval_individual qval_individual pval_joint qval_joint error model N N_not_zero
1 Escherichia coli Mul vs Control 7.629424946 0.131611244 0.18771331 0 0 0 0 NA abundance 23 8
3 uncultured Gemmiger sp. All vs Control 6.966441127 0.137195212 1.965344122 0.000852433 0.024392696 0.001704139 0.026741877 NA abundance 23 14
4 Floccifex porci All vs Control -4.631491308 0.137195212 1.638787302 0.00629218 0.156046061 0.012544768 0.170608848 NA abundance 23 12
5 Escherichia coli All vs Control 1.673184978 0.137195212 0.159000086 0.011208247 0.23163711 0.02229087 0.252629856 NA abundance 23 8
6 Paratractidigestivibacter faecalis Mul vs Control 5.479961949 0.131611244 2.038042088 0.013405374 0.249769021 0.026631045 0.272008957 NA abundance 23 10
10 Prevotella sp. P3-122 Mul vs Control 3.535368459 0.131611244 1.371536577 0.028621034 0.443626025 0.056422904 0.479594685 NA abundance 23 16
12 Limosilactobacillus reuteri Mul vs Control -3.051987987 0.131611244 1.333893365 0.036925876 0.528324066 0.072488231 0.568753812 NA abundance 23 14
17 Floccifex porci Mul vs Control -2.370745923 0.131611244 1.727433492 0.177839924 1 0.324052809 0.999991666 NA abundance 23 12
18 Limosilactobacillus reuteri All vs Control -0.821854061 0.137195212 1.333893365 0.491901061 1 0.565657478 0.999991666 NA abundance 23 14
19 Paratractidigestivibacter faecalis All vs Control 0.824710876 0.137195212 1.712963775 0.693539267 1 0.618858121 0.999991666 NA abundance 23 10
20 Prevotella sp. P3-122 All vs Control 1.663584533 0.137195212 2.038875704 0.459632814 1 0.708003305 0.999991666 NA abundance 23 16

I am a bit confused of which results should I report. Does it mean that if I increase the number of contrasts high enough (more than 6 in this case), the FDR might adjust (increase) q_value to reduce the false positives?

P/S: I also tried a different dataset (subset So instead of Ca). The coef and stdder are still the same, but null_hy is slightly difference between maaslin3 fit and constrast. Again, am I missing somethings here?

Thank you for your help.
Lenny

16

Hi Lenny,

The simple solution is the following:

  1. Fit MaAsLin 3 to the full data multiple times with the reference set to different groups on each run.
  2. Aggregate the comparisons you care about.
  3. FDR correct (p.adjust) over all the p-values of the comparisons you care about.

Regarding the actual question: Without the input files so I can run things myself, I’m still not sure why you’re seeing substantial differences in the p-values despite there being no difference in the coefficients, standard errors, and null hypotheses. We actually have a unit test checking whether the results from running the analysis in the two ways you’ve described are the same. This seems to work, so there might be something about your data that we haven’t considered. If I had to guess, my guess would be that when your p-values are as low as the ones you’re looking at, we might need to use more Monte Carlo draws to get the p-values closer. If you have a chunk of the data that reproduces this issue that I can run myself, I’d be happy to take a look.

Regarding the slightly different null hypothesis, the null hypothesis is chosen as the median of the coefficients across the features for each covariate. However, we do some filtering to make sure that coefficients that are likely wrong because of model misfits aren’t included in this calculation. It’s possible that the models that produced errors are slightly different in the original run and the contrast tests, so this might be why the null hypothesis are slightly different.

Will

17

Hi Will,
Thank you for your suggestions.
I tried to subset my input_data and compared between the results. It seemed to me that if I subset the input_data small enough, the differences between pvalue and qvalue (individual) would be minimised. Please see the attachment for the subset data for examples, which you can find the same issue (I tested with abundance model using maaslin3_fit ← maaslin3(input_data = input_data_maasslin3,input_metadata = input_metadata_test,output = ‘maaslin3’,formula = ‘~ Factor_2 + (1|Random_factor)’, reference = c(‘Factor_2,F2_1’), min_abundance = 0.05,cores = 1,small_random_effects=T, save_models = T)).
metadata_test.csv (1.7 KB)
input_data_test.csv (20.4 KB)

In the meantime, I will try to fit all the data to interaction model then change the reference group. It is difficult to grasp the correct comparison in interaction model output (~ Factor_1 + Facto_2 + Factor_1:Factor_2 + (1|Random_factor) with F1_1, and F2_1 as references. As I understand, the main factors in ‘name’ (F2_2, F2_3, F2_4) are comparison between them vs F2_1 only in the pool of F1_1 (ref of factor 1), but not across Factor_1. Similarly, F1_2 showed comparison between F1_2 vs F1_1 in the pool of F2_1 (ref of factor 2). However, in the interaction terms (for examples: F1_2:F2_2, F1_2:F2_3), they are comparison between F1_2 vs ref (F1_1) in the pool of F2_2, and F2_3, respectively. Are these the correct way to understand the interaction?

feature metadata value name coef null_hypothesis stderr pval_individual qval_individual pval_joint qval_joint error model N N_not_zero
92 Dorea longicatena Factor_1 F1_2 Factor_1F1_2 3.657541236 0.242062547 1.033260144 0.005915582 0.745363356 0.01179617 0.842246557 NA linear 46 30
93 Dorea longicatena Factor_1 F1_2:Factor_2F2_2 Factor_1F1_2:Factor_2F2_2 -4.844027341 0.92420314 1.633727736 0.004060544 0.745363356 0.0081046 0.842246557 NA linear 46 30
96 Dorea longicatena Factor_2 F2_2 Factor_2F2_2 3.819300532 -0.161349964 1.033260144 0.002382188 0.745363356 0.004758702 0.842246557 NA linear 46 30
98 Dorea longicatena Factor_2 F2_4 Factor_2F2_4 3.314044468 -0.162336157 0.980236641 0.005250315 0.745363356 0.010473064 0.842246557 NA linear 46 30
336 uncultured Gemmiger sp. Factor_2 F2_4 Factor_2F2_4 6.980533454 -0.162336157 2.050683751 0.003165618 0.745363356 0.006321216 0.842246557 NA linear 46 27
1 Agathobacter rectalis Factor_1 F1_2 Factor_1F1_2 1.089260192 0.242062547 1.637006757 0.609591992 1 0.582167147 1 NA linear 46 39
2 Agathobacter rectalis Factor_1 F1_2:Factor_2F2_2 Factor_1F1_2:Factor_2F2_2 -0.03530979 0.92420314 2.089074258 0.667829013 1 0.889662435 1 NA linear 46 39
3 Agathobacter rectalis Factor_1 F1_2:Factor_2F2_3 Factor_1F1_2:Factor_2F2_3 -0.768176935 -0.785596909 2.143344965 0.993580342 1 0.904816315 1 NA linear 46 39
4 Agathobacter rectalis Factor_1 F1_2:Factor_2F2_4 Factor_1F1_2:Factor_2F2_4 -0.813931695 -0.45620534 2.089074258 0.863397868 1 0.886099512 1 NA linear 46 39

Thank you for your help.
Lenny

18

Thanks for sending an example - I’ll look into it. Regarding the interaction model, if you want to compare every one of your 2*4=8 groups to each other, it’s probably easiest to turn them into an 8-group categorical variable and change the baseline each time. This will be equivalent to a fully interacted model. If you’re actually trying to do something else, let me know and I could advise on what to do instead.

1 Like