\section*{Aula 02: Experimentos hierárquicos}
\indent


\subsection*{Objetivo dessa aula}
\indent

Nesta aula veremos como analisar experimentos hierárquicos no \R.

Procure entender os comandos utilizados e também interpretar os resultados.

\subsection*{Trabalhando com o arquivo de dados}
\indent

Vamos considerar primeiro o exemplo da apostila retirado do livro de
Montgomery. 
\htmladdnormallink{Clique aqui}{../dados/exemplo06.txt} para ver e
copiar o arquivo com o conjunto de dados. 

Este experimento estuda a variabilidade de lotes e a diferença entre fornecedores com relação a
pureza 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
{\tt read.table}. 
\begin{verbatim}
> ex06 <- read.table("exemplo06.txt", header=T)
> ex06
\end{verbatim}

A seguir vamos examinar o objeto que contém os dados.

\begin{verbatim}
> 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  
\end{verbatim}

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

A seguir verificamos que {\tt forn} e {\tt lot} não foram lidas como
fatores. {\bf 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 {\tt 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êcia deveríamos fazer uma análise exploratória, alguns gráficos
descritivos etc, como na análise do experimento da aula anterior.
Vamos deixar isto por conta do leitor e passar
direto para a análise de variância.

\subsection*{Análise de variância}
\indent


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

\begin{verbatim}
> 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 
\end{verbatim}

Confira os graus de liberdade das causas de variação. Eles estão corretos?

Embora os elementos do quadro de análise de variância estejam corretos
o teste F para efeito dos fornecedores está {\bf 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 {\tt lotes} é um efeito aleatório deveríamos dividir 
o quadrado médio de to termo {\tt lot} pelo quadrado médio de
{\tt forn:lot}. Veja nas notas de aula o quadro com as esperanças dos quadrados médios dos efeitos para este modelo.

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. 

\begin{verbatim}
> 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 
\end{verbatim}

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 verificamos que o comando {\tt anova} produz
uma lista que tem entre seus elementos os graus de liberdade {\tt Df}
e os quadrados médios ({\tt Mean Sq}. A partir destes elementos podemos obter o valor da estatística F e o valor~P associado.

\begin{verbatim}
> 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 
\end{verbatim}


\subsubsection*{Usando o pacote NLME}

Uma outra possível e elegante solução no \R{} para este problema é utilizar a
função {\tt 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 todos os elementos relevantes da análise estão incluídos nos
resultados.

Vamos a seguir ver os comandos necessários e comentar os resultados.

Inicialmente temos que carregar o pacote \nlme{} com o comando {\tt  require}.

A seguir criamos uma variável para indicar o efeito aleatório que neste exemplo chamamos de {\tt  ex06\$fa} utilizando a função {\tt getGroups}.

Feito isto podemos rodar a função {\tt lme} que faz o ajuste do modelo.
  
\begin{verbatim}
> 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 
\end{verbatim}

Este modelo tem a variável {\tt forn} como efeito fixo e a variável
{\tt 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 {\tt 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 {\tt 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.
\begin{verbatim}
> 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 
\end{verbatim}

%{\bf NOTA:} Tiramos o intercepto do modelo com {\tt ~forn - 1}.
%Esta parametrização é mais adequada porque os coeficientes são
%diretamente interpretados como a média de cada fornecedor.
%Comparando com o modelo especificado apenas como {\tt ~forn} vemos que
%isto altera a interpretação dos coeficientes do modelo mas não afeta o
%teste do efeito de fornecedores na análise de variância.

Finalmente uma versão mais detalhada dos resultados pode ser obtida
com o comando {\tt summary}.
\begin{verbatim}
> 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 
\end{verbatim}

O próximo passo da seria fazer uma análise dos resíduos para
verificar os pressupostos, semalhante ao que foi feito no experimento
fatorial da aula anterior.
Vamos deixar isto por conta do leitor.


\section*{Exercícios}
\vspace{0.3in}
\begin{enumerate}
\item Faça uma análise completa do experimento hieráquico descrito
  acima.
\item Faça a análise do seguinte experimento descrito na apostila do curso, conduzido no arranjo hierárquico, com 3 repetições.

\begin{table}[ht]
\begin{center}
\caption{Diâmetro de peças produzidas por diferentes máquinas e operadores (mm).}
\begin{tabular}{l|cccccccc}\hline
        &  \multicolumn{2}{c}{Máquina 1}  & \multicolumn{2}{c} {Máquina 2} & \multicolumn{2}{c} {Máquina 3} & \multicolumn{2}{c} {Máquina 4} \\ \hline
Operador         & 1      &  2    &  1    &    2   &  1 & 2 &  1 &  2 \\ \hline
Amostra 1        & 30,66  & 30,85 & 30,97 &  30,67 &  31,03 & 30,76 & 30,77 & 31,03 \\ 
Amostra 2        & 30,78  & 31,00 & 30,97 &  30,75 & 31,01 & 30,69 &  30,69 & 31,09\\ 
Amostra 3        & 30,82  & 30,87 & 30,95 &  30,81 & 30,93 & 30,69 & 30,71 &  30,89 \\\hline
\end{tabular}
\end{center}
\end{table}

\end{enumerate}


