Estou pretendendo realizar o gráfico half-normal plot para um modelo de efeitos mistos que foi ajustado utilizando o pacote lme4
. Para visualizar a análise diagnóstica deste modelo, pretendo utilizar os resíduos de Pearson, no entanto, estou enfrentando algumas dificuldades.
Os dados utilizados podem ser encontrados aqui.
Na rotina computacional descrita abaixo, encontra-se o modelo ajustado.
library(lme4)
library(ggplot2)
library(gridExtra)
library(hnp)
lmm <- lmer(log(Var1)~log(Var2)+log(Var3)+
(Var4)+(Var5)+(Var7||Var6),
Dados2, REML=T)
Por default, o pacote hnp
informa que os resíduos deviance são considerados em casos que se realizou o ajuste de modelos de efeitos mistos pela função lmer
, no entanto, para utilizar outra opção disponível para esta classe de modelos dentro do pacote hnp
basta utilizar a função resid.type
. Sendo assim, abaixo encontra-se a construção do gráfico considerando o default (Ga) e também os resíduos de Pearson (Gb).
Grap5=hnp(lmm, xlab = 'Theoretical quantiles', ylab = 'Residuals')
G5 <- with(Grap5, data.frame(x, lower, upper, median, residuals))
Ga=ggplot(data = G5, aes(x)) +
geom_point(aes(y = residuals)) +
geom_line(aes(y = lower)) +
geom_line(aes(y = upper)) +
geom_line(aes(y = median), linetype = "dashed") +
ylab("Residuals") +
xlab("Theoretical quantiles")+ggtitle("Deviance")+
theme(legend.title = element_text(size = 17),
legend.text = element_text(size = 17),
axis.text.x = element_text(size = 20,color = "black"),
axis.text.y = element_text(size = 20,color = "black"),
axis.title = element_text(size = 20),
axis.text = element_text(size = 25),
plot.title = element_text(size = 20))
Grap6=hnp(lmm, xlab = 'Theoretical quantiles', ylab = 'Residuals',
resid.type="pearson")
G6 <- with(Grap6, data.frame(x, lower, upper, median, residuals))
Gb=ggplot(data = G6, aes(x)) +
geom_point(aes(y = residuals)) +
geom_line(aes(y = lower)) +
geom_line(aes(y = upper)) +
geom_line(aes(y = median), linetype = "dashed") +
ylab("Residuals") +
xlab("Theoretical quantiles")+ggtitle("Pearson")+
theme(legend.title = element_text(size = 17),
legend.text = element_text(size = 17),
axis.text.x = element_text(size = 20,color = "black"),
axis.text.y = element_text(size = 20,color = "black"),
axis.title = element_text(size = 20),
axis.text = element_text(size = 25),
plot.title = element_text(size = 20))
grid.arrange(Ga,Gb,ncol=2)
Vejam na figura abaixo que nenhuma diferença foi observada entre as opções desejadas, logo, me parece que algum problema deva estar acontecendo.
cor(fitted(lmm), resid(lmm))
ouhist(resid(lmm))
verá que o fit não é mau, Em vez de compara os resíduos"deviance"
com Pearson, compare com"response"
e vai ver diferenças (mas pequenas).