##-----------------------------------------------------------------------------
## 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_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"))
\[ \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.275 -1.644 0.063 1.056 5.650
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.550 1.714 11.99 7.4e-11 ***
## blocII 3.500 1.462 2.39 0.0261 *
## blocIII 2.275 1.462 1.56 0.1345
## blocIV 2.025 1.462 1.39 0.1805
## variB 1-52 -0.225 2.067 -0.11 0.9144
## variB 25-50 E -6.000 2.067 -2.90 0.0085 **
## variB 72-53 A 0.300 2.067 0.15 0.8860
## variBuena Vista -10.075 2.067 -4.87 8.1e-05 ***
## variHuinkul 2.550 2.067 1.23 0.2310
## variKennebec -11.800 2.067 -5.71 1.1e-05 ***
## variS. Rafalela 2.950 2.067 1.43 0.1682
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.92 on 21 degrees of freedom
## Multiple R-squared: 0.844, Adjusted R-squared: 0.77
## F-statistic: 11.4 on 10 and 21 DF, p-value: 2.15e-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 1
contrasts(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 51 16.8 1.97 0.15
## vari 7 920 131.4 15.37 5.7e-07 ***
## Residuals 21 179 8.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## 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.
## 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.50 1.462 21 15.393 6.524e-13 19.460 25.54
## 2 B 1-52 22.28 1.462 21 15.239 7.925e-13 19.235 25.31
## 3 B 25-50 E 16.50 1.462 21 11.288 2.229e-10 13.460 19.54
## 4 B 72-53 A 22.80 1.462 21 15.599 5.047e-13 19.760 25.84
## 5 Buena Vista 12.43 1.462 21 8.501 3.071e-08 9.385 15.46
## 6 Huinkul 25.05 1.462 21 17.138 8.028e-14 22.010 28.09
## 7 Kennebec 10.70 1.462 21 7.320 3.316e-07 7.660 13.74
## 8 S. Rafalela 25.45 1.462 21 17.412 5.878e-14 22.410 28.49
## 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.11 0.91436
## B 25-50 E - B 116-51 == 0 -6.000 2.067 -2.90 0.01703 *
## B 72-53 A - B 116-51 == 0 0.300 2.067 0.15 0.91436
## Buena Vista - B 116-51 == 0 -10.075 2.067 -4.87 0.00025 ***
## Huinkul - B 116-51 == 0 2.550 2.067 1.23 0.29398
## Kennebec - B 116-51 == 0 -11.800 2.067 -5.71 5.4e-05 ***
## S. Rafalela - B 116-51 == 0 2.950 2.067 1.43 0.24795
## B 25-50 E - B 1-52 == 0 -5.775 2.067 -2.79 0.01904 *
## B 72-53 A - B 1-52 == 0 0.525 2.067 0.25 0.89822
## Buena Vista - B 1-52 == 0 -9.850 2.067 -4.77 0.00029 ***
## Huinkul - B 1-52 == 0 2.775 2.067 1.34 0.27130
## Kennebec - B 1-52 == 0 -11.575 2.067 -5.60 5.9e-05 ***
## S. Rafalela - B 1-52 == 0 3.175 2.067 1.54 0.21697
## B 72-53 A - B 25-50 E == 0 6.300 2.067 3.05 0.01317 *
## Buena Vista - B 25-50 E == 0 -4.075 2.067 -1.97 0.10212
## Huinkul - B 25-50 E == 0 8.550 2.067 4.14 0.00109 **
## Kennebec - B 25-50 E == 0 -5.800 2.067 -2.81 0.01904 *
## S. Rafalela - B 25-50 E == 0 8.950 2.067 4.33 0.00075 ***
## Buena Vista - B 72-53 A == 0 -10.375 2.067 -5.02 0.00020 ***
## Huinkul - B 72-53 A == 0 2.250 2.067 1.09 0.35149
## Kennebec - B 72-53 A == 0 -12.100 2.067 -5.85 4.6e-05 ***
## S. Rafalela - B 72-53 A == 0 2.650 2.067 1.28 0.28510
## Huinkul - Buena Vista == 0 12.625 2.067 6.11 3.2e-05 ***
## Kennebec - Buena Vista == 0 -1.725 2.067 -0.83 0.48229
## S. Rafalela - Buena Vista == 0 13.025 2.067 6.30 2.8e-05 ***
## Kennebec - Huinkul == 0 -14.350 2.067 -6.94 1.0e-05 ***
## S. Rafalela - Huinkul == 0 0.400 2.067 0.19 0.91369
## S. Rafalela - Kennebec == 0 14.750 2.067 7.14 1.0e-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$let <- sprintf("%.2f %s", p0$estimate, l0$mcletters$Letters)
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)
})
