Análise de variância

A Análise de variância consiste em decompor a variação total das observações do experimento em partes que podem ser atribuídas a causas controladas (conhecidas) e em partes a causas não controladas e/ou não controláveis (desconhecidas), o erro ou resíduo. O erro ou resíduo pode ocorrer em função do material que se está trabalhando ou em função do ambiente em que o experimento é conduzido. Outra fonte de erro pode ser a maneira como o experimento é conduzido pelo experimentador.

VAR. TOTAL = VAR. CONTROLADA + VAR NÃO CONTROLADA

Pressupostos fundamentais da análise de variância

A análise de variância (ANOVA) é um processo de análise baseado na decomposição da variação total, existente entre uma série de observações, em partes que podem ser atribuídas a causas conhecidas (ex. tratamentos e blocos) e numa parte devido a causas desconhecidas (erro experimental ou resíduo).

Para que a análise de variância seja considerada válida, algumas pressuposições devem ser obedecidas:

  • Os efeitos principais devem ser aditivos.

Nos experimentos, cada observação segue um modelo matemático denominado modelo linear aditivo,

\(y_{ij}=\mu+\tau_{ij}+\epsilon_{ij}\)

ou seja, cada um dos efeitos que compõem o modelo devem se somar. Na maior parte das situações, essa exigência é satisfeita. Quando não, as principais conseqüências dizem respeito ao nível de significância que se pensa estar testando. Esse inconveniente pode ser contornado por meio de uma transformação da variável resposta, por exemplo, com o uso de logaritmos.

  • Os erros das observações devem ser independentes.

Cada observação possui um erro que deve ser independente dos demais, seja das observações referentes ao mesmo tratamento, seja de observações referentes a outros tratamentos. Isso implica que os efeitos de tratamentos devem ser independentes, ou seja, que não ocorra correlação (associação) entre eles ou ainda, que a

cov(\(\epsilon_{ij}\) ; \(\epsilon_{ij}\)’) = 0

A independência dos resíduos, pode ser avaliada através de um gráfico dos resíduos vs valores preditos. Se existir independência, haverá ausência de padrão neste gráfico.

Se existir o registro da ordem de obtenção dos dados, fazendo um gráfico dosresíduos vs valores das observações ordenadas, pode-se detectar dependência entre as observações. Essa dependência aparece sob a forma de algum padrão de distribuição dos erros.

  • Os erros devem ser normalmente distribuídos.

Os erros devem possuir uma distribuição normal de probabilidade. Para que a análise de variância seja considerada válida, os erros devem ser originários da mesma população. Isso implica que todos os erros referentes às observações possuam a mesma distribuição de probabilidade.

Também aqui, pode-se utilizar o recurso da transformação dos dados para corrigir esse inconveniente.

A normalidade dos resíduos pode ser verificada através das seguintes formas:

  1. Teste de Shapiro-Wilk

O objetivo do teste é fornecer uma estatística de teste para avaliar se uma amostra tem distribuição Normal. O teste pode ser utilizado para amostras de qualquer tamanho.

Hipóteses:

\(H_{0}\) : podemos afirmar que a variável (os resíduos) possuem distribuição Normal.

\(H_{1}\) : podemos afirmar que a variável (os resíduos) não tem distribuição Normal.

A estatística W do teste é avaliada através do P-valor. No caso de um valor significativo, isto indica falta de normalidade para a variável aleatória analisada (Rejeita \(H_{0}\)).

  1. Teste de Lilliefors

O teste de Lilliefors, é um teste não paramétrico, usado para testar a normalidade de variáveis onde \(\mu\) e \(\sigma²\) não são conhecidos mas, estimados a partir dos dados amostrais. Aqui, ele pode ser utilizado para avaliar se os resíduos de um experimento podem ser ou não considerados como provenientes de uma distribuição Normal.

Através do teste de Lilliefors estaremos testando a seguinte hipótese a respeito dos resíduos:

\(H_{0}\) : podemos afirmar que a variável (os resíduos) possuem distribuição Normal.

\(H_{1}\) : podemos afirmar que a variável (os resíduos) não tem distribuição Normal.

Se \(D_{calc}\) > \(D_{tab}\), rejeita-se \(H_{0}\),ou seja, podemos afirmar que os resíduos não têm distribuição Normal.

  1. Gráfico Normal de Probabilidade (Normal probability plot)

