##-----------------------------------------------------------------------------
## 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/volume.txt"
## Importa os dados a partir do arquivo na web.
da <- read.table(file=url, header=TRUE, sep="\t", nrow=27)
## Mostra a estrutura do objeto.
str(da)
## 'data.frame': 27 obs. of 4 variables:
## $ gen : Factor w/ 9 levels "ATF06B","ATF40B",..: 1 1 1 2 2 2 3 3 3 4 ...
## $ dose: int 0 0 0 0 0 0 0 0 0 0 ...
## $ rep : int 1 2 3 1 2 3 1 2 3 1 ...
## $ volu: num 1.19 1.17 1.12 1.27 1.85 ...
## Os dados no formato amplo (não empilhado).
unstack(da, form=volu~gen)
## ATF06B ATF40B ATF54B BR001B BR005B BR007B BR008B P9401 SC283
## 1 1.190 1.275 0.939 1.421 0.520 1.431 1.398 1.141 0.872
## 2 1.173 1.849 1.056 1.240 0.567 1.425 1.552 1.032 1.093
## 3 1.118 1.306 0.712 1.614 0.416 1.152 1.127 1.226 1.136
## Os dados são de um experimento para verificar o ganho de peso em
## função de tipos de rações. Foram usandos 4 níveis do fator ração e 5
## repetições.
##-----------------------------------------------------------------------------
## Análise exploratória.
## xyplot(volu~gen, data=da, type=c("p","a"))
xyplot(volu~reorder(gen, volu), data=da, type=c("p","a"))
\[ \begin{aligned} Y|\text{fontes de variação} &\sim \text{Normal}(\mu_i, \sigma^2) \newline \mu_i &= \mu+\alpha_i \end{aligned} \]
em que \(\alpha_i\) é o efeito do genótipo \(i\) sobre a variável resposta \(Y\), \(\mu\) é a média de \(Y\) na ausência do efeito dos genótipos, \(\mu_i\) é a média para as observações em cada nível de genótipo e \(\sigma^2\) é a variância das observações ao redor desse valor médio.
##-----------------------------------------------------------------------------
## Estimação.
m0 <- lm(volu~gen, data=da)
## Estimativas e erros padrões.
summary(m0)
##
## Call:
## lm(formula = volu ~ gen, data = da)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.232 -0.131 0.019 0.091 0.372
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.1603 0.1015 11.43 1.1e-09 ***
## genATF40B 0.3163 0.1436 2.20 0.04088 *
## genATF54B -0.2580 0.1436 -1.80 0.08920 .
## genBR001B 0.2647 0.1436 1.84 0.08185 .
## genBR005B -0.6593 0.1436 -4.59 0.00023 ***
## genBR007B 0.1757 0.1436 1.22 0.23700
## genBR008B 0.1987 0.1436 1.38 0.18345
## genP9401 -0.0273 0.1436 -0.19 0.85117
## genSC283 -0.1267 0.1436 -0.88 0.38937
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.176 on 18 degrees of freedom
## Multiple R-squared: 0.803, Adjusted R-squared: 0.716
## F-statistic: 9.18 on 8 and 18 DF, p-value: 5.46e-05
## Restrição paramétrica de zerer o primeiro nível.
contrasts(da$gen)
## ATF40B ATF54B BR001B BR005B BR007B BR008B P9401 SC283
## ATF06B 0 0 0 0 0 0 0 0
## ATF40B 1 0 0 0 0 0 0 0
## ATF54B 0 1 0 0 0 0 0 0
## BR001B 0 0 1 0 0 0 0 0
## BR005B 0 0 0 1 0 0 0 0
## BR007B 0 0 0 0 1 0 0 0
## BR008B 0 0 0 0 0 1 0 0
## P9401 0 0 0 0 0 0 1 0
## SC283 0 0 0 0 0 0 0 1
## model.matrix(m0)
##-----------------------------------------------------------------------------
## Existe efeito de genótipos? Ou seja, todos os \alpha_i podem ser
## considerados iguais à zero?
## Quadro de análise de variância.
anova(m0)
## Analysis of Variance Table
##
## Response: volu
## Df Sum Sq Mean Sq F value Pr(>F)
## gen 8 2.271 0.2839 9.18 5.5e-05 ***
## Residuals 18 0.557 0.0309
## ---
## 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)
##-----------------------------------------------------------------------------
## Ajuste com a variável transformada.
m1 <- lm(log(volu)~gen, data=da)
par(mfrow=c(2,2))
plot(m1); layout(1)
anova(m1)
## Analysis of Variance Table
##
## Response: log(volu)
## Df Sum Sq Mean Sq F value Pr(>F)
## gen 8 2.710 0.339 15.4 1.4e-06 ***
## Residuals 18 0.396 0.022
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Passa m1 para m0.
m0 <- m1
##-----------------------------------------------------------------------------
## Comparações múltiplas.
## Médias ajustadas (LSmeans).
p0 <- LSmeans(m0, effect="gen", level=0.95)
## Criando a tabela com as estimativas.
p0 <- cbind(p0$grid, as.data.frame(p0$coef))
p0 <- equallevels(p0, da)
p0
## gen estimate se df t.stat p.value lwr upr
## 1 ATF06B 0.14835 0.08566 18 1.7319 0.1003861 -0.03161 0.32831
## 2 ATF40B 0.37485 0.08566 18 4.3762 0.0003641 0.19489 0.55481
## 3 ATF54B -0.11604 0.08566 18 -1.3547 1.8077373 -0.29600 0.06392
## 4 BR001B 0.34840 0.08566 18 4.0673 0.0007229 0.16844 0.52836
## 5 BR005B -0.69946 0.08566 18 -8.1658 1.9999998 -0.87942 -0.51950
## 6 BR007B 0.28468 0.08566 18 3.3235 0.0037800 0.10472 0.46464
## 7 BR008B 0.29805 0.08566 18 3.4796 0.0026755 0.11809 0.47801
## 8 P9401 0.12239 0.08566 18 1.4288 0.1701868 -0.05757 0.30235
## 9 SC283 0.02649 0.08566 18 0.3093 0.7606674 -0.15347 0.20645
## Comparações múltiplas, contrastes de Tukey.
## Método de correção de p-valor: single-step.
g0 <- summary(glht(m0, linfct=mcp(gen="Tukey")),
test=adjusted(type="fdr"))
g0
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = log(volu) ~ gen, data = da)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## ATF40B - ATF06B == 0 0.2265 0.1211 1.87 0.13451
## ATF54B - ATF06B == 0 -0.2644 0.1211 -2.18 0.09575 .
## BR001B - ATF06B == 0 0.2000 0.1211 1.65 0.18982
## BR005B - ATF06B == 0 -0.8478 0.1211 -7.00 1.1e-05 ***
## BR007B - ATF06B == 0 0.1363 0.1211 1.13 0.36694
## BR008B - ATF06B == 0 0.1497 0.1211 1.24 0.33471
## P9401 - ATF06B == 0 -0.0260 0.1211 -0.21 0.85647
## SC283 - ATF06B == 0 -0.1219 0.1211 -1.01 0.42140
## ATF54B - ATF40B == 0 -0.4909 0.1211 -4.05 0.00299 **
## BR001B - ATF40B == 0 -0.0265 0.1211 -0.22 0.85647
## BR005B - ATF40B == 0 -1.0743 0.1211 -8.87 1.4e-06 ***
## BR007B - ATF40B == 0 -0.0902 0.1211 -0.74 0.55951
## BR008B - ATF40B == 0 -0.0768 0.1211 -0.63 0.62017
## P9401 - ATF40B == 0 -0.2525 0.1211 -2.08 0.10334
## SC283 - ATF40B == 0 -0.3484 0.1211 -2.88 0.02785 *
## BR001B - ATF54B == 0 0.4644 0.1211 3.83 0.00438 **
## BR005B - ATF54B == 0 -0.5834 0.1211 -4.82 0.00062 ***
## BR007B - ATF54B == 0 0.4007 0.1211 3.31 0.01173 *
## BR008B - ATF54B == 0 0.4141 0.1211 3.42 0.01003 *
## P9401 - ATF54B == 0 0.2384 0.1211 1.97 0.12247
## SC283 - ATF54B == 0 0.1425 0.1211 1.18 0.35262
## BR005B - BR001B == 0 -1.0479 0.1211 -8.65 1.4e-06 ***
## BR007B - BR001B == 0 -0.0637 0.1211 -0.53 0.68100
## BR008B - BR001B == 0 -0.0503 0.1211 -0.42 0.74466
## P9401 - BR001B == 0 -0.2260 0.1211 -1.87 0.13451
## SC283 - BR001B == 0 -0.3219 0.1211 -2.66 0.04124 *
## BR007B - BR005B == 0 0.9841 0.1211 8.12 1.8e-06 ***
## BR008B - BR005B == 0 0.9975 0.1211 8.23 1.8e-06 ***
## P9401 - BR005B == 0 0.8219 0.1211 6.78 1.4e-05 ***
## SC283 - BR005B == 0 0.7260 0.1211 5.99 5.9e-05 ***
## BR008B - BR007B == 0 0.0134 0.1211 0.11 0.91336
## P9401 - BR007B == 0 -0.1623 0.1211 -1.34 0.29549
## SC283 - BR007B == 0 -0.2582 0.1211 -2.13 0.09973 .
## P9401 - BR008B == 0 -0.1757 0.1211 -1.45 0.25706
## SC283 - BR008B == 0 -0.2716 0.1211 -2.24 0.09076 .
## SC283 - P9401 == 0 -0.0959 0.1211 -0.79 0.54483
## ---
## 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)
l0
## ATF06B ATF40B ATF54B BR001B BR005B BR007B BR008B P9401 SC283
## "bd" "d" "b" "d" "a" "cd" "cd" "bd" "bc"
## Formatando média seguida de letras.
p0$let <- sprintf("%.2f %s", p0$estimate, l0$mcletters$Letters)
p0$let
## [1] "0.15 bd" "0.37 d" "-0.12 b" "0.35 d" "-0.70 a" "0.28 cd" "0.30 cd" "0.12 bd"
## [9] "0.03 bc"
##-----------------------------------------------------------------------------
## Representação gráfica dos resultados.
xlab <- "Volume das raízes (log)"
ylab <- "Genótipo de sorgo"
sub <- list(
"Médias seguidas de mesma letra indicam diferença nula à 5%.",
font=1, cex=0.8)
p0 <- transform(p0, gen=reorder(gen, estimate))
levels(p0$gen)
## [1] "BR005B" "ATF54B" "SC283" "P9401" "ATF06B" "BR007B" "BR008B" "BR001B" "ATF40B"
p0 <- p0[order(p0$gen),]
segplot(gen~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)
})
##-----------------------------------------------------------------------------
## Representando na escala original.
sel <- c("estimate","lwr","upr")
p0[,sel] <- exp(p0[,sel])
xlab <- "Volume das raízes"
p0$let <- sprintf("%.2f %s", p0$estimate, l0$mcletters$Letters)
segplot(gen~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)
})
