Página Inicial > delineamento > Experimento fatorial duplo com um tratamento adicional

Experimento fatorial duplo com um tratamento adicional

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

#-----------------------------------------------------------------------------
About these ads
  1. Pablo Fernando
    03/05/2012 às 09:23

    Gostaria que montasse uma tabela com o exemplo. E demonstrasse passo a passo como fazia a ANAVA nesse caso. Iria ajudar muito. Grato

    • Walmes Zeviani
      07/05/2012 às 17:26

      Você não reproduziu o exemplo do post? Nele tem os passos para obter a análise de variância.

      • Pablo Fernando
        07/05/2012 às 23:21

        Eu sou estudante de graduação e não sou expert em estatistica e por isso não entendi essa codificação que foi usada. Precisa de algum software específico pra visualizar? Ou será que foi o meu navegador que tá desconfigurado? Tem como da um exemplo com uns valores em uma tabela, formato pdf ou doc. e passar o passo de como faz a anava para esse tipo de experimento. Desde já agradeço a atenção…Valeu

  2. Pablo Fernando
    07/05/2012 às 23:26

    De preferencia se pudesse detalhar melhor o exemplo do experimento com herbicidas (4×6)+2.

  3. Ricardo
    30/09/2012 às 18:46

    Olá Walmes,
    achei o seu post quando tentava analisar os dados de um delineamento fatorial. Estou tratando de um problema de desbalanceamento com células ausentes. Estou seguindo essa referência http://users.monash.edu.au/~murray/BDAR/
    capitulo 12.

    Pelo que li até agora recomenda se fazer uma anova simples, como se cada célula (combinação de fatores) fosse um tratamento, e então fazer os contrastes para descobrir os efeitos principais e interações. O capitulo do livro não explicou uma maneira sistematica de gerar os contrastes, e nem se ha restrições quanto a que interações eu posso testar. Gostaria de saber se você teria uma leitura para recomendar explicando como lidar com esse tipo de delineamento.

    Estou fazendo assim, seguindo o livro

    contrasts(fdata$ALL) &lt;- cbind(c(rep(1,46),rep(-1,45), rep(0,32)),c(rep(-1,46),rep(0,45), rep(1,32)))
    AnovaM(aov(resp ~ ALL, fdata), 
    split = list(ALL = list(&quot;PIE&quot; = 1:2, &quot; PIE1 vs PIE2&quot;=1, &quot; PIE1 vs PIE3&quot;=2)))
    

    dessa forma eu verifico o efeito principal de um fator, mas essa maneira de gerar os contrastes é lenta, e estou com dificuldades de extrapolar para a maneira como você fez, se você pudesse recomendar uma leitura ficaria muito grato.

    • Walmes Zeviani
      01/10/2012 às 10:38

      Ricardo, não tenho uma leitura pontual para lhe passar. Eu conheço esse procedimento de analisar o fatorial como se fosse apenas um único fator (tratando combinações como níveis individuais) e depois montando ons contrastes de interesse. Confesso que considero isso complicado porque você tem que especificar esses contrastes. Do jeito que eu fiz, parte do trabalho é evitado, os testes F são apresentados, você só precisa cuidar do confundimento gerado. Eu não considero complicado um experimento com perdas de caselas e acho suficiente usar uma anova de SQ sequencial (anova(modelo)). A complicação que existe é para testar efeitos principais na presença de interações significativas e isso não faz nenhum sentido, porém é praticado. Pra quê essa complicação desnecessária, eu me pergunto? Não sei se é o seu caso.

      • Ricardo
        01/10/2012 às 17:06

        Obrigado por responder Walmes,

        então, eu não tenho experiência com análises de variância, e as três referências que eu li ate agora sugeriram o procedimento da anova simples com contrastes. Graficamente eu consigo verificar que há interação, mas como são interações de 5 níveis eu estou com dificuldades em especificar os contrastes, vou dar uma olhada no procedimento sequêncial, se achar algo tento postar para continuar a discussão.

  1. No trackbacks yet.

Deixe um comentário

Preencha os seus dados abaixo ou clique em um ícone para log in:

Logotipo do WordPress.com

Você está comentando utilizando sua conta WordPress.com. Sair / Alterar )

Imagem do Twitter

Você está comentando utilizando sua conta Twitter. Sair / Alterar )

Foto do Facebook

Você está comentando utilizando sua conta Facebook. Sair / Alterar )

Foto do Google+

Você está comentando utilizando sua conta Google+. Sair / Alterar )

Conectando a %s

Seguir

Obtenha todo post novo entregue na sua caixa de entrada.

Junte-se a 52 outros seguidores

%d blogueiros gostam disto: