# No variance in binomial glmer yields problematic CI

I have a glmer with a binomial response variable and two binary predictor variables, where there are no observations for one of the four groups. I used emmeans to calculate confidence intervals, and for this group obviously the mean and variance are zero. The model predicts a CI from 0 to 1, and my question is whether the confidence interval and contrasts are calculated correctly when all observations and thus variance of a group are zero. Please see details below. I have a total of 804 observations, so a bit much to include here. Should someone want to reproduce my models, I can share the data somehow.

The data concern observations of two female sharks in an aquarium. The response variable is presence/absence in a shipwreck (in_wrak_inc). Presence was recorded by multiple people (nearly twenty observers) every five minutes between 9:00 and 17:00 for many days (a few in early January 2021, and then intermettently from 2021-12-14 to 2022-01-19). As random effects, the model includes observer and date. As fixed effects, the model includes individual (either the first or second female shark, jv or ov), presence/absence of a male shark (locatie_jm, where Kleine Oceaan == absence) and there interaction.

  individu locatie_jm    in_wrak_inc     n
1 jv       Kleine Oceaan nee            52
2 jv       Kleine Oceaan ja             27
3 jv       Oceaan        nee           213
4 ov       Kleine Oceaan nee            83
5 ov       Kleine Oceaan ja            165
6 ov       Oceaan        nee           177
7 ov       Oceaan        ja             87


As you can see there are no observations for jv when locatie_jm == Oceaan + in_wrak_inc == ja. My code is below.

model_wrak1 <- lme4::glmer(in_wrak_inc ~ locatie_jm * individu + (1|datum) + (1|observant),
data = scans_vrouwtjes_wrak,
family = binomial)

summary(model_wrak1)

# Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
# Family: binomial  ( logit )
# Formula: in_wrak_inc ~ locatie_jm * individu + (1 | datum) + (1 | observant)
# Data: scans_vrouwtjes_wrak
#
# AIC      BIC   logLik deviance df.resid
# 406.6    434.7   -197.3    394.6      798
#
# Scaled residuals:
#   Min      1Q  Median      3Q     Max
# -1.7776 -0.0625  0.0000  0.0980  8.6087
#
# Random effects:
#   Groups    Name        Variance Std.Dev.
# datum     (Intercept) 26.32    5.130
# observant (Intercept) 32.81    5.728
# Number of obs: 804, groups:  datum, 27; observant, 16
#
# Fixed effects:
#   Estimate Std. Error z value Pr(>|z|)
# (Intercept)                   1.9335     2.7759   0.697 0.486091
# locatie_jmOceaan            -31.3131    66.9586  -0.468 0.640036
# individuov                    1.1188     0.3358   3.331 0.000864 ***
#   locatie_jmOceaan:individuov  25.2350    66.9589   0.377 0.706268
# ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Correlation of Fixed Effects:
#   (Intr) lct_jO indvdv
# locat_jmOcn -0.013
# individuov  -0.115  0.000
# lct_jmOcn:n -0.012 -0.999  0.000

model_wrak_emm_loc_ind <- emmeans::emmeans(model_wrak1,
specs = pairwise ~ locatie_jm | individu,
type = "response")

model_wrak_emm_loc_ind$emmeans # individu = jv: # locatie_jm prob SE df asymp.LCL asymp.UCL # Kleine Oceaan 0.8736 0.306 Inf 0.029111 0.999 # Oceaan 0.0000 0.000 Inf 0.000000 1.000 # # individu = ov: # locatie_jm prob SE df asymp.LCL asymp.UCL # Kleine Oceaan 0.9549 0.119 Inf 0.086899 1.000 # Oceaan 0.0463 0.123 Inf 0.000208 0.919 # # Confidence level used: 0.95 # Intervals are back-transformed from the logit scale model_wrak_emm_loc_ind$contrasts |> confint()

# individu = jv:
#   contrast               odds.ratio       SE  df asymp.LCL asymp.UCL
# Kleine Oceaan / Oceaan   3.97e+13 2.66e+15 Inf      0.00  3.93e+70
#
# individu = ov:
#   contrast               odds.ratio       SE  df asymp.LCL asymp.UCL
# Kleine Oceaan / Oceaan   4.36e+02 1.34e+03 Inf      1.05  1.82e+05
#
# Confidence level used: 0.95
# Intervals are back-transformed from the log odds ratio scale


• I don't really understand why you consider the CI (0, 1), problematic. It is the most conservative CI possible (with guaranteed coverage), which seems appropriate if there is no information to derive a more narrow CI. Commented Oct 27, 2022 at 5:27
• If there are many observations and all are negative, surely that is an indication it is somewhat likely the CI should be closer to zero than to one? In this case, the probability of the other group is 0.87, and I expect a difference in the contrasts. Commented Oct 31, 2022 at 9:08
• Sure, but you would need to quantify "somewhat likely" to derive a more narrow CI. And you can't do that just based on the model. Commented Nov 1, 2022 at 6:01
• I understand the model cannot do that (because there is no variance?). Can it be done based on the number of observations? I have 213 observations where in_wrak_inc == 0. Surely the likelihood of in_wrak_inc == 1 is therefore lower than when I would have, say, only 15 observations. Commented Nov 2, 2022 at 9:19

You may find this paper helpful.

This problem, sparse data in one or more cells of the cross-tabulation of your predictor and outcome variables (as you present above) which causes extreme parameter estimates, is referred to as separation. The solution is usually one or another form of regularization

(Alternatively, you can omit the covariate(s), which cause the separation, but I'll ignore that here).

There are several such options discussed in the linked paper. Potentially the easiest to implement would be data augmentation. Another option would be to estimate the model in the Bayesian framework. Using even modestly informative priors on the log odds ratios would provide the regularization needed to get sane results (and narrower intervals).

You could try estimating the model with R package brms. It uses the same syntax as lme4 so you wouldn't really need modify your model specification, you'd just need to add priors for the model parameters.

• Thanks. That paper, and looking at separation in logistic regression generally, seems very useful. I'll check it out. Commented Oct 31, 2022 at 9:09

Trying to answer my own question. The lme4 package cannot deal with complete separation, but info how to fix it can be found here and here (based on a question posted here). Running a mixed model with priors with blme::bglmer() did the trick. The priors I chose (10 for the intercept and 2.5 for the various slopes) are default.

model_wrak1_bayes <- blme::bglmer(in_wrak_inc ~ locatie_jm * individu + (1|datum) + (1|observant),
data = scans_vrouwtjes_wrak,
family = binomial,
fixef.prior = normal(sd = c(10, 2.5)))


This gives the following emmeans, which I find very reasonable:

individu = jv:
locatie_jm        prob       SE  df asymp.LCL asymp.UCL
Kleine Oceaan 0.092783 0.216260 Inf  6.65e-04    0.9402
Oceaan        0.000234 0.000542 Inf  2.49e-06    0.0215

individu = ov:
locatie_jm        prob       SE  df asymp.LCL asymp.UCL
Kleine Oceaan 0.269436 0.505533 Inf  2.40e-03    0.9826
Oceaan        0.093029 0.198690 Inf  1.01e-03    0.9120

Confidence level used: 0.95
Intervals are back-transformed from the logit scale