Fatorial com desdobramento da interação, gráfico e tabela de médias

May 21, 2011
By

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")

#-----------------------------------------------------------------------------

Tags: , , , , , , , , ,

Comments are closed.