Experimento fatorial duplo com um tratamento adicional

August 9, 2011
By

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

Produção de cultivares em função do nível de fertilizante (95% de confiança). As linhas horizontais em cinza representam o desempenho da cultivar de referência. Linhas pretas representam o desempenho das demais como função do nível de fertilizante. A seta aponta o nível de fertilizante que confere a máxima produção da cultura. Análise baseada em dados simulados.

Em algumas áreas do conhecimento é quase regra fazer experimentos incluindo uma tratamento adicional. Um experimento fatorial é aquele que possui duas ou mais fontes de variação. Em geral, têm-se presente no experimento todos os níveis possíveis das combinações entre os níveis desses fatores, ou em outras palavras, todos os níveis possíveis do produto cartesiano. Os fatores podem ser qualitativos (ou categóricos) ou quantitativos (ou métricos). Um tratamento adicional é um nível que não corresponde à estrutura fatorial, ou seja, é um nível que não veio da combinação dos níveis dos fatores. Para ficar mais claro, segue alguns exemplos:

  • Para um determinado município, sabe-se que a melhor cultivar de soja é a A com produção máxima ao receber a aplicação de 100 kg de fertilizante. Uma empresa está estudando 5 novas cultivares para o município e precisa definir a dose de fertilizante para a máxima produção. Além de definir a dose ótima, deseja-se comparar a produção com a cultivar A. Para isso faz-se um experimento com os 5 níveis de cultivar de soja (B, C, D, E, F) e um número de níveis suficientes de fertilizante que permita determinar o ponto de máximo, diga-se, 5 níveis (50, 75, 100, 125, 150 kg). O experimento terá 25 níveis (5×5) vindos da porção fatorial e um nível de referência, chamado de testemunha ou controle, que é cultivar A com 100 kg de fertilizante. Esse experimento é um fatorial 5×5+1.
  • Em um experimento com herbicidas, deseja-se definir a dose mínima de aplicação para controlar plantas daninhas. São 4 níveis de novos herbicidas (A, B, C, D) e 6 níveis de concentração (0,1; 0,5; 1; 2; 5; 10%). A variável resposta de interesse é produtividade da cultura. Para verificar se os produtos causam toxidade à cultura, tem-se um nível que é capina manual (testemunha positiva), que é o controle de plantas daninhas sem produto químico, onde se espera a maior produtividade. Para ver a eficiência desses herbicidas, no estudo se incluí um nível que é sem a capina (testemunha negativa), onde se espera a menor produtividade. Este experimento é um fatorial 4×6+2.

Perceba que pelas descrições acima ficaram evidentes às hipóteses de cada experimento. No primeiro é definir a dose ótima de cada cultivar e depois comparar o desempenho das cultivares na melhor dose com a cultivar A. No segundo é verificar se os herbicidas causam toxides à cultura (comparar com a testemunha positiva), qual a perda de produção quando não se faz o controle de plantas daninhas (comparar com a testemunha negativa) e definir a dose mínima para controlar plantas daninhas para cada herbicida.

O modelo estatístico que representa o primeiro experimento é

Y_{ijk} = M+C_i + F_j + CF_{ij} + \epsilon_{ijk} \quad \text{e} \quad Y_{k} = M + A + \epsilon_{k},

em que Y são os valores observados, M um valor inerente à todas as observações, C representa o efeito dos níveis de cultivar, F representa o efeito dos níveis de fertilizante, CF representa o efeito da combinação dos níveis de cultivar e fertilizante, A é o efeito do tratamento adicional e \epsilon é um erro aleatório.

Veja que até aqui foi tudo muito comum, com exceção ao fato de ser necessário escrever o modelo com duas expressões, uma para descrever a parte fatorial e outra para a parte adicional. No R não conseguimos declarar o modelo com duas fórmulas e é justamente sobre isso que vou tratar nesse post.

Podemos representar o modelo acima como um fatorial incompleto. Pense assim: são 6 nível de cultivar e apenas 5 deles combinados com 5 níveis de fertilizante. É assim que vamos declarar o modelo no R. Com algumas manipulações, iremos repartir a soma de quadrados da porção fatorial da porção adicional. A repartição será ortogonal que nem sempre representa as hipóteses que nós estabelecemos. Então por que fazemos essa repartição? Só porque ela é ortogonal.

