SCRIPT “R” PARA ANOVA DE DELINEAMETNOS BÁSICOS (DIC, DBC, DQL),COM VERIFICAÇÃO DE PRESSUPOSTOS E TESTES DE MÉDIAS

# ======================================================================# ANÁLISE DE VARIÂNCIA - ENTRANDO DADOS PARA UMA ANOVA - "read.table ()"# ======================================================================help (read.table)exp1 <- read.table (file="C:/Users/Joao/Desktop/Biometria/2016/Des_Bas.csv", sep=",", header=TRUE) # Buscando dados de outro arquivo (neste caso, gerado como CSV no Excel)head(exp1)mean(PROD)mean(exp1$PROD)head(exp1) # lista o cabeçalho do conjunto de dados (seis primeiras linhas apenas)plot(exp1) # gráfico para inspeções iniciais (ex. correlações entre variáveis mensuradas nas parcelas)summary(exp1) # resumo estatísticos simples das variáveis do conjunto de dadosstr(exp1) # mostra a estrutura interna do objeto 'exp1'sort(TRAT) # um mero teste (que não funciona, pois o objeto em que está a variável TRAT não está sob "prioridade" no caminho de busca do R)sort(PROD)names(exp1) # lista as variáveis no objeto 'exp1'exp1$TRATexp1$PRODsort(exp1$TRAT) # outro teste (agora funciona, pois '$' remete a uma variável, no caso TRAT, no objeto específico, 'exp')sort(exp1$PROD)attach(exp1) # com este comando o problema se resolve , pois o objeto 'exp1' agora é anexado diretamente ao caminho de busca do Rhelp(attach)sort(TRAT) # constate agora que a função aplicada à variável funciona, sem indicar o objeto associadosort(PROD)exp1[2,] # listando linha escolhida (ex. linha 2)exp1[,4] # listando coluna escolhida (ex. coluna 4)exp1[PROD>500,] # listando observações com ALT>2 (todas as colunas)hist(PROD)boxplot(PROD)plot(PROD ~ TRAT, exp1) # observe como a mesma instrução gera gráficos diferentes (aqui, variável numérica vs. categórica)plot(PROD ~ ALT, exp1) # (aqui, variável numérica vs. variável numérica - gráfico de dispersão XY)str(exp1)

# PREPARANDO PARA UMA ANOVA - declarando Fatores (duas opções)# ============================================================exp1$TRAT <- factor(exp1$TRAT)exp1$LINHA <- factor(exp1$LINHA)exp1$COLUNA <- factor(exp1$COLUNA)exp1$BLOCOS1 <- factor(exp1$LINHA)exp1$BLOCOS2 <- factor(exp1$COLUNA)TRAT <- as.factor(TRAT) # Depende do comando 'attach' ter sido executado antesLINHA <- as.factor(LINHA)COLUNA <- as.factor(COLUNA)BLOCOS1 <- as.factor(LINHA)BLOCOS2 <- as.factor(COLUNA)

# ANOVA DE DADOS BALANCEADOS / MODELO FIXO (Delineamentos Básicos)# ================================================================

# 1) Delineamento Inteiramente Casualizado (D.I.C.)# ----------------------------------------------is.factor(TRAT) # checagem fundamental (pois, numa ANOVA, a variável TRAT deve ser um fator, não uma variável numérica)levels(TRAT)

# 1a. Opção: uso da função 'aov' (analysis of variance) - para ANOVA de dados balanceados modelo1 <- aov(PROD ~ TRAT, data=exp1) help(aov)modelo1names(modelo1)fitted.values(modelo1)residuals(modelo1)plot(residuals(modelo1)~fitted.values(modelo1))summary(modelo1)anova(modelo1)

shapiro.test(residuals(modelo1)) # teste de normalidade dos resíduosbartlett.test(residuals(modelo1)~exp1$TRAT) # teste de homogeneidade das variâncias de TRATbartlett.test(PROD ~ TRAT, data=exp1) # instrução equivalente (apenas em DIC)

# Opção mais geral: função 'lm' (linear model), que acomoda desbalanceamentosmodelo1b <- lm(PROD ~ TRAT, data=exp1)

help(lm)summary(modelo1b) # soluções do sistema de equações normais ("estimativas" dos parâmetros)anova(modelo1b) # quadro da análise de variância (nesta saída "Residual standard error" é o desvio padrão residual = Raiz(QMerro))names(modelo1b) # variáveis com resultados do ajuste do modelo definido no objeto 'modelo1b'par(mfrow=c(2,2)) # divide a matriz gráfica em duas linhas e duas colunas (recebe quatro gráficos)plot(modelo1b) # diagnóstico inicial de problemas relacionados às pressuposições da ANOVA clássica layout(1) # retorna à configuração gráfica padrão (um gráfico por página)

