Comment gérer la séparation quasi complète dans un GLMM logistique?

Mise à jour : Puisque je sais maintenant que mon problème est appelé séparation quasi-complète, j'ai mis à jour la question pour refléter cela (merci à Aaron).

J'ai un ensemble de données provenant d'une expérience dans laquelle 29 participants humains (facteur code) ont travaillé sur un ensemble d'essais et la valeur responseétait soit 1 soit 0. De plus, nous avons manipulé le matériel de manière à avoir trois facteurs croisés p.validity(valides versus invalides), type(affirmation contre déni), et counterexamples(peu contre beaucoup):

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

Dans l'ensemble, il n'y a qu'un petit nombre de 0:

mean(d.binom$response)
## [1] 0.9504

Une hypothèse est qu'il existe un effet de validity, cependant, une analyse préliminaire suggère qu'il pourrait y avoir un effet de counterexamples. Comme j'ai des données dépendantes (chaque participant a travaillé sur tous les essais), je voudrais utiliser un GLMM sur les données. Malheureusement, counterexamplesséparez presque complètement les données (au moins pour un niveau):

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

Cela se reflète également dans le modèle:

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

Les erreurs standard pour les paramètres sont tout simplement folles. Étant donné que mon objectif final est d'évaluer si certains effets sont importants ou non, les erreurs types ne sont pas totalement sans importance.

• Comment gérer la séparation quasi complète? Ce que je veux, c'est obtenir des estimations à partir desquelles je peux juger si oui ou non un certain effet est significatif ou non (par exemple, en utilisant le PRmodcomppackage pkrtest, mais c'est une autre étape non décrite ici).

Les approches utilisant d'autres packages sont également très bien.

J'ai bien peur qu'il y ait une faute de frappe dans votre titre: vous ne devriez pas essayer d'ajuster des modèles mixtes, et encore moins des modèles mixtes non linéaires, avec seulement 30 clusters. Sauf si vous pensez pouvoir ajuster une distribution normale à 30 points obstrués par une erreur de mesure, des non-linéarités et une séparation presque complète (aka prédiction parfaite).

Ce que je ferais ici, c'est d'exécuter cela comme une régression logistique régulière avec la correction de Firth :

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

La correction de Firth consiste à ajouter une pénalité à la probabilité et est une forme de rétrécissement. En termes bayésiens, les estimations résultantes sont les modes postérieurs du modèle avec un a priori de Jeffreys. En termes fréquentistes, la pénalité est le déterminant de la matrice d'information correspondant à une seule observation, et disparaît donc asymptotiquement.

Vous pouvez utiliser une approche bayésienne maximale a posteriori avec un faible a priori sur les effets fixes pour obtenir approximativement le même effet. En particulier, le paquet blme pour R (qui est une enveloppe mince autour du lme4paquet) fait cela, si vous spécifiez des priorités pour les effets fixes comme dans l'exemple ici (recherchez "séparation complète"):

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

Cet exemple provient d'une expérience où tttest un effet fixe catégorique avec 4 niveaux, donc le$\beta$ vecteur aura la longueur 4. La matrice de variance-covariance antérieure spécifiée est $\Sigma = 9 I$, c'est-à-dire que les paramètres d'effets fixes ont des $N(\mu=0,\sigma^2=9)$ (ou $\sigma$, c'est-à-dire l'écart type, $=3$) prieurs. Cela fonctionne assez bien, bien qu'il ne soit pas identique à la correction Firth (puisque Firth correspond à un Jeffreys antérieur , ce qui n'est pas tout à fait le même).

L'exemple lié montre que vous pouvez également le faire avec le MCMCglmmpackage, si vous voulez passer au Bayésien ...