Vou usar dados simulados para exemplificar o procedimento. O efeito dos níveis de fertilizante serão representados por um polinômio de segundo grau, apenas para que apresente ponto de máximo. No script abaixo eu gero os dados, faço os ajuste dos modelos, apresento os quadros de análise de variância, faço a repartição da soma de quadrados, a estimação da dose ótima de fertilizante para a resposta máxima e a comparação com a testemunha. O código foi comentado para orientar o leitor. Se você quer saber como é a análise do fatorial com dois tratamentos adicionais, deixe um comentário nesse post. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# gerando dados conforme o modelo estatístico 5x5+1, depois explico o fert=101

adi <- expand.grid(cult=gl(1,5,la=LETTERS[1]), fert=101)
fat <- expand.grid(cult=gl(5,5,la=LETTERS[2:6]), fert=seq(50,150,25))
da <- rbind(adi, fat)

theta <- c(c(-194.29, -197.26, -197.85, -203.03, -190.20, -190.45),
           c(9.1797, 8.2686, 8.6437, 9.3438, 8.8773, 8.1872),
           c(-0.03382, -0.03479, -0.03632, -0.03341, -0.03597, -0.03675))
X <- model.matrix(~-1+cult/(fert+I(fert^2)), data=da)
da$eta <- X%*%matrix(theta)
set.seed(2)
da$y <- da$eta+rnorm(nrow(da),0,30)

require(lattice)
xyplot(y~fert|cult, data=da, type=c("p","a"))

#-----------------------------------------------------------------------------
# usar o confundimento

da$Fert <- factor(da$fert) # esse passo explica o fert=101
levels(da$Fert) # o nível 101 ocorre nas mesmas observações onde
levels(da$cult) # ocorre o nível A, portanto estão confundidos

#-----------------------------------------------------------------------------
# ajuste do modelo considerando os fatores todos de níveis categórios
# poderia ir direto para a análise com polinômio mas esse passo ilustra
# o procedimento para 2 fatores de níveis categóricos

m0 <- aov(y~cult*Fert, data=da)
anova(m0) # compare os graus de liberdade com o número de níveis

# 6-1 gl para cult
# 5-1 gl para Fert
# (6-1)*(5-1) para cult:Fert? não, é (5-1)*(5-1)
# pois existe níveis de fert em apenas 5 níveis de cultivar

#-----------------------------------------------------------------------------
# os contrates entre níveis são os seguintes

contrasts(da$cult) # todo nível de cult contra o nível A
contrasts(da$Fert) # todo nível de Fert contra o nível 50

#-----------------------------------------------------------------------------
# para decompor a soma de quadrados em parte fatorial e parte adicional
# devemos mudar esses contrastes. Um deles deve ser o nível adicional vs todos
# os demais níveis daquele fator. Para isso podemos trocar a ordem dos níveis
# e o tipo de contraste para Helmert

contrasts(C(da$cult, contr=helmert)) # nível F contra os demais, trocar para A
contrasts(C(da$Fert, contr=helmert)) # nível 150 trocar para 101

da$cult <- factor(da$cult, levels=rev(levels(da$cult)))
da$Fert <- factor(da$Fert, levels=c(50,75,100,125,150,101))

levels(da$cult) # A é o último nível
levels(da$Fert) # 101 é o último nível

contrasts(C(da$cult, contr=helmert)) # A vs demais
contrasts(C(da$Fert, contr=helmert)) # 101 vs demias

m0 <- aov(y~cult*Fert, data=da, # passando o tipo de contraste
          contrasts=list(cult=contr.helmert, Fert=contr.helmert))
anova(m0) # mesmo quadro de anova

# agora vamos repartir a soma de quadrados. Note que existe uma hipótese em
# teste quando fazemos esse procedimento. Essa hipótese surge porque estamos
# fazendo um partição ortogonal específica.

summary(m0, expand=FALSE, split=list("cult"=list(fat=1:4, adi=5)))

