# How to deal with quasi-complete separation in a logistic GLMM?

Update: Since I now know that my problem is called quasi-complete separation I updated the question to reflect this (thanks to Aaron).

I have a dataset from an experiment in which 29 human participants (factor code) worked on a set of trials and the response was either 1 or 0. In addition, we manipulated the materials so that we had three crossed factors, p.validity (valid versus invalid), type (affirmation versus denial), and counterexamples (few versus many):

d.binom <- read.table("http://pastebin.com/raw.php?i=0yDpEri8")
str(d.binom)
## 'data.frame':   464 obs. of  5 variables:
##      $code : Factor w/ 29 levels "A04C","A14G",..: 1 1 1 1 1 1 1 1 1 1 ... ##$ response       : int  1 1 1 1 0 1 1 1 1 1 ...
##      $counterexamples: Factor w/ 2 levels "few","many": 2 2 1 1 2 2 2 2 1 1 ... ##$ type           : Factor w/ 2 levels "affirmation",..: 1 2 1 2 1 2 1 2 1 2 ...
##      $p.validity : Factor w/ 2 levels "invalid","valid": 1 1 2 2 1 1 2 2 1 1 ...  Overall there is only a small number of 0s: mean(d.binom$response)
##  0.9504


One hypothesis is that there is an effect of validity, however, preliminary analysis suggests there might be an effect of counterexamples. As I have dependent data (each participant worked on all trials) I would like to use a GLMM on the data. Unfortunately, counterexamplesquasi-completely separate the data (at least for one level):

with(d.binom, table(response, counterexamples))
##         counterexamples
## response few many
##        0   1   22
##        1 231  210


This is also reflected in the model:

require(lme4)
options(contrasts=c('contr.sum', 'contr.poly'))

m2 <- glmer(response ~ type * p.validity * counterexamples + (1|code),
data = d.binom, family = binomial)
summary(m2)
## [output truncated]
## Fixed effects:
##                                      Estimate Std. Error z value Pr(>|z|)
##   (Intercept)                            9.42     831.02    0.01     0.99
##   type1                                 -1.97     831.02    0.00     1.00
##   p.validity1                            1.78     831.02    0.00     1.00
##   counterexamples1                       7.02     831.02    0.01     0.99
##   type1:p.validity1                      1.97     831.02    0.00     1.00
##   type1:counterexamples1                -2.16     831.02    0.00     1.00
##   p.validity1:counterexamples1           2.35     831.02    0.00     1.00
##   type1:p.validity1:counterexamples1     2.16     831.02    0.00     1.00


The standard errors for the parameters are simply insane. As my final goal is to assess whether or not certain effects are significant, standard errors are not totally unimportant.

• How can I deal with the quasi complete separation? What I want is to obtain estimates from which I can judge whether or not a certain effect is significant or not (e.g., using PRmodcomp from package pkrtest, but this is another step not described here).

Approaches using other packages are fine as well.

• To get you started, try this: ats.ucla.edu/stat/mult_pkg/faq/general/… – Aaron left Stack Overflow Oct 2 '12 at 20:16
• @Aaron Looks great. Putting this into an answer would at least have brought you one upvote... – Henrik Oct 2 '12 at 20:25
• Not really an answer though, but thanks! – Aaron left Stack Overflow Oct 2 '12 at 20:29
• @Henrik You can upvote comments too. – Peter Flom Oct 2 '12 at 21:10
• See this paper by Paul Allison. Although he emphasizes SAS, the same points will be applicable in other languages. – Peter Flom Oct 2 '12 at 21:11

I am afraid there's a typo in your title: you should not attempt to fit mixed models, let alone nonlinear mixed models, with just 30 clusters. Not unless you believe you can fit a normal distribution to 30 points obstructed by measurement error, nonlinearities, and nearly complete separation (aka perfect prediction).

What I would do here is to run this as a regular logistic regression with Firth's correction:

library(logistf)
mf <- logistf(response ~ type * p.validity * counterexamples + as.factor(code),
data=d.binom)


Firth's correction consists of adding a penalty to the likelihood, and is a form of shrinkage. In Bayesian terms, the resulting estimates are the posterior modes of the model with a Jeffreys prior. In frequentist terms, the penalty is the determinant of the information matrix corresponding to a single observation, and hence disappears asymptotically.

• Actually I do believe in fitting mixed models with less than 30 clusters. But the analysis seems nevertheless promising (+1). Or is there Firth's method for GLMMs? – Henrik Oct 2 '12 at 22:13
• Right, you are the enthusiast of the minimal sample size requirements... Firth's correction works with i.i.d. data only. You can believe in anything, but you'd be better off running some simulations to see if any given belief is justified in a given data situation. With a perfectly balanced data set and continuous reponse, it MAY work OK. With a heavily imbalanced data set, in terms of response, you see only a far left tail of the normal distribution of the random effects, and you are willing to bet on this tail being well approximated by one-point Laplace in *lmer??? :-\ – StasK Oct 3 '12 at 13:58

You can use a Bayesian maximum a posteriori approach with a weak prior on the fixed effects to get approximately the same effect. In particular, the blme package for R (which is a thin wrapper around the lme4 package) does this, if you specify priors for the fixed effects as in the example here (search for "complete separation"):

cmod_blme_L2 <- bglmer(predation~ttt+(1|block),data=newdat,
family=binomial,
fixef.prior = normal(cov = diag(9,4)))


This example is from an experiment where ttt is a categorical fixed effect with 4 levels, so the $$\beta$$ vector will have length 4. The specified prior variance-covariance matrix is $$\Sigma = 9 I$$, i.e. the fixed effect parameters have independent $$N(\mu=0,\sigma^2=9)$$ (or $$\sigma$$, i.e. standard devation, $$=3$$) priors. This works pretty well, although it's not identical to Firth correction (since Firth corresponds to a Jeffreys prior, which is not quite the same).

The linked example shows you can also do it with the MCMCglmm package, if you want to go full-Bayesian ...