# CASO AS PRESSUPOSIÇÕES SEJAM ATENDIDAS => COMPARAÇÕES DE MÉDIAS (Trats. Qualitativos)TukeyHSD(aov(PROD ~ TRAT, data=exp1)) # Teste Tukey (alpha=5%) - função para dados balanceados (associada a 'aov')TukeyHSD(modelo1) # idem (Teste Tukey - somente associado à função 'aov' - sob balanceamento)result<- TukeyHSD(modelo1,"TRAT",ordered=T)plot(result)

# DMS (t-Student para comparação de médias duas a duas - necessidade de correção da taxa de erro tipo I)with(modelo1, pairwise.t.test(PROD, TRAT, p.adj="fdr")) # Comparação de médias - t-Student (alpha=5%) c/ ajust. FDR

# Teste Tukey c/ letras (script adapt. de Ribeiro Junior, 2009) - funciona bem para D.I.C.medias <- tapply(PROD, TRAT, mean) mediascompara <- TukeyHSD(modelo1)dms <- unname(0.5*diff(compara[[1]][1, 2:3]))medias.ord <- sort(medias, decreasing = TRUE)i <- pos <- letra <- 1letras <- character(nlevels(TRAT))while (i <= nlevels(TRAT)) {print(letters[letra])ind <- (medias.ord[i] - (medias.ord[-(1:i)])) < dmspos.i <- i + sum(ind)if (pos.i > pos) {letras.vec <- rep(" ", length(letras))letras.vec[i:pos.i] <- letters[letra]letras <- paste(letras, letras.vec, sep =" ")pos <- pos.iletra <- letra + 1}i <- i + 1

}data.frame(Medias = medias.ord, "Tukey" = letras)

# Testes de Comparações Múltiplas (TCM) - obtendo os resultados com o pacote 'agricolae'require(agricolae)xtabs(~TRAT, na.omit(exp1)) # repetições por nível de TRATglr <- df.residual(modelo1b) # graus de liberdade do Resíduoqmr <- deviance(modelo1b)/glr # Quadrado médio do Resíduo# teste de Tukeyc3<-with(exp1, HSD.test(PROD, TRAT, DFerror=glr, MSerror=qmr, alpha=0.05))c3

# Obtendo os resultados com o pacote ExpDes (inclui Testes Comparações Múltiplas)require(ExpDes)require(doBy)with(exp1, crd(treat=TRAT, resp=PROD, mcomp="tukey")) # Tukey (tradicional)with(exp1, crd(treat=TRAT, resp=PROD, mcomp="duncan")) # Duncan (tradicional)with(exp1, crd(treat=TRAT, resp=PROD, mcomp="snk")) # Student-Newman-Keulswith(exp1, crd(treat=TRAT, resp=PROD, mcomp="lsd")) # t-Student (LSD)with(exp1, crd(treat=TRAT, resp=PROD, mcomp="lsdb")) # t-Student / Bonferroniwith(exp1, crd(treat=TRAT, resp=PROD, mcomp="sk")) # Scott-Knott (sem ambiguidades) - sob balanceamentohelp(ExpDes) # vide 'index'

# Um teste livre de distribuições por Bootstrap (Ramos & Ferreira, 2009: Revista Ceres, v.56, p.140-149)crd(TRAT, PROD, quali = TRUE, mcomp='ccboot', sigF = 0.05)

#) 2) Delineamento em Blocos Completos Casualizados (D.B.C.)# ------------------------------------------------------# 1a Análise - tomando Linhas de 'exp1' como Blocosattach(exp1)help(str)is.factor(BLOCOS1)is.factor(TRAT)xyplot(PROD ~ TRAT, groups = BLOCOS1, data=exp1, jitter.x=TRUE)xtabs(PROD ~ TRAT + BLOCOS1) interaction.plot(TRAT, BLOCOS1, PROD)interaction.plot(BLOCOS1, TRAT, PROD)par(mfrow=c(1,2))

stripchart(PROD ~ TRAT, xlab="PROD", ylab="TRAT")stripchart(PROD ~ BLOCOS1, xlab="PROD", ylab="BLOCO")par(mfrow=c(1,1))

modelo2 <- aov(PROD ~ LINHA + TRAT, data=exp1) # comando p/ dados balanceados anova(modelo2)modelo2b <- lm(PROD ~ BLOCOS1 + TRAT, data=exp1) # comando mais geral (modelo linear)anova(modelo2b)par(mfrow=c(2,2)) # divide a matriz gráfica em duas linhas e duas colunas (recebe quatro gráficos)plot(modelo3) # diagnóstico inicial de problemas relacionados às pressuposições da ANOVA clássica layout(1) # retorna à configuração gráfica padrão (um gráfico por página)

