# ========================================================================================= #. Título : ANOVA - Diagnóstico de Pressupostos e Transformação de dados #. Autor : Laboratório de Estatística Computacional - LEC # Data : 2023-07-06 00:14:09 # Tutores: José C. Faria # Ivan B. Allaman # ========================================================================================= # Referências básicas # <<< # https://agronomiar.github.io/livroagro/transforma%C3%A7%C3%A3o-de-dados.html # http://www.leg.ufpr.br/Rpira/Rpira/node17.html # https://rpubs.com/gomes555/transformacao_dos_dados # >>> #. Pacotes necessários library(hnp) # Half-Normal Plots with Simulation Envelopes # hnp função library(car) # Companion to Applied Regression # Funções # - Boxplot # - leveneTest library(gplots) # Various R Programming Tools for Plotting Data # plotmeans função #. Transformações Box-Cox mais comuns # -------------------------- # lambda Transformação # -------------------------- # -3 1/y^3 # -2 1/y^2 # -1 1/y # -0.5 1/sqrt(y) # 0 log(y) # 0.5 sqrt(y) # 1 Nenhuma # 2 y^2 # 3 y^3 #----------------------- #. Exemplo: 1 # Dados adaptados de ZAMBÃO; SAMPAIO; BARBIN, 1982 (Livro Planejamento e Análise Estatística de Experimentos Agronômicos - Décio Barbin), # em que o pesquisador pretende comparar quatro cultivares de pêssego quanto ao enraizamento de estacas. # Foi utilizado cinco repetições por tratamento e o delineamento experimental foi inteiramente casualizado. #.. Dados Y <- c(02, 02, 01, 01, 00, 01, 00, 00, 01, 01, 12, 10, 14, 17, 11, 07, 09, 15, 08, 10) X <- gl(n=4, k=5, labels=LETTERS[1:4]) #.. Análise exploratória par(mfrow=c(2, 1)) car::Boxplot(Y ~ X, main='Boxplot c/ ident. de pontos') # Médias com intervalo de confiança plotmeans(Y ~ X, mean.labels=TRUE, digits=2, connect=FALSE, barcol=1, main='Médias, com I.C.', xlab='', pch=' ') #.. Análise de variância mod <- aov(Y ~ X) summary(mod) par(mfrow=c(2, 2)) plot(mod) #... Violação dos pressupostos! #... Interpretação dos gráficos # ------------- # | 1 | 3 | # ------------- # | 2 | 4 | # ------------- # 1: Falta de ajuste (Residuals vs. Fitted) # 2: Suposição de variância constante (Scale-Location) # 3: Gráfico de normalidade (Normal Q-Q) # 4: Observações influentes (Residuals vs Factor Levels) #.. Pressuposições #... Normalidade dos erros: Shapiro-Wilk shapiro.test(mod$residuals) # Violação: erros seguem não distribuição normal! par(mfrow=c(1, 1)) hist(mod$residuals, nclass=10) #... Homogeneidade das variâncias: Bartlett e fligner # Teste mais rigoroso bartlett.test(mod$residuals ~ X) # Violação: variâncias não são homogêneas! # Teste rigor intermediário fligner.test(Y ~ X) # OK: variâncias são homogêneas! # Teste menos rigoroso leveneTest(mod) # OK: variâncias são homogêneas! # Cálculo da razão entre a maior e a menor variância est_v <- tapply(Y, X, var) est_v # Observação da regra prática: V(Max)/V(Min) <= 4 max(est_v)/min(est_v) # Violação: variâncias não são homogêneas! #... Independência dos erros: Durbin-Watson lmtest::dwtest(mod) # OK! os erros são independentes! # As pressuposições: # - normalidade dos erros # - homogeneidade das variâncias # não foram atendidas. Vamos transformar os dados e conferir novamente as pressuposições! #... Gráfico de resíduos padronizados par(mfrow=c(1, 1)) hnp::hnp(mod, xlab="Quantis teóricos", pch=19) #.. Em busca da transformação ideal # A função boxcox (MASS) não aceita quando ocorre observações 0. # Para contornar, iremos adicionar uma constante com valor bem baixo. bc <- MASS::boxcox(aov(Y + 0.000001 ~ X)) #... Descobrindo o valor exato de lambda (lambda <- bc$x[which.max(bc$y)]) #.. Modelo transformado # A aproximação de lambda é 0.5 -> sqrt(Y) mod_t <- aov(sqrt(Y) ~ X) par(mfrow=c(2, 2)) plot(mod_t) #... OK! #.. Pressuposições após transformação #... Normalidade dos erros: Shapiro-Wilk shapiro.test(mod_t$residuals) # OK: os erros seguem distribuição normal! par(mfrow=c(1, 1)) hist(mod_t$residuals, nclass=10) #... Homogeneidade das variâncias: Bartlett # Teste mais rigoroso bartlett.test(mod_t$residuals ~ X) # OK: as variâncias são homogêneas! # Teste rigor intermediário fligner.test(Y^0.5 ~ X) # OK: as variâncias são homogêneas! # Teste menos rigoroso leveneTest(mod_t) # OK: as variâncias são homogêneas! # Cálculo da razão entre a maior e a menor variância est_v <- tapply(sqrt(Y), X, var) est_v # Observação da regra prática: V(Max)/V(Min) <= 4 max(est_v)/min(est_v) # OK: variâncias são homogêneas! #... Independência dos erros: Durbin-Watson lmtest::dwtest(mod_t) # OK: os erros são independentes! #... Gráfico de resíduos padronizados par(mfrow=c(1, 1)) hnp::hnp(mod_t, xlab="Quantis teóricos", pch=19) #. Exemplo: 2 # Experimento conduzido com o intuito de avaliar a inoculação de Trichoderma sp. (T4), Azospirillum sp. (T3) e associação de ambos (T2) em relação a testemunha, # quanto à altura de plantas de milho. O experimento foi conduzido em delineamento inteiramente casualizado com 8 repetições. #.. Dados Y <- c(124, 136, 124, 102, 112, 108, 102, 122, 130, 128, 118, 106, 126, 106, 128, 122, 132, 132, 190, 144, 090, 126, 142, 148, 140, 120, 118, 098, 110, 140, 104, 142) X <- gl(n=4, k=8, labels=paste0('T', 1:4)) dad <- data.frame(X, Y) head(dad, 10) #.. Análise exploratória par(mfrow=c(2, 1)) car::Boxplot(Y ~ X, dad, main='Boxplot c/ ident. de pontos') # Médias com intervalo de confiança plotmeans(Y ~ X, mean.labels=TRUE, digits=2, connect=FALSE, barcol=1, main='Médias, com I.C.', xlab='', pch=' ') #.. Análise de variância mod <- aov(Y ~ X, dad) summary(mod) par(mfrow=c(2, 2)) plot(mod) #... Violação! #.. Pressuposições #... Normalidade dos erros: Shapiro-Wilk shapiro.test(mod$residuals) # OK: erros seguem distribuição normal! #... Homogeneidade das variâncias: Bartlett # Teste mais rigoroso bartlett.test(mod$residuals ~ X) # Violação: variâncias não são homogêneas! # Teste rigor intermediário fligner.test(Y ~ X) # OK: variâncias são homogêneas! # Teste menos rigoroso leveneTest(mod) # OK: variâncias são homogêneas! # Cálculo da razão entre a maior e a menor variância est_v <- tapply(Y, X, var) est_v # Observação da regra prática: V(Max)/V(Min) <= 4 max(est_v)/min(est_v) # Violação: variâncias não são homogêneas! #... Independência dos erros: Durbin-Watson lmtest::dwtest(mod) # OK! os erros são independentes! #... Gráfico de resíduos padronizados par(mfrow=c(1, 1)) hnp::hnp(mod, xlab="Quantis teóricos", pch=19) # A pressuposiçõa homogeneidade das variâncias não foi atendida. # Dessa forma, vamos transformar os dados e conferir novamente as pressuposições! #.. Em busca da transformação ideal bc <- MASS::boxcox(mod) #... Descobrindo o valor exato de lambda (lambda <- bc$x[which.max(bc$y)]) # O valor de lambda para a transformação Box-Cox é -0,22. # Vamos usar a aproximação. Logo, iremos usar a transformação Log. #.. Modelo transformado mod_t <- aov(log(Y) ~ X) par(mfrow=c(2, 2)) plot(mod_t) #... OK! #.. Pressuposições após transformação #... Normalidade dos erros: Shapiro-Wilk shapiro.test(mod_t$residuals) # OK: os erros seguem distribuição normal! #... Homogeneidade das variâncias: Bartlett # Teste mais rigoroso bartlett.test(mod_t$residuals ~ X) # OK: as variâncias são homogêneas! # Teste rigor intermediário fligner.test(log(Y) ~ X) # OK: as variâncias são homogêneas! # Teste menos rigoroso leveneTest(mod_t) # OK: as variâncias são homogêneas! # Cálculo da razão entre a maior e a menor variância est_v <- tapply(log(Y), X, var) est_v # Observação da regra prática: V(Max)/V(Min) <= 4 max(est_v)/min(est_v) # Violação: variâncias não são homogêneas! #... Independência dos erros: Durbin-Watson lmtest::dwtest(mod_t) # OK: os erros são independentes! #... Gráfico de resíduos padronizados par(mfrow=c(1, 1)) hnp::hnp(mod_t, xlab="Quantis teóricos", pch=19) #. Exemplo: 3 # Controle de carrapatos #.. Dados dad <- read.table('https://lec.pro.br/download/R/dados/poly.txt', h=T) # Adequação dos dados str(dad) head(dad) dad$aplic <- factor(dad$aplic) dad$grupo <- factor(dad$grupo) str(dad) head(dad, 10) #.. Análise exploratória par(mfrow=c(2, 1)) car::Boxplot(ncarra ~ aplic*grupo, dad, main='Boxplot c/ ident. de pontos') plotmeans(ncarra ~ interaction(aplic, grupo, sep ="."), data=dad, mean.labels=TRUE, digits=2, connect=FALSE, barcol=1, main='Médias, com I.C.', xlab='aplic*grupo', pch=' ') #.. Análise de variância mod <- aov(ncarra ~ aplic*grupo, dad) summary(mod) par(mfrow=c(2, 2)) plot(mod) #... Violação! #.. Pressuposições #... Normalidade dos erros: Shapiro-Wilk shapiro.test(mod$residuals) # Violação: erros não seguem distribuição normal! #... Homogeneidade das variâncias: Bartlett # Definindo os tratamentos no fatorial (X <- with(dad, paste(aplic, grupo, sep='_'))) # 18 tratamentos no fatorial! length(levels(as.factor(X))) # Teste mais rigoroso bartlett.test(mod$residuals ~ X) # Violação: variâncias não são homogêneas! # 17 gl associado ao teste! # Teste rigor intermediário fligner.test(dad$ncarra ~ X) # Violação: variâncias não são homogêneas! # Teste menos rigoroso leveneTest(mod) # Violação: variâncias não são homogêneas! # Cálculo da razão entre a maior e a menor variância est_v <- with(dad, tapply(ncarra, list(aplic, grupo), var)) est_v # Observação da regra prática: V(Max)/V(Min) <= 4 max(est_v)/min(est_v) # Violação: variâncias não são homogêneas! #... Independência dos erros: Durbin-Watson lmtest::dwtest(mod) # Violação: os erros não são independentes! #... Gráfico de resíduos padronizados par(mfrow=c(1, 1)) hnp::hnp(mod, xlab="Quantis teóricos", pch=19) # A pressuposiçõa homogeneidade das variâncias não foi atendida. # Dessa forma, vamos transformar os dados e conferir novamente as pressuposições! #.. Em busca da transformação ideal par(mfrow=c(1, 1)) bc <- with(dad, boxcox((ncarra + 0.000001) ~ aplic*grupo, lam=seq(-0.5, 1.5, l=100), plotit=T)) #... Descobrindo o valor exato de lambda (lambda <- bc$x[which.max(bc$y)]) # O valor de lambda para a Transformação Box-Cox é 0.3282828. # Vamos usar a aproximação. Logo, iremos usar a transformação sqrt. #.. Modelo transformado # Transformação boxcox mod_t <- aov((ncarra^lambda - 1)/lambda ~ aplic*grupo, dad) par(mfrow=c(2, 2)) plot(mod_t) #... OK! #.. Análise exploratória par(mfrow=c(2, 1)) car::Boxplot((ncarra^lambda - 1)/lambda ~ aplic*grupo, dad, main='Boxplot c/ ident. de pontos') plotmeans((ncarra^lambda - 1)/lambda ~ interaction(aplic, grupo, sep ="."), data=dad, mean.labels=TRUE, digits=2, connect=FALSE, barcol=1, main='Médias, com I.C.', xlab='aplic*grupo', pch=' ') #.. Pressuposições após transformação #... Normalidade dos erros: Shapiro-Wilk shapiro.test(mod_t$residuals) # OK: os erros seguem distribuição normal! #... Homogeneidade das variâncias: Bartlett # Teste mais rigoroso bartlett.test(mod_t$residuals ~ X, dad) # Violação: as variâncias não são homogêneas! # Teste rigor intermediário fligner.test((dad$ncarra^lambda - 1)/lambda ~ X) # OK: variâncias são homogêneas! # Teste menos rigoroso leveneTest(mod_t) # OK: variâncias são homogêneas! # Cálculo da razão entre a maior e a menor variância est_v <- with(dad, tapply((ncarra^lambda - 1)/lambda, list(aplic, grupo), var)) est_v # Observação da regra prática: V(Max)/V(Min) <= 4 max(est_v)/min(est_v) # Violação: as variâncias não são homogêneas! #... Independência dos erros: Durbin-Watson lmtest::dwtest(mod_t) # Violação: os erros não são independentes! #... Gráfico de resíduos padronizados par(mfrow=c(1, 1)) hnp::hnp(mod_t, xlab="Quantis teóricos", pch=19) #... ------------------------------------------------------------------------ #... OBS: A tranformação proposta melhora os dados, mas ainda não é a ideal! #... ------------------------------------------------------------------------ #. Exemplo: 4 dad <- read.table('https://lec.pro.br/download/R/dados/cris.txt', h=T, dec=',', stringsAsFactors=T) str(dad) head(dad, 10) #.. Análise exploratória par(mfrow=c(2, 1)) car::Boxplot(hiperativo ~ gar*trat, dad, main='Boxplot c/ ident. de pontos') plotmeans(hiperativo ~ interaction(gar, trat, sep ="."), data=dad, mean.labels=TRUE, digits=2, connect=FALSE, barcol=1, main='Médias, com I.C.', xlab='gar*trat', pch=' ') #.. Análise de variância mod <- aov(hiperativo ~ gar*trat, dad) summary(mod) par(mfrow=c(2, 2)) plot(mod) #... Violação! #.. Pressuposições #... Normalidade dos erros: Shapiro-Wilk shapiro.test(mod$residuals) # Violação: erros não seguem distribuição normal! #... Homogeneidade das variâncias: Bartlett # Definindo os tratamentos no fatorial (X <- with(dad, paste(gar, trat, sep='_'))) # 16 tratamentos no fatorial! length(levels(as.factor(X))) # Teste mais rigoroso bartlett.test(mod$residuals ~ X) # Violação: variâncias não são homogêneas! # 15 gl associado ao teste! # Teste rigor intermediário fligner.test(dad$hiperativo ~ X) # OK: variâncias são homogêneas! # Teste menos rigoroso leveneTest(mod) # OK: variâncias são homogêneas! # Cálculo da razão entre a maior e a menor variância est_v <- with(dad, tapply(hiperativo, list(gar, trat), var)) est_v # Observação da regra prática: V(Max)/V(Min) <= 4 max(est_v)/min(est_v) # Violação: variâncias não são homogêneas! #... Independência dos erros: Durbin-Watson lmtest::dwtest(mod) # OK: os erros são independentes! #... Gráfico de resíduos padronizados par(mfrow=c(1, 1)) hnp::hnp(mod, xlab="Quantis teóricos", pch=19) # As pressuposições da normalidade dos erros e homogeneidade das variâncias não foi atendida. # Dessa forma, vamos transformar os dados e conferir novamente as pressuposições! #.. Em busca da transformação ideal par(mfrow=c(1, 1)) bc <- with(dad, boxcox((hiperativo + 0.000001) ~ gar*trat, lam=seq(-0.5, 1.5, l=100), plotit=T)) #... Descobrindo o valor exato de lambda (lambda <- bc$x[which.max(bc$y)]) #.. Modelo transformado # Transformação boxcox mod_t <- aov((hiperativo^lambda - 1)/lambda ~ gar*trat, dad) par(mfrow=c(2, 2)) plot(mod_t) #... OK! #.. Análise exploratória par(mfrow=c(2, 1)) car::Boxplot((hiperativo^lambda - 1)/lambda ~ gar*trat, dad, main='Boxplot c/ ident. de pontos') plotmeans((hiperativo^lambda - 1)/lambda ~ interaction(gar, trat, sep ="."), data=dad, mean.labels=TRUE, digits=2, connect=FALSE, barcol=1, main='Médias, com I.C.', xlab='gar*trat', pch=' ') #.. Pressuposições após transformação #... Normalidade dos erros: Shapiro-Wilk shapiro.test(mod_t$residuals) # OK: os erros seguem distribuição normal! #... Homogeneidade das variâncias: Bartlett # Teste mais rigoroso bartlett.test(mod_t$residuals ~ X, dad) # Violação: as variâncias não são homogêneas! # Teste rigor intermediário fligner.test((dad$hiperativo^lambda - 1)/lambda ~ X) # OK: variâncias são homogêneas! # Teste menos rigoroso leveneTest(mod_t) # OK: variâncias são homogêneas! # Cálculo da razão entre a maior e a menor variância est_v <- with(dad, tapply((hiperativo^lambda - 1)/lambda, list(gar, trat), var)) est_v # Observação da regra prática: V(Max)/V(Min) <= 4 max(est_v)/min(est_v) # Violação: as variâncias não são homogêneas! #... Independência dos erros: Durbin-Watson lmtest::dwtest(mod_t) # OK: os erros são independentes! #... Gráfico de resíduos padronizados par(mfrow=c(1, 1)) hnp::hnp(mod_t, xlab="Quantis teóricos", pch=19) #... ------------------------------------------------------------------------ #... OBS: A tranformação proposta melhora os dados, mas ainda não é a ideal! #... ------------------------------------------------------------------------ #. Exemplo: 5 dad <- read.table('https://lec.pro.br/download/R/dados/marc.txt', h=T, colClasses=c('factor', 'factor', 'numeric')) str(dad) head(dad, 15) #.. Análise exploratória par(mfrow=c(2, 1)) car::Boxplot(cortisol ~ tempo*den, dad, main='Boxplot c/ ident. de pontos') plotmeans(cortisol ~ interaction(tempo, den, sep ="."), data=dad, mean.labels=TRUE, digits=2, connect=FALSE, barcol=1, main='Médias, com I.C.', xlab='tempo*den', pch=' ') #.. Análise de variância mod <- aov(cortisol ~ tempo*den, dad) summary(mod) par(mfrow=c(2, 2)) plot(mod) #... Violação! #.. Pressuposições #... Normalidade dos erros: Shapiro-Wilk shapiro.test(mod$residuals) # Violação: erros não seguem distribuição normal! #... Homogeneidade das variâncias: Bartlett # Definindo os tratamentos no fatorial (X <- with(dad, paste(tempo, den, sep='_'))) # 16 tratamentos no fatorial! length(levels(as.factor(X))) # Teste mais rigoroso bartlett.test(mod$residuals ~ X) # Violação: variâncias não são homogêneas! # 11 gl associado ao teste! # Teste rigor intermediário fligner.test(dad$cortisol ~ X) # OK: variâncias são homogêneas! # Teste menos rigoroso (geralmente!) leveneTest(mod) # OK: variâncias são homogêneas! # Cálculo da razão entre a maior e a menor variância est_v <- with(dad, tapply(cortisol, list(tempo, den), var)) est_v # Observação da regra prática: V(Max)/V(Min) <= 4 max(est_v)/min(est_v) # Violação: variâncias não são homogêneas! #... Independência dos erros: Durbin-Watson lmtest::dwtest(mod) # OK: os erros são independentes! #... Gráfico de resíduos padronizados par(mfrow=c(1, 1)) hnp::hnp(mod, xlab="Quantis teóricos", pch=19) # As pressuposições da normalidade dos erros e homogeneidade das variâncias não foi atendida. # Dessa forma, vamos transformar os dados e conferir novamente as pressuposições! #.. Em busca da transformação ideal par(mfrow=c(1, 1)) bc <- with(dad, boxcox((cortisol) ~ tempo*den, lam=seq(-0.5, 1.5, l=100), plotit=T)) #... Descobrindo o valor exato de lambda (lambda <- bc$x[which.max(bc$y)]) #.. Modelo transformado # Transformação boxcox mod_t <- aov((cortisol^lambda - 1)/lambda ~ tempo*den, dad) ## Pode aproximar para 0.5, uma vez que 0.5 está no intervalo de confiança de lambda. #mod_t <- aov(sqrt(cortisol) ~ tempo*den, # dad) par(mfrow=c(2, 2)) plot(mod_t) #... OK! #.. Análise exploratória par(mfrow=c(2, 1)) car::Boxplot((cortisol^lambda - 1)/lambda ~ tempo*den, dad, main='Boxplot c/ ident. de pontos') plotmeans((cortisol^lambda - 1)/lambda ~ interaction(tempo, den, sep ="."), data=dad, mean.labels=TRUE, digits=2, connect=FALSE, barcol=1, main='Médias, com I.C.', xlab='tempo*den', pch=' ') #.. Pressuposições após transformação #... Normalidade dos erros: Shapiro-Wilk shapiro.test(mod_t$residuals) # OK: os erros seguem distribuição normal! #... Homogeneidade das variâncias: Bartlett # Teste mais rigoroso bartlett.test(mod_t$residuals ~ X, dad) # Violação: as variâncias não são homogêneas! # Teste rigor intermediário fligner.test((dad$cortisol^lambda - 1)/lambda ~ X) # OK: variâncias são homogêneas! # Teste menos rigoroso (geralmente) leveneTest(mod_t) # OK: variâncias são homogêneas! # Cálculo da razão entre a maior e a menor variância est_v <- with(dad, tapply((cortisol^lambda - 1)/lambda, list(tempo, den), var)) est_v # Observação da regra prática: V(Max)/V(Min) <= 4 max(est_v)/min(est_v) # Violação: as variâncias não são homogêneas! #... Independência dos erros: Durbin-Watson lmtest::dwtest(mod_t) # OK: os erros são independentes! #... Gráfico de resíduos padronizados par(mfrow=c(1, 1)) hnp::hnp(mod_t, xlab="Quantis teóricos", pch=19) #... ------------------------------------------------------------------------ #... OBS: A tranformação proposta melhora os dados, mas ainda não é a ideal! #... ------------------------------------------------------------------------