Hello,
I am attempting to construct a Bernoulli mixed-effects model and a binomial mixed-effects model for a dataset using various R packages, including glmmTMB, glmmML, and lme4::glmer().
However, I am skeptical about the output shown by the R packages. Sometimes, I receive a warning about a singular fit or model that is not converging. However, at other times, even when I am not receiving any warnings, I still encounter them when diagnosing a model, for example, using glmmTMB::diagnose().
Here are some details regarding the Problem Description, model formulation, and experimental design setup.
Experimental Setup: https://www.dropbox.com/scl/fi/f4q6r4tnkyuw7rzop8puq/Design_Tree_Screenshot-2025-12-05-134907.png?rlkey=u7fe7yyfmzyq0azr5am1haf5w&st=1art9674&dl=0
Problem Description: https://www.dropbox.com/scl/fi/yjxo0ggf4opymeodcy43e/Exp_Design_Screenshot-2025-12-05-134839.png?rlkey=38zvrh1rwmibk3zq81dwv53m7&st=anppw5xu&dl=0
Binary Response Variable, Factor-1, Factor-2, and a Continuous Variable (x). Factor-2 Levels nested within levels of Factor-1. Unbalanced Design. A New Derived variable is used by combining Factor1 with Factor2 and is represented as "Factor1_Factor2". To represent nested effect with random intercept, I can either do this (1 | Factor1/Factor2) or as this (1 | Factor1_Factor2). The responses y_{i,j,k} as seen in the Tree Diagram above are repeated measures. The variance term for the Random effects is expected to be separate for each treatment effect. And, "x_scaled" is the scaled version of "x", i.e. x_scaled = (x - mean(x))/sd(x).
Bernoulli GLMM Model Formulation: https://www.dropbox.com/scl/fi/e0qu3wpfysjeagldzxs5w/Model_formulation_Bernoulli_Screenshot-2025-12-05-135101.png?rlkey=r516o40eia3zu78hz7tkmblgo&st=0qbhtcye&dl=0
Binomial GLMM Model Formulation
https://www.dropbox.com/scl/fi/3aqcq22q030zcroagq1x4/Binomial_GLMM_Screenshot-2025-12-05-140353.png?rlkey=ofcyz11l74eul38g410opdlsx&st=3b2f8n99&dl=0
Initially, I attempted modeling using the lme4::glmer() function. It was for Bernoulli GLMM. But, lme4 produces common variance for all random effect terms. Later, I switched to the glmmTMB package; I also attempted the glmmML package recently, but things did not work out as expected.
****************Bernoulli GLMM using lme4::glmer() : ******************
lme4::glmer(formula = bernoulli_y ~ Factor_1 + (1 | Factor_1/Factor_2), data = mydata, family = binomial(link = "logit")) # For this, I got Boundary (Singular) Fit.
lme4::glmer(formula = bernoulli_y ~ Factor_1 + (1 | Factor_1:Factor_2), data = mydata, family = binomial(link = "logit")) # For this, I got Boundary (Singular) Fit.
lme4::glmer(formula = bernoulli_y ~ (1 | Factor_1_Factor_2), data = mydata, family = binomial(link = "logit")) # For this, no fixed Effect; it came with no error or warning. (Factor1_Factor2) It is a derived factor variable.
me4::glmer(formula = bernoulli_y ~ x + (1 | Factor_1_Factor_2), data = mydata, family = binomial(link = "logit")) # Here, 'x' is a continuous variable; # Got No error or warnings.
*********** Bernoulli GLMM using glmmTMB::glmmTMB(): This package was used to get separate variances for random effects for ith treatment effects (i.e. \sigma^2_i, where I = 1,2,3,4,5). ****************
summary(glmmTMB::glmmTMB(formula = bernoulli_y ~ factor1 + diag(0+factor1|Factor1_Factor2), family=binomial(link = "logit"),data=mydata)) # This came without error/warning. However, when I diagnosed the model using glmmTMB::diagnose(), I received numerous warnings about unusually high coefficient values.
********* Binomial GLMM using lme4::glmer() : ********************
summary(lme4::glmer(formula = cbind(success,total-success) ~ Factor1+ (1 | Factor1/Factor2),data = mydata,family = binomial(link = "logit"))) # Singular Fit
summary(lme4::glmer(formula = cbind(success,total-success) ~ (1 | Factor1_Factor2), data = mydata, family = binomial(link = "logit"))) # No fixed effects, but worked without error/warning.
************** Binomial GLMM using glmmTMB(): ********************
glmmTMB::glmmTMB(formula = cbind(success,total-success) ~ factor_1 + (1|Factor1_Factor2), family=binomial(link = "logit"),data=mydata) # Initially, this model appeared to run successfully without any error/warning. For this case, I got a common variance. However, upon further inspection, I found that glmmTMB:check_singularity(model) returned TRUE, and glmmTMB::diagnose(model) also issued some warnings.
Another attempt using glmmTMB::glmmTMB was made to obtain the variances of random effects corresponding to each level of the Treatment effect (i.e., Factor 1).
glmmTMB::glmmTMB(formula = cbind(success,total-success) ~ Factor1 + diag(0+Factor1| Factor1_Factor2), family=binomial(link = "logit"),data=mydata) # Initially, this model appeared to run successfully without any error/warning. For this case, I got a common variance. However, upon further inspection, I found that glmmTMB:check_singularity(model) returned TRUE, and glmmTMB::diagnose(model) also issued some warnings.
************** Using glmmML package() for binomial GLMM *******************
summary(glmmML::glmmML(formula = cbind(success,total-success) ~ (1|Factor1/Factor2), family=binomial(link = "logit"),data=mydata,cluster = factor1_factor2)) # This resulted in error/warning.
summary(glmmML::glmmboot(formula = cbind(success,total-success) ~ (1|Factor1/Factor2), family=binomial(link = "logit"),data=mydata,cluster = factor1_factor2)) # This resulted in error/warning.
************** Using glmmML package() for bernoulli GLMM *******************
# warning/error
glmmML::glmmML(formula = bernoulli_y ~ x_scaled + Factor_1, family=binomial(link = "logit"),data=mydata, cluster = Factor1_Factor2)
# warning/error
glmmML::glmmML(formula = bernoulli_y ~ x_scaled, family=binomial(link = "logit"),data=mydata, cluster = Factor1_Factor2)
Anyway, if possible, could you please suggest the type of modelling approach that should be undertaken for this modelling formulation? Which of the above R packages perform as per the modelling requirements in this context? Or, will it be better to write explicitly a code for computing the Marginal log-likelihood of the Mixed Effects Model, and maximizing it using optim() or Laplace approximation?
Please provide any helpful suggestions in this regard.
????
With Best Regards,
Sami
(Pronouns: He/Him/His)
Iowa State University,
Ames, Iowa,
USA.