Large variance estimates in glmer

Question

I am running a piecewise glmer and the variance estimates we are getting for the random effects are very large:

Random effects:
 Groups Name        Variance  Std.Dev. Corr       
 id     (Intercept) 5.947e-02  0.2439             
        T1          2.813e+03 53.0372   0.09      
        T2          2.817e+03 53.0762  -0.07 -1.00
Number of obs: 196, groups:  id, 50

I have tried to narrow down the issue, and this only seems to occur when I include a specific covariate (variable c2). Without this covariate, the random effects seem more what I would expect for this model:

Random effects:
 Groups Name        Variance Std.Dev. Corr       
 id     (Intercept) 0.6838   0.8269              
        T1          0.3065   0.5536   -0.76      
        T2          9.6504   3.1065    0.25 -0.82
Number of obs: 247, groups:  id, 63

Keeping this covariate in is desirable as it is conceptually important and used in many similar models in existing literature.

This question was somewhat addressed in a similar question but it doesn't seem to be answered satisfactorily (in my opinion) as the answers pertain to the unit of measurement rather than the reason for such large variance estimates.

My model is:

model <-glmer(y ~ T1 + T2 +          
                 # predictors
                 x1 + x2 + x3 +
                 # covariates
                 c1 + 
                 c2 + # covariate affecting random effects
                 c3 +
                 c4 +
                 c5 +
                 (T1 + T2| id), 
               data = d, family="binomial" (link="logit"),               
               glmerControl(optimizer = "nlminbwrap"))

This issue is somewhat fixed by using nAGQ=0 in the model, though according to answers here it is not preferable to use this option if there is another way to fix the model. So, I would like to understand why this occurs, and only with the inclusion of a specific covariate.

The reproducible dataframe is here (values have been scaled otherwise warnings were thrown). I have been using lmerTest for analysis.

d <- read.table("https://pastebin.com/raw.php?i=mz4c7qBp")
colnames(d) <-c("id","y","T1","T2","x1","x2","x3","c1","c2","c3","c4","c5")

Any help is greatly appreciated!

Answer 1

I noticed that you are comparing two models fitted to (partially) different datasets: one has 196 observations and the other 247. You've helpfully provided a link to the data, so I took a look.

The summary

The random slopes in your models have high variance and high negative correlation. The unusual results are due to the unusual structure of your data. You should take the glmer warnings and results as a hint that the mixed-effects model you are considering (with or without the c2 covariate) might not be appropriate for your data.

The details

There are 63 IDs and 13 have missing c2 values. This explains why the number of observations is different depending on whether c2 is included or not. It's hard to compare models fitted to different data, so I'll drop those 13 IDs (about 20% of the data).

In this analysis there are 50 IDs and each ID has 4 observations. (There are 4 missing responses, so there are 196 observations instead of 50×4 = 200.)

The T1 and T2 variables are binary and for each ID we have one observation at (T1=0,T2=0), one observation at (T1=1,T2=0) and two observations at (T1=1,T2=1). The predictors x1,x2,x3 (numeric) and the covariates c1,...,c5 (one binary, four numeric) are ID-level variables: they don't vary within ID. Aside: What do you mean by the distinction between predictors and covariates?

d %>% count(T1,T2)
#> # A tibble: 3 × 3
#>      T1    T2     n
#>   <dbl> <dbl> <int>
#> 1     0     0    50
#> 2     1     0    50
#> 3     1     1   100

We now have some idea of your experiment: For each participant/item you recorded their x1,x2,x3 and c1,...,c5; then you measured the outcome under the three different settings for T1 & T2.

Let's write down the mixed-effects logistic regression as well:

$$ \begin{aligned} \operatorname{Y}_{i} &\sim \operatorname{Binomial}(n = 1, p_i) \\ \log\left(\frac{p_i}{1 - p_i} \right) &= \alpha_{j[i]} + \beta_{1j[i]}(\operatorname{T1}) + \beta_{2j[i]}(\operatorname{T2}) \end{aligned} $$ where $$ \begin{aligned} \operatorname{E}\alpha_j &= \gamma_{0} + \gamma_1\operatorname{x1}_j + \gamma_2\operatorname{x2}_j + \gamma_3\operatorname{x3}_j + \gamma_4\operatorname{c1}_j + \gamma_5\operatorname{c2}_j + \gamma_6\operatorname{c3}_j + \gamma_7\operatorname{c4}_j + \gamma_8\operatorname{c5}_j \end{aligned} $$

Since T1 and T2 are binary, there are three random intercepts: for ID $j$ the log odds are $\alpha_j$ when (T1=0,T2=0), $\alpha_j + \beta_{1i}$ when (T1=1,T2=0) and $\alpha_j + \beta_{1j} + \beta_{2j}$ when (T1=1,T2=1). The $\alpha_j$s are a function of the predictors $x_k$s and $c_k$s.

Now that we understand the data and the model, let's fit two logistic regressions: m0 without c2 and m1 with c2. There are warnings for both models, which I ignore; we already suspect that the models don't fit.

Here is the output about the random effect estimates. Notice the perfect correlations.

m0
#> Formula: 
#> y ~ T1 + T2 + x1 + x2 + x3 + c1 + c3 + c4 + c5 + (T1 + T2 | id)
#>
#> Random effects:
#>  Groups Name        Std.Dev. Corr       
#>  id     (Intercept) 0.3770              
#>         T1          0.7296   -1.00      
#>         T2          2.2582    1.00 -1.00
#> Number of obs: 196, groups:  id, 50
m1
#> Formula: 
#> y ~ T1 + T2 + x1 + x2 + x3 + c1 + c2 + c3 + c4 + c5 + (T1 + T2 | id)
#>
#> Random effects:
#>  Groups Name        Std.Dev. Corr

#>  id     (Intercept)  0.2439

#>         T1          53.0374   0.09

#>         T2          53.0764  -0.07 -1.00
#> Number of obs: 196, groups:  id, 50

To see what these perfect correlation mean, I plot the random effects $(\hat{\alpha}_j, \hat{\beta}_{1j}, \hat{\beta}_{2j})$ for the 50 ids. I've highlighted the three ids that drive the high variance we observed when c2 is included in the model.

What does this all mean? I don't know since I don't know anything about your experiment. It seems clear, however, that you should think more about what model to use to analyze the data.