Let's say we have a model
mod <- Y ~ X*Condition + (X*Condition|subject) # Y = logit variable # X = continuous variable # Condition = values A and B, dummy coded; the design is repeated # so all participants go through both Conditions # subject = random effects for different subjects summary(model) Random effects: Groups Name Variance Std.Dev. Corr subject (Intercept) 0.85052 0.9222 X 0.08427 0.2903 -1.00 ConditionB 0.54367 0.7373 -0.37 0.37 X:ConditionB 0.14812 0.3849 0.26 -0.26 -0.56 Number of obs: 39401, groups: subject, 219 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) 2.49686 0.06909 36.14 < 2e-16 *** X -1.03854 0.03812 -27.24 < 2e-16 *** ConditionB -0.19707 0.06382 -3.09 0.00202 ** X:ConditionB 0.22809 0.05356 4.26 2.06e-05 ***
Here we observe a singular fit, because the correlation between intercept and x random effects is -1. Now, according to this helpful link one way to deal with this model is to remove higher-order random effects (e.g., X:ConditionB) and see whether that makes a difference when testing for singularity. The other is to use the Bayesian approach, e.g., the
blme package to avoid singularity.
What is the preffered method and why?
I am asking this because using the first or the second one leads to different results - in the first case, I will remove the X:ConditionB random effect and won't be able to estimate the correlation between X and X:ConditionB random effects. On the other hand, using
blme allows me to keep X:ConditionB and to estimate the given correlation. I see no reason why I should even use the non-bayesian estimations and remove random effects when singular fits occur when I can estimate everything with the Bayesian approach.
Can someone explain to me the benefits and problems using either method to deal with singular fits?