# 2a Análise - tomando Colunas de 'exp1' como Blocosxtabs(PROD ~ TRAT + BLOCOS2) is.factor(BLOCOS2)modelo3 <- lm(PROD ~ BLOCOS2 + TRAT, data=exp1) anova(modelo3)par(mfrow=c(2,2)) # divide a matriz gráfica em duas linhas e duas colunas (recebe quatro gráficos)plot(modelo3) # diagnóstico inicial de problemas relacionados às pressuposições da ANOVA clássica layout(1) # retorna à configuração gráfica padrão (um gráfico por página) names(modelo3)x<-fitted.values(modelo3)y<-residuals(modelo3)plot(y~x)

# Testando Pressuposições da Análise:# Gráfico de resíduos:par(mfrow=c(2,2)) # divide a matriz gráfica em duas linhas e duas colunas (recebe quatro gráficos)plot(modelo3) # diagnóstico inicial de problemas relacionados às pressuposições da ANOVA clássica layout(1) # retorna à configuração gráfica padrão (um gráfico por página)

# Homocedasticia:bartlett.test(residuals(modelo3)~exp1$TRAT) # teste de homogeneidade das variâncias de TRATbartlett.test(PROD ~ BLOCOS2 + TRAT, data=exp1) # NÃO USARbartlett.test(PROD ~ TRAT + BLOCOS2, data=exp1) # NÃO USAR

#Normalidade:shapiro.test(residuals(modelo3)) # teste de normalidade dos resíduos

qqnorm(residuals(modelo3))qqline(residuals(modelo3))hist(residuals(modelo3))

# Teste dos desvios de aditividade (Tukey, 1949; citado por Parente, 1984, p. 17-21): # ------------(script J. B. Duarte / UFG)------------------------------------------- qij <- modelo3$fitted.values^2mod_NA <- lm(qij ~ BLOCOS2 + TRAT, data=exp1) # atualizar o modelo em caso de outro delineamentoA <- sum(mod_NA$residuals^2)B <- sum(modelo3$residuals*qij)SQ_NA <- (B^2)/A # SQ de não-aditividade, com um grau de liberdade (GL=1)QM_NA <- SQ_NA/1 # QM de não-aditividade (GL=1)SQ_Res <- sum(modelo3$residuals^2)-SQ_NA # SQ do Resíduo (descontando-se os desvios de não-aditividade)GL_Res <- modelo3$df.residual-1 # GL do Resíduo (descontando-se os desvios de não-aditividade)QM_Res <- SQ_Res/GL_ResF_NA <- QM_NA/QM_Res Pr_NA <- 1-pf(F_NA,1,GL_Res) result1 <- matrix(c(1,SQ_NA,QM_NA,F_NA,Pr_NA,GL_Res,SQ_Res,QM_Res,1,1),nrow=2,ncol=5,byrow=T,dimnames=list(c("N_Aditiv.","Residuals"),c("Df","Sum Sq","Mean Sq","F value","Pr(>F)")))result1 # Significância da fonte de variação "N_Aditiv.", com "1" grau de liberdade, indica desvios de pressuposição de aditividade do modelo.

# Usando o pacote "asbio": require(asbio)tukey.add.test(PROD, TRAT, BLOCOS2) # teste de aditividade

# Independência:library(lmtest)dwtest(PROD~BLOCOS2+TRAT,alternative = c("two.sided")) # teste de independência dos resíduos

#================================================================# CASO AS PRESSUPOSIÇÕES SEJAM ATENDIDAS => COMPARAÇÕES DE MÉDIAS

# Comparação de médias duas a duas (Testes de Comparações Múltiplas)library(gregmisc)TukeyHSD(aov(PROD ~ BLOCOS2 + TRAT, data=exp1), TRAT) # Teste Tukey (alpha=5%) - função para dados balanceadoscompara <- TukeyHSD(aov(PROD ~ BLOCOS2 + TRAT, data=exp1))

compara$TRAT

# Testes de comparação de médias pelo pacote 'agricolae'require(agricolae)glr <- df.residual(modelo3) # graus de liberdade do Resíduoqmr <- deviance(modelo3)/glr # quadrado médio do Resíduo# teste de Tukeycomp <-with(exp1, HSD.test(PROD, TRAT, DFerror=glr, MSerror=qmr, alpha=0.05))comp

# Obtendo os resultados com o pacote ExpDes (inclui Testes Comparações Múltiplas)require(ExpDes)require(doBy)with(exp1, rbd(treat=TRAT, block=BLOCOS2, resp=PROD, mcomp="tukey")) # Tukey (tradicional)with(exp1, rbd(treat=TRAT, block=BLOCOS2, resp=PROD, mcomp="duncan")) # Duncan (tradicional)with(exp1, rbd(treat=TRAT, block=BLOCOS2, resp=PROD, mcomp="snk")) # Student-Newman-Keulswith(exp1, rbd(treat=TRAT, block=BLOCOS2, resp=PROD, mcomp="lsd")) # t-Student (LSD)with(exp1, rbd(treat=TRAT, block=BLOCOS2, resp=PROD, mcomp="lsdb")) # t-Student / Bonferroniwith(exp1, rbd(block=BLOCOS2, treat=TRAT, resp=PROD, mcomp="sk")) # Scott-Knott (sem ambiguidades) - sob balanceamentowith(exp1, rbd(block=BLOCOS2, treat=TRAT, resp=PROD, mcomp="ccboot")) # TCM bootstrap

