====== Exemplo: Introdução a análise de experimentos (cont.) ======

Neste outro exemplo ainda mostramos exemplos simplificados de recursos para análise de experimentos.\\
A ênfase aqui continua sendo  ilustrar o uso da linguagem. Mais detalhes sobre recursos para análise serão fornecidos em outros tutoriais.


===== Experimento em blocos casualizados =====
Vamos agora analisar o experimento  em blocos ao acaso.
Os dados são fornecidos na tabela abaixo.

**Conteúdo de óleo de //S. linicola//, em percentagem, em vários estágios de crescimento (Steel \& Torrie, 1980, p.199).**
^Estágios ^  Blocos ^^^^ 
|                |       I       | II    |  III  |  IV |
|Estágio 1       |       4,4     | 5,9   | 6,0   | 4,1 |
|Estágio 2       |       3,3     | 1,9   | 4,9   | 7,1 |
|Estágio 3       |       4,4     | 4,0   | 4,5   | 3,1 |
|Estágio 4       |       6,8     | 6,6   | 7,0   | 6,4 |
|Estágio 5       |       6,3     | 4,9   | 5,9   | 7,1 |
|Estágio 6       |       6,4     | 7,3   | 7,7   | 6,7 |

Inicialmente vamos entrar com os dados no **R**.
Há várias possíveis maneiras de fazer isto.
Vamos aqui usar a função ''scan()'' e entrar com os dados por linha da tabela.
Digitamos o comando abaixo e e função ''scan()'' recebe os dados. Depois de digitar o último dado digitamos 
ENTER em um campo em branco e a função encerra a entrada de dados retornando para o //prompt// do programa.

**OBS: ** Note que, sendo um programa escrito na língua inglesa, os decimais devem ser indicados por ''.'' e não por vírgulas.

<code R>
y <- scan()
#!1: 4.4
#!2: 5.9
#!3: 6.0
#!...
#!24: 6.7
#!25: 
#!Read 24 items
</code>

Agora vamos montar um //data.frame// com os dados e os indicadores de blocos e tratamentos.
<code R>
ex02 <- data.frame(estag = factor(rep(1:6, each=4)), bloco=factor(rep(1:4, 6)), resp=y)
</code>

Note que usamos a função ''factor''  para indicar que as variáveis ''blocos'' e ''estag'' são níveis de fatores e não
valores numéricos.

Vamos agora explorar um pouco os dados.

<code R>
names(ex02)
summary(ex02)

attach(ex02)

plot(resp ~ estag + bloco)

interaction.plot(estag, bloco, resp)
interaction.plot(bloco, estag, resp)

ex02.mt <- tapply(resp, estag, mean)
ex02.mt
ex02.mb <- tapply(resp, bloco, mean)
ex02.mb

plot.default(estag, resp)
points(ex02.mt, pch="x", col=2, cex=1.5)

plot.default(bloco, resp)
points(ex02.mb, pch="x", col=2, cex=1.5)
</code>

Nos gráficos e resultados acima procuramos captar os principais aspectos dos dados bem como verificar se não há interação entre blocos e tratamentos, o que não deve acontecer neste tipo de experimento.

A seguir vamos ajustar o modelo e obter outros resultados, incluindo a análise de resíduos e testes para verificar a validades dos pressupostos do modelo.
<code R>
ex02.av <- aov(resp ~ bloco + estag)
anova(ex02.av)

names(ex02.av)

par(mfrow=c(2,2))
plot(ex02.av)

par(mfrow=c(2,1))
residuos <- (ex02.av$residuals)

plot(ex02$bloco,residuos)
title("Resíduos vs Blocos")

plot(ex02$estag,residuos)
title("Resíduos vs Estágios")

par(mfrow=c(2,2))
preditos <- (ex02.av$fitted.values)
plot(residuos,preditos)
title("Resíduos vs Preditos")
respad <- (residuos/sqrt(anova(ex02.av)$"Mean Sq"[2]))
boxplot(respad)
title("Resíduos Padronizados")
qqnorm(residuos,ylab="Residuos", main=NULL)
qqline(residuos)
title("Grafico Normal de \n Probabilidade dos Resíduos")
</code>

Teste para normalidade
<code R>
shapiro.test(residuos)
</code>
Testando a não aditividade:
Primeiro vamos extrair coeficientes de tratamentos e blocos
<code R>
ex02.av$coeff
bl <- c(0, ex02.av$coeff[2:4])
tr <- c(0, ex02.av$coeff[5:9])
bl
tr
</code>

e agora criar um novo termo e testar sua significancia na ANOVA
<code R>
bltr <- rep(bl, 6) * rep(tr, rep(4,6))

ttna <- update(ex02.av, .~. + bltr)
anova(ttna)
</code>

Os resultados acima indicam que os pressupostos estão obedecidos para este conjunto de dados e a análise de variância é válida.
Como foi detectado efeito de tratamentos vamos proceder fazendo um teste de comparações múltiplas e encerrar as análises desanexando o objeto do caminho de procura.

<code R>
ex02.tk <- TukeyHSD(ex02.av, "estag", ord=T)
ex02.tk
plot(ex02.tk)

detach(ex02)
</code>