I am a biologist how is working with relative abundance data of the microbiome derived from 16S rRNA sequencing. I'm comparing the changes of abundance from different bacterial taxa in people who undergo a dietary intervention. The data is derived from two different cohorts, one being normal weight and the other obese subjects. In each cohort subjects can further be stratified into high and low conversion types. The research questions are:
- 1.) Is the abundance of taxa changing on both dietary interventions in the same way (Pre vs Post)?
- 2.) Is the abundance of taxa depending on the conversion type (high vs low vs not classified)?
- 3.) Does the conversion type affect the dietary response?
All questions should be tested if they are independent on the cohort.
As I have already analysed both cohorts separately, I know that both cohort show the same taxa alterations, but I need to proof this in a combined approach as well, which controls for the cohorts. Therefore, I wanted to use Linear Mixed Models in R using
As the more I dig into these models, they appear to be more complicated and my working group and I are very new to these models, so I find myself overwhelmed and came here for advise. After thorough self-education (mainly http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html) I came up with the following:
I want to determine the variance of the dietary effect (Time) and the conversion type (converter). Both are categorial with Time having two levels (PRE vs POST) and converter having three levels (low vs high vs not classified). Both of them are my fixed effects. For my random effects I added the two cohorts (Study) and Subject (Proband) as I have repeated measurements of the individual.
The two fixed effects are fully crossed with every conversion type being present at each time point, also in each cohort. But the subjects are nested with the cohort.
Since I'm using relative abundance data, the data was clr transformed before using the LMM.
The data looks like this:
df taxon Proband Study Time converter 1 8.337073e-02 1 CFS PRE high 2 3.556427e-02 1 CFS POST high 3 -3.065652e-03 2 CFS PRE high 4 -3.076171e-03 2 CFS POST nc 5 3.001781e-03 3 CFS PRE nc 6 1.363824e-03 3 CFS POST high 7 2.986941e-03 4 CFS PRE high 8 3.264113e-03 4 CFS POST high 9 1.458944e-01 5 CFS PRE high 10 5.375029e-03 5 CFS POST high 11 1.616321e-02 6 KD PRE high 12 5.263631e-02 6 KD POST high 13 9.498285e-02 7 CFS PRE high 14 1.290295e-01 8 CFS PRE high 15 5.002668e-02 10 CFS PRE high 16 -3.064218e-03 11 CFS PRE high 17 4.524587e-02 12 KD PRE low 18 6.753473e-02 12 KD POST low ...
So I was first thinking about using the following model:
model1 <- lmerTest::lmer(data = df, taxon ~ Time + converter + Time*converter+ (1|Proband/Study) + (1|Study))
But this resulted in the error messages:
boundary (singular) fit: see ?isSingular Model failed to converge with 1 negative eigenvalue: -1.1e-03
I think that this might be due to the low levels of the Study effect, with only two. I read that a random effect should at least contain 5-6 levels for optimal performance.
So I thought using the cohorts (Study) as either fixed effect (model2) or only as nested effect of the Subject (model3)
model2 <- lmerTest::lmer(data = df, taxon ~ Time + converter + Study + Time:converter+ (1|Proband/Study)) model3 <- lmerTest::lmer(data = df, taxon ~ Time + converter + Time*converter+ (1|Proband/Study))
Again, model2 resulted in the error that the model failed to converge. But both results of model2 and model 3 are very similar.
summary(model2) REML criterion at convergence: -935.6 Scaled residuals: Min 1Q Median 3Q Max -1.2881 -0.5684 -0.2194 0.3822 4.2572 Random effects: Groups Name Variance Std.Dev. Study:Proband (Intercept) 0.0002481 0.01575 Proband (Intercept) 0.0003285 0.01813 Residual 0.0013602 0.03688 Number of obs: 286, groups: Study:Proband, 176; Proband, 176 Fixed effects: Estimate Std. Error df t value Pr(>|t|) (Intercept) 0.047016 0.004455 269.333078 10.553 < 2e-16 *** TimePOST -0.019315 0.005881 152.070853 -3.284 0.00127 ** converterlow 0.015898 0.008302 278.684763 1.915 0.05652 . converternc 0.005332 0.008789 272.074328 0.607 0.54461 StudyKD -0.006141 0.007816 154.956926 -0.786 0.43325 TimePOST:converterlow 0.023732 0.014015 178.192218 1.693 0.09213 . TimePOST:converternc -0.009934 0.012860 234.676008 -0.772 0.44063 summary(model3) REML criterion at convergence: -942.8 Scaled residuals: Min 1Q Median 3Q Max -1.2626 -0.5612 -0.2137 0.3981 4.2835 Random effects: Groups Name Variance Std.Dev. Study:Proband (Intercept) 3.122e-05 0.005587 Proband (Intercept) 5.396e-04 0.023230 Study (Intercept) 0.000e+00 0.000000 Residual 1.362e-03 0.036908 Number of obs: 286, groups: Study:Proband, 176; Proband, 176; Study, 2 Fixed effects: Estimate Std. Error df t value Pr(>|t|) (Intercept) 0.046211 0.004329 277.001787 10.675 < 2e-16 *** TimePOST -0.019753 0.005858 153.624537 -3.372 0.000944 *** converterlow 0.015647 0.008292 279.734337 1.887 0.060190 . converternc 0.005276 0.008785 272.890143 0.601 0.548624 TimePOST:converterlow 0.024858 0.013936 181.891132 1.784 0.076142 . TimePOST:converternc -0.010498 0.012837 236.197672 -0.818 0.414291
Which model is the right for my approach? Is the error of the model failed to converge critical? I'm unsure as it does not seem to change the results much. I rather tend to use model 3 as it does not result in error messages but I'm not sure if the cohort effect is accounted for in model3 correctly if it is only added as nested effect of Proband.
I hope that anyone with LMM experience can give me a recommendation. Thank you in advance!!