#) 3) Delineamento em Quadrado Latino (Q.L.)# --------------------------------------attach(exp1)xtabs(PROD ~ LINHA + COLUNA)matrix(TRAT,5,5)LINHA <- as.factor(LINHA)COLUNA <- as.factor(COLUNA)is.factor(LINHA)is.factor(COLUNA)is.factor(TRAT)modelo4 <- lm(PROD ~ LINHA + COLUNA + TRAT, data=exp1) # comando "lm" (neste caso poderia ser "aov" - dados balanceados)aov(PROD ~ LINHA + COLUNA + TRAT, data=exp1)anova(modelo4) names(modelo4)desv_padr_res <- sqrt(sum(modelo4$residuals^2)/modelo4$df.residual)desv_padr_res boxplot(residuals(modelo4))

hist(residuals(modelo4))plot(fitted(modelo4), residuals(modelo4)/desv_padr_res, xlab="Ajustados", ylab="Resíduos")abline(h=0)par(mfrow=c(2,2))plot(modelo4)layout(1)qqnorm(residuals(modelo4)/desv_padr_res)qqline(residuals(modelo4)/desv_padr_res)library(lmtest)dwtest(PROD~LINHA+COLUNA+TRAT,alternative = c("two.sided")) # teste de independência dos resíduosshapiro.test(residuals(modelo4)) # teste de normalidade dos resíduosshapiro.test(modelo4$residuals) # idem anteriorbartlett.test(residuals(modelo4)~exp1$TRAT) # teste de homogeneidade das variâncias de TRATbartlett.test(PROD ~ LINHA + COLUNA + TRAT, data=exp1) # NÃO UTILIZARbartlett.test(PROD ~ TRAT + LINHA +COLUNA, data=exp1) # NÃO UTILIZARrequire(asbio)tukey.add.test(PROD, TRAT, LINHA) # teste de aditividadetukey.add.test(PROD, TRAT, COLUNA) # teste de aditividade

#================================================================# CASO AS PRESSUPOSIÇÕES SEJAM ATENDIDAS => COMPARAÇÕES DE MÉDIAS

TukeyHSD(aov(PROD ~ LINHA + COLUNA + TRAT, data=exp1)) # Teste Tukey (alpha=5%) - para dados balanceadoscompara <- TukeyHSD(aov(PROD ~ LINHA + COLUNA + TRAT, data=exp1))compara$TRAT

# Testes de comparação de médias pelo pacote 'agricolae'require(agricolae)glr <- df.residual(modelo4) # graus de liberdade do Resíduoqmr <- deviance(modelo4)/glr # quadrado médio do Resíduo# teste de Tukeycomp<- with(exp1, HSD.test(PROD, TRAT, DFerror=glr, MSerror=qmr, alpha=0.05))comp

# Obtendo os resultados com o pacote ExpDes (inclui Testes Comparações Múltiplas)require(ExpDes)require(doBy)with(exp1,latsd(treat=TRAT,row=LINHA,col=COLUNA,resp=PROD, mcomp="tukey")) # Tukey (tradicional)

with(exp1,latsd(treat=TRAT,row=LINHA,col=COLUNA,resp=PROD, mcomp="duncan")) # Duncan (tradicional)with(exp1,latsd(treat=TRAT,row=LINHA,col=COLUNA,resp=PROD, mcomp="snk")) # Student-Newman-Keulswith(exp1,latsd(treat=TRAT,row=LINHA,col=COLUNA,resp=PROD, mcomp="lsd")) # t-Student (LSD)with(exp1,latsd(treat=TRAT,row=LINHA,col=COLUNA,resp=PROD, mcomp="lsdb")) # t-Student / Bonferroniwith(exp1,latsd(treat=TRAT,row=LINHA,col=COLUNA,resp=PROD, mcomp="sk")) # Scott-Knott (sem ambiguidades) - sob balanceamento

# ===========================================================# Exportando resultados para arquivos texto (explorar melhor)# ----------------------------------------------------------help(write.table)tabela<- matrix(compara)write(compara,"arquivo.txt")write.table(tabela, file="tabela1.txt", sep=",", row.names=FALSE)write.table(t(tabela), file="tabela2.txt",sep=",", row.names=FALSE)sink('tabela3.txt')tabelasink()help(sink)sink(file="teste.txt", split=T)

# Outros pacotes úteis neste contexto (instalar e carregar)require(contrast)require(multicomp)require(lattice)require(latticeExtra)

# --o--FIM--o--

SCRIPT “R” – ANOVA (ex. DBC), VERIFICAÇÃO DE PRESSUPOSTOS E TRANSFORMAÇÃO DE DADOS

