O experimento fatorial descrito na apostila do curso de Planejamento de Experimentos II comparou o crescimento de mudas de eucalipto considerando diferentes recipientes e espécies.
Vamos considerar agora que os dados já esteajm digitados em um arquivo texto. Clique aqui para ver e copiar o arquivo com conjunto de dados para o seu diretório de trabalho.
A seguir vamos ler (importar) os dados para R com o comando read.table:
> ex04 <- read.table("exemplo04.txt", header=T)
> ex04
Antes de começar as análise vamos inspecionar o objeto que contém os dados para saber quantas observações e variáveis há no arquivo, bem como o nome das variáveis. Vamos tembém pedir o R que exiba um rápido resumo dos dados.
> dim(ex04) [1] 24 3 > names(ex04) [1] "rec" "esp" "resp" > attach(ex04) > is.factor(rec) [1] TRUE > is.factor(esp) [1] TRUE > is.factor(resp) [1] FALSE > is.numeric(resp) [1] TRUE
Nos resultados acima vemos que o objeto ex04 que contém os dados tem 24 linhas (observações) e 3 colunas (variáveis). As variáveis tem nomes rec, esp e resp, sendo que as duas primeiras são fatores enquanto resp é uma variável numérica, que neste caso é a variável resposta. O objeto ex04 foi incluído no caminho de procura usando o comando attach para facilitar a digitação.
Inicialmente vamos obter um resumo de nosso conjunto de dados usando a função summary.
> summary(ex04)
rec esp resp
r1:8 e1:12 Min. :18.60
r2:8 e2:12 1st Qu.:19.75
r3:8 Median :23.70
Mean :22.97
3rd Qu.:25.48
Max. :26.70
Note que para os fatores são exibidos o número de dados em cada nível do fator. Já para a variável numérica são mostrados algumas medidas estatísticas. Vamos explorar um pouco mais os dados
> ex04.m <- tapply(resp, list(rec,esp), mean)
> ex04.m
e1 e2
r1 25.650 25.325
r2 25.875 19.575
r3 20.050 21.325
> ex04.mr <- tapply(resp, rec, mean)
> ex04.mr
r1 r2 r3
25.4875 22.7250 20.6875
> ex04.me <- tapply(resp, esp, mean)
> ex04.me
e1 e2
23.85833 22.07500
Nos comandos acima calculamos as médias para cada fator, assim como para os cruzamentos entre os fatores. Note que podemos calcular outros resumos além da média. Experimente nos comandos acima substituir mean por var para calcular a variância de cada grupo, e por summary para obter um outro resumo dos dados.
Em experimentos fatoriais é importante verificar se existe interação entre os fatores. Inicialmente vamos fazer isto graficamente e mais a frente faremos um teste formal para presença de interação. Os comandos a seguir são usados para produzir os gráficos exibidos na Figura 1.
> par(mfrow=c(1,2)) > interaction.plot(rec, esp, resp) > interaction.plot(esp, rec, resp)
Seguindo o modelo adequado, o análise de variância para este experimento inteiramente casualizado em esquema fatorial pode ser obtida com o comando:
> ex04.av <- aov(resp ~ rec + esp + rec * esp)
Entretanto o comando acima pode ser simplificado produzindo os mesmos resultados com o comando
![\includegraphics[width=0.45\textwidth]{figuras/ex04inter1.ps}](node10%20Arquivos/img2.png)
![\includegraphics[width=0.45\textwidth]{figuras/ex04inter2.ps}](node10%20Arquivos/img3.png)
> ex04.av <- aov(resp ~ rec * esp)
> summary(ex04.av)
Df Sum Sq Mean Sq F value Pr(>F)
rec 2 92.861 46.430 36.195 4.924e-07 ***
esp 1 19.082 19.082 14.875 0.001155 **
rec:esp 2 63.761 31.880 24.853 6.635e-06 ***
Residuals 18 23.090 1.283
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Isto significa que no R, ao colocar uma interação no modelo, os efeitos principais são incluídos automaticamente. Note no quadro de análise de variância que a interação é denotada por rec:esp. A análise acima ostra que este efeito é significativo, confirmando o que verificamos nos gráficos de interação vistos anteriormente.
O objeto ex04.av guarda todos os resultados da análise e pode ser explorado por diversos comandos. Por exemplo a função model.tables aplicada a este objeto produz tabelas das médias definidas pelo modelo. O resultado mostra a média geral, médias de cada nível fatores e das combinações dos níveis dos fatores. Note que no resultado está incluído também o número de dados que gerou cada média.
> ex04.mt <- model.tables(ex04.av, ty="means")
> ex04.mt
Tables of means
Grand mean
22.96667
rec
r1 r2 r3
25.49 22.73 20.69
rep 8.00 8.00 8.00
esp
e1 e2
23.86 22.07
rep 12.00 12.00
rec:esp
esp
rec e1 e2
r1 25.650 25.325
rep 4.000 4.000
r2 25.875 19.575
rep 4.000 4.000
r3 20.050 21.325
rep 4.000 4.000
Mas isto não é tudo! O objeto ex04.av possui vários elementos que guardam informações sobre o ajuste.
> names(ex04.av) [1] "coefficients" "residuals" "effects" "rank" [5] "fitted.values" "assign" "qr" "df.residual" [9] "contrasts" "xlevels" "call" "terms" [13] "model" > class(ex04.av) [1] "aov" "lm"
O comando class mostra que o objeto ex04.av pertence às classes aov e lm. Isto significa que devem haver métodos associados a este objeto que tornam a exploração do resultado mais fácil. Na verdade já usamos este fato acima quando digitamos o comando summary(ex04.av). Existe uma função chamada summary.aov que foi utilizada já que o objeto é da classe aov. Iremos usar mais este mecanismo no próximo passo da análise.
Após ajustar o modelo devemos proceder a análise dos resíduos para verificar os pressupostos. O R produz automaticamente 4 gráficos básicos de resíduos conforme a Figura 2 com o comando plot.
> par(mfrow=c(2,2)) > plot(ex04.av)
Os gráficos permitem uma análise dos resíduos que auxiliam no julgamento da adequacidade do modelo. Evidentemente voce não precisa se limitar os gráficos produzidos automaticamente pelo R - voce pode criar os seus próprios gráficos muito facilmente. Neste gráficos voce pode usar outras variáveis, mudar texto de eixos e títulos, etc, etc, etc. Examine os comandos abaixo e os gráficos por eles produzidos.
![\includegraphics[width=\textwidth]{figuras/ex04res1.ps}](node10%20Arquivos/img4.png)
> par(mfrow=c(2,1))
> residuos <- resid(ex04.av)
> plot(ex04$rec, residuos)
> title("Resíduos vs Recipientes")
> plot(ex04$esp, residuos)
> title("Resíduos vs Espécies")
> par(mfrow=c(2,2))
> preditos <- (ex04.av$fitted.values)
> plot(residuos, preditos)
> title("Resíduos vs Preditos")
> s2 <- sum(resid(ex04.av)^2)/ex04.av$df.res
> respad <- residuos/sqrt(s2)
> boxplot(respad)
> title("Resíduos Padronizados")
> qqnorm(residuos,ylab="Residuos", main=NULL)
> qqline(residuos)
> title("Grafico Normal de \n Probabilidade dos Resíduos")
Além disto há alguns testes já programados. Como exemplo vejamos e teste de Shapiro-Wilk para testar a normalidade dos resíduos.
> shapiro.test(residuos)
Shapiro-Wilk normality test
data: residuos
W = 0.9293, p-value = 0.09402
Conforma visto na apostila do curso, quando a interação entre os
fatores é significativa podemos desdobrar os graus de liberdade de um
fator dentro de cada nível do outro.
A forma de fazer isto no R é reajustar o modelo utilizando
a notação / que indica efeitos aninhados.
Desta forma podemos desdobrar os efeitos de espécie dentro de cada
recipiente e vice versa conforme mostrado a seguir.
> ex04.avr <- aov(resp ~ rec/esp)
> summary(ex04.avr, split=list("rec:esp"=list(r1=1, r2=2, r3=3)))
Df Sum Sq Mean Sq F value Pr(>F)
rec 2 92.861 46.430 36.1952 4.924e-07 ***
rec:esp 3 82.842 27.614 21.5269 3.509e-06 ***
rec:esp: r1 1 0.211 0.211 0.1647 0.6897
rec:esp: r2 1 79.380 79.380 61.8813 3.112e-07 ***
rec:esp: r3 1 3.251 3.251 2.5345 0.1288
Residuals 18 23.090 1.283
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
> ex04.ave <- aov(resp ~ esp/rec)
> summary(ex04.ave, split=list("esp:rec"=list(e1=c(1,3), e2=c(2,4))))
Df Sum Sq Mean Sq F value Pr(>F)
esp 1 19.082 19.082 14.875 0.001155 **
esp:rec 4 156.622 39.155 30.524 8.438e-08 ***
esp:rec: e1 2 87.122 43.561 33.958 7.776e-07 ***
esp:rec: e2 2 69.500 34.750 27.090 3.730e-06 ***
Residuals 18 23.090 1.283
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Há vários testes de comparações múltiplas disponíveis na literatura, e muitos deles implementados no R. Os que não estão implementados podem ser facilmente calculados utilizando os recursos do R.
Vejamos por exemplo duas formas de usar o Teste de Tukey, a primeira usando uma implementação com a função TukeyHSD e uma segunda fazendo ops cálculos necessários com o R.
Poderíamos simplesmente digitar:
> ex04.tk <- TukeyHSD(ex04.av) > plot(ex04.tk) > ex04.tk
e obter diversos resultados. Entretanto nem todos nos interessam. Como a interação foi significativa na análise deste experimento a comparação dos níveis fatores principais não nos interessa.
Podemos então pedir a função que somente mostre a comparação de médias entre as combinações dos nívies dos fatores.
> ex04.tk <- TukeyHSD(ex04.ave, "esp:rec")
> plot(ex04.tk)
> ex04.tk
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = resp ~ esp/rec)
$"esp:rec"
diff lwr upr
[1,] -0.325 -2.8701851 2.220185
[2,] 0.225 -2.3201851 2.770185
[3,] -6.075 -8.6201851 -3.529815
[4,] -5.600 -8.1451851 -3.054815
[5,] -4.325 -6.8701851 -1.779815
[6,] 0.550 -1.9951851 3.095185
[7,] -5.750 -8.2951851 -3.204815
[8,] -5.275 -7.8201851 -2.729815
[9,] -4.000 -6.5451851 -1.454815
[10,] -6.300 -8.8451851 -3.754815
[11,] -5.825 -8.3701851 -3.279815
[12,] -4.550 -7.0951851 -2.004815
[13,] 0.475 -2.0701851 3.020185
[14,] 1.750 -0.7951851 4.295185
[15,] 1.275 -1.2701851 3.820185
Mas ainda assim temos resultados que não interessam. Mais especificamente estamos intessados nas comparações dos níveis de um fator dentro dos nívies de outro. Por exemplo, vamos fazer as comparações dos recipientes para cada uma das espécies.
Primeiro vamos obter
> s2 <- sum(resid(ex04.av)^2)/ex04.av$df.res
> dt <- qtukey(0.95, 3, 18) * sqrt(s2/4)
> dt
[1] 2.043945
>
> ex04.m
e1 e2
r1 25.650 25.325
r2 25.875 19.575
r3 20.050 21.325
>
> m1 <- ex04.m[,1]
> m1
r1 r2 r3
25.650 25.875 20.050
> m1d <- outer(m1,m1,"-")
> m1d
r1 r2 r3
r1 0.000 -0.225 5.600
r2 0.225 0.000 5.825
r3 -5.600 -5.825 0.000
> m1d <- m1d[lower.tri(m1d)]
> m1d
r2 r3 <NA>
0.225 -5.600 -5.825
>
> m1n <- outer(names(m1),names(m1),paste, sep="-")
> names(m1d) <- m1n[lower.tri(m1n)]
> m1d
r2-r1 r3-r1 r3-r2
0.225 -5.600 -5.825
>
> data.frame(dif = m1d, sig = ifelse(abs(m1d) > dt, "*", "ns"))
dif sig
r2-r1 0.225 ns
r3-r1 -5.600 *
r3-r2 -5.825 *
>
> m2 <- ex04.m[,2]
> m2d <- outer(m2,m2,"-")
> m2d <- m2d[lower.tri(m2d)]
> m2n <- outer(names(m2),names(m2),paste, sep="-")
> names(m2d) <- m2n[lower.tri(m2n)]
> data.frame(dif = m2d, sig = ifelse(abs(m2d) > dt, "*", "ns"))
dif sig
r2-r1 -5.75 *
r3-r1 -4.00 *
r3-r2 1.75 ns