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
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:
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.
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 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:
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}\)).
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.
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 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:
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.
Alguns testes para a verificação da homocedasticidade:
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.
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.
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.
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
##
##