##-----------------------------------------------------------------------------
## 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
## TRUE
source("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] splines stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] wzRfun_0.1 plyr_1.8.1 reshape_0.8.5 multcomp_1.3-3
## [5] TH.data_1.0-3 mvtnorm_0.9-99992 doBy_4.5-10 MASS_7.3-33
## [9] survival_2.37-7 latticeExtra_0.6-26 RColorBrewer_1.0-5 lattice_0.20-29
## [13] knitr_1.5
##
## loaded via a namespace (and not attached):
## [1] evaluate_0.5.5 formatR_0.10 grid_3.1.1 lme4_1.1-6
## [5] Matrix_1.1-4 methods_3.1.1 minqa_1.2.3 nlme_3.1-117
## [9] Rcpp_0.11.1 RcppEigen_0.3.2.1.2 sandwich_2.3-0 stringr_0.6.2
## [13] tools_3.1.1 zoo_1.7-11
## obs: Para instalar um pacote faça:
## install.packages("nome_do_pacote", dependencies=TRUE)
##-----------------------------------------------------------------------------
## 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)))
\[ \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 25 4.89 0.0019 **
## A 2 447 223 43.55 4.6e-12 ***
## K 4 2592 648 126.39 < 2e-16 ***
## A:K 8 286 36 6.97 2.6e-06 ***
## Residuals 55 282 5
## ---
## 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 25.1 4.89 0.0019 **
## A 2 447 223.3 43.55 4.6e-12 ***
## A:K 12 2878 239.8 46.78 < 2e-16 ***
## Residuals 55 282 5.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ef <- 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 25 4.89 0.0019 **
## A 2 447 223 43.55 4.6e-12 ***
## A:K 12 2878 240 46.78 < 2e-16 ***
## A:K: 37.5 4 502 125 24.47 1.1e-11 ***
## A:K: 50 4 714 178 34.81 1.8e-14 ***
## A:K: 62.5 4 1662 416 81.06 < 2e-16 ***
## Residuals 55 282 5
## ---
## 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.47 1.149e-11
## 50 34.81 1.765e-14
## 62.5 81.06 0.000e+00
##-----------------------------------------------------------------------------
## 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.52 1.013 55 13.35 6.766e-19 11.49 15.55 37.5 0
## 2 15.71 1.013 55 15.52 9.605e-22 13.68 17.74 50 0
## 3 16.09 1.013 55 15.89 3.284e-22 14.06 18.12 62.5 0
## 4 20.33 1.013 55 20.08 5.458e-27 18.30 22.36 37.5 30
## 5 22.57 1.013 55 22.29 3.261e-29 20.54 24.60 50 30
## 6 20.99 1.013 55 20.73 1.168e-27 18.96 23.02 62.5 30
## 7 23.93 1.013 55 23.63 1.775e-30 21.90 25.96 37.5 60
## 8 28.69 1.013 55 28.33 1.719e-34 26.66 30.72 50 60
## 9 29.83 1.013 55 29.46 2.313e-35 27.80 31.86 62.5 60
## 10 25.31 1.013 55 24.99 1.048e-31 23.28 27.34 37.5 120
## 11 29.79 1.013 55 29.41 2.488e-35 27.76 31.82 50 120
## 12 36.34 1.140 55 31.87 3.843e-37 34.05 38.62 62.5 120
## 13 25.39 1.013 55 25.07 8.894e-32 23.36 27.42 37.5 180
## 14 28.76 1.013 55 28.40 1.522e-34 26.73 30.79 50 180
## 15 37.15 1.013 55 36.68 2.368e-40 35.12 39.18 62.5 180
X <- 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.52 11.49 15.55 c
## 2 37.5 30 20.33 18.30 22.36 b
## 3 37.5 60 23.93 21.90 25.96 a
## 4 37.5 120 25.31 23.28 27.34 a
## 5 37.5 180 25.39 23.36 27.42 a
## 6 50 0 15.71 13.68 17.74 c
## 7 50 30 22.57 20.54 24.60 b
## 8 50 60 28.69 26.66 30.72 a
## 9 50 120 29.79 27.76 31.82 a
## 10 50 180 28.76 26.73 30.79 a
## 11 62.5 0 16.09 14.06 18.12 d
## 12 62.5 30 20.99 18.96 23.02 c
## 13 62.5 60 29.83 27.80 31.86 b
## 14 62.5 120 36.34 34.05 38.62 a
## 15 62.5 180 37.15 35.12 39.18 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))
})
