| This post was kindly contributed by Ridículas - go there to comment and to read the full post. |
Rendimento de grãos em função dos níveis de potássio e umidade do solo. Letras minúsculas comparam níveis de potassio em um mesmo nível de umidade, as letras maisculas fazem o inverso.
O R não é um aplicativo específico de análise de experimentos. Mesmo assim, ele oferece funções que permitem a análise dos mais diversos tipos de experimentos, ficando a cargo do usuário saber usar corretamente as funções/pacotes que o R oferece.
Pode-se dizer que os experimentos fatoriais são os mais comuns na ciência. Nas ciências agrárias são amplamente usados. Têm-se uma tradição de representar a interação entre os fatores (qualitativos) empregando testes de médias para níveis de um fator fixando níveis dos demais.
Nessa ridícula, apresento procedimentos para obter o desdobramento de interação em um fatorial duplo. Os dados são do rendimento de grãos de soja em função dos fatores umidade do solo e dose de fertilizante potássico. O experimento foi instalado em casa de vegetação com as unidades experimentais (vasos) agrupadas em blocos que controlaram variações de luminosidade/umidade/ventilação.
Os procedimentos incluem o ajuste do modelo, o teste sobre os efeitos dos fatores pela análise de variância, a checagem das pressuposições do modelo, o fatiamento das somas de quadrados, o teste de médias. No final são apresentados procedimentos para representar os resultados em gráfico e tabela de dupla entrada. É importante dizer que os fatores deveriam ser analisados como quantitativos, onde espera-se uma interpretação mais interessante sobre o fenômeno. Entretanto, o método abaixo é ilustrativo do procedimento para fatores qualitativos. Em uma futura oportunidade analisarei dos dados considerando os fatores como quantitativos.
#-----------------------------------------------------------------------------
# importa e prepara os dados
rend <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/rendimento.txt",
header=TRUE, sep="\t")
rend <- transform(rend, K=factor(K), A=factor(A), bloc=factor(bloc))
str(rend)
#-----------------------------------------------------------------------------
# análise gráfica
require(lattice)
xyplot(rg~K|A, groups=bloc, data=rend, type="b", auto.key=TRUE)
#-----------------------------------------------------------------------------
# ajuste do modelo e faz checagem
m0 <- aov(rg~bloc+A*K, data=rend)
summary(m0)
par(mfrow=c(2,2)); plot(m0); layout(1)
#-----------------------------------------------------------------------------
# desdobrando somas de quadrados para a variação de K dentro de A e vice-versa
m1 <- aov(rg~bloc+K/A, data=rend)
desAinK <- sapply(paste("K", levels(rend$K), sep=""), simplify=FALSE,
grep, x=names(coef(m1)[m1$assign==3]))
summary(m1, split=list("K:A"=desAinK))
m2 <- aov(rg~bloc+A/K, data=rend)
desKinA <- sapply(paste("A", levels(rend$A), sep=""), simplify=FALSE,
grep, x=names(coef(m2)[m2$assign==3]))
summary(m2, split=list("A:K"=desKinA))
#-----------------------------------------------------------------------------
# desdobrando a interação em testes de Tukey
require(agricolae)
require(plyr)
KinA <- sapply(levels(rend$A), simplify=FALSE,
function(a){
with(subset(rend, A==a),
HSD.test(rg, K,
DFerror=df.residual(m0),
MSerror=deviance(m0)/df.residual(m0)))
})
KinA <- ldply(KinA, NULL)
KinA$M <- gsub(" ", "", KinA$M, fixed=TRUE)
KinA$trt <- as.factor(as.numeric(as.character(KinA$trt)))
str(KinA)
AinK <- sapply(levels(rend$K), simplify=FALSE,
function(k){
with(subset(rend, K==k),
HSD.test(rg, A,
DFerror=df.residual(m0),
MSerror=deviance(m0)/df.residual(m0)))
})
AinK <- ldply(AinK, NULL)
AinK$M <- toupper(gsub(" ", "", AinK$M, fixed=TRUE))
AinK$trt <- as.factor(as.numeric(as.character(AinK$trt)))
str(AinK)
#-----------------------------------------------------------------------------
# combinando os data.frames para obter o gráfico e a tabela
KinA$comb <- paste(KinA$trt, KinA$.id)
AinK$comb <- paste(AinK$.id, AinK$trt)
desd <- merge(KinA, AinK, by=c("comb","comb"))
desd$let <- with(desd, paste(M.y, M.x, sep=""))
desd <- desd[,c("trt.x","trt.y","means.x","let")]
names(desd) <- c("K","A","means","let")
str(desd)
#-----------------------------------------------------------------------------
# gráfico das médias com as letras sobre as barras
#png("f005.png", w=500, h=250)
barchart(means~K|A, data=desd, horiz=FALSE, layout=c(3,1),
ylab="Rendimento de grãos (g/vaso)", ylim=c(NA,45),
xlab="Dose de potássio (mg/dm³)",
key=list(text=list("Níveis de umidade do solo (%)")),
sub=list("*letras minúsculas comparam médias em um mesmo nível de umidade.", cex=0.8, font=3),
strip=strip.custom(bg="gray90"), col=colors()[51],
panel=function(x, y, subscripts, ...){
panel.barchart(x, y, subscripts=subscripts, ...)
panel.text(x, y, label=desd[subscripts,"let"], pos=3)
})
#dev.off()
#-----------------------------------------------------------------------------
# tabela de dupla entrada com as médias das combinações A e K
require(reshape)
desd$means.let <- with(desd, paste(means, let))
cast(desd, K~A, value="means.let")
#-----------------------------------------------------------------------------
