O experimento fatorial descrito em Banzato & Kronka (1989) comparou o crescimento de mudas de eucalipto considerando como fatores diferentes tipos de recipientes e espécies.
Vamos considerar agora que os dados já estejam digitados em um arquivo texto. Clique aqui para ver e/ou copiar o arquivo com conjunto de dados para o seu diretório de trabalho. A seguir deve-se ler ("importar") os dados para R com o comando read.table(): Se você não tiver restrições de acesso (firewall, etc) pode importar o arquivo diretamente fornecendo a URL (endereço web) do arquivo.
> ex04 <- read.table("http://www.leg.ufpr.br/~paulojus/dados/exemplo04.txt", head=T, stringsAsFactors=TRUE) > # ex04 <- read.table("DEST/aulasR/dados/exemplo04.txt", head=T)
Antes de começar as análise vamos usar alguns comandos para 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 também pedir o R que exiba um rápido resumo dos dados e verificar se cada variável possui o "tipo"correto.
> head(ex04)
rec esp resp 1 r1 e1 26.2 2 r1 e1 26.0 3 r1 e1 25.0 4 r1 e1 25.4 5 r1 e2 24.8 6 r1 e2 24.6
> dim(ex04)
[1] 24 3
> names(ex04)
[1] "rec" "esp" "resp"
> is.factor(ex04$rec)
[1] TRUE
> is.factor(ex04$esp)
[1] TRUE
> is.factor(ex04$resp)
[1] FALSE
> is.numeric(ex04$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 no caso deste experimento é a variável resposta.
Inicialmente vamos obter um resumo de nosso conjunto de dados usando a função summary(). 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.
> 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
Vamos explorar um pouco mais os dados calculando as médias para cada nível de cada fator e também para as combinações dos níveis dos fatores.
> ex04.mr <- with(ex04,tapply(resp, rec, mean)) > ex04.mr
r1 r2 r3 25.4875 22.7250 20.6875
> ex04.me <- with(ex04,tapply(resp, esp, mean)) > ex04.me
e1 e2 23.85833 22.07500
> ex04.m <- with(ex04,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
As combinações dos níveis dos fatores podem ainda ser obtidas com interaction() que produz uma saída na forma de um vetor com nomes que combinam os níveis dos fatores envolvidos.
> with(ex04, tapply(resp, interaction(rec,esp), mean))
r1.e1 r2.e1 r3.e1 r1.e2 r2.e2 r3.e2 25.650 25.875 20.050 25.325 19.575 21.325
Nos comandos mostrados anteriormente a função mean() pode ser substituída por qualquer outra função de interesse seja pré definida ou definida pelo usuário. Nos exemplos a seguir ilustramos ambas situações onde são obtidas as medianas com a função pré-definida median e o número de observações acima de 22 para cada combinação dos fatores, com uma função definida por nós.
> with(ex04, tapply(resp, interaction(rec,esp), median))
r1.e1 r2.e1 r3.e1 r1.e2 r2.e2 r3.e2 25.70 26.00 19.30 25.00 19.30 21.35
> with(ex04, tapply(resp, interaction(rec,esp), function(x) sum(x > 22)))
r1.e1 r2.e1 r3.e1 r1.e2 r2.e2 r3.e2 4 4 1 4 0 1
As médias para os fatores e suas combinações também poderiam ser obtidas com o comando model.tables() o que será mostrado mais adiante. Entretanto neste estágio de análise descritiva, preferimos o mecanismo mais geral de tapply() que permite o cálculo de 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 62.
> with(ex04, interaction.plot(rec, esp, resp, ylab="médias", xlab="recipiente", xpd=F)) > with(ex04, interaction.plot(esp, rec, resp, ylab="médias", xlab="espécie", xpd=F))
Pode-se usar o R para obter outros tipos de gráficos de acordo com o interesse de quem está analisando os dados. Os comandos a seguir ilustram alguns outros tipos de gráficos que podemos produzir. Na figura 25.2 são mostrados gráficos semelhantes aos mostrados anteriormente, porém com pontos referentes às observações o que permite visualizar a variabilidade em cada grupo definido pelas combinações dos níveis dos fatores.
> with(ex04, plot.default(rec, resp, ty="n", ylab="Resposta", xlab="recipiente", )) > with(ex04, points(rec[esp=="e1"], resp[esp=="e1"], col=1)) > points(ex04.m[,1], pch="x", col=1, cex=1.5) > with(ex04, points(rec[esp=="e2"], resp[esp=="e2"], col=2)) > points(ex04.m[,2], pch="x", col=2, cex=1.5) > with(ex04, interaction.plot(rec, esp, resp, xpd=F, lty=1, col=1:2)) > > with(ex04, plot.default(esp, resp, ty="n", ylab="Resposta", xlab="espécie")) > with(ex04, points(esp[rec=="r1"], resp[rec=="r1"], col=1)) > points(ex04.m[1,], pch="x", col=1, cex=1.5) > with(ex04, points(esp[rec=="r2"], resp[rec=="r2"], col=2)) > points(ex04.m[2,], pch="x", col=2, cex=1.5) > with(ex04, points(esp[rec=="r3"], resp[rec=="r3"], col=3)) > points(ex04.m[3,], pch="x", col=3, cex=1.5) > with(ex04, interaction.plot(esp, rec, resp, xpd=F, lty=1, col=1:3))
Além destes gráficos produzidos pelo sitema básico de gráficos do R pode-se usar comandos fornecidos pelo pacote lattice que implementam um poderoso conjunto alternativo de gráficos mas não serão abordados aqui.
Seguindo o modelo adequado, a análise de variância para este experimento inteiramente casualizado em esquema fatorial pode ser obtida com as funções aov() ("analysis of variance") ou lm() ("linear model"). A primeira usa a segunda internamente visto que o modelo é linear, porém ajusta os resultados em um formato em geral mais adequado para análise de experimentos. Nestas funções os modelos são declarados por "fórmulas". A seguir vemos duas fórmulas que especificam o mesmo modelo.
> ex04.av <- aov(resp ~ rec + esp + rec:esp, data=ex04) > ex04.av <- aov(resp ~ rec * esp, data=ex04) > summary(ex04.av)
Df Sum Sq Mean Sq F value Pr(>F) rec 2 92.86 46.43 36.20 4.92e-07 *** esp 1 19.08 19.08 14.88 0.00116 ** rec:esp 2 63.76 31.88 24.85 6.64e-06 *** Residuals 18 23.09 1.28 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Isto significa que dentro de uma fórmula no R, o símbolo ":" define o termo de interação e "*" indica da inclusão dos efeitos principais e interações. A análise acima mostra que neste caso o efeito de interação é significativo, confirmando o que foi indicado nos gráficos exploratórios do efeito 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 da classe aov produz tabelas dos efeitos (se type="effects") ou das médias (se type="means") definidas pelo modelo. O resultado mostra a média geral, médias de cada nível dos fatores e das combinações dos níveis dos fatores. No resultado está incluído também o número de dados que gerou cada média.
> model.tables(ex04.av, type="means")
Tables of means Grand mean 22.96667 rec rec r1 r2 r3 25.487 22.725 20.687 esp esp e1 e2 23.858 22.075 rec:esp esp rec e1 e2 r1 25.650 25.325 r2 25.875 19.575 r3 20.050 21.325
Mas isto ainda não é tudo que se pode extrair da análise! O objeto ex04.av possui vários elementos que guardam diversas outras informações sobre o ajuste do modelo e que podem ser exploradas subsequentemente por métodos de funções para as classes aov e lm ou por requisições definidas pelo usuário. A seguir veremos alguns exemplos.
> 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"
A chamada 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, a análise de resíduos.
A análise de resíduos é útil para verificar os pressupostos do modelo. Usando o mecanismos de classes, o comando plot(ex04.av) aplicado sobre o objeto que contém o ajuste do modelo produz uma figura com quatro gráficos básicos para análise dos resíduos conforme mostrado na Figura 64.
Os gráficos permitem uma análise dos resíduos que auxilia no julgamento da adequação do modelo. Evidentemente não é necessário limitar-se aos gráficos produzidos automaticamente pelo R – você pode criar os seus próprios gráficos. Nestes gráficos pode-se usar outras variáveis, tipos de gráficos, mudar texto de eixos e títulos, etc, etc, etc. Os comandos a seguir mostram como obter os gráficos boxplot dos resíduos para os níveis de cada um dos fatores como mostrado na Figura 65.
> residuos <- resid(ex04.av) > plot(ex04$rec, residuos) > title("Resíduos vs Recipientes") > plot(ex04$esp, residuos) > title("Resíduos vs Espécies")
A Figura 66 mostra outros gráficos definidos pelo usuário: resíduos versus valores preditos, um boxplot dos resíduos padronizados, e um qqplot dos resíduos do modelo. Note que o objeto que contém o ajuste foi utilizado para extrair resíduos, valores preditos e a estimativa s2 da variância dos resíduos.
> preditos <- fitted(ex04.av) > plot(residuos, preditos) > title("Resíduos vs Preditos") > s2 <- sum(residuos^2)/ex04.av$df.res > respad <- residuos/sqrt(s2) > boxplot(respad) > title("Resíduos Padronizados") > qqnorm(residuos, ylab="Resíduos", main=NULL) > qqline(residuos) > title("Gráfico Normal de \n Probabilidade dos Resíduos")
Além da análise gráfica de resíduos há alguns testes já programados em funções. Como exemplo vejamos o teste de Shapiro-Wilks para testar a normalidade dos resíduos.
> shapiro.test(residuos)
Shapiro-Wilk normality test
data: residuos
W = 0.9293, p-value = 0.09402
Quando a interação entre os fatores é significativa pode-se adotar como estratégia de análise o
desdobramento dos graus de liberdade de um fator dentro de cada nível do outro fator. Uma forma de
obter tal desdobramento 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, data=ex04) > 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.86 46.43 36.195 4.92e-07 *** rec:esp 3 82.84 27.61 21.527 3.51e-06 *** rec:esp: r1 1 0.21 0.21 0.165 0.690 rec:esp: r2 1 79.38 79.38 61.881 3.11e-07 *** rec:esp: r3 1 3.25 3.25 2.535 0.129 Residuals 18 23.09 1.28 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> ex04.ave <- aov(resp ~ esp/rec, data=ex04) > 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.08 19.08 14.88 0.00116 ** esp:rec 4 156.62 39.16 30.52 8.44e-08 *** esp:rec: e1 2 87.12 43.56 33.96 7.78e-07 *** esp:rec: e2 2 69.50 34.75 27.09 3.73e-06 *** Residuals 18 23.09 1.28 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Neste ponto vale uma explicação mais detalhada sobre como obter os desdobramentos da interação usando o argumento split, em particular como definir os elementos da lista que, no caso de ex04.ave foram e1=c(1,3) e e2=c(2,4). Iniciamente vamos extrair usando effects() os efeitos ajustados pelo modelo.
> effects(ex04.ave)[1:6]
(Intercept) espe2 espe1:recr2 espe2:recr2 espe1:recr3 espe2:recr3 -112.513229 -4.368257 4.939804 6.123724 -7.919596 -5.656854
Os efeitos que temos interesse no desdobramento são os da interação, que são: espe1:recr2, espe2:recr2, espe1:recr3 e espe2:recr3. Portanto temos que localizar no vetor de efeitos as posições desses efeitos de interação que são: \(1^o:\) espe1:recr2, \(2^o:\) espe2:recr2, \(3^o:\) espe1:recr3 e \(4^o:\) espe2:recr3. Isto mostra que a posição dos efeitos que contém a espécie1 (e1) são 1 e 3, e especie2 (e2) são 2 e 4 o que define os valores nos vetores indicados no argumento split.
Há vários testes de comparações múltiplas disponíveis na literatura, e muitos deles são implementados nos pacotes básicos do R e/ou em pacotes contribuídos. Por exemplo, o pacote multcomp é inteiramente dedicado à implementação de diversos procedimentos de comparações múltiplas no R. Além disto, procedimentos que não estejam implementados podem ser calculados utilizando os recursos usuais do R utilizando os objetos com o ajuste dos modelos. Como ilustração mostramos a seguir duas formas de obter resultados para o Teste de Tukey, a primeira usando uma implementação já disponível com a função TukeyHSD() e uma segunda fazendo os cálculos necessários passo a passo com operações básicas do R. Para a função já disponível simplesmente digitamos os comandos a seguir e os resultados podem ser mostrados na forma texto ou gráfica como na Figura 67 que é produzida com o comando plot(ex04.tk1).
> ex04.tk1 <- TukeyHSD(ex04.av) > ex04.tk1
Tukey multiple comparisons of means 95% family-wise confidence level Fit: aov(formula = resp ~ rec * esp, data = ex04) $rec diff lwr upr p adj r2-r1 -2.7625 -4.207787 -1.3172128 0.0003395 r3-r1 -4.8000 -6.245287 -3.3547128 0.0000003 r3-r2 -2.0375 -3.482787 -0.5922128 0.0055472 $esp diff lwr upr p adj e2-e1 -1.783333 -2.75476 -0.8119067 0.0011553 $`rec:esp` diff lwr upr p adj r2:e1-r1:e1 0.225 -2.3201851 2.770185 0.9997185 r3:e1-r1:e1 -5.600 -8.1451851 -3.054815 0.0000204 r1:e2-r1:e1 -0.325 -2.8701851 2.220185 0.9983324 r2:e2-r1:e1 -6.075 -8.6201851 -3.529815 0.0000068 r3:e2-r1:e1 -4.325 -6.8701851 -1.779815 0.0004825 r3:e1-r2:e1 -5.825 -8.3701851 -3.279815 0.0000120 r1:e2-r2:e1 -0.550 -3.0951851 1.995185 0.9811892 r2:e2-r2:e1 -6.300 -8.8451851 -3.754815 0.0000041 r3:e2-r2:e1 -4.550 -7.0951851 -2.004815 0.0002705 r1:e2-r3:e1 5.275 2.7298149 7.820185 0.0000444 r2:e2-r3:e1 -0.475 -3.0201851 2.070185 0.9902110 r3:e2-r3:e1 1.275 -1.2701851 3.820185 0.6135909 r2:e2-r1:e2 -5.750 -8.2951851 -3.204815 0.0000143 r3:e2-r1:e2 -4.000 -6.5451851 -1.454815 0.0011258 r3:e2-r2:e2 1.750 -0.7951851 4.295185 0.2914242
Esta saída fornece resultados detalhados de várias comparações possíveis entre os níveis dos fatores e suas combinações. Entretanto, neste caso, nem todos os resultados mostrados 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íveis dos fatores e o gráfico com tais resultados pode ser obtido com plot(ex04.tk2).
> ex04.tk2 <- TukeyHSD(ex04.ave, "esp:rec") > ex04.tk2
Tukey multiple comparisons of means 95% family-wise confidence level Fit: aov(formula = resp ~ esp/rec, data = ex04) $`esp:rec` diff lwr upr p adj e2:r1-e1:r1 -0.325 -2.8701851 2.220185 0.9983324 e1:r2-e1:r1 0.225 -2.3201851 2.770185 0.9997185 e2:r2-e1:r1 -6.075 -8.6201851 -3.529815 0.0000068 e1:r3-e1:r1 -5.600 -8.1451851 -3.054815 0.0000204 e2:r3-e1:r1 -4.325 -6.8701851 -1.779815 0.0004825 e1:r2-e2:r1 0.550 -1.9951851 3.095185 0.9811892 e2:r2-e2:r1 -5.750 -8.2951851 -3.204815 0.0000143 e1:r3-e2:r1 -5.275 -7.8201851 -2.729815 0.0000444 e2:r3-e2:r1 -4.000 -6.5451851 -1.454815 0.0011258 e2:r2-e1:r2 -6.300 -8.8451851 -3.754815 0.0000041 e1:r3-e1:r2 -5.825 -8.3701851 -3.279815 0.0000120 e2:r3-e1:r2 -4.550 -7.0951851 -2.004815 0.0002705 e1:r3-e2:r2 0.475 -2.0701851 3.020185 0.9902110 e2:r3-e2:r2 1.750 -0.7951851 4.295185 0.2914242 e2:r3-e1:r3 1.275 -1.2701851 3.820185 0.6135909
Mas ainda assim temos resultados que podem não interessar. Mais especificamente, considere que estamos interessados nas comparações dos níveis de um fator dentro de cada um dos níveis do outro fator. Neste ponto, vamos fazer as comparações dos recipientes para cada uma das espécies, fazendo os cálculos passo a passo. Primeiro vamos obter a estimativa da variância dos resíduos, que é usada junto com o valor da amplitude estudantizada fornecida por qtukey() para obter o valor da diferença mínima significativa que no código a seguir armazenamos no objeto dt.
> s2 <- sum(resid(ex04.av)^2)/ex04.av$df.res > dt <- qtukey(0.95, 3, 18) * sqrt(s2/4) > dt
[1] 2.043945
Este valor é então usado para comparar as médias de interesse. Anteriormente armazenamos as médias para as combinações de todos os níveis dos fatores no objeto ex04.m onde as linhas se referem aos recipientes e colunas às espécies. No objeto m1 armazenamos as médias para espécie1 e na sequência são feitos cálculos para verificar a significância da diferença entre as médias dos recipientes para esta espécie.
> # comparação de médias de recipientes para espécie 1 : > 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
[1] 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 *
> # comparação de médias de recipientes para espécie 2 : > 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
No código mostrado anteriormente fazemos alguma manipulação dos objetos para formatar a saída. Esta sequência pode ser usada para definir uma função o que evitaria a digitação de todos estes comandos a cada comparação de médias desejada. Procedimento análogo pode ser adotado para fazer outras comparações de interesse.