Fatorial duplo com dois tratamentos adicionais

December 3, 2012
By

This post was kindly contributed by Ridículas - go there to comment and to read the full post.

Número de sementes em 40 infectadas por Fusarium (tranformado) em função da aplicação de fungicida em combinação com polímero.

Experimentos fatoriais são aqueles em que dois ou mais fatores estão presentes. Nos fatoriais incompletos nem todas as combinações possíveis estão presentes. Algumas das combinações podem não ser de interesse e por isso não são realizadas, alguma combinação é acidentalmente perdida durante o experimento, certas combinações não identificam novos níveis, como por exemplo aplicação do produto A e B na dose 0 que, apesar dos rótulos diferentes, são a mesmo tratamento.

Um caso comum de fatorial incompleto são os fatoriais (completos) com a presença de tratamentos adicionais. Experimentos fatorais são realizados para verificar se os fatores interagem e daí escolher combinação de níveis mais promissora para o processo estudado. Adiciona-se tratamentos adicionais ao experimento quando se deseja comparar o promissor aos de referência ou extremos.

Esse é um tipo de ensaio típico na avaliação do efeito de insumos agrícolas (fungicidas, insecticidas e herbicidas). Os ensaios testam produtos e suas doses e precisam ter as chamadas testemunhas para avaliar a eficiência das combinações geradas. Devido não ser um delineamento de níveis completamente cruzados a análise requer mais atenção.

Nessa matéria faço análise de um fatorial com dois tratamentos adicionais. Os dados analisados foram retirados do Zimmermann (2004) de um experimento que estudou três fungicidas (benlate, captam, derosal) em três formas de aplicação (puro, antes de polímero, misturado com polímero). Adicionou-se dois tratamentos que foram a testemunha e aplicação do polímero puro. A resposta
avaliada foi o número de sementes infectadas com Fusarium em um total de 40 sementes.

O código está comentado para informar o leitor. São realizados a importação e reordenação dos níveis, a análise de variância com decomposição das somas de quadrados, contrastes entre efeitos, e testes sobre as médias ajustadas. Foi dada prioridade para os procedimentos matriciais de análise que podem ser adaptados para situações de desbalanceamento. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# definições da sessão

require(lattice)
require(latticeExtra)
require(grid)
require(MASS)
require(doBy)
require(gmodels)
require(multcomp)

#-----------------------------------------------------------------------------
# leitura

fu <- read.table(
  "http://www.leg.ufpr.br/~walmes/data/zimmermann_fusarium.txt",
                 header=TRUE, sep="\t")
str(fu)

#-----------------------------------------------------------------------------
# modificar ordem dos níveis, primeiro níveis fatoriais, depois adicionais

levels(fu$fungicida)               # ordem correta
levels(fu$aplicacao)               # precisa alterar
levels(fu$aplicacao)[c language="(4,1,2,3,5)"][/c] # para ficar assim

fu$aplicacao <- factor(fu$aplicacao,
                       levels=levels(fu$aplicacao)[c language="(4,1,2,3,5)"][/c])
levels(fu$aplicacao) # agora ordem correta

#-----------------------------------------------------------------------------
# análise exploratória

xtabs(~fungicida+aplicacao, fu)

# os níveis não se cruzam completamente, portanto é um fatorial incompleto
# no caso 3x3+2, fatorial com 2 tratamentos adicionais

xyplot(infec40~tratamentos, fu) # resposta do tipo contagem de sucessos
xyplot(infec40~tratamentos, fu, jitter.x=TRUE)

#-----------------------------------------------------------------------------
# Vamos primeito analisar fora da estrutura fatorial, ou seja, considerando
# 3x3+2 = 11 níveis de um único fator.
# Vamos considerar distribuição para reproduzir os resultados do livro
# do Zimmermann (2004). Todavia, não é a análise mais correta devido os dados
# serem discretos. A distribuição binomial parece ser mais indicada.

m0 <- lm(infec40~tratamentos, data=fu)
par(mfrow=c(2,2)); plot(m0); layout(1) # não tão ruim, tentar boxcox?

m0 <- lm(infec40+1~tratamentos, data=fu) # soma 1 para eliminar os 0
boxcox(m0) # aponta log

m1 <- lm(log(infec40+1)~tratamentos, data=fu)
par(mfrow=c(2,2)); plot(m1); layout(1) # alguma melhoria?

# Zimmermann usa asin(sqrt(x/40)).
m2 <- lm(asin(sqrt(infec40/40))~tratamentos, data=fu)
par(mfrow=c(2,2)); plot(m2); layout(1) # alguma melhoria?

#-----------------------------------------------------------------------------
# Visualmente todos os resíduos são semelhantes, apresentam fugas.

lapply(list(m0, m1, m2), function(i) shapiro.test(residuals(i)))
lapply(list(m0, m1, m2), function(i) ks.test(scale(residuals(i)), "pnorm"))

