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 lme4::lmer().
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!!
