\section*{Aula 01: Experimentos em esquema fatorial}
\indent

\subsection*{Objetivo dessa aula}
\indent

O objetivo dessa aula é realizar a análise de variância de experimentos fatoriais. Será utilizado o Software R. 

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


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

Nesta aula iremos detalhar e discutir a análise do experimento
fatorial $2 \times 3$.

Este experimento, descrito na apostila do curso de Planejamento de
Experimentos II, comparou a resposta (altura) de mudas a diferentes
recipientes e espécies de eucalipto.

No restante destas notas as linhas que começam com 
o símbolo {\tt >} são comandos a serem digitados no \R{}.
Outros textos com a font {\tt typewriter como esta} são saídas
produzidas pelo programa.

\begin{enumerate}
\item{\bf Lendo os dados}

  \htmladdnormallink{Clique aqui}{../dados/exemplo04.txt} para ver e
  copiar o arquivo com conjunto de dados. 
  
  A seguir vamos ler (importar) os dados para \R{} com o comando {\tt read.table}:
\begin{verbatim}
> ex04 <- read.table("exemplo04.txt", header=T)
> ex04
\end{verbatim}

Lembre-se de fornecer o caminho onde você salvou os dados.


Antes de começar as análises 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 para o \R{} que exiba um
rápido resumo dos dados.

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

Nos resultados acima vemos que o objeto {\tt ex04} que contém os dados
tem 24 linhas (observações) e 3 colunas (variáveis).
As variáveis tem nomes {\tt rec}, {\tt esp} e {resp}, sendo que as
duas primeiras são {\it fatores} enquanto {\tt resp} é uma variável
numérica, que neste caso é a variável resposta.
O objeto {\tt ex04}  foi incluído no caminho de procura usando o comando {\tt attach} para facilitar a digitação.


\subsection*{Análise descritiva}
\indent

