-
Faça uma análise de variância e avalie os pressupostos.
Resposta
Fazendo a análise de variância tem-se:
mod1 <- aov(GPD ~ Tratamentos,dad1)
Antes de vermos os resultados vamos avaliar os pressupostos.
erros1 <- residuals(mod1)
par(mfrow=c(1,2))
qqnorm(erros1)
qqline(erros1)
plot(erros1 ~ fitted(mod1))
abline(h=0,col='red')
Os pressupostos foram atendidos. Vamos as resultados.
summary(mod1)
## Df Sum Sq Mean Sq F value Pr(>F)
## Tratamentos 3 0.02369 0.007897 4.048 0.0335 *
## Residuals 12 0.02341 0.001951
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Considerando um \(\alpha=5\%\) rejeita-se \(H_0\), ou seja, pelo menos um dos tratamentos difere dos demais.
-
Utilize um teste de comparação múltipla e justifique o uso de tal teste.
Resposta
Vamos calcular o coeficiente de variação e ver como está a dispersão dos dados em relação a média geral.
library(TukeyC)
## Loading required package: doBy
cv1 <- cv(mod1)
cv1
## cv
## 2.19
Como o experimento tem um baixo coeficiente de variação, mostrando deste modo que as médias são estimadas com uma boa precisão, vamos utilizar o teste de Tukey, que é um teste rigoroso.
tk1 <- TukeyC(mod1, which = 'Tratamentos')
summary(tk1)
## Goups of means at sig.level = 0.05
## Means G1
## 1 2.07 a
## 3 2.05 a
## 2 1.99 a
## 4 1.98 a
##
## Matrix of the difference of means above diagonal and
## respective p-values of the Tukey test below diagonal values
## 1 3 2 4
## 1 0.000 0.017 0.080 0.089
## 3 0.946 0.000 0.063 0.072
## 2 0.101 0.238 0.000 0.009
## 4 0.062 0.153 0.991 0.000
Uma coisa interessante aconteceu neste exemplo. O teste de Tukey não concordou com o teste de F, ou seja, o teste de F atesta que existe pelo menos um tratamento que diferen dos demais e o teste de Tukey não aponta nenhuma diferença entre os tratamentos, isso ao um nível de 5% de significância. Não há um concenso sobre tal discordância, e de fato algumas vezes isto pode acontecer. No nível que esta disciplina está sendo apresentada, não iremos adentrar nas vastas discussões teóricas que há na literatura e vamos adotar outro teste de comparações múltiplas. Neste caso, vamos adotar o ScottKnott.
library(ScottKnott)
sk1 <- SK(mod1,which='Tratamentos')
summary(sk1)
## Levels Means SK(5%)
## 1 2.06725 a
## 3 2.05025 a
## 2 1.98750 b
## 4 1.97850 b
em termos gráficos:
plot(sk1)
Portanto, podemos concluir que os tratamentos 1 e 3 são semelhantes entre si e superiores aos tratamentos 2 e 4. Os tratamentos 2 e 4 são semelhantes entre si.
As hipóteses são:
\[
H_0: \mu_{lab1} = \mu_{lab2} = \mu_{lab3} = \mu_{lab4}\\
H_a: \mu_{lab_i} \neq \mu_{lab_j}
\]
Vamos fazer uma anova e em seguida uma análise de resíduos.
dados <- data.frame(labs = factor(rep(1:4,rep(3,4))),
almet = c(l1,l2,l3,l4))
mod <- aov(almet ~ labs,dados)
erros <- residuals(mod)
par(mfrow=c(1,2))
qqnorm(erros)
qqline(erros)
plot(erros ~ fitted(mod))
abline(h=0,col='red')
Quanto a normalidade podemos afirmar que os erros são normais. Quanto a homocedasticidade, vamos aplicar o teste F-máximo de Hartley pois pela análise gráfica houve dúvidas.
variancias <- with(dados, tapply(almet,labs,var))
Fmax <- max(variancias)/min(variancias)
pf(Fmax,2,2,lower.tail=F)
## [1] 0.1982064
Logo podemos afirmar que as variâncias são homocedásticas. Vamos ao resultado da anova.
summary(mod)
## Df Sum Sq Mean Sq F value Pr(>F)
## labs 3 1.0559 0.3520 3.958 0.0531 .
## Residuals 8 0.7114 0.0889
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Portanto, considerando um nível de significância de 10% podemos afirmar que há diferenças significativas entre os latoratórios quanto a porcentagem de álcool metílico. Vamos aplicar o teste de Tukey.
Vamos avaliar o CV.
cv2 <- cv(mod)
cv2
## cv
## 0.35
Um baixo CV. Vamos aplicar o teste de Tukey.
tk2 <- TukeyC(mod, which='labs',sig.level=0.1)
summary(tk2)
## Goups of means at sig.level = 0.1
## Means G1 G2
## 1 85.06 a
## 3 84.77 a b
## 2 84.72 a b
## 4 84.23 b
##
## Matrix of the difference of means above diagonal and
## respective p-values of the Tukey test below diagonal values
## 1 3 2 4
## 1 0.000 0.293 0.343 0.827
## 3 0.641 0.000 0.050 0.533
## 2 0.527 0.997 0.000 0.483
## 4 0.038 0.205 0.269 0.000
Graficamente tem-se:
plot(tk2,dispersion='cip',ylim=c(83,87))
Então podemos concluir o seguinte: O laboratório 1 apresentou em média uma porcentagem de álcool metílico superior ao laboratório 4 e semelhantes aos laboratórios 3 e 2. Os laboratórios 3, 2 e 4 foram semelhantes entre si.
Famos fazer primeiramente uma ANOVA e apresentar primeiramente uma análise de resíduos.
mod3 <- aov(resis ~ Trat, dad4)
erros <- residuals(mod3)
par(mfrow=c(1,2))
qqnorm(erros)
qqline(erros)
plot(erros ~ fitted(mod3))
abline(h=0,col='red')
Quanto a normalidade tudo certo. No entanto, quanto a homocedasticidade, parece que está heterogênea. Vamos ao teste F-máximo de Hartley.
variancias <- with(dad4, tapply(resis,Trat,var))
Fmax <- max(variancias)/min(variancias)
pf(Fmax,2,2,lower.tail=FALSE)
## [1] 0.03425614
Considerando um \(\alpha=0,05\) podemos rejeitar a hipótese de homocedasticidade. Neste caso precisamos fazer uma transformação nos dados. Vamos utilizar a transformação logarítmica.
dad4$logres <- log(dad4$resis)
mod31 <- aov(logres ~ Trat, dad4)
erros1 <- residuals(mod31)
par(mfrow=c(1,2))
qqnorm(erros1)
qqline(erros1)
plot(erros1 ~ fitted(mod31))
abline(h=0,col='red')
Parece que não houve melhora quanto a homocedasticidade. Vamos verificar por meio do teste F-máximo de Hartley.
varianciaslog <- with(dad4, tapply(logres,Trat,var))
Fmaxlog <- max(varianciaslog)/min(varianciaslog)
pf(Fmaxlog,2,2,lower.tail=FALSE)
## [1] 0.08679278
Como p-valor é maior do que \(\alpha\), vamos admitir que a transformação resolveu a o problema da homocedasticidade. Vamos ao teste de comparações múltiplas. Avaliando a precisão do experimento.
cv1 <- cv(mod31)
cv1
## cv
## 0.86
Vamos então ao teste de Tukey.
tk4 <- TukeyC(mod31, which='Trat')
summary(tk4)
## Goups of means at sig.level = 0.05
## Means G1 G2 G3 G4
## AIP 4.86 a
## AR 4.75 b
## ACC 3.52 c
## ACP 3.40 d
##
## Matrix of the difference of means above diagonal and
## respective p-values of the Tukey test below diagonal values
## AIP AR ACC ACP
## AIP 0.000 0.11 1.337 1.463
## AR 0.022 0.00 1.228 1.354
## ACC 0.000 0.00 0.000 0.126
## ACP 0.000 0.00 0.010 0.000
Conforme vimos em aula, não é adequado a apresentação dos dados na escala transformada. Neste caso, iremos fazer a transformação de volta das médias apresentadas para a escala original.
medlog <- tk4$info$Means$means
mediorig <- exp(medlog)
names(mediorig) <- tk4$info$Means$Trat
mediorig
## AIP AR ACC ACP
## 129.24577 115.82730 33.93132 29.91109
Vamos apresentar uma tabela sugestiva de como se deve colocar nos artigos, do mesmo modo que os exercícios anteriores. Neste caso, precisamos calcular o erro padrão da média. Não tem como fazer a transformação de volta com o erro padrão da anova. Neste caso, vamos utilizar uma aproximação razoável como apresentado nos slides de aula.
EPM <- round(sqrt((sum(variancias)/4)/3),2)
EPM
## [1] 1.57
Logo, a tabela juntamente com o p-valor da anova é:
|
|
AIP
|
AR
|
ACC
|
ACP
|
EPM
|
P-valor
|
|
Resistência
|
129.25 a
|
115.83 b
|
33.93 c
|
29.91 d
|
1.57
|
2.73130910803365e-11
|
Fazendo uma anova e testando os pressupostos graficamente.
mod = aov(peso_final ~ bloco + tratamentos,daddfm)
par(mfrow=c(2,2))
plot(mod)
Inferencialmente.
library(nortest)
resi = residuals(mod)
resi
## 1 2 3 4 5 6 7 8
## -0.21275 0.17725 -0.28275 -0.05275 0.23725 -0.00925 0.01075 -0.01925
## 9 10 11 12 13 14 15 16
## 0.04075 0.11075 0.21725 -0.13275 -0.18275 -0.02275 -0.18275 0.14075
## 17 18 19 20 21 22 23 24
## 0.06075 -0.21925 0.20075 0.12075 -0.01975 -0.01975 0.07025 -0.05975
## 25 26 27 28 29 30 31 32
## -0.03975 -0.07625 0.22375 0.23375 -0.13625 -0.17625 -0.06675 0.07325
## 33 34 35 36 37 38 39 40
## 0.26325 0.09325 0.14325 0.19675 -0.13325 -0.31325 -0.27325 0.01675
lillie.test(resi)# testando normalidade
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: resi
## D = 0.071728, p-value = 0.8704
varis = with(daddfm, tapply(peso_final, tratamentos,var))
varis
## sorgo sorgoMaisCE milheto milhetoMaisCE
## 0.02756000 0.03247111 0.02031222 0.03291111
Fmax = max(varis)/min(varis)
Fmax
## [1] 1.620261
pf(Fmax,9,9,lower.tail=FALSE)#2blocos * 5 repetições = 10 repetições no geral
## [1] 0.2417089
Os pressupostos estão ok! Vamos aos resultados da anova.
summary(mod)
## Df Sum Sq Mean Sq F value Pr(>F)
## bloco 1 0.0216 0.02162 0.759 0.38972
## tratamentos 3 0.5652 0.18842 6.610 0.00118 **
## Residuals 35 0.9977 0.02850
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Houve diferenças significativas entre os tratamentos (P=0,00118). Vamos então elaborar os contrastes de interesse do pesquisador.
cmat <- rbind('sorgo vs milho' = c(+1, 0, -1, 0),
'sorgo vs sorgoCE' = c(+1,-1, 0, 0),
'milho vs milhoCE' = c( 0, 0, +1,-1))
cmat
## [,1] [,2] [,3] [,4]
## sorgo vs milho 1 0 -1 0
## sorgo vs sorgoCE 1 -1 0 0
## milho vs milhoCE 0 0 1 -1
Percebam que os contrastes não são mutuamente ortogonais, pois o primeiro contraste não é ortogonal ao segundo e nem ao terceiro. Apenas os contrastes dois e três são ortogonais. Isso não impede o andamento das análises, pois cada contraste é testado individualemente. No entanto, não podemos desdobrar os contrastes no quadro da anova, pois, como não há ortogonalidade entre os contrastes, a soma de quadrados dos contrastes não irão ser iguais a soma de quadrados dos tratamentos.
Fazendo os cálculos manualmente, tem-se:
# Estimando as médias
medias = with(daddfm, tapply(peso_final,tratamentos,mean))
medias
## sorgo sorgoMaisCE milheto milhetoMaisCE
## 2.346 2.306 2.573 2.550
# Calculando os contrastes
Y1 = +1*2.346 +0*2.306 -1*2.573 +0*2.550
Y1
## [1] -0.227
Y2 = +1*2.346 -1*2.306
Y2
## [1] 0.04
Y3= -1*2.573 +1*2.550
Y3
## [1] -0.023
# Calculando a soma de quadrados de cada contraste
SQ1 = (Y1^2*10)/2# o valor 10 é porque as médias são estimadas de 10 repetições (2 blocos * 5 repetições)
SQ1
## [1] 0.257645
SQ2 = (Y2^2*10)/2
SQ2
## [1] 0.008
SQ3 = (Y3^2*10)/2
SQ3
## [1] 0.002645
# Calculando a estatística F para cada constraste
QME = summary(mod)[[1]][3,3]#quadrado médio do erro da anova
QME
## [1] 0.02850479
Fy1 = SQ1/QME
Fy2 = SQ2/QME
Fy3 = SQ3/QME
# Calculando o p-valor para cada contraste
pfy1 = pf(Fy1, 1, 35, lower.tail=FALSE)
pfy1
## [1] 0.004866588
pfy2 = pf(Fy2, 1, 35, lower.tail=FALSE)
pfy2
## [1] 0.5996172
pfy3 = pf(Fy3, 1, 35, lower.tail=FALSE)
pfy3
## [1] 0.762462
Fazendo os contrastes de forma automatizada, tem-se:
# Fazendo de modo automático
library(gmodels)
# Como os contrastes não são ortogonais, iremos utilizar a função fit.contrast do pacote gmodels
tratamentos=daddfm$tratamentos #este objeto deve ter o mesmo nome do objeto do modelo da anova.
fit.contrast(mod,
tratamentos,
cmat) # fit.contrast (gmodels): testa constraste(s)
## Estimate Std. Error t value Pr(>|t|)
## tratamentossorgo vs milho -0.227 0.07550468 -3.0064360 0.004866588
## tratamentossorgo vs sorgoCE 0.040 0.07550468 0.5297685 0.599617235
## tratamentosmilho vs milhoCE 0.023 0.07550468 0.3046169 0.762461958
## attr(,"class")
## [1] "fit_contrast"
Neste caso a estatística \(t^2 = F\).