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:

enter image description here

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:


      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 +
      (1|Proband/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 +

model3 <- lmerTest::lmer(data = df, 
        taxon ~ 
          Time + 
          converter +

Again, model2 resulted in the error that the model failed to converge. But both results of model2 and model 3 are very similar.


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 


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!!


1 Answer 1


Regarding model1:

taxon ~ Time * converter + (1|Proband/Study) + (1|Study)

..note that (1|Proband/Study) means that Study is nested within Proband which is not the case according to the description and the diagram. Also note that (1|Study) does not makes sense because there are only 2 levels of Study. That may explain the singular fit.

Regarding model2:

 taxon ~ Time * converter + Study + (1|Proband/Study)

...this also has the misspecified (1|Proband/Study), but addiionally fits Study as a fixed effect. This may explain the convergence problem.

Regarding model3:

taxon ~ Time * converter + (1|Proband/Study)

...something is not right with the model formula or the output, since the random effect term (1|Proband/Study) is not reflected correctly in the random effects estimates - the output shows Study:Proband, Study and Proband

 Study:Proband (Intercept) 3.122e-05 0.005587
 Proband       (Intercept) 5.396e-04 0.023230
 Study         (Intercept) 0.000e+00 0.000000

whereas (1|Proband/Study) should only resuly in estimates for Study:Proband and Proband.

I would suggest the following model:

taxon ~ Time * converter + Study + (1|Proband)
  • 1
    $\begingroup$ Thank you Robert for your fast and thorough answer! I tried the suggested model and it worked well. The singular fit problem is solved. And I totally see the syntax faults and now can understand why your suggestion fits. Learning from mistakes is always best! Thanks for helping me discussing the models! $\endgroup$
    – AMF
    Commented Jun 18, 2021 at 6:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.