# ==============================================================================# ANÁLISE DE VARIÂNCIA - D.B.C (Exercício) c/ Verificação de Suposições da ANOVA# ==============================================================================

sink(file="Saida.txt", split=T) # armazenará os resultados das operações no R, exceto gráficos, num arquivo texto; neste caso, no mesmo diretório)# Buscando dados de outro arquivo (neste caso, gerado como CSV no Excel)exp1 <- read.table (file="C:/Users/Joao/Desktop/Biometria/2016/Verif_sup_ANOVA_soja_nemat.csv", sep=";", header=TRUE) exp1head(exp1)plot(exp1)plot(Y ~ Trat, exp1)summary(exp1)str(exp1)

# PREPARANDO PARA UMA ANOVA - declarando Fatores# ==============================================attach(exp1)exp1$Trat <- factor(exp1$Trat)exp1$Bloco <- factor(exp1$Bloco)Trat <- as.factor(Trat)Bloco <- as.factor(Bloco)

# ANOVA D.B.C. / MODELO FIXO# ==========================#) Delineamento em Blocos Completos Casualizados (D.B.C.)# ------------------------------------------------------is.factor(Bloco)is.factor(Trat)class(Bloco)class(Trat)xtabs(Y ~ Trat + Bloco) modelo2 <- lm(Y ~ Bloco + Trat, data=exp1) modelo2names(modelo2)summary(modelo2)anova(modelo2)

library(nortest) # Pacote para testes de normalidadelibrary(normtest) # Pacote para testes de normalidadelibrary(asbio) # Pacote útil para o teste de aditividadelibrary(dae) # Pacote útil para o teste de aditividade (exclusivo p/ "aov" -> problemas)library(lmtest) # Pacote útil para o teste de independêncialibrary(Rcmdr) # Pacote útil para teste de homocedasticidadelibrary(MASS) library(agricolae)library(ExpDes)require(doBy)

# Verificação das suposições da ANOVA#------------------------------------modelo2$fitted.values # valores ajustados (preditos pelo modelo)xtabs(modelo2$fitted.values ~ Trat + Bloco) modelo2$residualsxtabs(modelo2$residuals ~ Trat + Bloco) residuos_padronizados <- modelo2$residuals/sqrt(sum(modelo2$residuals^2)/modelo2$df.residual)xtabs(residuos_padronizados ~ Trat + Bloco) xtabs(rstandard(modelo2) ~ Trat + Bloco) plot(modelo2$fitted.values,residuos_padronizados) # gráfico de resíduos (permite avaliar independência e homocedasticidade)hist(modelo2$residuals)qqnorm(modelo2$residuals)qqline(modelo2$residuals)qqnorm(rstandard(modelo2))qqline(rstandard(modelo2))par(mfrow=c(2,2))plot(modelo2) layout(1)shapiro.test(residuals(modelo2)) # teste de normalidade dos resíduosrequire(nortest) # Pacote para testes de normalidadelillie.test(modelo2$residuals) ad.test(modelo2$residuals)cvm.test(modelo2$residuals)pearson.test(modelo2$residuals)sf.test(modelo2$residuals)

bartlett.test(residuals(modelo2)~exp1$Trat) # teste de homogeneidade das variâncias de TRAT

help(dwtest)library(lmtest) # Pacote útil para o teste de independênciadwtest(Y ~ Bloco + Trat, data=exp1) # teste de independência dos resíduos (default: autocorrelação nula (H0) vs. positiva (H1) - unilateral)dwtest(Y ~ Bloco + Trat, alternative = c("greater")) # teste de independência - idem anterior: H0 (autocorrelação nula) vs. H1 (autocorrelação positiva)dwtest(Y ~ Bloco + Trat, alternative = c("less")) # teste de independência - H0 (autocorrelação nula) vs. H1 (autocorrelação negativa)dwtest(Y ~ Bloco + Trat, alternative = c("two.sided")) # teste de independência - H0 (autocorrelação nula) vs. H1 (autocorrelação não nula)=> PREFERÍVEL AQUI.

library(asbio) # Pacote útil para o teste de aditividadetukey.add.test(Y, Trat, Bloco) # teste de aditividade (1a. opção)

# Teste dos desvios de aditividade (Tukey, 1949; citado por Parente, 1984, p. 17-21): # ------------(script J. B. Duarte / UFG)------------------------------------------- qij <- modelo2$fitted.values**2mod_NA <- lm(qij ~ Bloco + Trat, data=exp1) # atualizar o modelo em caso de outro delineamentoA <- sum(mod_NA$residuals^2)B <- sum(modelo2$residuals*qij)SQ_NA <- (B^2)/A # SQ de não-aditividade, com um grau de liberdade (GL=1)QM_NA <- SQ_NA/1 # QM de não-aditividade, com um grau de liberdade (GL=1)SQ_Res <- sum(modelo2$residuals^2)-SQ_NA # SQ do Resíduo (descontando-se os desvios de não-aditividade)GL_Res <- modelo2$df.residual-1 # GL do Resíduo (descontando-se os desvios de não-aditividade)QM_Res <- SQ_Res/GL_ResF_NA <- QM_NA/QM_Res Pr_NA <- 1-pf(F_NA,1,GL_Res) result1 <- matrix(c(1,SQ_NA,QM_NA,F_NA,Pr_NA,GL_Res,SQ_Res,QM_Res,1,1),nrow=2,ncol=5,byrow=T,dimnames=list(c("N_Aditiv.","Residuals"),c("Df","Sum Sq","Mean Sq","F value","Pr(>F)")))result1 # Significância da fonte de variação "N_Aditiv.", com "1" grau de liberdade, indica desvios de pressuposição de aditividade do modelo.

