25 Análise de experimentos em esquema fatorial

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.

25.1 Lendo os dados

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 voce 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)

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.

25.2 Análise exploratória

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ívies 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 mediam 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 so 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 25.2.

  > 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))


PIC

Figura 62: Gráficos de interação entre os fatores.


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))


PIC

Figura 63: Gráficos de pontos examinando a interação entre os fatores.


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.

25.3 Análise de variância

Seguindo o modelo adequado, o 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.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 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 for 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="effects") 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.688
  
   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 residuos.

25.4 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 25.4.


PIC

Figura 64: Gráficos de resíduos produzidos para objetos da classe lm.


Os gráficos permitem uma análise dos resíduos que auxilia no julgamento da adequacidade do modelo. Evidentemente não é necessario limitar-se aos gráficos produzidos automaticamente pelo R – voce pode criar os seus próprios gráficos. Neste 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 25.4.

  > residuos <- resid(ex04.av)
  > plot(ex04$rec, residuos)
  > title("Resíduos vs Recipientes")
  > plot(ex04$esp, residuos)
  > title("Resíduos vs Espécies")


PIC

Figura 65: Gráficos de resíduos para cada um dos fatores.


A Figura 25.4 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")


PIC

Figura 66: Alguns gráficos de resíduos definidos pelo usuário.


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

25.5 Desdobrando interações

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.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, 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.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

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: 1o : espe1:recr2, 2o : espe2:recr2, 3o : espe1:recr3 e 4o : 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.

25.6 Teste de Tukey para comparações múltiplas

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 sem fazendo os cálculos necessários passo a passo com operações básicas do R. Para 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 25.6 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


PIC

Figura 67: Visualização dos resultados do teste de Tukey de comparações múltiplas.


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 intessados 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.