Inicialmente vamos obter um resumo de nosso conjunto de dados usando a
função {\tt summary}.

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

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

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 {\tt mean} por {\tt var}
para calcular a variância de cada grupo, e por {\tt 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~\ref{fig:ex04inter.eps}.  

\begin{verbatim}
> par(mfrow=c(1,2))
> interaction.plot(rec, esp, resp)
> interaction.plot(esp, rec, resp)
\end{verbatim}

\begin{figure}[ht!]
\centerline{
\mbox{
\includegraphics[width=8.0cm,height=8.0cm]{../CE214/internet/figuras/ex04inter1.eps}
\includegraphics[width=8.0cm,height=8.0cm]{../CE214/internet/figuras/ex04inter2.eps}
}}
\vspace{-0.3in}
\caption{Gráficos de interação entre os fatores.}
\label{fig:ex04inter.eps}
\end{figure}

Pode-se usar o \R{} para obter outros tipos de gráficos de acordo com
o interesse de quem está analisando os dados.
Por exemplo, os comandos abaixo ilustram outros tipos de gráficos.
Experimente estes comandos, verifique os gráficos produzidos e
certifique-se que você entendeu cada comando.

\begin{verbatim}
> plot.default(rec, resp, ty="n")
> points(rec[esp=="e1"], resp[esp=="e1"], col=1)
> points(ex04.m[,1], pch="x", col=1, cex=1.5)
> points(rec[esp=="e2"], resp[esp=="e2"], col=2)
> points(ex04.m[,2], pch="x", col=2, cex=1.5)

> plot.default(esp, resp, ty="n")
> points(esp[rec=="r1"], resp[rec=="r1"], col=1)
> points(ex04.m[1,], pch="x", col=1, cex=1.5)
> points(esp[rec=="r2"], resp[rec=="r2"], col=2)
> points(ex04.m[2,], pch="x", col=2, cex=1.5)
> points(esp[rec=="r3"], resp[rec=="r3"], col=3)
> points(ex04.m[3,], pch="x", col=3, cex=1.5)

> coplot(resp ~ rec|esp)
> coplot(resp ~ esp|rec)
\end{verbatim}

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

Seguindo o modelo adequado, 
a análise de variância para esse experimento inteiramente casualizado
em esquema fatorial pode ser obtida com o comando:
\begin{verbatim}
> ex04.av <- aov(resp ~ rec + esp + rec * esp)
\end{verbatim}

Entretanto o comando acima pode ser simplificado produzindo os mesmos
resultados com o comando
\begin{verbatim}
> 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 
\end{verbatim}

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
{\tt rec:esp}.
A análise acima mostra que este efeito é significativo, confirmando o
que verificamos nos gráficos de interação vistos anteriormente.

O objeto {\tt ex04.av} guarda todos os resultados da análise e pode
ser explorado por diversos comandos.
Por exemplo a função {\tt 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. 
\begin{verbatim}  
> 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
\end{verbatim}

O objeto {\tt ex04.av} possui vários elementos que guardam
informações sobre o ajuste.
\begin{verbatim}
> 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" 
\end{verbatim}

O comando {\tt class} mostra que o objeto {\tt ex04.av} pertence às
classes {\tt aov} e {\tt lm}.
Isto significa que devem haver {\it 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 
{\tt summary(ex04.av)}.
Existe uma função chamada {\tt summary.aov} que foi utilizada já que o
objeto é da classe {\tt aov}.
Iremos usar mais este mecanismo no próximo passo da análise.

\subsection*{Análise de resíduos}
\indent

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~\ref{fig:ex04res1.eps} com o comando {\tt plot}.
\begin{verbatim}
> par(mfrow=c(2,2))
> plot(ex04.av)
\end{verbatim}

\begin{figure}[htbp]
\centerline{\includegraphics[width=12.0cm,height=12.0cm]{../CE214/internet/figuras/ex04res1.eps}}
\caption{Gráficos de resíduos produzidos automaticamente pelo \R{}.}
\label{fig:ex04res1.eps}
\end{figure}

Os gráficos permitem uma análise dos resíduos que auxiliam no
julgamento da adequacidade do modelo.
Evidentemente você não precisa se limitar aos gráficos produzidos
automaticamente pelo \R{} -- você pode criar os seus próprios gráficos
muito facilmente. Neste gráficos você 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.
\begin{verbatim}
> 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")
\end{verbatim}

Além disto há alguns testes já programados. Como exemplo vejamos e
teste de Shapiro-Wilk para testar a normalidade dos resíduos.
\begin{verbatim}
> shapiro.test(residuos)

        Shapiro-Wilk normality test

data:  residuos 
W = 0.9293, p-value = 0.09402
\end{verbatim}

\item{\bf Desdobrando interações}

Conforme 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 \verb+/+ que indica efeitos aninhados.
Desta forma podemos desdobrar os efeitos de espécie dentro de cada
recipiente e vice-versa conforme mostrado a seguir.
\begin{verbatim}
> 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 
\end{verbatim}


\subsection*{Teste de Tukey para comparações múltiplas}
\indent

Vejamos por exemplo duas formas de usar o {\it Teste de Tukey}, a
primeira usando uma implementação com a função {\tt TukeyHSD} e uma
segunda fazendo ops cálculos necessários com o \R.

Poderíamos simplesmente digitar:
\begin{verbatim}
> ex04.tk <- TukeyHSD(ex04.av)
> plot(ex04.tk)
> ex04.tk
\end{verbatim}

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íveis dos fatores.
\begin{verbatim} 
> 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
\end{verbatim}

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

\end{enumerate}

{\large EXERCÍCIO}
\begin{enumerate}
  \itemsep 1ex
\item Analise o experimento de adubação NPK da apostila de
  Planejamento de Experimentos II.\\
  \htmladdnormallink{Clique aqui}{../dados/exemplo05.txt} para ver e
  copiar o arquivo de dados.
{\bf sugere-se realizar essa análise seguindo os comandos da aula 1.}
\end{enumerate}


\section*{Exercício: Experimento no esquema fatorial com três fatores}
\indent

\subsection*{Objetivo do exercício}

O objetivo desse exercício é realizar a análise de variância de um experimento fatorial com três fatores, utilizando o software R.


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


Nesse exercício será analisado o experimento fatorial $2^3$ conduzido em um delineamento em blocos completos casualizados apresentado nas notas de aula. Esse experimento refere-se a adubação de pés de café, utilizando-se os nutrientes N (nitrogênio), P (fósforo) e K (potássio) com dois níveis cada. 

O arquivo de dados pode ser obtido \htmladdnormallink{aqui.}{../dados/exemplo05.txt}.


\begin{verbatim}
ex05<-read.table("exemplo05.txt",header=T)
\end{verbatim}


Lembre-se que cada coluna indicadora dos fatores e blocos deve ser um fator no R. Por essa razão, os vetores n, p, k e bloc precisam ser transformados para fatores.

\begin{verbatim}
> is.factor(ex05$n)
[1] FALSE
> is.factor(ex05$p)
[1] FALSE
> is.factor(ex05$k)
[1] FALSE
> is.factor(ex05$bloc)
[1] FALSE
\end{verbatim}


Para facilitar a digitação, pode-se incluir no caminho de procura o objeto {\tt ex05} com o comando {\tt attach(ex05)}.

Agora, deve-se transformar as colunas em fatores

\begin{verbatim}
 ex05$n<-as.factor(ex05$n)
 ex05$p<-as.factor(ex05$p)
 ex05$k<-as.factor(ex05$k)
 ex05$bloc<-as.factor(ex05$bloc)
\end{verbatim}



Com o comando {\tt summary} obtém-se um resumo dos dados do objeto {\tt ex05}.

\begin{verbatim}
> summary(ex05)
 n      p      k      bloc       resp     
 0:24   0:24   0:24   1:8   Min.   :1876  
 1:24   1:24   1:24   2:8   1st Qu.:3171  
                      3:8   Median :3743  
                      4:8   Mean   :3863  
                      5:8   3rd Qu.:4521  
                      6:8   Max.   :6829  
\end{verbatim}


Explorando a variável resposta com o comando {\tt tapply}

\begin{verbatim}

> ex05.mn <- tapply(resp, n, mean)
> ex05.mn

> ex05.mp <- tapply(resp, p, mean)
> ex05.mp

> ex05.mk <- tapply(resp, k, mean)
> ex05.mk
\end{verbatim}
 


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


A análise de variância para o experimento fatorial conduzido no delineamento em blocos completos casualizados pode ser obtida com o comando:

\begin{verbatim}
> ex05.av <- aov(resp ~ bloc + n + p + k + n*p + n*k + p*k + n*p*k)
\end{verbatim}

ou

\begin{verbatim}
> ex05.av <- aov(resp ~ bloc + n*p*k)
> summary(ex05.av)
            Df   Sum Sq  Mean Sq F value    Pr(>F)    
bloc         5  2134332   426866  0.7127 0.6180190    
n            1 10229610 10229610 17.0797 0.0002123 ***
p            1   194438   194438  0.3246 0.5724714    
k            1  6279257  6279257 10.4841 0.0026366 ** 
n:p          1   553196   553196  0.9236 0.3431151    
n:k          1  8728749  8728749 14.5738 0.0005274 ***
p:k          1   474218   474218  0.7918 0.3796409    
n:p:k        1   288765   288765  0.4821 0.4920416    
Residuals   35 20962662   598933                      
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 
\end{verbatim}


O uso da função {\tt model.tables} permite a obtenção de estatísticas para cada um dos fatores do modelo e a combinação entre eles.


\begin{verbatim}
> ex05.mt <- model.tables(ex05.av, ty="means")
> ex05.mt
\end{verbatim}


\subsection{Análise de resíduos}
\indent

Uma forma de analisar os resíduos é atraves do comando {\tt plot(ex05.av)}.  Esta função realiza a análise de resíduos do objeto ex05.av. 

\begin{verbatim}
>par(mfrow=c(2,2))
>plot(ex05.av)
\end{verbatim}


Uma outra maneira de realizar a análise de resíduos é gerando os gráficos individualmente. Desse modo existe uma maior liberdade para 

\begin{verbatim}
> par(mfrow=c(1,1))
>  residuos <- resid(ex05.av)

>  plot.default(ex05$n, residuos)
> title("Resíduos vs N")

>  plot.default(ex05$p, residuos)
> title("Resíduos vs P")

>  plot.default(ex05$k, residuos)
> title("Resíduos vs K")


> par(mfrow=c(2,2))
> preditos <- (ex05.av$fitted.values)
> plot(residuos, preditos)
> title("Resíduos vs Preditos")

s2 <- sum(resid(ex05.av)^2)/ex05.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")
\end{verbatim}



O teste de Shapiro-Wilk também pode ser utilizado para avaliar a normalidade dos resíduos.

\begin{verbatim}
shapiro.test(residuos)

        Shapiro-Wilk normality test

data:  residuos 
W = 0.9619, p-value = 0.1202
\end{verbatim}


Uma vez verificado os pressupostos do modelo de análise de variância, pode-se continuar a análise. Nesse exemplo, observa-se que o efeito da interação NK foi significativo. Por esse motivo, faz-se necessário efetuar o desdobramento dos graus de liberdade da interação.


\subsection*{Análise da Interação}
\indent


Inicialmente, faz-se o gráfico da interação para visualizar-se o comportamento entre as combinações dos diferentes níveis dos diferentes fatores.

\begin{verbatim}
par(mfrow=c(2,2))

interaction.plot(n, p, resp)
title("Interação N vs P")
interaction.plot(n, k, resp)
title("Interação N vs K")
interaction.plot(p, k, resp)
title("Interação P vs K")
\end{verbatim}


\subsection*{Desdobramento da interação NK:}
\indent

Estudaremos os níveis de N dentro de K.

Inicialmente, deve-se utilizar o modelo com o efeito de N aninhado em K.

\begin{verbatim}
ex05.avn<-aov(resp~bloc+k/n+n*p*k-n,data=ex05)
\end{verbatim}

Observe que apenas nos interessa interpretar os efeitos do desdobramento.

\begin{verbatim}
summary(ex05.avn,split=list("k:n"=list(N.d.K0=1, N.d.K1=2)))
\end{verbatim}


Para estudar os níveis de K dentro de N, basta considerar $n/k$.
