Você tem que estabelecer um critério de parada porque você pode encontrar outliers se continuar o ciclo de buscar por eles.
A sua estratégia de remoção acaba sendo modelo dependente, ou seja, você ajusta, remove outliers, ajusta de novo, vê que o modelo mudou (e deveria mesmo), mas surgem outros outliers pois na presença dos dados restantes, novas observações se tornam vilãs (influentes). Se o seu modelo estiver mau especificado, os danos são ainda maiores porque os "outliers" (que não seriam) são apenas observações injustiçadas pela má especificação do modelo. Se você olhar a curva de cada indivíduo, vai concluir que as observações dentro dos indivíduos estão coesas, ou seja, desvios da curva imaginária (uso panel.smoother()) para me auxiliar com a curva imaginária).
# Olhando para cada a curva de cada ID.
xyplot(PESO ~ SEMANA | ID,
data = DADOS,
strip = FALSE,
type = "o",
as.table = TRUE,
jitter.x = TRUE) +
layer(panel.smoother(..., col = 1))
# Olhando por tratamento.
xyplot(PESO ~ SEMANA | TRAT,
groups = ID,
data = DADOS,
type = "o",
jitter.x = TRUE) +
layer(panel.smoother(..., col = 1, lwd = 3))
Acredito que o modelo mais apropriado para os seus dados seja um modelo de efeitos aleatórios, coisa que você não está fazendo. A variação entre ID você está estruturando já de início com função de variância e correlação entre observações. Acredito que primeiro deveria acomodar a variação do estrato superior (ID) para depois fazer um ajuste mais fino (com função de variância e correlação se forem necessários) especificando elementos para o estrato inferior (residual). Você está modelando a camada inferior sem dar a devida a atenção para o estrato de cima. Em muitos casos, a inclusão dos efeitos aleatórios vais dispensar modelagem da variância porque a variância marginal aumenta com as semanas mas não é porque ela cresce dentro dos ID mas sim entre ID. E a correlação pode desaparecer também, porque muitas vezes ela resulta de um modelo que apresenta falta de ajuste. Abaixo eu uso um polinômio de 4 grau, que é bem flexível, ou seja, o modelo "corre para os dados". Nos gráficos eu não fiquei convencido de que existem observações influentes.
# Criando objeto de classe groupedData.
da <- groupedData(PESO ~ SEMANA | ID,
data = DADOS,
order.groups = FALSE)
# Ajuste de um modelo suficientemente flexível.
m0 <- lme(PESO ~ SEMANA + I(SEMANA^2) + I(SEMANA^3) + I(SEMANA^4),
random = ~ 1 | ID,
na.action = na.omit,
data = da)
summary(m0)
plot(augPred(m0), strip = FALSE)
plot(residuals(m0))
plot(residuals(m0, type = "normalized"))
xyplot(residuals(m0) ~ fitted(m0) | TRAT, data = da)
qqmath(~residuals(m0) | TRAT, data = da)
qqmath(~residuals(m0) | ID, data = da, strip = FALSE)