##-----------------------------------------------------------------------------
## 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/pimentel_canadeacucar2.txt"
## Importa os dados a partir do arquivo na web.
da <- read.table(file=url, header=TRUE, sep="\t")
## Para facilitar, truncar os nomes.
names(da) <- substr(names(da), 1, 4)
## Converte linha e coluna de inteiro para fator.
da <- transform(da, linh=factor(linh), colu=factor(colu))
## Mostra a estrutura do objeto.
str(da)
## 'data.frame': 25 obs. of 4 variables:
## $ linh: Factor w/ 5 levels "1","2","3","4",..: 1 1 1 1 1 2 2 2 2 2 ...
## $ colu: Factor w/ 5 levels "1","2","3","4",..: 1 2 3 4 5 1 2 3 4 5 ...
## $ vari: Factor w/ 5 levels "A","B","C","D",..: 4 3 5 2 1 1 5 2 4 3 ...
## $ prod: int 432 724 489 494 515 518 478 384 500 660 ...
## Vendo o quadrado no plano.
cast(da, linh~colu, value="vari")
## linh 1 2 3 4 5
## 1 1 D C E B A
## 2 2 A E B D C
## 3 3 B A C E D
## 4 4 C B D A E
## 5 5 E D A C B
cast(da, linh~colu, value="prod")
## linh 1 2 3 4 5
## 1 1 432 724 489 494 515
## 2 2 518 478 384 500 660
## 3 3 458 524 556 313 438
## 4 4 583 550 297 486 394
## 5 5 331 400 420 501 318
## Tem registros perdidos?
sum(is.na(da$prod))
## [1] 0
##-----------------------------------------------------------------------------
## Análise exploratória.
xyplot(prod~vari, data=da, type=c("p","a"))
## Gráfico do layout.
levelplot(prod~linh+colu,
data=da, aspect="iso")+
layer(with(da,
panel.text(x=linh, y=colu,
label=paste(vari, prod))))
\[ \begin{aligned} Y|\text{fontes de variação} &\sim \text{Normal}(\mu_{ijk}, \sigma^2) \newline \mu_{ijk} &= \mu+\alpha_i+\beta_j+\tau_k \end{aligned} \]
em que \(\alpha_i\) é o efeito da linha \(i\), \(\beta_j\) é o efeito da coluna \(j\) e \(\tau_j\) o efeito da variedade \(j\) sobre a variável resposta \(Y\), \(\mu\) é a média de \(Y\) na ausência do efeito dos termos mencionados, \(\mu_{ijk}\) é a média para as observações de procedencia \(ijk\) e \(\sigma^2\) é a variância das observações ao redor desse valor médio.
##-----------------------------------------------------------------------------
## Estimação.
m0 <- lm(prod~linh+colu+vari, data=da)
##-----------------------------------------------------------------------------
## Diganóstico dos resíduos.
par(mfrow=c(2,2)); plot(m0, which=1:3); layout(1)
##-----------------------------------------------------------------------------
## Quadro de análise de variância.
anova(m0)
## Analysis of Variance Table
##
## Response: prod
## Df Sum Sq Mean Sq F value Pr(>F)
## linh 4 55641 13910 4.89 0.01423 *
## colu 4 30481 7620 2.68 0.08313 .
## vari 4 137488 34372 12.09 0.00036 ***
## Residuals 12 34115 2843
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Comparações múltiplas.
## Médias ajustadas (LSmeans).
p0 <- LSmeans(m0, effect="vari", level=0.95)
## Criando a tabela com as estimativas.
p0 <- cbind(p0$grid, as.data.frame(p0$coef))
p0 <- equallevels(p0, da)
p0
## vari estimate se df t.stat p.value lwr upr
## 1 A 492.6 23.84 12 20.66 9.552e-11 440.6 544.6
## 2 B 440.8 23.84 12 18.49 3.489e-10 388.8 492.8
## 3 C 604.8 23.84 12 25.36 8.572e-12 552.8 656.8
## 4 D 413.4 23.84 12 17.34 7.344e-10 361.4 465.4
## 5 E 401.0 23.84 12 16.82 1.044e-09 349.0 453.0
## Comparações múltiplas, contrastes de Tukey.
## Método de correção de p-valor: single-step.
g0 <- summary(glht(m0, linfct=mcp(vari="Tukey")),
test=adjusted(type="fdr"))
g0
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = prod ~ linh + colu + vari, data = da)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## B - A == 0 -51.8 33.7 -1.54 0.21493
## C - A == 0 112.2 33.7 3.33 0.01507 *
## D - A == 0 -79.2 33.7 -2.35 0.06134 .
## E - A == 0 -91.6 33.7 -2.72 0.03747 *
## C - B == 0 164.0 33.7 4.86 0.00130 **
## D - B == 0 -27.4 33.7 -0.81 0.48035
## E - B == 0 -39.8 33.7 -1.18 0.32596
## D - C == 0 -191.4 33.7 -5.68 0.00051 ***
## E - C == 0 -203.8 33.7 -6.04 0.00051 ***
## E - D == 0 -12.4 33.7 -0.37 0.71949
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- fdr method)
## Resumo compacto com letras.
l0 <- cld(g0, decreasing=TRUE)
l0
## A B C D E
## "b" "bc" "a" "bc" "c"
## Formatando média seguida de letras.
p0$let <- sprintf("%.1f %s", p0$estimate, l0$mcletters$Letters)
p0$let
## [1] "492.6 b" "440.8 bc" "604.8 a" "413.4 bc" "401.0 c"
##-----------------------------------------------------------------------------
## Representação gráfica dos resultados.
xlab <- "Produção"
ylab <- "Variedades de cana-de-açúcar"
sub <- list(
"Médias seguidas de mesma letra indicam diferença nula à 5%.",
font=1, cex=0.8)
p0 <- transform(p0, vari=reorder(vari, estimate))
levels(p0$vari)
## [1] "E" "D" "B" "A" "C"
p0 <- p0[order(p0$vari),]
segplot(vari~lwr+upr, centers=estimate,
data=p0, draw=FALSE,
xlab=xlab, ylab=ylab, sub=sub, letras=p0$let,
panel=function(x, y, z, centers, letras, ...){
panel.segplot(x=x, y=y, z=z, centers=centers, ...)
panel.text(y=as.numeric(z), x=centers,
label=letras, pos=3)
})
