Experimentos com fatores hierárquicos

Vamos considerar o exemplo da apostila retirado do livro de Montgomery. Clique aqui para ver e copiar o arquivo com conjunto de dados para sua área de trabalho.

O experimento estuda a variabilidade de lotes e fornecedores na puraza da matéria prima. A análise assume que os fornecedores são um efeito fixo enquanto que lotes são efeitos aleatórios.

Inicialmente vamos ler (importar) os dados para R com o comando read.table (certifique-se que o arquivo exemplo06.txt está na sua área de trabalho ou coloque o caminho do arquivo no comando abaixo). A seguir vamos examinar o objeto que contém os dados.

> ex06 <- read.table("exemplo06.txt", header=T)
> ex06

> dim(ex06)
[1] 36  3

> names(ex06)
[1] "forn" "lot"  "resp"

> is.factor(ex06$forn)
[1] FALSE
> is.factor(ex06$lot)
[1] FALSE

> ex06$forn <- as.factor(ex06$forn)
> ex06$lot <- as.factor(ex06$lot)

> is.factor(ex06$resp)
[1] FALSE
> is.numeric(ex06$resp)
[1] TRUE

> summary(ex06)
 forn   lot        resp        
 1:12   1:9   Min.   :-4.0000  
 2:12   2:9   1st Qu.:-1.0000  
 3:12   3:9   Median : 0.0000  
        4:9   Mean   : 0.3611  
              3rd Qu.: 2.0000  
              Max.   : 4.0000

Nos comandos acima verificamos que o objeto ex06 possui 36 linhas correspondentes às observações e 3 colunas que correspondem às variáveis forn (fornecedor), lot (lote) e resp (a variável resposta).

A seguir verificamos que forn e lot não foram lidas como fatores. NÃO podemos seguir as análise desta forma pois o R leria os valores destas variáveis como quantidades numéricas e não como indicadores dos nívies dos fatores. Para corrigir isto usamos o comando as.factor para indicar ao R que estas variáveis são fatores.

Finalmente verificamos que a variável resposta é numérica e produzimos um rápido resumo dos dados.

Na sequência deveríamos fazer uma análise exploratória, alguns gráficos descritivos etc, como na análise dos experimentos mostrados anteriormente. Vamos deixar isto por conta do leitor e passar direto para a análise de variância.

A notação para indicar efeitos aninhados no modelo é /. Desta forma poderíamos ajustar o modelo da seguinte forma:

> ex06.av <- aov(resp ~ forn/lot, data=ex06)
> summary(ex06.av)
            Df Sum Sq Mean Sq F value  Pr(>F)  
forn         2 15.056   7.528  2.8526 0.07736 .
forn:lot     9 69.917   7.769  2.9439 0.01667 *
Residuals   24 63.333   2.639                  
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

Embora os elementos do quadro de análise de variância estejam corretos o teste F para efeito dos fornecedores está ERRADO. A análise acima considerou todos os efeitos como fixos e portanto dividiu os quadrados médios dos efeitos pelo quadrado médio do resíduo. Como lotes é um efeito aleatório deveríamos dividir o quadrado médio de to termo lot pelo quadrado médio de forn:lot

Uma forma de indicar a estrutura hierárquica ao R é especificar o modelo de forma que o termo de resíduo seja dividido de maneira adequada. Veja o resultado abaixo.

> ex06.av1 <- aov(resp ~ forn/lot + Error(forn) , data=ex06)
> summary(ex06.av1)

Error: forn
     Df  Sum Sq Mean Sq
forn  2 15.0556  7.5278

Error: Within
          Df Sum Sq Mean Sq F value  Pr(>F)  
forn:lot   9 69.917   7.769  2.9439 0.01667 *
Residuals 24 63.333   2.639                  
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