# Outra opção para testar não-aditividade (exclusiva para D.B.C): # (Fonte: Prof. Ribeiro Jr, P.J.- UFPR. <http://leg.ufpr.br/Rpira/Rpira/node15.html>) modelo2$coeffbl <- c(0, modelo2$coeff[2:5]) # atualizar [2:5] para selecionar os efeitos de blocos (listados primeiramente no modelo)tr <- c(0, modelo2$coeff[6:10]) # atualizar [6:10] para selecionar os efeitos de tratamentos (listados após blocos no modelo)bltr <- rep(tr,5) * rep(bl, rep(6,5)) # atualizar (tr,5) para o número de blocos e (6,5) para tratamentos e blocos, respectivamente.ttna <- update(modelo2, .~. + bltr) # acrescenta a F.V. associada à não-aditividadeanova(ttna) # Significância da F.V. "bltr" (não-aditividade), com "1" GL, indica desvios de aditividade.

# Caso necessite: busque uma transformação - ex. Box & Cox (1964):# ---------------------------------------------------------------require(MASS) boxcox(Y ~ Bloco+Trat, data=exp1, plotit=T)boxcox(Y ~ Bloco+Trat, data=exp1, lam=seq(-1, 1, 1/10))boxcox(Y ~ Bloco+Trat, data=exp1, lam=seq(-0.75, 0.25, 1/20)) # lâmbida=-0,20 (Transf. Box & Cox)MASS::boxcox(modelo2) # dispensa o carregamento do pacote todo (MASS) e executa apenas a sua função 'boxcox'MASS::boxcox(modelo2, lam=seq(-0.75, 0.25, 1/20)) # como o IC(lamb) inclui o "zero", a transformação logaritmica também é aplicável

# ANOVA com dados transformados - ex. (Box & Cox, 1964):# -----------------------------------------------------head(exp1)Yt <- ((Y**(-0.2))-1)/(-0.2) # Neste caso, a potência escolhida, que maximiza a verossimilhança, é -0,2.Yt <- Y**(-0.3)Yt <- log(Y); Yt3 <- log10(Y); Ytxtabs(Yt ~ Trat + Bloco) modelo2 <- lm(Yt ~ Bloco + Trat, data=exp1) modelo2summary(modelo2)anova(modelo2)

# Verificação das suposições da ANOVA#------------------------------------modelo2$fitted.values # valores ajustados (preditos pelo modelo)modelo2$residualsresiduos_padronizados <- modelo2$residuals/sqrt(sum(modelo2$residuals^2)/modelo2$df.residual)plot(modelo2$fitted.values,residuos_padronizados) # gráfico de resíduos (permite avaliar independência e homocedasticidade)hist(modelo2$residuals)par(mfrow=c(2,2))plot(modelo2)layout(1)shapiro.test(residuals(modelo2)) # teste de normalidade dos resíduosbartlett.test(residuals(modelo2)~exp1$Trat) # teste de homogeneidade das variâncias de TRATdwtest(Yt~Bloco+Trat,alternative = c("two.sided")) # teste de independência dos resíduostukey.add.test(Yt, Trat, Bloco) # teste de aditividade

#NOTA: No exemplo em estudo (contagem de nematóides em soja), a transformação Box & Cox melhora simultaneamente as condições # de normalidade, homocedasticidade e independência dos resíduos, além da pressuposição de aditividade. # É importante salientar que, neste caso, transformações como Y'= Log(Y) ou Y'=Y**(-0,30), esta última resultante# da busca de transformação "estabilizadora da variância", via regressão de Ln(Var(i))= a + b.Ln(Med(i)), também # resolveriam bem os problemas identificados de violação das pressuposicões da ANOVA clássica.

# =======================================================================================================# Quando os resultados indicarem que os pressupostos estão obedecidos, a análise de variância será válida. # =======================================================================================================# Se detectado efeito significativo de tratamentos e estes forem qualitativos, deve-se seguir a análise com um teste de # comparação múltipla - TCM (ex. teste Tukey, Duncan, entre outros). Se os tratamentos forem quantitativos, não se # recomenda TCM; mas sim, a análise de regressão.

