lme4 models with random effects with suppressed intercepts, (0 + f | g), fail

[Ax3man]

Hi, this package is a really good idea and already very useful.

Just wanted to report that I can't get any lme4 models that have random coefficients (other than intercept) to work. Examples:

library(equatiomatic)
library(lme4)
mtcars2 <- mtcars
mtcars2$carb <- factor(mtcars$carb)
mtcars2$am <- factor(mtcars$am)
extract_eq(lmer(hp ~ 1 + (1 | carb), mtcars2)) 
# works
extract_eq(lmer(hp ~ 1 + (1 + am | carb), mtcars2)) 
# Error in `[<-`(`*tmp*`, v[1], v[2], value = greek_vcov) : 
#  subscript out of bounds
extract_eq(lmer(hp ~ 1 + (0 + am | carb), mtcars2)) 
# Error in dimnames(mat) <- list(unique(unlist(rand_lev$vcov_greek[sd_rows])),  : 
#  length of 'dimnames' [1] not equal to array extent
extract_eq(lmer(hp ~ 1 + (1 + qsec | carb), mtcars2))
# Error in `[<-`(`*tmp*`, v[1], v[2], value = greek_vcov) : 
#  subscript out of bounds

I'm running a fresh github install, v0.2.0.9000.

[datalorax]

Thanks for the report - I believe in this case it's because you are specifying a model with a random effect that does not have a corresponding fixed effect. Extending your example you can see they work if you add in the fixed effect.

library(equatiomatic)
library(lme4)
#> Loading required package: Matrix
mtcars2 <- mtcars
mtcars2$carb <- factor(mtcars$carb)
mtcars2$am <- factor(mtcars$am)

# this fails
# extract_eq(lmer(hp ~ 1 + (1 + am | carb), mtcars2))

# But this works
extract_eq(lmer(hp ~ 1 + am + (1 + am | carb), mtcars2))
#> boundary (singular) fit: see ?isSingular

# This fails
# extract_eq(lmer(hp ~ 1 + (1 + qsec | carb), mtcars2))

# But this works
extract_eq(lmer(hp ~ 1 + qsec + (1 + qsec | carb), mtcars2))
#> boundary (singular) fit: see ?isSingular

^{Created on 2021-02-18 by the reprex package (v0.3.0)}

[Ax3man]

Wonderful, my actual model does have a corresponding fixed effect so this is very helpful. I really appreciate the rapid response!

The following model, including the suppression of the intercept in the random coefficients, still fails however. Consider these two examples:

library(equatiomatic)
library(lme4)
mtcars2 <- mtcars
mtcars2$am <- factor(mtcars$am)

extract_eq(lmer(hp ~ am + (1 + am | carb), mtcars2)) 
# works
extract_eq(lmer(hp ~ am + (0 + am | carb), mtcars2)) 
# still errors

Note that am needs to be coded as a factor here, because the correlation between the two levels should be modeled.

[datalorax]

Okay... yeah that is not a model type that I ever tested so it doesn't surprise me that it's failing. I will work on that in the next little bit.

As an aside - are there situations where a model like lmer(hp ~ 1 + (1 + am | carb), mtcars2) makes sense? I've never encountered a model like that and it doesn't make intuitive sense to me, so I'd be interested if you could share an application where a model like that is preferred.

[Ax3man]

Yeah, I wrote it that way just to give the most minimal example. My actual model has a lot more parameters, but the problem was with random effects so I just removed all fixed effects. I can't give you a definitive answer on whether people use such a model, but I haven't.

The (0 + f | g) construction is useful, however.

[datalorax]

Great. Yes, the lmer(hp ~ am + (0 + am | carb), mtcars2) model definitely makes sense. I'll work on trying to get that solved but not worry about the other. Thanks again for the report!

[datalorax]

So just to be totally clear, in this example

mtcars2 <- mtcars
mtcars2$carb <- factor(mtcars$carb)
mtcars2$am <- factor(mtcars$am)

equatiomatic::extract_eq( lme4::lmer(hp ~ am + (1 + am | carb), mtcars2))

Would you expect it to look like this?

Or something different? I guess I haven't really thought through these types of models.

[Ax3man]

I actually don't know, sorry! I'm quite comfortable with the lme4 syntax and what coefficients are modeled, but not at all with the mathematical notation. I was hoping to use equatiomatic to give me the correct notation.