# de 1:4 são estimativas dos contrastes entre níveis da porção fatorial
# 5 é a estimativa do contraste entre o nível adicional contra os outros
# e essa é a hipótese que a partição ortogonal testa
# mas que nem sempre é a hipótese que temos interesse
# agora é só aplicar os testes para as hipótese de interesse, sejam elas
# comparações múltiplas, contrastes planejados, etc.

#-----------------------------------------------------------------------------
# ajustando modelo polinomio de 2 grau para o efeito de fertilizante

m1 <- lm(y~-1+cult/(fert+I(fert^2)), data=da)
summary(m1)

# usamos essa fórmula para termos interpretação direta das estimativas
# perceba que para o nível A obtivemos NA, lógico, não temos como estimar
# efeito de fertilizante porque não temos níveis para essa cultivar

#-----------------------------------------------------------------------------
# estimando as doses que conferem máxima resposta em cada cultivar
# solução é derivar o modelo quadrático e igualar a zero

D(expression(b0+b1*x+b2*x^2), "x") # igualar a zero e isolar x
fert.max <- -0.5*coef(m1)[6+1:6]/coef(m1)[12+1:6]
fert.max <- fert.max[-6]
fert.max # as doses de fertilizante que conferem o máximo para cada cultivar

#-----------------------------------------------------------------------------
# vamos ajustar com a aov() que omite os NA para então usar a estimable()

m2 <- aov(y~-1+cult/(fert+I(fert^2)), data=da)
summary.lm(m2) # não contém NA

#-----------------------------------------------------------------------------
# agora é só estimar a produção máxima para cada nível de cultivar passando a
# para estimable() a matrix que identifica as funções lineares envolvidas

con <- cbind(diag(5), 0, diag(5)*fert.max, diag(5)*fert.max^2)
con <- rbind(con, 0); con[6,6] <- 1
rownames(con) <- paste(levels(da$cult),
                       c(round(fert.max,1),"100.0"), sep="-")
con

require(gmodels)
e0 <- estimable(m2, cm=con, conf.int=0.95)
e0 # quadro de estimativas da produção nas doses de máximo

#-----------------------------------------------------------------------------
# fazer um gráfico da produção máxima com os intervalos

xlim <- extendrange(e0[,6:7])
barchart(rownames(e0)~Estimate, data=e0, xlim=xlim,
         xlab="Produção", ylab="Cultivar-dose",
         panel=function(x,y,...){
           panel.barchart(x,y,...)
           panel.arrows(e0$Upper.CI, y, e0$Lower.CI, y,
                        code=3, angle=90, length=0.1)
         })

#-----------------------------------------------------------------------------
# fazendo as comparações das produções máximas com a testemunha

contr <- t(sapply(1:5, function(l) con[l,]-con[6,]))
rownames(contr) <- paste(rownames(con)[-6], "vs A")
e1 <- estimable(m2, cm=contr, conf.int=0.95)
e1

# segundo os contrastes, apenas o cultivar D é superior ao A quando recebe
# a dose 158 de fertilizante

#-----------------------------------------------------------------------------
# confeccionando o gráfico final

pred <- expand.grid(cult=levels(da$cult)[-6], fert=seq(50,170,1))
aux <- predict(m2, newdata=pred, interval="confidence")
aux <- stack(as.data.frame(aux))
pred <- cbind(pred, aux)
str(pred)

require(latticeExtra)

#png("f020.png", w=600, h=300)
xyplot(y~fert|cult, subset(da, cult!="A"), col="green4",
       xlim=extendrange(pred$fert), ylim=c(0,525), layout=c(5,1),
       xlab="Nível de fertilizante", ylab="Produtividade da cultura",
       strip=strip.custom(bg="gray90"),
       panel=function(...){
         panel.xyplot(...)
         panel.abline(h=e0[6,c(1,6,7)], lty=c(1,2,2), col="gray50")
         panel.arrows(fert.max[which.packet()], e0[which.packet(),1],
                      fert.max[which.packet()], 0, length=0.1)
       })+
  as.layer(xyplot(values~fert|cult, groups=ind, data=pred,
                  type="l", col=1, lty=c(1,2,2)))
trellis.focus("panel", 1, 1)
panel.text(50, e0[6,1], label="Cultivar A", pos=4, font=3, cex=0.8)
trellis.unfocus()
#dev.off()

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

Tags: , , , , , ,

Comments are closed.