## Objetivos:
## a) Permitir a modelagem da ANOVA com número variado de grupos e ##repetições
## b) Exibir resultados numéricos e gráficos
## c) Comparar os resultados por 2 métodos diferentes de ##discriminação dos grupos
## Pacotes necessários:
## - plotrix
## - TukeyC
## - ScottKnott
##===============================================================================
## Dados para modelagem
alfa <- .08 # Probabilidade do erro tipo I na inferência da ANOVA
r <- 7 # Número de repetições de cada grupo (tratamento)
M <- c(60,
60,
# 60,
# 63,
65) # Médias dos grupos
SD <- 6 # Desvio padrão dos grupos
##===============================================================================
## Opção de cores
ocol <- rainbow(length(M))
## Funções
plotar_info_Ft <- function()
{
polygon(x=c(Ft,
seq(Ft,
20,
length=1e3)),
y=c(0,
df(seq(Ft,
20,
length=1e3),
length(M) - 1,
length(M) * (r-1))),
col='darkred')
segments(x0=Ft,
x1=Ft,
y0=0,
y1=1,
lty=3,
col='darkred')
text(x=Ft,
y=0.4,
pos=4,
offset=0,
bquote(alpha == .(alfa)),
cex=1.4,
col='darkred')
points(x=Ft,
y=-.06,
pch=17,
cex=1.5,
col='darkred')
}
plotar_info_Fc <- function()
{
polygon(x=c(Fc,
seq(Fc,
20,
length=1e3)),
y=c(0,
df(seq(Fc,
20,
length=1e3),
length(M) - 1,
length(M) * (r-1))),
col=gray(.5))
if (Fc <= 15)
{
segments(x0=Fc,
x1=Fc,
y0=0,
y1=.15,
lty=3)
}
if (Fc <= 15)
{
text(x=Fc,
y=0.12,
pos=4,
offset=0,
bquote(P == .(tab[1, 5])),
cex=1.4)
} else {
text(x=15,
y=0.12,
pos=4,
offset=0,
bquote(P == .(tab[1, 5])),
cex=1.4)
}
}
## Início da modelagem
## Gerando grupos com r repetições com distribuição normal em um data.frame
dad <- data.frame(tr=gl(length(M),
r,
labels=LETTERS[1:length(M)]),
y=rnorm(length(M) * r,
rep(M,
each=r),
sd=SD))
## Análise de variância
av <- aov(y ~ tr,
data=dad)
sav <- summary(av)
tab <- as.data.frame(round(sav[[1]][, 1:5],
2))
cv <- 100 * sqrt(tab[2, 3]) / mean(dad$y)
colnames(tab)=c('Gl',
'SQ',
'QM',
'Fcal',
'P(>Fcal)')
rownames(tab)=c('Grupos',
'Resíduo')
## Calculando médias e desvios dos grupos
medias <- with(dad,
tapply(y,
tr,
mean))
desvios <- with(dad,
tapply(y,
tr,
sd))
## O ordenamento irá facilitar a padronização de cores nos plots
ord <- order(medias)
medias <- sort(medias)
desvios <- desvios[ord]
## Valores F
Ft <- qf(alfa,
length(M) - 1,
length(M) * (r-1),
lower.tail=FALSE)
Fc <- tab[1,
4]
## Verificar dispositivo gráfico
if (is.null(dev.list()))
X11(w=8,
h=8)
## Layout do dispositivo gráfico
layout(matrix(1:6,
ncol=2,
byrow=T),
heights=c(3,
4,
4))
## Plot da tabela
plot(0.5,
type = 'n',
axes=F,
xlab='',
ylab='')
## Tabela da ANOVA no plot
suppressMessages(require(plotrix))
addtable2plot(.5,
.4,
tab,
display.rownames=TRUE,
bg=gray(.8),
cex=1.8)
mtext('ANOVA')
mtext(paste('cv = ',
round(cv,
2),
'%',
sep=''),
side=1)
## Distribuição F e decisão da inferência
curve(df(x,
length(M) - 1,
length(M) * (r-1)),
from=0,
to=20,
n=1e3,
xlab='F',
ylab='f(F)',
ylim=c(-.09,
1),
col=gray(.4))
text(x=Ft,
y=0.8,
lab=c('RAHo',
'RRHo'),
pos=c(2,
4),
offset=.5,
cex=1.4,
col=c('darkgreen',
'darkred'))
## Decisão da ANOVA
if(Fc < Ft){
points(Fc,
y=-.08,
pch=17,
cex=2,
col=gray(.5))
mtext(text = 'Não rejeita Ho',
cex=1.5,
col='darkgreen')
} else {
mtext('Rejeita Ho',
cex=1.5,
col='darkred')
if(Fc <= 15){
points(Fc,
y=-.08,
pch=17,
cex=2,
col=gray(.5))
} else {
text(15,
-.08,
pos=4,
offset=0,
expression(Delta > 15),
cex=1.3,
col=gray(.5))
}
}
## Plotando áreas sob a distribuição F
if(Fc >= Ft){
plotar_info_Ft()
plotar_info_Fc()
} else {
plotar_info_Fc()
plotar_info_Ft()
}
legend('topright',
pch=17,
col=c('darkred',
gray(.5)),
legend=c('Ftab',
'Fcal'),
cex=1.5)
means_dec <- with(dad,
reorder(tr,
-y,
mean))
boxplot(y ~ means_dec,
data=dad,
xlab='Grupos',
ylab='Y',
col=ocol)
## Plotando as densidades dos grupos
plot(1,
type='n',
ylab='Densidade',
xlab='Médias dos grupos',
xlim=c( .3 * min(medias),
1.7 * max(medias)),
ylim=c(-.005,
max(c(dnorm(medias[1:length(M)],
medias[1:length(M)],
desvios[1:length(M)])))))
for(i in 1:length(M)){
curve(dnorm(x,
medias[i],
desvios[i]),
.2 * min(medias[i]),
1.8 * max(medias[i]),
n=1e3,
add=TRUE,
col=rev(ocol)[i],
lw=1.5)
}
for(i in 1:length(M)){
segments(medias[i],
max(dnorm(medias[i],
medias[i],
desvios[i])),
medias[i],
0,
col=rev(ocol)[i],
lty=3)
}
for(i in 1:length(M)){
text(medias[i],
-.005,
substitute(mu[i],
list(i = LETTERS[ord][i])),
col=rev(ocol)[i])
}
legend('topright',
rev(LETTERS[ord]),
lty=rep(1,
length(M)),
lw=2,
col=ocol)
## Tukey
suppressMessages(require(TukeyC))
tk <- TukeyC(av,
sig.level=alfa)
plot(tk,
xlab='Grupos',
ylab='Y',
col=ocol,
title=paste('Tukey',
' (',
100*alfa,
'%',
')',
sep=''),
cex=1.2)
## ScottKnott
suppressMessages(require(ScottKnott))
sk <- SK(av,
sig.level=alfa)
plot(sk,
xlab='Grupos',
ylab='Y',
title=paste('Scott & Knott',
' (',
100*alfa,
'%',
')',
sep=''),
cex=1.2,
col=gray(seq(0,
.8,
length=length(unique(sk$groups)))))
# Coeficiente de Variação
# Arguments:
# av : an lm, aov or aovlist object
# round : rounds the values of cv to the specified number of decimal places (default 2)
cv <- function(av, round=2)
{
if(is.null(av))
stop('Please, check the parameter!')
if(inherits(av,
'lm')) {
qmee <- with(av,
sum(residuals^2) / df.residual)
m <- mean(av$fitted.values)
cv <- 100 * sqrt(qmee) / m
names(cv) <- 'cv'
return(round(cv,
round))
}
if(inherits(av,
'aovlist')) {
sm.av <- summary(av)
qmee <- sapply(sm.av,
function(i){
i[[1]]["Residuals", 3]
})
qmee <- qmee[!is.na(qmee)]
m <- av$'(Intercept)'$coefficients[[1]]
cv <- 100 * sqrt(qmee) / m
names(cv) <- paste('cv(',
letters[1:length(qmee)],
')',
sep='')
return(round(cv,
round))
}
}
