detecção de outliers modelos mistos não lineares
Caros colegas, Nos comandos abaixo, realizo a análise de um modelo não linear misto. Após selecionar o modelo, realizei um teste outliers do modelo. Identifiquei os outliers e os removi. Ajustei novamente o modelo, porém a nova análise gráfica dos resíduos desses modelo indica novos outliers. Porque isso ocorre? Como lido com isso? Removê-los não é uma opção, uma vez que estarei perdendo muitas observações. A minha outra dúvida é: Qual o melhor momento para realizar o teste de outliers? Eu realizo após a seleção do melhor modelo, pois com a correta definição dos efeitos aleatórios e da matriz R, o número de outliers normalmente reduz. Porém a remoção de outliers pode alterar o modelo e as estruturas até então selecionadas para as matrizes D e R podem não ser mais adequadas. Qual procedimento mais indicado? Att library(nlme) install.packages("tidyr") library(tidyr) DADOS.CURTO <- read.csv(" https://www.dropbox.com/s/b4cckybgrcg86me/galinhas.csv?raw=1",head = TRUE) head(DADOS.CURTO) colnames(DADOS.CURTO) <- c("0","1","2","3","4","5","6","TRAT") DADOS<- DADOS.CURTO %>% gather('0','1','2','3','4','5','6',key ="SEMANA",value = "PESO", -TRAT) head(DADOS) str(DADOS) ID <-rep(seq(1,370), times=7) DADOS<-cbind(ID,DADOS) DADOS <- transform(DADOS, SEMANA = as.numeric(SEMANA),ID=factor(ID)) #-----------------------------------ajuste dos modelos---------------------------------------------------- MNG8<-gnls(PESO ~ a * exp(-b * exp( -c * SEMANA)),params=list(a~TRAT,b~1,c~1),start = list(c(a = c(4.10,4.79,4.75),b = 4.53,c = 0.37)),na.action=na.omit,data =DADOS) MNG8.4 <- update(MNG8,weights = varExp(form = ~fitted(.)|SEMANA),correlation = corARMA(form = ~SEMANA |ID,q = 6)) #-----------------------------------outliers---------------------------------------------------------------- (GRAF8 <- plot(MNG8final,resid(.,type = "pearson")~fitted(.),id=0.05,xlab="Valores ajustados",ylab="Residuos Padronizados",abline = 0,main = "MNG8.4")) #------------------------------------ Identificação de outliers--------------------------------------- DADOS <- DADOS[complete.cases(DADOS),] ifelse(abs(residuals(MNG8,type="pearson"))>qnorm(0.975),1,0)-> DADOS$cookMNG #-----------------------------------novo modelo---------------------------------------------- MNG8final <- update(MNG8,data=subset(DADOS,cookMNG==0)) #------------------------------- (GRAF8 <- plot(MNG8final,resid(.,type = "pearson")~fitted(.),id=0.05,xlab="Valores ajustados",ylab="Residuos Padronizados",abline = 0,main = "MNG8.4")) -- ========================================= Fernando Souza Zootecnista, DSc. Produção e Alimentação Animal Celular: (31)99796-8781 (Vivo) E-mail:nandodesouza@gmail.com <e-mail%3Anandodesouza@gmail.com> Lattes: http://lattes.cnpq.br/6519538815038307 Blog: https://producaoanimalcomr.wordpress.com/ ==========================================
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) À disposição. Walmes.
participantes (3)
-
Fernando Antonio de souza -
Fernando Souza
-
Walmes Zeviani