[datalorax]

Okay no worries, I think I actually got it figured out. I played with a model with more levels and it made more sense. Still need to implement it but I at least think I know what needs to be implemented now.

[datalorax]

Okay so I'm going to leave this for now, but I did play around with it a bit I'm not exactly sure what the best way forward is. The issue is that with typical models all your variance components have a corresponding fixed effect. When you have dummy-coded variables, for example, you leave one level out, but that goes into your intercept. If you allow the intercept to randomly vary, then all is good. But in this case you have a fixed effect for your intercept, which you're not allowing to vary randomly, but you also have a random effect for your slope. Because you've dropped the intercept, you end up with variance components for all levels of the categorical variable, which I suspect is what you were after, but that also means you don't have a corresponding fixed effect to refer back to.

I'll illustrate with an example. Here's a standard model, and all is good.

mtcars2 <- mtcars
mtcars2$carb <- factor(mtcars$carb)
mtcars2$am <- factor(mtcars$am)

m1 <- lme4::lmer(hp ~ 1 + carb + am +  (1 + carb | gear), mtcars2)
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
#> unable to evaluate scaled gradient
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
#> Model failed to converge: degenerate Hessian with 2 negative eigenvalues
equatiomatic::extract_eq(m1)
#> Registered S3 method overwritten by 'broom.mixed':
#>   method      from 
#>   tidy.gamlss broom
#> NULL

We can also drop the intercept from both the fixed effect and random effects, and we're all good, because now all levels of the categorical variable are estimated.

m2 <- lme4::lmer(hp ~ 0 + carb + am +  (0 + carb | gear), mtcars2)
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
#> unable to evaluate scaled gradient
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
#> Model failed to converge: degenerate Hessian with 4 negative eigenvalues
equatiomatic::extract_eq(m2)

Although the subscripts should probably go 1 to n instead of 0 to n, but this is minor.

The problem comes when we're estimating the variance components for all the level of the categorical variable, but the fixed effect for the first level is part of the intercept.

m3 <- lme4::lmer(hp ~ 1 + carb + am +  (0 + carb | gear), mtcars2)
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
#> unable to evaluate scaled gradient
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
#> Model failed to converge: degenerate Hessian with 2 negative eigenvalues
broom.mixed::tidy(m3)
#> Registered S3 method overwritten by 'broom.mixed':
#>   method      from 
#>   tidy.gamlss broom
#> # A tibble: 29 x 6
#>    effect   group term             estimate std.error statistic
#>    <chr>    <chr> <chr>               <dbl>     <dbl>     <dbl>
#>  1 fixed    <NA>  (Intercept)        93.7       10.9      8.59 
#>  2 fixed    <NA>  carb2              25.1       18.3      1.37 
#>  3 fixed    <NA>  carb3             105.        56.2      1.87 
#>  4 fixed    <NA>  carb4             112.        39.1      2.86 
#>  5 fixed    <NA>  carb6              91.2       20.8      4.39 
#>  6 fixed    <NA>  carb8             211.        67.4      3.13 
#>  7 fixed    <NA>  am1                -7.33       8.57    -0.855
#>  8 ran_pars gear  sd__carb1          12.8       NA       NA    
#>  9 ran_pars gear  cor__carb1.carb2    0.919     NA       NA    
#> 10 ran_pars gear  cor__carb1.carb3   -0.767     NA       NA    
#> # … with 19 more rows

In Line 8 above you can see that there is a variance component for carb1, but the corresponding fixed effect for this is really part of the intercept.

So I'm not sure how the notation for this should work. The variance for carb1 is not the same as intercept variance, but it doesn't have a fixed effect to refer to.

As I mentioned, I'll keep this open for now, but unless there's a clear direction forward, this may just be an edge case that we end up not supporting.

[Ax3man]

I see, this really clarifies the problem, thanks. I think it would make sense to just use m2 instead of m3 in most cases.

My actual model has a random effects structure something like (0 + f | g1) + (1 | g2), so dropping the fixed intercept doesn't solve the problem, but I now understand that it's not really clear what the notation for such a model should look like.

Thanks again for all your help.

[datalorax]

I appreciate you raising the issue. I think I will close this for now because the direction forward is not clear to me. If it becomes a more common problem then perhaps we can figure out a standard way to communicate these types of models clearly.