Como lidar com a separação quase completa em um GLMM logístico?

8

Atualização : Como agora sei que meu problema se chama separação quase completa , atualizei a pergunta para refletir isso (graças a Aaron).


Eu tenho um conjunto de dados de um experimento no qual 29 participantes humanos (fator code) trabalharam em um conjunto de ensaios e responsefoi 1 ou 0. Além disso, manipulamos os materiais para termos três fatores cruzados p.validity(válido versus inválido), type(afirmação versus negação) e counterexamples(poucos versus muitos):

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

No geral, há apenas um pequeno número de 0s:

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

Uma hipótese é que existe um efeito de validity, no entanto, análises preliminares sugerem que pode haver um efeito de counterexamples. Como tenho dados dependentes (cada participante trabalhou em todos os ensaios), gostaria de usar um GLMM nos dados. Infelizmente, counterexamplesquase que completamente separe os dados (pelo menos para um nível):

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

Isso também se reflete no modelo:

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

Os erros padrão para os parâmetros são simplesmente insanos. Como meu objetivo final é avaliar se certos efeitos são significativos ou não, os erros padrão não são totalmente sem importância.

  • Como posso lidar com a quase completa separação? O que eu quero é obter estimativas a partir das quais eu possa julgar se um determinado efeito é significativo ou não (por exemplo, usando PRmodcompfrom package pkrtest, mas essa é outra etapa não descrita aqui).

Abordagens usando outros pacotes também são boas.

Henrik
fonte
2
Para começar, tente o seguinte: ats.ucla.edu/stat/mult_pkg/faq/general/…
Aaron deixou o Stack Overflow em
@ Aaron Parece ótimo. Colocar isso em uma resposta teria pelo menos você trouxe um upvote ...
Henrik
Não é realmente uma resposta, mas obrigado!
Aaron saiu de Stack Overflow
@Henrik Você também pode votar nos comentários.
Peter Flom
Veja este artigo de Paul Allison. Embora ele enfatize o SAS, os mesmos pontos serão aplicáveis ​​em outros idiomas.
Peter Flom

Respostas:

8

Receio que haja um erro de digitação no seu título: você não deve tentar ajustar modelos mistos, muito menos modelos mistos não lineares, com apenas 30 clusters. A menos que você acredite que pode ajustar uma distribuição normal a 30 pontos obstruídos por erro de medição, não linearidades e separação quase completa (também conhecida como previsão perfeita).

O que eu faria aqui é executar isso como uma regressão logística regular com a correção de Firth :

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

A correção de Firth consiste em adicionar uma penalidade à probabilidade e é uma forma de retração. Em termos bayesianos, as estimativas resultantes são os modos posteriores do modelo com um Jeffreys anterior. Em termos freqüentes, a penalidade é o determinante da matriz de informações correspondente a uma única observação e, portanto, desaparece assintoticamente.

StasK
fonte
4
Na verdade , acredito na montagem de modelos mistos com menos de 30 clusters. Mas a análise parece, no entanto, promissora (+1). Ou existe o método de Firth para GLMMs?
Henrik
2
Certo, você é o entusiasta dos requisitos mínimos de tamanho de amostra ... A correção do Firth funciona apenas com dados iid. Você pode acreditar em qualquer coisa, mas seria melhor executar algumas simulações para ver se alguma crença é justificada em uma determinada situação de dados. Com um conjunto de dados perfeitamente equilibrado e resposta contínua, pode funcionar bem. Com um conjunto de dados fortemente desequilibrado, em termos de resposta, você vê apenas uma extremidade esquerda da distribuição normal dos efeitos aleatórios, e está disposto a apostar que essa cauda seja bem aproximada por um ponto de Laplace em *lmer??? : - \
StasK
5

Você pode usar uma abordagem bayesiana máxima a posteriori com um fraco antes dos efeitos fixos para obter aproximadamente o mesmo efeito. Em particular, o pacote blme para R (que é um invólucro fino ao redor do lme4pacote) faz isso, se você especificar anteriores para os efeitos fixos, como no exemplo aqui (procure por "separação completa"):

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

tttβΣ=9EuN(μ=0 0,σ2=9)σ=3

O exemplo vinculado mostra que você também pode fazê-lo com o MCMCglmmpacote, se quiser ficar totalmente bayesiano ...

Ben Bolker
fonte