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

Question

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)
## [1] 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.

Answer 1

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.

Answer 2

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 ...