Agora o teste F errado não é mais mostrado, mas o teste correto também não foi feito! Isto não é problema porque podemos extrair os elementos que nos interessam e fazer o teste desejado. Primeiro vamos guardar verificamos que o comando anova produz uma lista que tem entre seus elementos os graus de liberdade Df e os quadrados médios (Mean Sq. A partir destes elementos podemos obtemos o valor da estatística F e o valor P associado.

> ex06.anova <- anova(ex06.av)
> is.list(ex06.anova)
[1] TRUE

> names(ex06.anova)
[1] "Df"      "Sum Sq"  "Mean Sq" "F value" "Pr(>F)" 

> ex06.anova$Df
 1  2    
 2  9 24 
> ex06.anova$Mean
       1        2          
7.527778 7.768519 2.638889 

> Fcalc <- ex06.anova$Mean[1]/ex06.anova$Mean[2]
> Fcalc
        1 
0.9690107 
> pvalor <- 1 - pf(Fcalc, ex06.anova$Df[1], ex06.anova$Df[2])
> pvalor
        1 
0.4157831

USANDO O PACOTE NLME
Uma outra possível e elegante solução no R para este problema é utilizar a função lme do pacote nlme. Note que a abordagem do problema por este pacote é um pouco diferente da forma apresentada no curso por se tratar de uma ferramente geral para modelos com efeitos aleatórios. Entretanto os todos os elementos relevantes da análise estão incluídos nos resultados. Vamos a seguir ver os comandos necessários comentar os resultados.

Inicialmente temos que carregar o pacote nlme com o comando require.
A seguir criamos uma variável para indicar o efeito aleatório que neste exemplo chamamos de ex06$fa utilizando a função getGroups.
Feito isto podemos rodar a função lme que faz o ajuste do modelo.

> require(nlme)
[1] TRUE
> ex06$fa <- getGroups(ex06, ~ 1|forn/lot, level=2)
> ex06.av <- lme(resp ~ forn, random=~1|fa, data=ex06)
> ex06.av
Linear mixed-effects model fit by REML
  Data: ex06 
  Log-restricted-likelihood: -71.42198
  Fixed: resp ~ forn 
(Intercept)       forn2       forn3 
 -0.4166667   0.7500000   1.5833333 

Random effects:
 Formula: ~1 | fa
        (Intercept) Residual
StdDev:    1.307561 1.624483

Number of Observations: 36
Number of Groups: 12

Este modelo tem a variável forn como efeito fixo e a variável lot como efeito aleatório com o componente de variância $\sigma^2_{lote}$. Além disto temos a variância residual $\sigma^2$. A saída acima mostra as estimativas destes componentes da variância como sendo $\hat{\sigma}^2_{lote} = (1.307)^2=1.71$ e $\hat{\sigma}^2 = (1.624)^2=2.64$.

O comando anova vai mostrar a análise de variância com apenas os efeitos principais. O fato do programa não incluir o efeito aleatório de lotes na saída não causa problema algum. O comando intervals mostra os intervalos de confiança para os componentes de variância. Portanto para verificar a significância do efeito de lotes basta ver se o intervalo para este componente de variância exclui o valor $0$, o que é o caso neste exemplo conforme vamos abaixo.

> anova(ex06.av)
            numDF denDF   F-value p-value
(Intercept)     1    24 0.6043242  0.4445
forn            2     9 0.9690643  0.4158
> intervals(ex06.av)
Approximate 95% confidence intervals

 Fixed effects:
                 lower       est.    upper
(Intercept) -2.0772277 -0.4166667 1.243894
forn2       -1.8239746  0.7500000 3.323975
forn3       -0.9906412  1.5833333 4.157308

 Random Effects:
  Level: fa 
                    lower     est.    upper
sd((Intercept)) 0.6397003 1.307561 2.672682

 Within-group standard error:
   lower     est.    upper 
1.224202 1.624483 2.155644

Finalmente uma versão mais detalhada dos resultados pode ser obtida com o comando summary.

> summary(ex06.av)
Linear mixed-effects model fit by REML
 Data: ex06 
       AIC      BIC    logLik
  152.8440 160.3265 -71.42198

Random effects:
 Formula: ~1 | fa
        (Intercept) Residual
StdDev:    1.307561 1.624483

Fixed effects: resp ~ forn 
                 Value Std.Error DF    t-value p-value
(Intercept) -0.4166667 0.8045749 24 -0.5178718  0.6093
forn2        0.7500000 1.1378407  9  0.6591432  0.5263
forn3        1.5833333 1.1378407  9  1.3915246  0.1975
 Correlation: 
      (Intr) forn2 
forn2 -0.707       
forn3 -0.707  0.500

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.4751376 -0.7500844  0.0812409  0.7060895  1.8720268 

Number of Observations: 36
Number of Groups: 12

O próximo passo da seria fazer uma análise dos resíduos para verificar os pressupostos, semelhante ao que foi feito nos experimentos anteriormente analisados. Vamos deixar isto por conta do leitor.