#======================# COMPARAÇÃO DE MÉDIAS modelo2b<-aov(Yt~Bloco+Trat)anova(modelo2b)TukeyHSD(modelo2b)result<-TukeyHSD(modelo2b)plot(result)result<-TukeyHSD(modelo2b,"Trat",ordered=T)plot(result)

require(agricolae)xtabs(~Trat, na.omit(exp1)) # repetições por nível de TRATgle <- df.residual(modelo2) # graus de liberdade do Resíduoqme <- deviance(modelo2)/gle # Quadrado médio do Resíduocomp<-with(exp1, HSD.test(Yt, Trat, DFerror=gle, MSerror=qme, alpha=0.05))comphelp(agricolae) # vide 'index'

require(ExpDes)with(exp1,rbd(block=Bloco,treat=Trat,resp=Yt,mcomp="tukey")) # Tukey (tradicional)with(exp1,rbd(block=Bloco,treat=Trat,resp=Yt,mcomp="snk")) # Student-Newman-Keulswith(exp1,rbd(block=Bloco,treat=Trat,resp=Yt,mcomp="sk")) # Scott-Knott (sem ambiguidades) - sob balanceamento

# Opções de testes livre de distribuições (bootstrap (Ramos & Ferreira, 2009: Revista Ceres, v.56, p.140-149)ccboot(y=Yt,trt=Trat,DFerror=gle,SSerror=deviance(modelo2),alpha=0.05,group=TRUE,main=NULL,B=1000)

rbd(Trat, Bloco, Y, quali = TRUE, mcomp='ccboot', sigF = 0.05)help(ExpDes) # vide 'index'

#======================detach(exp1) # Encerra a análises, desanexando o objeto do caminho de procura - "detach()".

#===============================================# Consulte agora o arquivo "Saída.txt", # salvo no mesmo diretório de trabalho, # que contém resultados das análises realizadas. #===============================================

# --o--FIM--o--

ROTINA “SAS” PARA ANOVA, VERIFICAÇÃO DE PRESSUPOSTOS E BUSCA DE TRANSFORMAÇÃO

options nodate ls=80;data NEMAT_SOJA; input TRAT BLOCO Y; Yt1=((Y**(-0.2))-1)/(-0.2); /* Transf. Box & Cox */ Yt2=Y**(1-2.6/2); /* Transf. Estabilizadora da Variância */ Yt3=log10(Y); Yt4=log(Y); /* Transf. Logaritmica (Yt3: base 10 e Yt4: base 'e' - neperiano) */datalines;1 1 503.232 1 164.503 1 1272.414 1 724.925 1 192.006 1 742.571 2 1338.982 2 283.333 2 200.004 2 851.715 2 140.466 2 954.971 3 1850.132 3 581.823 3 516.924 3 7466.675 3 444.446 3 1647.941 4 2013.162 4 345.763 4 658.074 4 937.245 4 561.806 4 562.701 5 2236.422 5 1168.453 5 4326.044 5 815.535 5 345.716 5 235.60;run;

title 'VERIFICAÇÃO DE PRESSUPOSIÇÕES ASSOCIADAS À ANÁLISE DE VARIÂNCIA';

title2 'ANOVA - DELINEAMENTO BLOCOS CASUALIZADOS';proc glm data= NEMAT_SOJA outstat=FV_ANOVA; class BLOCO TRAT; model Y= BLOCO TRAT /p solution; output out=RESIDUOS predicted=PRED residual=RESID student=RES_PADR;run; quit;/* Uma vez gerados os valores preditos e os respectivos resíduos (opção 'p' do comando 'model'), os quais são armazenados num arquivo de saída (comando 'output ... out=NOME'), a sequência a seguir pode ser implementada para quaisquer outros delineamentos. */ title2 'Teste de Aditividade do modelo (Tukey, 1949)';data ADITIV1; set RESIDUOS; qij=PRED*PRED; bij=RESID*qij;run;proc glm data=ADITIV1 noprint; class BLOCO TRAT; model qij= BLOCO TRAT /p; output out=RES_ADIT predicted=P_NA residual=R_NA;run; quit;data ADITIV2; set RES_ADIT; aij=R_NA*R_NA;run;proc means noprint data=ADITIV2; var aij bij; output out=N_ADIT;run;data N_ADIT (keep=_SOURCE_ DF SS MS); set N_ADIT; if _STAT_^='MEAN' then delete; _SOURCE_='Non-Additivity'; A=aij*_FREQ_; B=bij*_FREQ_; DF=1; SS=B**2/A; MS=SS/1;run;data FV_ANOVA2 (keep=_SOURCE_ DF SS MS); set FV_ANOVA; MS=SS/DF; if _SOURCE_^='ERROR' then delete;run;data ADITIVIDADE (drop=_SOURCE_); set FV_ANOVA2 N_ADIT;run;proc iml; use work.ADITIVIDADE; setin work.ADITIVIDADE; read all var {DF SS MS} into FV; df_res=FV[1,1]-FV[2,1]; SS_res=FV[1,2]-FV[2,2]; MS_res=SS_res/df_res; F_na=FV[2,2]/MS_res; p_na=1-probf(F_na,1,df_res); x1={. .}; x2=F_na||p_na; x=x1//x2;

