{\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*{Experimento no esquema fatorial com três fatores}
\indent

\subsection*{Objetivo da aula}

O objetivo dessa aula é 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


Nessa aula 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}



ex05.mnk <- tapply(resp, list(n,k), mean)

plot.default(n, resp, ty="n")

points(n[k=="0"], resp[k=="0"], col=1)
points(ex05.mnk[,1], pch="x", col=1, cex=1.5)

points(n[k=="1"], resp[k=="1"], col=1)
points(ex05.mnk[,2], pch="x", col=2, cex=1.5)


plot.default(k, resp, ty="n")

points(k[n=="0"], resp[n=="0"], col=1)
points(ex05.mnk[1,], pch="x", col=1, cex=1.5)

points(k[n=="1"], resp[n=="1"], col=1)
points(ex05.mnk[2,], pch="x", col=2, cex=1.5)



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

Não definido.





ex05.mnk <- tapply(resp, list(n,k), mean)

plot.default(n, resp, ty="n")

points(n[k=="0"], resp[k=="0"], col=1)
points(ex05.mnk[,1], pch="x", col=1, cex=1.5)

points(n[k=="1"], resp[k=="1"], col=1)
points(ex05.mnk[,2], pch="x", col=2, cex=1.5)


plot.default(k, resp, ty="n")

points(k[n=="0"], resp[n=="0"], col=1)
points(ex05.mnk[1,], pch="x", col=1, cex=1.5)

points(k[n=="1"], resp[n=="1"], col=1)
points(ex05.mnk[2,], pch="x", col=2, cex=1.5)



ex05.avn <- aov(resp ~ n/k)
summary(ex05.avn, split=list("n:k"=list(0=0, 1=1)))




ex05.avn <- aov(resp ~ n/k + bloc + p + n*p + p*k + n*p*k -k )

ex05.avk <- aov(resp ~ k/n + bloc + p + n*p + p*k + n*p*k -n )


summary(ex05.avn)


summary(ex05.avk, split=list("k:n"=list(0=c(1,2,3), 1=c(4,6,7))))

summary(ex05.avk, split=list("k:n"=list(0=c(1,3,4), 1=c(2,5,6), "p"=list() )))

summary(ex05.avn, split=list("k:n"=list(0=c(1,2,3), 1=c(4,6,7))))
summary(ex05.avn, split=list("k:n"=list(0=c(1,3,4), 1=c(2,5,6))))


summary(ex05.avn, split=list("n:k"=list(0=c(1,3,4), 1=c(2,5,6))))
summary(ex04.ave, split=list("esp:rec"=list(e1=c(1,3), e2=c(2,4))))



summary(ex05.avn, split=list("n:k"=list(0=c(1,2,3), 1=c(4,6,7))))

summary(ex05.avk, split=list("n:k"=list(0=c(1,3,4), 1=c(2,5,6))))
summary(ex05.avk, split=list("n:k"=list(0=c(1,2,3), 1=c(4,6,7))))



sapply(split(list(ex05$n[0]=c(1,3,4), ex05$n[1]=c(2,5,6)),mean)



