##===============================================================================
## Título: Análise didática de experimento em parcela subdividida
## Versão: v4
## Objetivos:
##===============================================================================
## a) Mostrar diferentes alternativas de análise dos experimentos em parcelas
## subdivididas
## b) Testes de comparação de médias e análise de contrastes
##===============================================================================
# Leitura dos dados
# https://sites.google.com/site/agrobiostat/files
#=======================================
dados<-read.table('./dados/eps.txt',
header=T)
head(dados, 13)
str(dados)
is.factor(dados$bl)
is.factor(dados$adubo)
is.factor(dados$dose)
is.numeric(dados$y)
summary(dados)
#=======================================
#============== ANOVA ==================
#=======================================
av <-aov(y ~ bl + adubo * dose + Error(bl/adubo),
data=dados)
summary(av) # Notar que a interação entre adubo e dose
# é significativa
library(TukeyC)
cv(av)
#=======================================
#==== Ftab. GL adubo 3; GL resíduo 9 ===
#=======================================
(Ftab1 <- qf(5 / 100,
3,
9,
lower.tail=FALSE))
#=======================================
#Ftab. GL adubo:dose 6 ; GL resíduo 24
#=======================================
(Ftab2 <- qf(5 / 100,
6,
24,
lower.tail=FALSE))
#=======================================
#===== Ftab.GL dose 2; GL resíduo 24 ===
#=======================================
(Ftab3 <- qf(5 / 100,
2,
24,
lower.tail=FALSE))
#=======================================
#======== GRÁFICO ADUBO ================
#=======================================
curve(expr=df(x, 3, 9),
main="ADUBO - F calculado F(3, 9) gl",
xlab="F",
ylab="Densidade de probabilidades, f(F)",
xlim=c(0, 10),
n=1e3)
abline(v=Ftab1,
col="red")
abline(h=0,
lty=2)
xf <- seq(Ftab1,
10,
0.01)
ydf <- df(xf,
3,
9)
polygon(c(Ftab1, xf),
c(0, ydf),
col="red")
text(x=6,
y=.1,
paste("pf(>=",
round(Ftab1, 2),
") =",
round(pf(Ftab1,
3,
9,
lower.tail=F),
4)),
cex=0.8,
col="black")
#=======================================
#============ GRÁFICO DOSE =============
#=======================================
curve(expr=df(x, 2, 24),
main="DOSE F calculado F(2, 24) gl",
xlab="Valor F",
ylab="Densidade Probabilística F(f)",
xlim=c(0, 10),
n=1e3)
abline(v=Ftab3,
col="red")
abline(h=0,
lty=2)
xf <- seq(Ftab3,
10,
0.001)
ydf <- df(xf,
2,
24)
polygon(c(Ftab3, xf),
c(0, ydf),
col="red")
text(x=5,
y=.1,
paste("pf(>=",
round(Ftab3, 2),
") =",
round(pf(Ftab3,
2,
24,
lower.tail=F),
4)),
cex=0.8,
col="black")
#=======================================
#======== GRÁFICO ADUBO:DOSE ===========
#=======================================
curve(expr=df(x, 6, 24),
main="ADUBO:DOSE - F calculado F(6, 24) gl",
xlab="F",
ylab="Densidade Probabilística F(f)",
xlim=c(0, 10),
n=1e3)
abline(v=Ftab2,
col="red")
abline(h=0,
lty=2)
xf <- seq(Ftab2,
10,
0.001)
ydf <- df(xf,
6,
24)
polygon(c(Ftab2, xf),
c(0, ydf),
col="red")
text(x=4,
y=.1,
paste("pf(>=",
round(Ftab2, 2),
") =",
round(pf(Ftab2,
6,
24,
lower.tail=F),
4)),
cex=0.8,
col="black")
#=======================================
#============ INTERAÇÃO ==============
#=======================================
par(mfrow = c(2, 1))
with(dados,
interaction.plot(adubo,
dose,
y,
col = 1:4,
lwd = 2,
ylab = "Kg/ha"))
with(dados,
interaction.plot(dose,
adubo,
y,
col = 1:4,
lwd = 2,
ylab = "Kg/ha"))
#=======================================
#============== TUKEY ==================
#=======================================
# Testando efeito principal das parcelas (adubo)
# (finalidade didática pois a interação é significativa: necessário desdobrar!)
ta <- TukeyC(av,
which='adubo',
error='bl:adubo',
sig.level=.05)
summary(ta)
# Testando efeito principal da subparcela (dose)
# (finalidade didática pois a interação é significativa: necessário desdobrar!)
td <- TukeyC(av,
which='dose',
error='Within',
sig.level=0.05)
summary(td)
#=======================================================
# adubo/dose: estudando dose dentro dos níveis de adubo
#=======================================================
#===: nível 1
tad1 <- TukeyC.nest(av,
which='adubo:dose',
error='Within',
sig.level=0.05,
fl1=1) # 1: SALITRE C
summary(tad1)
#===: nível 2
tad2 <- TukeyC.nest(av,
which='adubo:dose',
error='Within',
sig.level=0.05,
fl1=2) # 2: SULFATO
summary(tad2)
#===: nível 3
tad3 <- TukeyC.nest(av,
which='adubo:dose',
error='Within',
sig.level=0.05,
fl1=3) # 3: URÉIA
summary(tad3)
#===: nível 4
tad4 <- TukeyC.nest(av,
which='adubo:dose',
error='Within',
sig.level=0.05,
fl1=4) # 4: CALNITRO
summary(tad4)
#=======================================
#============ CONTRASTES ===============
#=======================================
# É bem trabalhoso no caso dos experimentos em parcelas subdividadas
# devido a mais de uma estimativa do resíduo!
# Matriz de contrastes dos adubos
#A1 A2 A3 A4
mca <- rbind('A2 vs. (A1, A3, A4)' = c(-1, 3, -1, -1),
'A3 vs. (A1, A4)' = c(-1, 0, 2, -1),
'A1 vs. A4' = c( 1, 0, 0, -1))
# Matriz de contrastes das doses
#D1 D2 D3
mcd <- rbind('D3 vs. (D1, D2)' = c(-1, -1, 2),
'D1 vs. D2' = c( 1, -1, 0))
library(gmodels) # Necessário para a função: make.contrasts
## O modelo usado para a ANOVA
## não fornece todos os resultados corretos:
# av2 <-aov(y ~ bl + adubo * dose + Error(bl/adubo),
# data=dados,
# contrasts=list('adubo'=make.contrasts(mca),
# 'dose'=make.contrasts(mcd)))
##-----------------------------------------------------------------------------------
## Alternativa 1:
## Usar o modelo fatorial
## e corrigir manualmente as estimativas dos contrastes associados a subparcela
## (um pouco trabalhoso)
##-----------------------------------------------------------------------------------
## Necessário aninhar os efeitos!
av2 <-aov(y ~ bl + adubo + adubo/dose,
data=dados,
contrasts=list('adubo'=make.contrasts(mca),
'dose'=make.contrasts(mcd)))
effects(av2)[1:15]
sm <- summary(av2,
split=list('adubo'=list('A2 vs. (A1, A3, A4)' = 1,
'A3 vs. (A1, A4)' = 2,
'A1 vs. A4' = 3),
'adubo:dose'=list('A1/dose*' = c(1, 5),
'A1/d3 vs. (d1, d2)'= 1,
'A1/d1 vs. d2' = 5,
'A2/dose*' = c(2, 6),
'A2/d3 vs. (d1, d2)'= 2,
'A2/d1 vs. d2' = 6,
'A3/dose*' = c(3, 7),
'A3/d3 vs. (d1, d2)'= 3,
'A3/d1 vs. d2' = 7,
'A4/dose*' = c(4, 8),
'A4/d3 vs. (d1, d2)'= 4,
'A4/d1 vs. d2' = 8
)))
sm # Como o modelo informado não é correto,
# foi usado apenas para estimar os contrastes, observar que n erro b (subparcela):
# O qmd (141717) está errada, o correto é 74480!
# os graus de liberdade (33) estão errados, o correto é 24!
# Corrigindo os testes:
# Efeitos de dose em cada nível de adubo
round(pf(sm[[1]][c(7, 10, 13, 16), 3] / 74480, # 74480 erro associado a subparcela
2,
24, # Graus de liberdade associados a subparcela
lower.tail=FALSE),
3)
# Contrastes de dose em cada nível de adubo
round(pf(sm[[1]][c(8, 9, 11, 12, 14, 15, 17, 18), 3] / 74480,
1,
24,
lower.tail=FALSE),
3)
##-----------------------------------------------------------------------------------
## Alternativa 2:
## Com a participação do prof. Ivan Bezerra Allaman
## foi possível fazer a decomposição da SQD desse experimento corretamente!
##-----------------------------------------------------------------------------------
av3 <-aov(y ~ bl + adubo + adubo/dose + Error(bl/adubo),
data=dados,
contrasts=list('adubo'=make.contrasts(mca),
'dose'=make.contrasts(mcd)))
effects(av3) # Mensagem de erro -> método não aplicado a classe!
coef(av3) # Dá para pegar a mesma informação com a função coef!
summary(av3)
summary(av3,
split=list('adubo'=list('A2 vs. (A1, A3, A4)' = 1,
'A3 vs. (A1, A4)' = 2,
'A1 vs. A4' = 3),
'adubo:dose'=list('A1/dose*' = c(1, 5),
'A1/d3 vs. (d1, d2)'= 1,
'A1/d1 vs. d2' = 5,
'A2/dose*' = c(2, 6),
'A2/d3 vs. (d1, d2)'= 2,
'A2/d1 vs. d2' = 6,
'A3/dose*' = c(3, 7),
'A3/d3 vs. (d1, d2)'= 3,
'A3/d1 vs. d2' = 7,
'A4/dose*' = c(4, 8),
'A4/d3 vs. (d1, d2)'= 4,
'A4/d1 vs. d2' = 8
)))
