No help da função hnp
(páginas 11 e 12 deste documento), você pode ver a lista de modelos que podem ser utilizado. Modelos feitos com o pacote glmer
não constam nessa lista.
Pelo help, você pode ver também que é possível construir este gráfico para um modelo não implementado. Neste caso você precisará fornecer três funções:
- simfun: simula números aleatórios com a distribuição do modelo
- fitfun: ajusta um modelo novamente
- diagfun: calcula a medida de diagnóstico para a qual você quer fazer o envelope.
Por exemplo, considere o seguinte banco de dados:
df <- data_frame(
x1 = runif(100, min = 0, max = 0.5),
x2 = runif(100, min = 0, max = 0.5),
y = as.factor(rbinom(n = 100, size = 1, prob = x1 + x2))
)
df
# A tibble: 100 × 3
x1 x2 y
<dbl> <dbl> <fctr>
1 0.05263752 0.44304477 1
2 0.23866105 0.38094325 1
3 0.16474934 0.12541530 0
4 0.47944698 0.46145428 1
5 0.31256457 0.01484219 1
6 0.28411475 0.17338685 1
7 0.34083505 0.35194770 0
8 0.30603019 0.40298573 1
9 0.22193441 0.44697690 1
10 0.36469374 0.34312517 1
# ... with 90 more rows
Neste exemplo y
é diretamente relacionado a x1
e x2
.
Ajustando um modelo misto neste caso, e rodando a funçãpo hnp
obtenho o mesmo erro que você:
fit_glmer <- glmer(y ~ x1 + x2 + (1 | x1), data = df, family = "binomial")
hnp(fit_glmer)
Binomial-normal model
Error in `[.default`(simulate(object)[, 1], , 1) :
incorrect number of dimensions
Isso acontecxe porque nenhum método está implementado para lidar com objetos do tipo glmer
. Vamos então implementar as funções necessárias para a construção do gráfico:
diagfun <- function(obj) {
resid(obj, type="pearson")
}
simfun <- function(n, obj) {
as.factor(rbinom(n, size = 1, prob = predict(obj, type = "response")))
}
fitfun <- function(y.) {
glmer(y. ~ x1 + x2 + (1 | x1), data = df, family = "binomial")
}
Agora você pode chamar a função hnp
desta forma:
hnp(fit_glmer, newclass=TRUE, diagfun=diagfun, simfun=simfun, fitfun=fitfun)
Ela produzirá o seguinte gráfico:

Claro, este é apenas um exemplo, você precisará adaptar este código para o seu problema específico. A princípio, acredito que o único lugar que você precisa mexer é na função fitfun
. Talvez você precise mexer na função diagfun
, se quiser que o tipo de resíduo não seja Pearson por exemplo (vi na legenda do seu gráfico que você usa resíduos deviance).
hnp
trabalha de maneira correta com os objetos da classe glmerMod.dput(dadosesc)
? Assim, o teu conjunto de dados seria partilhado conosco e seria possível, para a gente, reproduzir o teu problema.