# Pelos p-valores, a tranformação Box-Cox foi melhor, embora os resíduos ainda
# não apresentem um comportamento satísfatório.
# O foco não é testar transformações, então vamos seguir o que Zimmermann fez
# e trabalhar com asin(sqrt(x/40)).

#-----------------------------------------------------------------------------
# análise exploratória com a variável transformada

fu$y <- asin(sqrt(fu$infec40/40)) # cria a transformada

xyplot(y~tratamentos, fu, jitter.x=TRUE)

xyplot(y~fungicida, groups=aplicacao, fu, jitter.x=TRUE,
       type=c("p","a"), auto.key=TRUE)

# Apresenta indicativos de efeito de interação e que os adicionais
# são piores que o os fatoriais.

#-----------------------------------------------------------------------------
# análise na estrutura fatorial 3x3+2

m3 <- lm(y~fungicida*aplicacao, data=fu)
anova(m3)

# Os graus de liberdade são o que são pois:
# fungicida 4: são 5 niveis (3 fatoriais, 2 adicionais), 5-1=4;
# aplicacao 4: são 5, mas 2 estão confundidos (os adicionais), 5-2-1=2;
# fun:apli  4: (3-1)*(3-1)=4;

# O que queremos é a anova com somas de quadrados particionadas da seguinte
# forma
# fonte de variação     gl
# fungicidas             2
# aplicacao              2
# fun:apli               4
# fatorial vs adicionais 1
# polimero vs testemunha 1

#-----------------------------------------------------------------------------
# Na minha opinião, obter esse quadro de análise de variância é desperdício
# de tempo pois:
# * não tem sentido prático comparar o efeito médio dos fatoriais contra o
#   efeito médio dos adicionais, principalmente sabendo que a interação
#   fungicida:aplicacao foi significativa;
# * a interação fungicida:aplicacao significativa torna não interpretável
#   o teste para os efeitos principais;
# * polímero vs testemunha tem só um grau de liberdade, posso comprar os
#   efeitos por um testes t;

coef(m3) # tem NA para os efeitos não estimados (não presentes)
efstm <- colnames(vcov(m3)) # nomes dos efeitos que foram estimados
b <- coef(m3)[efstm] # só os estimados
b

# matriz de coeficientes para estimação das médias ajustadas
X <- popMatrix(m3, effect=c("fungicida", "aplicacao"))
str(X)

Xn <- attr(X, "grid") # 25 linhas, 5x5=25 (faz pensando que é completo)
Xo <- unique(fu[,c("fungicida","aplicacao")]) # que ocorrem (3x3+2)
Xo

#-----------------------------------------------------------------------------
# calculando as médias ajustadas matricialmente

rownames(X) <- apply(Xn, 1, paste, collapse=":")
on <- apply(Xo, 1, paste, collapse=":")
X <- X[on,efstm] # só linhas e colunas que correspondem à efeitos que ocorrem

X%*%b                                             # são as médias ajustadas
with(fu,
     tapply(y, list(fungicida, aplicacao), mean)) # são médias amostrais

# elas coincidem por ser um experimento regular (balanceado/ortogonal)

#-----------------------------------------------------------------------------
# fazendo contrastes
# não funcionam com presenta de NA
# contrast::contrast()
# multcomp::glht()
# para isso usar a gmodels::estimable()

m4 <- aov(formula(m3), data=fu) # aov não NA em coef

# polímero contra a testemunha
p.vs.t <- matrix(X[10,]-X[11,], nrow=1)
rownames(p.vs.t) <- "p.vs.t"
estimable(m4, cm=p.vs.t)

# fatorial contra adicional
f.vs.a <- matrix(apply(X[1:9,], 2, mean)-apply(X[10:11,], 2, mean), 1)
rownames(f.vs.a) <- "f.vs.a"
estimable(m4, cm=f.vs.a)

estimable(m4, cm=rbind(p.vs.t, f.vs.a)) # os dois de uma só vez

# E aí, complicado fazer tudo isso, né? Qual a sensação em saber de que
# pouco disso vai ser útil? Se deu interação teremos que estudar níveis
# de um fator fixando os níveis do outro.

#-----------------------------------------------------------------------------
# fazendo a decomposição das somas de quadrados dentro da anova

cf <- contrasts(C(fu$fungicida, helmert))
cf # o segredo está em montar os contrastes corretamente
cf[,3] <- c(-1,-1,-1,1.5,1.5) # fatorial contra adicional
cf[,4] <- c(0,0,0,-1,1)       # polímero vs testemunha
cf
contrasts(fu$fungicida) <- cf

m5 <- aov(y~fungicida*aplicacao, data=fu)
summary(m5, expand=FALSE,
        split=list("fungicida"=list(
                     "fungicida"=1:2,
                     "fat.vs.adi"=3,
                     "pol.vs.tes"=4)))

# Esse procedimento depende de balanceamento pois as somas de quadrados são
# sequênciais.