O objetivo deste gráfico é avaliar o ajuste dos resíduos à distribuição Normal.

Esta análise gráfica consiste em se colocar em papel de normalidade (logNormal) os \(\epsilon_{ij}\)(\(P_{k}\) 100), onde \(P_{k}=\frac(k-{(\frac{1}{2})}{N}\).

Em que,

k = número de ordem do erro; N = número total de observações; \(P_{k}\)=probabilidade acumulada dos resíduos

No gráfico os valores devem formar uma linha reta se estes seguem uma distribuição normal. A maioria dos dados deve estar concentrada no meio da reta para que possamos considerar que os dados possuam uma distribuição Normal.

Os valores das caudas da distribuição não devem ser considerados com rigor.Mas, pontos extremos podem indicar a presença de outliers (pontos discrepantes).

  • Os erros devem ter variância comum (homocedasticidade).

Os componentes do erro devem ser todos estimados de uma mesma população. Isso implica que cada tratamento deve ter aproximadamente a mesma variância para que os testes da análise de variância tenham validade. Essa pressuposição pode ser testada por meio do teste de Bartley, de F máximo ou F de Hartley. Ou também pode ser verificada por um boxplot.

Quando as variâncias não são homogêneas, diz-se que existe heterocedasticidade.

A variabilidade entre repetições de um mesmo tratamento deve ser semelhante aos outros tratamentos. Isso pode ser verificado:

  1. Por meio de uma análise gráfica (mais usual):
  1. Fazer um Box-Plot para os tratamentos vs resíduos.
    • Se existe homocedasticidade, esperamos que os Box-plots sejam semelhantes, ou seja, a variabilidade é a mesma em todas as caixas.
    • Se existe heterocedasticidade, a variabilidade é diferente entre as caixas. As vezes, a heterocedasticidade pode ser também um indício de falta de normalidade.
  2. Através de um gráfico de dispersão para os tratamentos vs resíduos.
    • Se existe homocedasticidade, espera-se que os desvios variem de forma homogênea, dentro de uma mesma amplitude.
    • Se existe heterocedasticidade, os desvios não variam dentro da mesma amplitude. Isso pode ocorrer de forma homogênea, que é chamada de heterocedasticidade regular ou de forma heterogênea, heterocedasticidade irregular. Quando ocorre heterocedasticidade regular, há também uma associação entre as médias dos tratamentos e os resíduos, pode acontecer que quando aumentam as médias, os erros também aumentam e vice-versa.

O que fazer quando ocorre heterocedasticidade ?

Uma das saidas mais simples é através da transformação de dados do tipo seno, cosseno ou raiz.

O uso de uma transformação, um artifício matemático, é aplicável quando existe uma relação conhecidade entre média e variância (heterocedasticidade regular). Quando existe uma heterocedasticidade irregular, as transformações não resolvem o problema.

Veremos mais adiante como escolher uma transformação adequada.

  1. Através de testes

Alguns testes para a verificação da homocedasticidade:

  • Hartley: exige um número igual de repetições entre os tratamentos.

O teste de Hartley também é conhecido como F máximo. A estatística do teste é dada por:

\(F_{max}=\frac{\sigma^2_{max}}{\sigma^2_{min}}\)

Em que:

\(sigma^2_{max}\)=tratamento que apresentou a maior variância

\(sigma^2_{min}\)=tratamento que apresentou a menor variância

Se a razão for muito alta, tem tratamentos com variâncias muito diferentes.

Hipótese testada:

\(H_{0}= Há\ homocedasticidade\)

\(H_{1}= Não\ há\ homocedasticidade\)

O valor calculado de \(F_{max}\) é comparado com o valor tabelado para \(H_{(g, r-1)}\) da tabela de PEARSON & HARTLEY (g=número de grupos ou tratamentos e r=número de repetições).

Se \(F{max} > H_{(g,r-1) \alpha}\) rejeita-se a hipótese de homocedasticidade e conclui-se que não existe homogeneidade de variâncias entre os tratamentos. Caso contrário,aceitamos \(H_{0}\).

  • Cochran: pode ser usado com número de repetições diferentes por tratamento.

  • Bartlett: pode ser usado com número de repetições diferentes por tratamento.

O teste de Bartlett testa a seguinte hipótese:

\(H_{0}: \sigma²_{1}=\sigma²_{2}=...=\sigma²_{i}\)

\(H_{1}: \sigma²_{1}\neq \sigma²_{i'}\ para\ pelo\ menos\ um\ par\ i \neq i'\)

A hipótese de variãncias iguais será rejeitada se \(U>\chi_{(\alpha, t-1)}\)

  • Levene: Anova para os resíduos. Pode levar a resultados conflitantes.

    As pressuposições em relação aos erros podem ser resumidas na seguinte expressão:

\(\epsilon_{ij}\) \(\stackrel{\rm IID}{=}\) N(0,\(\sigma²\)), significa erros independente e identicamente, distribuídos, com distribuição normal

Todas essas pressuposições podem ser verificadas, na prática, realizando-se a análise de resíduos.

Identificação de Outliers ou Dados Discrepantes

A identificação de outilers ou dados discrepantes faz parte da análise exploratória de dados através de elementos da estatística descritiva.

Em planejamento de experimentos, em função da sua natureza, existem alguns procedimentos específicos que auxiliam na verificação da existência de outliers. Um deles, através da análise de resíduos, podemos construir um gráfico dos valores preditos vs resíduos padronizados.

Se os erros têm distribuição N(0;\(\sigma²\)), pode-se esperar que a média \(\pm\) 1 contém 68% dos dados, a média \(\pm\) 2 contém 95% dos dados e a média \(\pm\) 3 contém 99% dos dados.

Desta forma, podem ser considerados outliers, valores que ultrapassem \(\pm 3 \sigma\). Na verdade, estatisticamente, esta conclusão é válida mas, na prática, é o pesquisador quem determina se o outlier pode realmente ser assim considerado. Pois, os outliers podem fornecer informações importantes sobre o experimento, como problemas de condução e execução do experimento; novos fatos relevantes e não explorados pelo pesquisador e estatisticamente, podem revelar que outra distribuição possa explicar melhor o comportamento dos dados.

Delineamento Inteiramente casualizado

Modelo Matemático

As observaçõeses oriundas deste delineamento se adequam ao modelo matemático:

\(y_{ij}=\mu+\tau_{ij}+\epsilon_{ij}\)

onde:

\(y_{ij}\) = valor de uma observação correspondente ao i-ésimo tratamento na jésima parcela.

\(\tau_{ij}\)= efeito do i-ésimo tratamento.

\(\epsilon_{ij}\) = erro experimental associado ao iésimo tratamento na j-ésima parcela.

O esquema do quadro de análise de variâcia é dado da seguinte forma:

Fator de Variação GL Soma Quadrados Quadrado Médio Teste F
Tratamentos (k-1) SQTr QMTr QMTr/QMR
Resíduo (n-k) SQR QMR
Total (n-1) SQT

** Atividade Prática do Delineamento Inteiramente Casualizado

Nesta sessão iremos usar os dados do exemplo abaixo, para analizar um experimento em delineamento inteiramente casualizado.

Exemplo de análise de experiemnto em del. Int. ao acaso

Considere o seguinte experimento que foi conduzido considerando um delineamento inteiramente casualizado. Foram comparadas 9 linhagens de fungos medindo-se as taxas de crescimento em micras/hora.

R e p e t i
Linhagens I II III IV V VI Total
L1 385 323 417 370 437 340 2272
L2 406 385 444 443 474 437 2589
L3 354 292 389 312 432 299 2078
L4 271 208 347 302 370 264 1762
L5 344 292 354 354 401 306 2051
L6 354 354 410 453 448 417 2436
L7 167 115 194 130 240 139 985
L8 344 385 410 437 437 410 2423
L9 385 385 396 453 458 417 2494
Total 19090
ex01<-read.table("/home/roberval/Documentos/UFAM-2016-Mestrado_Estat_Exp/Notas_de_Aulas-2016/Aula5-del_Int_acaso/dados/exemplo01.txt", header=TRUE,)
is.data.frame(ex01)
## [1] TRUE
names(ex01)
## [1] "trat" "resp"
ex01$resp
##  [1] 385 323 417 370 437 340 406 385 444 443 474 437 354 292 389 312 432
## [18] 299 271 208 347 302 370 264 344 292 354 354 401 306 354 354 410 453
## [35] 448 417 167 115 194 130 240 139 344 385 410 437 437 410 385 385 396
## [52] 453 458 417
ex01$trat
##  [1] t1 t1 t1 t1 t1 t1 t2 t2 t2 t2 t2 t2 t3 t3 t3 t3 t3 t3 t4 t4 t4 t4 t4
## [24] t4 t5 t5 t5 t5 t5 t5 t6 t6 t6 t6 t6 t6 t7 t7 t7 t7 t7 t7 t8 t8 t8 t8
## [47] t8 t8 t9 t9 t9 t9 t9 t9
## Levels: t1 t2 t3 t4 t5 t6 t7 t8 t9
is.factor(ex01$trat)
## [1] TRUE
is.numeric(ex01$resp)
## [1] TRUE

Portanto, o objeto é um data.frame com duas variáveis, sendo uma delas um fator (a variáavel trat) e a outra uma variável numérica.

Nota: na ANOVA, as variáveis independentes precisam possuir a característica “factor”. Caso contrário, o R realizará uma análise de regressão entre as variáveis.

Vamos agora fazer uma rápida análise descritiva:

summary(ex01)
##       trat         resp      
##  t1     : 6   Min.   :115.0  
##  t2     : 6   1st Qu.:307.5  
##  t3     : 6   Median :377.5  
##  t4     : 6   Mean   :353.5  
##  t5     : 6   3rd Qu.:417.0  
##  t6     : 6   Max.   :474.0  
##  (Other):18
tapply(ex01$resp, ex01$trat, mean)
##       t1       t2       t3       t4       t5       t6       t7       t8 
## 378.6667 431.5000 346.3333 293.6667 341.8333 406.0000 164.1667 403.8333 
##       t9 
## 415.6667

Há um mecanismo no R de “anexar” objetos ao caminho de procura que permite economizar um pouco de digitação.

attach(ex01)
search()
##  [1] ".GlobalEnv"        "ex01"              "package:stats"    
##  [4] "package:graphics"  "package:grDevices" "package:utils"    
##  [7] "package:datasets"  "package:methods"   "Autoloads"        
## [10] "package:base"
tapply(resp, trat, mean)
##       t1       t2       t3       t4       t5       t6       t7       t8 
## 378.6667 431.5000 346.3333 293.6667 341.8333 406.0000 164.1667 403.8333 
##       t9 
## 415.6667

Pode-se “desanexar” o objeto com os dados (embora isto não seja obrigatório) com o comando.

detach(ex01)

Prosseguindo com a análise exploratória:

ex01.m <- tapply(resp, trat, mean)
ex01.m
##       t1       t2       t3       t4       t5       t6       t7       t8 
## 378.6667 431.5000 346.3333 293.6667 341.8333 406.0000 164.1667 403.8333 
##       t9 
## 415.6667
ex01.v <- tapply(resp, trat, var)
ex01.v
##       t1       t2       t3       t4       t5       t6       t7       t8 
## 1916.267  987.500 3117.867 3494.667 1513.767 1903.600 2173.367 1242.167 
##       t9 
## 1091.067
plot(ex01) # ou plot(trat, resp)
points(ex01.m, pch="x", col=2, cex=1.5)

O comando aov, realiza a análise dos dados do data.frame. Compare a saída do comando aov com o comando ANOVA.

ex01.av <- aov(resp ~ trat, data = ex01)
ex01.av
## Call:
##    aov(formula = resp ~ trat, data = ex01)
## 
## Terms:
##                     trat Residuals
## Sum of Squares  332918.1   87201.3
## Deg. of Freedom        8        45
## 
## Residual standard error: 44.02053
## Estimated effects may be unbalanced
summary(ex01.av)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## trat         8 332918   41615   21.48 5.44e-13 ***
## Residuals   45  87201    1938                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(ex01.av)
## Analysis of Variance Table
## 
## Response: resp
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## trat       8 332918   41615  21.475 5.445e-13 ***
## Residuals 45  87201    1938                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Portanto o objeto ex01.av guarda os resultados da análise de variância e outras informações.

names(ex01.av)
##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "contrasts"     "xlevels"       "call"          "terms"        
## [13] "model"
ex01.av$coef
## (Intercept)      tratt2      tratt3      tratt4      tratt5      tratt6 
##   378.66667    52.83333   -32.33333   -85.00000   -36.83333    27.33333 
##      tratt7      tratt8      tratt9 
##  -214.50000    25.16667    37.00000
ex01.av$res
##          1          2          3          4          5          6 
##   6.333333 -55.666667  38.333333  -8.666667  58.333333 -38.666667 
##          7          8          9         10         11         12 
## -25.500000 -46.500000  12.500000  11.500000  42.500000   5.500000 
##         13         14         15         16         17         18 
##   7.666667 -54.333333  42.666667 -34.333333  85.666667 -47.333333 
##         19         20         21         22         23         24 
## -22.666667 -85.666667  53.333333   8.333333  76.333333 -29.666667 
##         25         26         27         28         29         30 
##   2.166667 -49.833333  12.166667  12.166667  59.166667 -35.833333 
##         31         32         33         34         35         36 
## -52.000000 -52.000000   4.000000  47.000000  42.000000  11.000000 
##         37         38         39         40         41         42 
##   2.833333 -49.166667  29.833333 -34.166667  75.833333 -25.166667 
##         43         44         45         46         47         48 
## -59.833333 -18.833333   6.166667  33.166667  33.166667   6.166667 
##         49         50         51         52         53         54 
## -30.666667 -30.666667 -19.666667  37.333333  42.333333   1.333333
residuals(ex01.av)
##          1          2          3          4          5          6 
##   6.333333 -55.666667  38.333333  -8.666667  58.333333 -38.666667 
##          7          8          9         10         11         12 
## -25.500000 -46.500000  12.500000  11.500000  42.500000   5.500000 
##         13         14         15         16         17         18 
##   7.666667 -54.333333  42.666667 -34.333333  85.666667 -47.333333 
##         19         20         21         22         23         24 
## -22.666667 -85.666667  53.333333   8.333333  76.333333 -29.666667 
##         25         26         27         28         29         30 
##   2.166667 -49.833333  12.166667  12.166667  59.166667 -35.833333 
##         31         32         33         34         35         36 
## -52.000000 -52.000000   4.000000  47.000000  42.000000  11.000000 
##         37         38         39         40         41         42 
##   2.833333 -49.166667  29.833333 -34.166667  75.833333 -25.166667 
##         43         44         45         46         47         48 
## -59.833333 -18.833333   6.166667  33.166667  33.166667   6.166667 
##         49         50         51         52         53         54 
## -30.666667 -30.666667 -19.666667  37.333333  42.333333   1.333333

Apóss realizar a ANOVA, é necessário verificar os pressupostos do modelo:

Homocedasticidade:

Graficamente, podemos analisar a homocedasticidade através de um box-plot,

boxplot(ex01.av$res ~ trat)

Através de um gráfico dos resíduos vs tratamentos,

plot.default(ex01.av$res,trat)

Ainda, pode-se avaliar através de um teste, como por exemplo o teste de Bartlett:

bartlett.test(ex01.av$res, trat)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  ex01.av$res and trat
## Bartlett's K-squared = 3.6738, df = 8, p-value = 0.8853

Graficamente, pode-se avaliar a normalidade dos resíduos fazendo

hist(ex01.av$res, main=NULL)
title("Histograma dos resíduos padronizados")

stem(ex01.av$res)
## 
##   The decimal point is 1 digit(s) to the right of the |
## 
##    -8 | 6
##    -6 | 0
##    -4 | 64220977
##    -2 | 96441106530
##    -0 | 99
##     0 | 123466668812222
##     2 | 03378
##     4 | 22337389
##     6 | 66
##     8 | 6
qqnorm(ex01.av$res,ylab="Residuos", xlab="Resíduos", main=NULL)
qqline(ex01.av$res)
title("Grafico Normal de Probabilidade dos Resíduos")

Através do teste de Shapiro-Wilk

shapiro.test(ex01.av$res) #teste para normalidade
## 
##  Shapiro-Wilk normality test
## 
## data:  ex01.av$res
## W = 0.97159, p-value = 0.2263

A independência, com algumas restrições, pode ser analisada graficamente, através de

plot(ex01.av$fit, ex01.av$res, xlab="valores ajustados", ylab="resíduos")
title("resíduos vs Preditos")

Utilizando o critério de +3 ou -3 desvios padronizados, pode-se avaliar a existência de candidatos à outlier utilizando os seguintes comandos:

plot(ex01.av) # pressione a tecla enter para mudar o gráfico

par(mfrow=c(2,2))
plot(ex01.av)

par(mfrow=c(1,1))
names(anova(ex01.av))
## [1] "Df"      "Sum Sq"  "Mean Sq" "F value" "Pr(>F)"
s2 <- anova(ex01.av)$Mean[2] # estimativa da variância
res <- ex01$res # extraindo resíduos
respad <- (res/sqrt(s2)) # resíduos padronizados
boxplot(respad)
title("Resíduos Padronizados" )

plot.default(ex01$trat,respad)
title("Resíduos Padronizados" )

No R, o teste de Tukey é apresentado através de intervalos de confiança. A interpretação é: se o intervalo de confiança para a diferença entre duas médias não incluir o valor zero, siginifica que rejeita-se a hipótese nula, caso contrário, não rejeita-se.

O resultado pode ser visto através de uma tabela e/ou graficamente:

ex01.tu <- TukeyHSD(ex01.av)
ex01.tu
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = resp ~ trat, data = ex01)
## 
## $trat
##              diff         lwr         upr     p adj
## t2-t1   52.833333  -29.947633  135.614299 0.4998060
## t3-t1  -32.333333 -115.114299   50.447633 0.9342210
## t4-t1  -85.000000 -167.780966   -2.219034 0.0401018
## t5-t1  -36.833333 -119.614299   45.947633 0.8721075
## t6-t1   27.333333  -55.447633  110.114299 0.9749062
## t7-t1 -214.500000 -297.280966 -131.719034 0.0000000
## t8-t1   25.166667  -57.614299  107.947633 0.9849417
## t9-t1   37.000000  -45.780966  119.780966 0.8693183
## t3-t2  -85.166667 -167.947633   -2.385701 0.0394343
## t4-t2 -137.833333 -220.614299  -55.052367 0.0000730
## t5-t2  -89.666667 -172.447633   -6.885701 0.0247945
## t6-t2  -25.500000 -108.280966   57.280966 0.9836416
## t7-t2 -267.333333 -350.114299 -184.552367 0.0000000
## t8-t2  -27.666667 -110.447633   55.114299 0.9730043
## t9-t2  -15.833333  -98.614299   66.947633 0.9993743
## t4-t3  -52.666667 -135.447633   30.114299 0.5040619
## t5-t3   -4.500000  -87.280966   78.280966 1.0000000
## t6-t3   59.666667  -23.114299  142.447633 0.3369467
## t7-t3 -182.166667 -264.947633  -99.385701 0.0000002
## t8-t3   57.500000  -25.280966  140.280966 0.3855262
## t9-t3   69.333333  -13.447633  152.114299 0.1671352
## t5-t4   48.166667  -34.614299  130.947633 0.6203900
## t6-t4  112.333333   29.552367  195.114299 0.0018566
## t7-t4 -129.500000 -212.280966  -46.719034 0.0002153
## t8-t4  110.166667   27.385701  192.947633 0.0024139
## t9-t4  122.000000   39.219034  204.780966 0.0005599
## t6-t5   64.166667  -18.614299  146.947633 0.2479215
## t7-t5 -177.666667 -260.447633  -94.885701 0.0000004
## t8-t5   62.000000  -20.780966  144.780966 0.2886707
## t9-t5   73.833333   -8.947633  156.614299 0.1146645
## t7-t6 -241.833333 -324.614299 -159.052367 0.0000000
## t8-t6   -2.166667  -84.947633   80.614299 1.0000000
## t9-t6    9.666667  -73.114299   92.447633 0.9999849
## t8-t7  239.666667  156.885701  322.447633 0.0000000
## t9-t7  251.500000  168.719034  334.280966 0.0000000
## t9-t8   11.833333  -70.947633   94.614299 0.9999286
plot(ex01.tu)

Uma outra maneira é utilizar a biblioteca Agricolae ou laercio.

require(laercio)
## Loading required package: laercio
LTukey(ex01.av, which = "", conf.level = 0.95)
## 
##  TUKEY TEST TO COMPARE MEANS 
##  
##  Confidence level:  0.95 
##  Dependent variable:  resp
##  Variation Coefficient:  12.45212 % 
##  
## Independent variable:  trat 
##   Factors Means                
##   t2      431.5            a   
##   t9      415.666666666667 ab  
##   t6      406              ab  
##   t8      403.833333333333 ab  
##   t1      378.666666666667 ab  
##   t3      346.333333333333  bc 
##   t5      341.833333333333  bc 
##   t4      293.666666666667   c 
##   t7      164.166666666667    d
## 
##