Experimento em delineamento de quadrado latino
##-----------------------------------------------------------------------------
## 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/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 Bcast(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))))
Especificação e estimação do modelo
\[ \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_k\) o efeito da variedade \(k\) 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); 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.8930 0.0142293 *
## colu 4 30481 7620 2.6804 0.0831343 .
## vari 4 137488 34372 12.0905 0.0003585 ***
## Residuals 12 34115 2843
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Comparações múltiplas
##-----------------------------------------------------------------------------
## 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.84489 12 20.65852 9.552243e-11 440.6465 544.5535
## 2 B 440.8 23.84489 12 18.48614 3.489183e-10 388.8465 492.7535
## 3 C 604.8 23.84489 12 25.36393 8.571958e-12 552.8465 656.7535
## 4 D 413.4 23.84489 12 17.33705 7.344429e-10 361.4465 465.3535
## 5 E 401.0 23.84489 12 16.81702 1.044341e-09 349.0465 452.9535## 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.80 33.72 -1.536 0.214930
## C - A == 0 112.20 33.72 3.327 0.015072 *
## D - A == 0 -79.20 33.72 -2.349 0.061340 .
## E - A == 0 -91.60 33.72 -2.716 0.037468 *
## C - B == 0 164.00 33.72 4.863 0.001298 **
## D - B == 0 -27.40 33.72 -0.813 0.480346
## E - B == 0 -39.80 33.72 -1.180 0.325962
## D - C == 0 -191.40 33.72 -5.676 0.000515 ***
## E - C == 0 -203.80 33.72 -6.044 0.000515 ***
## E - D == 0 -12.40 33.72 -0.368 0.719490
## ---
## 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$cld <- l0$mcletters$Letters
p0$let <- sprintf("%.1f %s", p0$estimate, p0$cld)
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)
})
##-----------------------------------------------------------------------------
## Análise com o pacote ExpDes.
require(ExpDes)## Loading required package: ExpDes
##
## Attaching package: 'ExpDes'
##
## The following object is masked from 'package:MASS':
##
## ginv## Funções do pacote ExpDes.
ls("package:ExpDes")## [1] "ccboot" "crd" "duncan" "fat2.ad.crd" "fat2.ad.rbd"
## [6] "fat2.crd" "fat2.rbd" "fat3.ad.crd" "fat3.ad.rbd" "fat3.crd"
## [11] "fat3.rbd" "ginv" "lastC" "latsd" "lsd"
## [16] "lsdb" "order.group" "order.stat.SNK" "rbd" "reg.poly"
## [21] "scottknott" "snk" "split2.crd" "split2.rbd" "tapply.stat"
## [26] "tukey"## Latin squared design.
with(da, latsd(treat=vari, row=linh, column=colu,
resp=prod, quali=TRUE, mcomp="tukey"))## ------------------------------------------------------------------------
## Analysis of Variance Table
## ------------------------------------------------------------------------
## DF SS MS Fc Pr>Fc
## Treatament 4 137488 34372 12.0905 0.000358
## Row 4 55641 13910 4.8930 0.014229
## Column 4 30481 7620 2.6804 0.083134
## Residuals 12 34115 2843
## Total 24 257724
## ------------------------------------------------------------------------
## CV = 11.33 %
##
## ------------------------------------------------------------------------
## Shapiro-Wilk normality test
## p-value: 0.8201574
## According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
## ------------------------------------------------------------------------
##
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a C 604.8
## b A 492.6
## b B 440.8
## b D 413.4
## b E 401
## ------------------------------------------------------------------------