#-----------------------------------------------------------------------------
# Uma vez que deu interação é de prática fazer o desdobramento que consiste
# em estudar o efeito de aplicação dentro dos níveis de fungicida e
# vice-versa. Além do mais, por ter tratamentos adicionais, pode-se comparar
# com os níveis dos fatoriais. Vamos listas o número de hipóteses que teremos
# que fazer:
# * 3 níveis de aplicação em cada fungicida: 3*choose(3,2)=9
# * 3 níveis de fungicida em cada aplicação: 3*choose(3,2)=9
# * testemunha contra cada um dos 9 níveis do 3x3: 9
# * polimero contra cada um dos 9 níveis do 3x3: 9
# Com isso são 4*9=36 comparações. Já que é assim, vou fazer todos
# contra todos, choose(11, 2)=55, que o usuário pode olhar a que lhe
# interessa, afinal de contas, de 36 para 55 não têm diferença

mm <- model.matrix(m4)
dim(mm)
mm <- mm[,efstm]

# modelo de classe lm equivalente à m5, sem NA
m6 <- lm(y~-1+mm, data=fu)
anova(m6)

cbind(coef(m4), coef(m6)) # são o mesmo modelo

#-----------------------------------------------------------------------------
# médias ajustadas (ma)

cima <- confint(glht(m6, linfct=X)) # intervalos de confiança para ma
cima # tem cobertura *global* de 95% (e não indivídual)

str(cima)
cima$confint

result <- cbind(Xo, cima$confint)
str(result)

segplot(fungicida:aplicacao~lwr+upr, data=result,
        xlab="Plantas infectadas (asin(sqrt(x/40))",
        ylab="Fungicida:aplicação",
        centers=Estimate, draw.bands=FALSE,
        segments.fun=panel.arrows, ends="both",
        angle=90, length=1, unit="mm")

#-----------------------------------------------------------------------------
# carrega algumas funções que serão úteis para comparar médias e fazer
# gráficos com intervalos de confiança

objs <- ls()
objs <- c(objs, "apc", "prepanel.cbH", "panel.cbH")
source("http://www.leg.ufpr.br/~walmes/ridiculas/ridiculas_functionsi.R",
       encoding="latin1")

rm(list=ls()[!(ls()%in%objs)])
ls()

#-----------------------------------------------------------------------------
# comparando as médias ajustadas, choose(11, 2), contrastes de Tukey

str(X)
Xc <- apc(X)
str(Xc)

#-----------------------------------------------------------------------------
# fazendo as comparações com a glht, método para correção de p-valor
# fdr, o single-step demora muito

c0 <- summary(glht(m6, linfct=Xc), test=adjusted(type="fdr"))
c0

c0$focus <- "comparacoes" # para poder usar a cld()
cld(c0)

result$cld <- cld(c0)$mcletters$Letters

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

segplot(fungicida:aplicacao~lwr+upr, data=result,
        xlab="Plantas infectadas (asin(sqrt(x/40))",
        ylab="Fungicida:aplicação",
        centers=Estimate, draw.bands=FALSE,
        segments.fun=panel.arrows, ends="both",
        angle=90, length=1, unit="mm",
        panel=function(x, y, z, centers, subscripts, ...){
          panel.segplot(x, y, z, centers=centers, subscripts=subscripts, ...)
          panel.text(centers, z,
                     paste(format(centers, digits=2),
                           result$cld[subscripts], sep=" "), pos=3)
        })

#-----------------------------------------------------------------------------
# gráfico final

result$desloc <- 0.05*c(0,-1,1,0,0)[match(result$aplicacao,
                                          levels(fu$aplicacao))]

#png("f048.png", 500, 400)
xyplot(y~fungicida, groups=aplicacao, data=fu, jitter.x=TRUE, amount=0.05,
       xlab="Fungicida", ylab=expression(asin(sqrt(x/40))),
       prepanel=prepanel.cbH, ly=result$lwr, uy=result$upr,
       auto.key=list(title="Aplicação", cex.title=1.2, columns=2,
         type="o", divide=1, lines=TRUE, points=FALSE))+
  as.layer(xyplot(Estimate~fungicida, groups=aplicacao, data=result, type="l",
                  ly=result$lwr, uy=result$upr,
                  desloc=result$desloc,
                  cty="bars",
                  panel.groups=panel.cbH,
                  panel=panel.superpose))
trellis.focus("panel", 1, 1, h=FALSE)
x <- as.numeric(result$fungicida)+result$desloc
y <- result$Estimate
lab <- paste(format(y, digits=2), result$cld, sep=" ")
a <- paste(rep("0", max(nchar(lab))), collapse="")
grid.rect(x=unit(x, "native"), y=unit(y, "native"),
          width=unit(1, data=a, "strwidth"),
          height=unit(1.2, data=a, "strheight"),
          hjust=-0.1, gp=gpar(fill="gray90", col=NA, fontsize=10))
grid.text(lab, x=unit(x, "native"), y=unit(y, "native"),
          hjust=-0.2, gp=gpar(col="black", fontsize=10))
trellis.unfocus()
#dev.off()

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


Tags: , , ,

Comments are closed.