FV1=FV||x; FV2=df_res||SS_res||MS_res||x1; Source=FV1//FV2; OBS={'Error' 'Non-Additivity' 'Residual'}; Item={1,2,3}; VAR={'DF' 'SS' 'MS' 'F' 'Pr>F'}; print ' ' Source [rowname=OBS colname=VAR format=12.4]; quit;

title2 'Avaliação de Normalidade dos resíduos'; proc univariate data=RESIDUOS normal plot; var RESID; run;proc chart data=RESIDUOS; vbar RESID;run;

title2 'Gráfico de Resíduos - Independência e Homocedasticidade'; proc plot data=RESIDUOS; plot RESID*PRED;run;/* Estas verificações também podem ser feitas utilizando-se o módulo interativo do SAS, por meio dos seguintes passos (a partir do 'Menu' superior): Solutions->Analysis->Interactive Data Analysis->WORK->RESIDUOS->Open->Analyze-> ->Distribution(Y)->...->(Normal QQ Plot / Tests for Normality), para avaliação da normalidade dos erros. Seguindo-se, no mesmo módulo: Analyze->Scatter Plot(YX), para a construção de 'Gráficos de Resíduos' e avaliação de independência e homocedasticidade. */

title2 'Independência dos resídudos - teste de Durbin-Watson';proc autoreg data=RESIDUOS; model RESID = / dw=4 dwprob;run;

title2 'Homocedasticidade - F máx (teste de Hartley)';proc sort data=RESiDUOS; by TRAT;run;

proc means noprint data=RESIDUOS; var Y; /* atualizar esta variável para avaliar eficiência de transformações */ by TRAT; output out=VAR_MED mean=MED_T var=VAR_T;run;data F_MAX; set VAR_MED;run;proc print data=F_MAX; run;proc means noprint data=F_MAX; var VAR_T; output out=MAX_MIN min=V_min max=V_max;run;data HARTLEY; set MAX_MIN; F_max=V_max/V_min;run; proc print data=HARTLEY; var V_min V_max F_max;run; /* O valor de F_max obtido deve ser avaliado na distribuição de Pearson & Hartley: se o valor F_max > F(alpha,n_trat,r-1) => rejeita-se H0 (hipótese de nulidade), de que as variâncias são homogêneas; caso contrário, não se rejeita H0, assumindo-se que os erros são homocedásticos. */

title2 'Homocedasticidade - Teste de Levene';data RES_ABS; set RESIDUOS; RESID_AB=abs(RESID);run;proc glm data=RES_ABS; class TRAT; model RESID_AB= TRAT;run; title;/* No teste de Levene, a signficância do teste F para tratamentos, numa ANOVA aplicada aos valores absolutos dos resíduos estimados (êij), indica que os erros não possuem homogeneidade de variâncias nos diferentes tratamentos; caso contrário, assume-se homocedasticidade. */

title2 'Buscando Transformação Estabilizadora da Variância';data TRANSF; set VAR_MED; LMED=log(MED_T); LVAR=log(VAR_T);run;proc reg data=TRANSF; model LVAR=LMED; run; title;/* A transformação que estabiliza a variância é dada por Y^(lambda), sendo lambda = 1-(b/2), em que 'b' é o coeficiente de regressão estimado da função: Log(Var)= a + b.Log(Med), i.e. Y = a + b.X. As estimativas dos parâmetros 'a' e 'b' são fornecidos como saída da última rotina ('proc reg').

Outra opção é buscar a transformação potência ótima, de Box & Cox (1964). Neste caso, pode-se usar: Solutions->Analysis->Guided Data Analysis->(dados, modelo)->Analyze->Assumptions->Response scaling,para escolher a melhor transformação (ex. Optimal Power Transformation).

A geração de variável transformada (Yt) no arquivo de dados pode ser feita diretamente na rotina de leitura dos dados,incluindo-se o(s) comando(s) relativo(s) à transformação escolhida (isto implica em rodar novamente essa rotina, após inserir a instrução de transformação); ou pode ser feita criando-se um novo arquivo a partir do arquivo original, conforme instruções a seguir: */

data NEMAT_SOJA; set NEMAT_SOJA; Yr=sqrt(Y); /* Ex. Transf. Raiz quadrada' <=> Yt=Y**0.5 */ Ya=arsin(sqrt(P/100)); /* Ex. Transf. Angular = Arco_seno[raiz(P/100)], com P em %. */run;/* Feito isto, reanalisa-se os dados (agora transformados) e avalia-se novamente as pressuposições. Tal verificação deve ser feita até que se encontre uma transformação apropriada (caso não seja possível, deve-se buscar procedimento alternativo à ANOVA clássica). */run;