Experimento em DBC com arranjo fatorial duplo desbalanceado
##-----------------------------------------------------------------------------
## Definições da sessão, pacotes a serem usados.
pkg <- c("lattice", "latticeExtra", "doBy", "multcomp", "reshape",
"plyr", "wzRfun")
sapply(pkg, library, character.only=TRUE, logical.return=TRUE)## lattice latticeExtra doBy multcomp reshape plyr
## TRUE TRUE TRUE TRUE TRUE TRUE
## wzRfun
## TRUEsource("lattice_setup.R")
##-----------------------------------------------------------------------------
## Informações sobre as versões dos pacotes.
sessionInfo()## R version 3.1.1 (2014-07-10)
## Platform: i686-pc-linux-gnu (32-bit)
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=pt_BR.UTF-8
## [4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=pt_BR.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=pt_BR.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] methods splines stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] wzRfun_0.3 plyr_1.8.1 reshape_0.8.5 multcomp_1.3-7
## [5] TH.data_1.0-3 mvtnorm_0.9-9997 doBy_4.5-11 MASS_7.3-34
## [9] survival_2.37-7 latticeExtra_0.6-26 RColorBrewer_1.0-5 lattice_0.20-29
## [13] rmarkdown_0.3.3 knitr_1.7
##
## loaded via a namespace (and not attached):
## [1] digest_0.6.4 evaluate_0.5.5 formatR_1.0 grid_3.1.1 htmltools_0.2.6
## [6] Matrix_1.1-4 Rcpp_0.11.3 sandwich_2.3-0 stringr_0.6.2 tools_3.1.1
## [11] yaml_2.1.13 zoo_1.7-10## obs: Para instalar um pacote faça:
## install.packages("nome_do_pacote", dependencies=TRUE)Importação dos dados
##-----------------------------------------------------------------------------
## Lendo arquivos de dados.
## Url de um arquivo com dados.
url <- "http://www.leg.ufpr.br/~walmes/data/soja.txt"
## Importa os dados a partir do arquivo na web.
da <- read.table(file=url, header=TRUE, sep="\t", dec=",")
## Criando cópias de númericas como fator.
da <- transform(da, A=factor(agua), K=factor(potassio))
## Mostra a estrutura do objeto.
str(da)## 'data.frame': 75 obs. of 12 variables:
## $ potassio: int 0 30 60 120 180 0 30 60 120 180 ...
## $ agua : num 37.5 37.5 37.5 37.5 37.5 50 50 50 50 50 ...
## $ bloco : Factor w/ 5 levels "I","II","III",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ rengrao : num 14.6 21.5 24.6 21.9 28.1 ...
## $ pesograo: num 10.7 13.5 15.8 12.8 14.8 ...
## $ kgrao : num 15.1 17.1 19.1 18.1 19.1 ...
## $ pgrao : num 1.18 0.99 0.82 0.85 0.88 1.05 1.08 0.74 1.01 1.01 ...
## $ ts : int 136 159 156 171 190 140 193 200 208 237 ...
## $ nvi : int 22 2 0 2 0 20 6 6 7 10 ...
## $ nv : int 56 62 66 68 82 63 86 94 86 97 ...
## $ A : Factor w/ 3 levels "37.5","50","62.5": 1 1 1 1 1 2 2 2 2 2 ...
## $ K : Factor w/ 5 levels "0","30","60",..: 1 2 3 4 5 1 2 3 4 5 ...## Os dados são de um experimento para verificar o efeito a adubação com
## potássio nos componentes de produção de soja sob diferentes níveis de
## umidade.
##-----------------------------------------------------------------------------
## Análise exploratória.
xyplot(rengrao~potassio, groups=A, data=da,
type=c("p","a"), auto.key=list(columns=nlevels(da$A)))
Especificação e estimação do modelo
\[latex \begin{aligned} Y|\text{fontes de variação} &\sim \text{Normal}(\mu_{ijk}, \sigma^2) \newline \mu_i &= \mu+\beta_i+\alpha_j+\tau_k \end{aligned} \]
em que \(\beta_i\) é o efeito do bloco \(i\), \(\alpha_j\) o efeito do nível de água \(j\) e \(\tau_k\) o efeito do nível de potássio \(k\) sobre a variável resposta \(Y\). \(\mu\) é a média de \(Y\) na ausência do efeito dos termos mencioandos. \(\mu_{ijk}\) é a média para as observações na combinação \(ijk\) e \(\sigma^2\) é a variância das observações ao redor desse valor médio.
##-----------------------------------------------------------------------------
## Por uma questão didática, os fatores que são métricos (quatitativos)
## serão considerados como categóricos (qualitativos) para demonstrar
## como proceder a análise. Em outro script eles serão analisados como
## métricos que é mais recomendado.
## Espeficicação do modelo.
m0 <- lm(rengrao~bloco+A*K, data=da)
## Resíduos.
par(mfrow=c(2,2)); plot(m0); layout(1)
## Remove obs 74 (uma planta que foi comprimetida).
m0 <- lm(rengrao~bloco+A*K, data=da[-74,])
par(mfrow=c(2,2)); plot(m0); layout(1)
## Quadro de anova.
anova(m0)## Analysis of Variance Table
##
## Response: rengrao
## Df Sum Sq Mean Sq F value Pr(>F)
## bloco 4 100.39 25.10 4.8948 0.001896 **
## A 2 446.53 223.26 43.5456 4.617e-12 ***
## K 4 2592.13 648.03 126.3925 < 2.2e-16 ***
## A:K 8 286.00 35.75 6.9726 2.629e-06 ***
## Residuals 55 281.99 5.13
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Não inferir sobre os efeitos de menor ordem se os termos estiverem
## envolvidos em alguma interação (termo de ordem maior). A interação
## mascara os efeitos principais. Não havendo interação pode-se inferir
## sobre os efeitos principais.
##-----------------------------------------------------------------------------
## Partição da soma de quadrados.
## Essa etapa não é obrigatória porém muitos aplicativos reportam a
## particação do quadro de análise de variância como parte dos
## resultados do desdobramente da interação. Aqui será feito à titulo de
## demonstração mas é válido lembrar que ao fazer a remoção da
## observação 74 não se tem mais um experimento com efeito ortogonais,
## logo a partição da soma de quadrados depende da ordem dos termos.
## Partição para efeito de K dentro dos níveis de A.
m1 <- aov(rengrao~bloco+A/K, data=da[-74,])
summary(m1)## Df Sum Sq Mean Sq F value Pr(>F)
## bloco 4 100.4 25.10 4.895 0.0019 **
## A 2 446.5 223.26 43.546 4.62e-12 ***
## A:K 12 2878.1 239.84 46.779 < 2e-16 ***
## Residuals 55 282.0 5.13
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ef <- data.frame(fv=m1$assign, par=names(coef(m1)),
stringsAsFactors=FALSE)
ef3 <- subset(ef, fv==3)
i <- paste0("^A", levels(da$A), ":K")
KinA <- lapply(i, grep, x=ef3$par)
names(KinA) <- levels(da$A)
lapply(KinA, function(i) ef3[i,])## $`37.5`
## fv par
## 8 3 A37.5:K30
## 11 3 A37.5:K60
## 14 3 A37.5:K120
## 17 3 A37.5:K180
##
## $`50`
## fv par
## 9 3 A50:K30
## 12 3 A50:K60
## 15 3 A50:K120
## 18 3 A50:K180
##
## $`62.5`
## fv par
## 10 3 A62.5:K30
## 13 3 A62.5:K60
## 16 3 A62.5:K120
## 19 3 A62.5:K180## Quadro de análise de variância particionado.
summary(m1, split=list("A:K"=KinA))## Df Sum Sq Mean Sq F value Pr(>F)
## bloco 4 100.4 25.1 4.895 0.0019 **
## A 2 446.5 223.3 43.546 4.62e-12 ***
## A:K 12 2878.1 239.8 46.779 < 2e-16 ***
## A:K: 37.5 4 501.8 125.5 24.470 1.15e-11 ***
## A:K: 50 4 714.0 178.5 34.813 1.76e-14 ***
## A:K: 62.5 4 1662.3 415.6 81.055 < 2e-16 ***
## Residuals 55 282.0 5.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##-----------------------------------------------------------------------------
## Outra forma de avaliar a mesma hipótese é pela estatística F do teste
## de Wald. Esse mecanismo de teste de hipótese é mais geral e pode ser
## usado em outras classes de modelo, como glm, modelos de sobrevivência
## e modelos mistos.
i <- paste0("^A", levels(da$A), ":K")
KinA <- lapply(i, grep, x=ef$par)
names(KinA) <- levels(da$A)
lapply(KinA, function(i) ef[i,])## $`37.5`
## fv par
## 8 3 A37.5:K30
## 11 3 A37.5:K60
## 14 3 A37.5:K120
## 17 3 A37.5:K180
##
## $`50`
## fv par
## 9 3 A50:K30
## 12 3 A50:K60
## 15 3 A50:K120
## 18 3 A50:K180
##
## $`62.5`
## fv par
## 10 3 A62.5:K30
## 13 3 A62.5:K60
## 16 3 A62.5:K120
## 19 3 A62.5:K180## Matriz de covariância dos efeitos e estimativas. É um teste baseado
## em aproximação quadrática.
Sigma <- vcov(m1)
b <- coef(m1)
## Requer que o pacote 'aod' esteja instalado.
Fwald <- sapply(KinA, simplify=FALSE,
function(i){
wt <- aod::wald.test(
Sigma=Sigma,
b=b, Terms=i,
df=df.residual(m1))$result$Ftest[c(1,4)]
data.frame(as.list(wt))
})
## str(Fwald)
## Valores de F e correspondente p-valor.
do.call(rbind, Fwald)## Fstat P
## 37.5 24.46954 1.148526e-11
## 50 34.81275 1.765255e-14
## 62.5 81.05541 0.000000e+00Comparações múltiplas
##-----------------------------------------------------------------------------
## Comparações múltiplas para o desdobramento.
LSmeans(m0, effect=c("A","K"))## estimate se df t.stat p.value lwr upr A K
## 1 13.52000 1.012635 55 13.35131 6.766315e-19 11.49063 15.54937 37.5 0
## 2 15.71200 1.012635 55 15.51596 9.604700e-22 13.68263 17.74137 50 0
## 3 16.09000 1.012635 55 15.88924 3.284043e-22 14.06063 18.11937 62.5 0
## 4 20.33400 1.012635 55 20.08029 5.457627e-27 18.30463 22.36337 37.5 30
## 5 22.57000 1.012635 55 22.28839 3.260762e-29 20.54063 24.59937 50 30
## 6 20.98800 1.012635 55 20.72613 1.168374e-27 18.95863 23.01737 62.5 30
## 7 23.92600 1.012635 55 23.62747 1.775150e-30 21.89663 25.95537 37.5 60
## 8 28.69200 1.012635 55 28.33401 1.719045e-34 26.66263 30.72137 50 60
## 9 29.82800 1.012635 55 29.45584 2.312889e-35 27.79863 31.85737 62.5 60
## 10 25.30800 1.012635 55 24.99223 1.047716e-31 23.27863 27.33737 37.5 120
## 11 29.78600 1.012635 55 29.41436 2.487930e-35 27.75663 31.81537 50 120
## 12 36.33929 1.140218 55 31.87047 3.842691e-37 34.05424 38.62433 62.5 120
## 13 25.39000 1.012635 55 25.07321 8.894077e-32 23.36063 27.41937 37.5 180
## 14 28.76000 1.012635 55 28.40116 1.521591e-34 26.73063 30.78937 50 180
## 15 37.14600 1.012635 55 36.68253 2.368324e-40 35.11663 39.17537 62.5 180X <- LSmatrix(m0, effect=c("A","K"))
grid <- equallevels(attr(X, "grid"), da)
## Divide a matriz de acordo com os níveis de A.
L <- by(X, INDICES=grid$A, FUN=as.matrix)
L <- lapply(L, "rownames<-", levels(da$K))
t(L[[1]])## 0 30 60 120 180
## (Intercept) 1.0 1.0 1.0 1.0 1.0
## blocoII 0.2 0.2 0.2 0.2 0.2
## blocoIII 0.2 0.2 0.2 0.2 0.2
## blocoIV 0.2 0.2 0.2 0.2 0.2
## blocoV 0.2 0.2 0.2 0.2 0.2
## A50 0.0 0.0 0.0 0.0 0.0
## A62.5 0.0 0.0 0.0 0.0 0.0
## K30 0.0 1.0 0.0 0.0 0.0
## K60 0.0 0.0 1.0 0.0 0.0
## K120 0.0 0.0 0.0 1.0 0.0
## K180 0.0 0.0 0.0 0.0 1.0
## A50:K30 0.0 0.0 0.0 0.0 0.0
## A62.5:K30 0.0 0.0 0.0 0.0 0.0
## A50:K60 0.0 0.0 0.0 0.0 0.0
## A62.5:K60 0.0 0.0 0.0 0.0 0.0
## A50:K120 0.0 0.0 0.0 0.0 0.0
## A62.5:K120 0.0 0.0 0.0 0.0 0.0
## A50:K180 0.0 0.0 0.0 0.0 0.0
## A62.5:K180 0.0 0.0 0.0 0.0 0.0## Faz comparações múltiplas de médias entre os níveis de K dentro de
## cada nível de A.
grid <- lapply(L, apmc, model=m0, focus="K", test="fdr")
grid <- ldply(grid); names(grid)[1] <- "A"
grid <- equallevels(grid, da)
grid## A K estimate lwr upr cld
## 1 37.5 0 13.52000 11.49063 15.54937 c
## 2 37.5 30 20.33400 18.30463 22.36337 b
## 3 37.5 60 23.92600 21.89663 25.95537 a
## 4 37.5 120 25.30800 23.27863 27.33737 a
## 5 37.5 180 25.39000 23.36063 27.41937 a
## 6 50 0 15.71200 13.68263 17.74137 c
## 7 50 30 22.57000 20.54063 24.59937 b
## 8 50 60 28.69200 26.66263 30.72137 a
## 9 50 120 29.78600 27.75663 31.81537 a
## 10 50 180 28.76000 26.73063 30.78937 a
## 11 62.5 0 16.09000 14.06063 18.11937 d
## 12 62.5 30 20.98800 18.95863 23.01737 c
## 13 62.5 60 29.82800 27.79863 31.85737 b
## 14 62.5 120 36.33929 34.05424 38.62433 a
## 15 62.5 180 37.14600 35.11663 39.17537 a##-----------------------------------------------------------------------------
## Representação gráfica dos resultados.
ylab <- expression("Rendimento de grãos"~(g~vaso^{-1}))
xlab <- expression("Nível de potássio"~(mg~dm^{-3}))
sub <- list(
"Médias seguidas de mesma letra indicam diferença nula à 5%.",
font=1, cex=0.8)
## Gráfico final.
segplot(K~lwr+upr|A, centers=estimate, horizontal=FALSE,
data=grid, draw=FALSE, layout=c(NA,1),
xlab=xlab, ylab=ylab, sub=sub,
strip=strip.custom(
factor.levels=sprintf("Umidade: %s %s", levels(da$A), "%")),
letras=sprintf("%0.2f %s", grid$estimate, grid$cld),
panel=function(x, y, z, subscripts, centers, letras, ...){
panel.segplot(x=x, y=y, z=z, centers=centers,
subscripts=subscripts, ...)
panel.text(x=as.numeric(z)[subscripts],
y=centers[subscripts],
label=letras[subscripts],
srt=90, adj=c(0.5,-0.5))
})