0

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.

inserir a descrição da imagem aqui

1
  • Não há nada de errado no código, os gráficos são assim mesmo. Se tentar cor(fitted(lmm), resid(lmm)) ou hist(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). 16/04/2021 às 17:14

1 Resposta 1

2

Eis uma solução para ter gráficos com resíduos "deviance" e "response", que me parece ser a mais útil.

O código foi simplificado também para mostrar como se pode usar as facetas diretamente, sem ter de traçar dois gráficos e depois juntá-los com grid.arrange.

1. Definir um novo theme.

Em primeiro lugar, defini um theme custom.

theme_BrenoG <- function(){ 
  theme_bw() %+replace%    #replace elements we want to change
    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),
          strip.text.x = element_text(size = 20))
}

2. O modelo.

library(lme4)
library(ggplot2)
library(gridExtra)
library(hnp)

lmm <- lmer(log(Var1) ~ log(Var2) + log(Var3) + Var4 + Var5 + (Var7||Var6),
            data = Dados2,
            REML = TRUE)

Grap5 <- hnp(lmm, xlab = 'Theoretical quantiles', ylab = 'Residuals',
             resid.type = "deviance")
Grap6 <- hnp(lmm, xlab = 'Theoretical quantiles', ylab = 'Residuals',
             resid.type = "response")

G5 <- with(Grap5, data.frame(x, lower, upper, median, residuals))
G6 <- with(Grap6, data.frame(x, lower, upper, median, residuals))

3. O gráfico.

Ao criar uma coluna com o tipo de resíduos que os modelos têm, pode-se depois facetar por essa variável.

G5$type <- "Deviance"
G6$type <- "Response"


ggplot(data = rbind(G5, 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") +
  facet_wrap(~ type) +
  theme_BrenoG()

inserir a descrição da imagem aqui

4. Leitura dos dados.

Os dados são lidos com read.csv porque os separadores são ";" e as decimais são marcadas com ",".

google_id <- "1q4XuSPnZXx7TsrT3JVoucFyxeOSegWFz"
google_file <- sprintf("https://docs.google.com/uc?id=%s&export=download", google_id)
Dados2 <- read.csv2(google_file)
0

Você deve fazer log-in para responder a esta pergunta.

Esta não é a resposta que você está procurando? Pesquise outras perguntas com a tag .