Experimento em delineamento de blocos casualizados
##-----------------------------------------------------------------------------
## 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_batatinha.txt"
## Importa os dados a partir do arquivo na web.
da <- read.table(file=url, header=TRUE, sep="\t")
## Para facilitar, truncar nomes para 4 caracteres.
names(da) <- substr(names(da), start=1, stop=4)
## Mostra a estrutura do objeto.
str(da)## 'data.frame': 32 obs. of 3 variables:
## $ bloc: Factor w/ 4 levels "I","II","III",..: 1 1 1 1 1 1 1 1 2 2 ...
## $ vari: Factor w/ 8 levels "B 116-51","B 1-52",..: 7 6 8 5 3 2 1 4 7 6 ...
## $ prod: num 9.2 21.1 22.6 15.4 12.7 20 23.1 18 13.4 27 ...## Tabela de frequência.
xtabs(~vari+bloc, data=da)## bloc
## vari I II III IV
## B 116-51 1 1 1 1
## B 1-52 1 1 1 1
## B 25-50 E 1 1 1 1
## B 72-53 A 1 1 1 1
## Buena Vista 1 1 1 1
## Huinkul 1 1 1 1
## Kennebec 1 1 1 1
## S. Rafalela 1 1 1 1##-----------------------------------------------------------------------------
## Análise exploratória.
xyplot(prod~vari, groups=bloc, data=da, type=c("p","a"))
Especificação e estimação do modelo
\[ \begin{aligned} Y|\text{fontes de variação} &\sim \text{Normal}(\mu_{ij}, \sigma^2) \newline \mu_{ij} &= \mu+\alpha_i+\tau_j \end{aligned} \]
em que \(\alpha_i\) é o efeito do bloco \(i\) 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 de bloco e variedade, \(\mu_{ij}\) é a média para as observações em cada combinação de níveis de bloco e variedade e \(\sigma^2\) é a variância das observações ao redor desse valor médio.
##-----------------------------------------------------------------------------
## Estimação.
m0 <- lm(prod~bloc+vari, data=da)
## Estimativas e erros padrões.
summary(m0)##
## Call:
## lm(formula = prod ~ bloc + vari, data = da)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.2750 -1.6438 0.0625 1.0563 5.6500
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.550 1.714 11.990 7.40e-11 ***
## blocII 3.500 1.462 2.395 0.02605 *
## blocIII 2.275 1.462 1.556 0.13455
## blocIV 2.025 1.462 1.385 0.18047
## variB 1-52 -0.225 2.067 -0.109 0.91436
## variB 25-50 E -6.000 2.067 -2.903 0.00851 **
## variB 72-53 A 0.300 2.067 0.145 0.88599
## variBuena Vista -10.075 2.067 -4.874 8.08e-05 ***
## variHuinkul 2.550 2.067 1.234 0.23098
## variKennebec -11.800 2.067 -5.708 1.15e-05 ***
## variS. Rafalela 2.950 2.067 1.427 0.16825
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.923 on 21 degrees of freedom
## Multiple R-squared: 0.8439, Adjusted R-squared: 0.7696
## F-statistic: 11.35 on 10 and 21 DF, p-value: 2.153e-06## Restrição paramétrica de zerar o primeiro nível.
contrasts(da$bloc)## II III IV
## I 0 0 0
## II 1 0 0
## III 0 1 0
## IV 0 0 1contrasts(da$vari)## B 1-52 B 25-50 E B 72-53 A Buena Vista Huinkul Kennebec S. Rafalela
## B 116-51 0 0 0 0 0 0 0
## B 1-52 1 0 0 0 0 0 0
## B 25-50 E 0 1 0 0 0 0 0
## B 72-53 A 0 0 1 0 0 0 0
## Buena Vista 0 0 0 1 0 0 0
## Huinkul 0 0 0 0 1 0 0
## Kennebec 0 0 0 0 0 1 0
## S. Rafalela 0 0 0 0 0 0 1##-----------------------------------------------------------------------------
## Existe efeito das fontes de variação? Ou seja, todos os \alpha_i
## podem ser considerados iguais à zero? Todos os \tau_j podem ser
## considerados à zero.
## Quadro de análise de variância.
anova(m0)## Analysis of Variance Table
##
## Response: prod
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 3 50.53 16.843 1.9709 0.1493
## vari 7 919.72 131.389 15.3744 5.723e-07 ***
## Residuals 21 179.46 8.546
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Diagnóstico do modelo
##-----------------------------------------------------------------------------
## Análise de diagnóstico dos pressupostos via resíduos.
## As inferências acima só passam ser valiadas após confirmação de não
## haver afastamento dos pressupostos.
par(mfrow=c(2,2)); plot(m0); layout(1)
##-----------------------------------------------------------------------------
## Transformar os dados pode diminuir a fuga dos pressupostos?
MASS::boxcox(m0)
Comparaçõ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 B 116-51 22.500 1.461673 21 15.393319 6.524482e-13 19.460284 25.53972
## 2 B 1-52 22.275 1.461673 21 15.239386 7.925028e-13 19.235284 25.31472
## 3 B 25-50 E 16.500 1.461673 21 11.288434 2.228887e-10 13.460284 19.53972
## 4 B 72-53 A 22.800 1.461673 21 15.598564 5.047112e-13 19.760284 25.83972
## 5 Buena Vista 12.425 1.461673 21 8.500533 3.070928e-08 9.385284 15.46472
## 6 Huinkul 25.050 1.461673 21 17.137896 8.027793e-14 22.010284 28.08972
## 7 Kennebec 10.700 1.461673 21 7.320379 3.315981e-07 7.660284 13.73972
## 8 S. Rafalela 25.450 1.461673 21 17.411555 5.877653e-14 22.410284 28.48972## 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 ~ bloc + vari, data = da)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## B 1-52 - B 116-51 == 0 -0.225 2.067 -0.109 0.914357
## B 25-50 E - B 116-51 == 0 -6.000 2.067 -2.903 0.017029 *
## B 72-53 A - B 116-51 == 0 0.300 2.067 0.145 0.914357
## Buena Vista - B 116-51 == 0 -10.075 2.067 -4.874 0.000251 ***
## Huinkul - B 116-51 == 0 2.550 2.067 1.234 0.293976
## Kennebec - B 116-51 == 0 -11.800 2.067 -5.708 5.36e-05 ***
## S. Rafalela - B 116-51 == 0 2.950 2.067 1.427 0.247947
## B 25-50 E - B 1-52 == 0 -5.775 2.067 -2.794 0.019042 *
## B 72-53 A - B 1-52 == 0 0.525 2.067 0.254 0.898222
## Buena Vista - B 1-52 == 0 -9.850 2.067 -4.765 0.000293 ***
## Huinkul - B 1-52 == 0 2.775 2.067 1.342 0.271297
## Kennebec - B 1-52 == 0 -11.575 2.067 -5.600 5.91e-05 ***
## S. Rafalela - B 1-52 == 0 3.175 2.067 1.536 0.216969
## B 72-53 A - B 25-50 E == 0 6.300 2.067 3.048 0.013172 *
## Buena Vista - B 25-50 E == 0 -4.075 2.067 -1.971 0.102125
## Huinkul - B 25-50 E == 0 8.550 2.067 4.136 0.001095 **
## Kennebec - B 25-50 E == 0 -5.800 2.067 -2.806 0.019042 *
## S. Rafalela - B 25-50 E == 0 8.950 2.067 4.330 0.000752 ***
## Buena Vista - B 72-53 A == 0 -10.375 2.067 -5.019 0.000201 ***
## Huinkul - B 72-53 A == 0 2.250 2.067 1.088 0.351487
## Kennebec - B 72-53 A == 0 -12.100 2.067 -5.854 4.62e-05 ***
## S. Rafalela - B 72-53 A == 0 2.650 2.067 1.282 0.285098
## Huinkul - Buena Vista == 0 12.625 2.067 6.108 3.25e-05 ***
## Kennebec - Buena Vista == 0 -1.725 2.067 -0.834 0.482294
## S. Rafalela - Buena Vista == 0 13.025 2.067 6.301 2.81e-05 ***
## Kennebec - Huinkul == 0 -14.350 2.067 -6.942 1.04e-05 ***
## S. Rafalela - Huinkul == 0 0.400 2.067 0.194 0.913685
## S. Rafalela - Kennebec == 0 14.750 2.067 7.136 1.04e-05 ***
## ---
## 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## B 116-51 B 1-52 B 25-50 E B 72-53 A Buena Vista Huinkul Kennebec
## "a" "a" "b" "a" "bc" "a" "c"
## S. Rafalela
## "a"## Formatando média seguida de letras.
p0$cld <- l0$mcletters$Letters
p0$let <- sprintf("%.2f %s", p0$estimate, p0$cld)
p0$let## [1] "22.50 a" "22.28 a" "16.50 b" "22.80 a" "12.43 bc" "25.05 a" "10.70 c"
## [8] "25.45 a"##-----------------------------------------------------------------------------
## Representação gráfica dos resultados.
xlab <- "Produção"
ylab <- "Variedades de batata"
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] "Kennebec" "Buena Vista" "B 25-50 E" "B 1-52" "B 116-51" "B 72-53 A"
## [7] "Huinkul" "S. Rafalela"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)
})
Por se tratar de um experimento regular (não desbalanceado), pode-se fazer essa análise com o pacote ExpDes. Deve ser mecionado que é aconselhável antes declarar o modelo na função lm() e fazer a análise dos resíduos. Apenas depois de confirmada a adequação dos pressupostos, feita as transformações, caso necessárias, fazer a análise com o pacote ExpDes.
##-----------------------------------------------------------------------------
## 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"## Complete randomized design.
with(da, rbd(treat=vari, block=bloc, resp=prod,
quali=TRUE, mcomp="tukey"))## ------------------------------------------------------------------------
## Analysis of Variance Table
## ------------------------------------------------------------------------
## DF SS MS Fc Pr>Fc
## Treatament 7 919.72 131.389 15.3744 0.000001
## Block 3 50.53 16.843 1.9709 0.149251
## Residuals 21 179.46 8.546
## Total 31 1149.72
## ------------------------------------------------------------------------
## CV = 14.83 %
##
## ------------------------------------------------------------------------
## Shapiro-Wilk normality test
## p-value: 0.4709591
## According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
## ------------------------------------------------------------------------
##
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a S. Rafalela 25.45
## a Huinkul 25.05
## ab B 72-53 A 22.8
## ab B 116-51 22.5
## ab B 1-52 22.275
## bc B 25-50 E 16.5
## c Buena Vista 12.425
## c Kennebec 10.7
## ------------------------------------------------------------------------## Tem um possível bug.
## with(da, rbd(treat=vari2, block=bloc, resp=prod,
## quali=TRUE, mcomp="sk"))
## Error in scan(file, what, nmax, sep, dec, quote, skip, nlines, na.strings, :
## line 1 did not have 3 elements
## Opções para fazer testes de comparação múltiplas são os pacotes:
## * agricolae: http://cran.r-project.org/web/packages/agricolae/
## * tukeyC: http://cran.r-project.org/web/packages/TukeyC/
## * Skott-Knott: http://cran.r-project.org/web/packages/ScottKnott/