Glmm binomial com uma variável categórica com sucessos totais

11

Estou executando um glmm com uma variável de resposta binomial e um preditor categórico. O efeito aleatório é dado pelo design aninhado usado para a coleta de dados. Os dados são assim:

m.gen1$treatment
 [1] sucrose      control      protein      control      no_injection .....
Levels: no_injection control sucrose protein
m.gen1$emergence 
 [1]  1  0  0  1  0  1  1  1  1  1  1  0  0....
> m.gen1$nest
 [1] 1  1  1  2  2  3  3  3  3  4  4  4  .....
Levels: 1 2 3 4 5 6 8 10 11 13 15 16 17 18 20 22 24

O primeiro modelo que eu corro fica assim

m.glmm.em.<-glmer(emergence~treatment + (1|nest),family=binomial,data=m.gen1)

Recebo dois avisos parecidos com este:

Warning messages:
1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model failed to converge with max|grad| = 0.0240654 (tol = 0.001, component 4)
2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model is nearly unidentifiable: large eigenvalue ratio
 - Rescale variables?

O resumo do modelo mostra que um dos tratamentos tem um erro padrão incomumente grande, que você pode ver aqui:

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)  
(Intercept)         2.565      1.038   2.472   0.0134 *
treatmentcontrol   -1.718      1.246  -1.378   0.1681  
treatmentsucrose   16.863   2048.000   0.008   0.9934  
treatmentprotein   -1.718      1.246  -1.378   0.1681 

Eu tentei os diferentes otimizadores do controle glmer e funções de outros pacotes e recebo uma saída semelhante. Eu executei o modelo usando glm ignorando o efeito aleatório, e o problema persiste. Ao explorar os dados, percebi que o tratamento com um alto padrão. O erro tem apenas sucessos na variável de resposta. Apenas para verificar se isso pode estar causando o problema, adicionei um ponto de dados falso com uma "falha" para esse tratamento e o modelo funciona sem problemas e gera um erro padrão razoável. Você pode ver isso aqui:

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)  
(Intercept)        3.4090     1.6712   2.040   0.0414 *
treatmentcontrol  -1.8405     1.4290  -1.288   0.1978  
treatmentsucrose  -0.2582     1.6263  -0.159   0.8738  
treatmentprotein  -2.6530     1.5904  -1.668   0.0953 .

Fiquei me perguntando se minha intuição está certa sobre a falta de falhas para esse tratamento, impedindo uma boa estimativa, e como posso solucionar esse problema.

Desde já, obrigado!

ameba diz Restabelecer Monica
fonte

Respostas:

15

Sua intuição está exatamente certa. Esse fenômeno é chamado de separação completa . Você pode encontrar bastante (agora que você sabe o nome) pesquisando no Google ... É bastante discutido aqui em um contexto geral e aqui no contexto dos GLMMs . A solução padrão para esse problema é adicionar um pequeno termo que empurre os parâmetros de volta para zero - em contextos freqüentistas, isso é chamado de método penalizado ou corrigido por viés . O algoritmo padrão é devido a Firth (1993, "Redução de polarização das estimativas de máxima verossimilhança" Biometrika 80, 27-38) e é implementado no pacote logistfno CRAN. Nos contextos bayesianos, isso é descrito como adicionando um valor fraco antes dos parâmetros de efeito fixo.

Que eu saiba, o algoritmo de Firth não foi estendido aos GLMMs, mas você pode usar o truque bayesiano usando o pacote blme , que coloca uma fina camada bayesiana sobre a parte superior do lme4pacote. Aqui está um exemplo da discussão do GLMM acima vinculada:

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

As duas primeiras linhas neste exemplo são exatamente as mesmas que usamos no glmermodelo padrão ; o último especifica que o anterior para os efeitos fixos é uma distribuição normal multivariada com uma matriz diagonal de variância-covariância. A matriz é 4x4 (porque temos 4 parâmetros de efeito fixo neste exemplo) e a variação anterior de cada parâmetro é 9 (correspondendo a um desvio padrão de 3, que é bastante fraco - isso significa que +/- 2 SD é ( -6,6), que é uma faixa muito grande na escala logit).

Os erros padrão muito grandes dos parâmetros em seu exemplo são um exemplo de um fenômeno intimamente relacionado à separação completa (ocorre sempre que obtemos valores extremos de parâmetros em um modelo logístico) chamado efeito Hauck-Donner .

Mais duas referências potencialmente úteis (eu não as procurei ainda):

  • Gelman A, Jakulin A, Pittau MG e Su TS (2008) Uma distribuição anterior padrão pouco informativa para modelos de regressão logística e outros. Annals of Applied Statistics , 2, 1360-383.
  • José Cortiñas Abrahantes e Marc Aerts (2012) Uma solução para separação para dados binários agrupados Modelagem Estatística 12 (1): 3–27 doi: 10.1177 / 1471082X1001200102

Uma pesquisa mais recente no Google por "bglmer 'complete separação'" mostra:

  • Quiñones, AE e WT Wcislo. “Cuidado críptico estendido para ninhadas na abelha suor facultativamente eusocial Megalopta genalis .” Insectes Sociaux 62.3 (2015): 307-313.
Ben Bolker
fonte
uau muito obrigado !! Isso faz todo o sentido, e o modelo agora funciona sem problemas com o bglmer. Gostaria apenas de mais uma pergunta. Posso usar os métodos como no lme4 para avaliar os efeitos aleatórios e fixos, em outras palavras, para comparar modelos diferentes?
2
Eu diria que sim, mas eu não sei se há qualquer apoio formal e / ou peer-reviewed a minha opinião ...
Ben Bolker
Obrigado! Este é exatamente o meu problema também. Um acompanhamento rápido: em contraste com o seu exemplo, que possui um fator com 4 níveis, eu tenho um design 2 x 2 em que cada fator possui 2 níveis (portanto, o total ainda é de 4 níveis). Também posso usar diag (9,4) para o meu modelo? Eu não sou bem versado em matrizes, então eu queria checar novamente. De maneira semelhante, para justificar esta solução em meu artigo, devo citar Firth (1993) ou existe um artigo mais diretamente relevante, que implementou sua solução usando bglmer ()?
Sol
2
veja resposta atualizada.
Ben Bolker
2
Acho que sim - só importa quantos parâmetros de efeito fixo existem no total.
precisa