##-----------------------------------------------------------------------------
## 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)
Os dados são de um experimento feito para avaliar o comprimento de conídios para 6 isolados do fungo Colletrotrichum lindemuthianum. O Experimento foi organizado em delineamento inteiramente casualizado com no máximo 10 repetições por isolado. Houveram parcelas perdidas.
##-----------------------------------------------------------------------------
## Lendo arquivos de dados.
## Url de um arquivo com dados.
url <- "http://www.leg.ufpr.br/~walmes/data/ExpGen_p66.txt"
## Importa os dados a partir do arquivo na web.
da <- read.table(file=url, header=TRUE, sep="\t", dec=",")
## Nomes das variáveis em minúsculo.
names(da) <- tolower(names(da))
## Mostra a estrutura do objeto.
str(da)
## 'data.frame': 60 obs. of 2 variables:
## $ iso : Factor w/ 6 levels "i1","i2","i3",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ comp: num 12.5 13.3 13.5 15.3 15.3 ...
##-----------------------------------------------------------------------------
## Tabela com o valores observados.
do.call(cbind, split(da$comp, da$iso))
## i1 i2 i3 i4 i5 i6
## [1,] 12.47 15.08 15.30 15.75 13.2 13.50
## [2,] 13.27 12.60 15.98 16.20 11.4 13.28
## [3,] 13.50 14.40 13.50 15.75 12.4 12.15
## [4,] 15.30 13.05 12.15 16.20 12.6 11.48
## [5,] 15.30 13.72 16.42 16.42 12.6 13.28
## [6,] 16.42 12.80 14.18 17.55 13.4 13.78
## [7,] 14.18 13.05 16.42 14.40 11.2 12.60
## [8,] 16.20 15.08 15.25 15.75 13.4 12.38
## [9,] 11.70 NA 13.72 20.25 17.4 12.82
## [10,] 14.62 NA 15.52 15.75 NA NA
##-----------------------------------------------------------------------------
## Análise exploratória.
xyplot(comp~iso, data=da, jitter.x=TRUE,
xlab="Isolado", ylab="Comprimento do conídio")+
layer(panel.abline(v=1:nlevels(da$iso), col="gray50"))
\[ \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 isolado \(i\) sobre a variável resposta \(Y\), \(\mu\) é a média de \(Y\) na ausência do efeito dos isolados, \(\mu_i\) é a média para as observações em cada nível de isolado e \(\sigma^2\) é a variância das observações ao redor desse valor médio.
##-----------------------------------------------------------------------------
## Estimação.
m0 <- lm(comp~iso, data=da)
## Avaliação dos pressupostos.
par(mfrow=c(2,2))
plot(m0); layout(1)
## Estimativas e erros padrões.
summary(m0)
##
## Call:
## lm(formula = comp ~ iso, data = da)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.694 -0.672 -0.159 0.681 4.333
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.296 0.446 32.06 <2e-16 ***
## isoi2 -0.574 0.669 -0.86 0.3953
## isoi3 0.548 0.631 0.87 0.3890
## isoi4 2.106 0.631 3.34 0.0016 **
## isoi5 -1.229 0.648 -1.90 0.0635 .
## isoi6 -1.488 0.648 -2.30 0.0258 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.41 on 50 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.456, Adjusted R-squared: 0.402
## F-statistic: 8.39 on 5 and 50 DF, p-value: 8.1e-06
## Quadro de análise de variância.
anova(m0)
## Analysis of Variance Table
##
## Response: comp
## Df Sum Sq Mean Sq F value Pr(>F)
## iso 5 83.4 16.68 8.39 8.1e-06 ***
## Residuals 50 99.4 1.99
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Comparações múltiplas.
## Médias ajustadas (LSmeans).
p0 <- LSmeans(m0, effect="iso", level=0.95)
## Criando a tabela com as estimativas.
p0 <- cbind(p0$grid, as.data.frame(p0$coef))
p0 <- equallevels(p0, da)
p0
## iso estimate se df t.stat p.value lwr upr
## 1 i1 14.30 0.4459 50 32.06 5.249e-35 13.40 15.19
## 2 i2 13.72 0.4985 50 27.53 7.196e-32 12.72 14.72
## 3 i3 14.84 0.4459 50 33.29 8.691e-36 13.95 15.74
## 4 i4 16.40 0.4459 50 36.78 7.166e-38 15.51 17.30
## 5 i5 13.07 0.4700 50 27.80 4.515e-32 12.12 14.01
## 6 i6 12.81 0.4700 50 27.25 1.156e-31 11.86 13.75
## Comparações múltiplas, contrastes de Tukey.
## Método de correção de p-valor: single-step.
g0 <- summary(glht(m0, linfct=mcp(iso="Tukey")),
test=adjusted(type="fdr"))
g0
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = comp ~ iso, data = da)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## i2 - i1 == 0 -0.574 0.669 -0.86 0.4235
## i3 - i1 == 0 0.548 0.631 0.87 0.4235
## i4 - i1 == 0 2.106 0.631 3.34 0.0060 **
## i5 - i1 == 0 -1.229 0.648 -1.90 0.1059
## i6 - i1 == 0 -1.488 0.648 -2.30 0.0484 *
## i3 - i2 == 0 1.121 0.669 1.68 0.1498
## i4 - i2 == 0 2.680 0.669 4.01 0.0010 **
## i5 - i2 == 0 -0.656 0.685 -0.96 0.4235
## i6 - i2 == 0 -0.915 0.685 -1.34 0.2562
## i4 - i3 == 0 1.558 0.631 2.47 0.0363 *
## i5 - i3 == 0 -1.777 0.648 -2.74 0.0211 *
## i6 - i3 == 0 -2.036 0.648 -3.14 0.0084 **
## i5 - i4 == 0 -3.335 0.648 -5.15 3.3e-05 ***
## i6 - i4 == 0 -3.594 0.648 -5.55 1.6e-05 ***
## i6 - i5 == 0 -0.259 0.665 -0.39 0.6986
## ---
## 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)
## Formatando média seguida de letras.
p0$let <- sprintf("%.2f %s", p0$estimate, l0$mcletters$Letters)
p0$let
## [1] "14.30 bc" "13.72 ac" "14.84 c" "16.40 d" "13.07 ab" "12.81 a"
##-----------------------------------------------------------------------------
## Representação gráfica dos resultados.
xlab <- expression("Comprimento de conídios"~(mu*m))
ylab <- expression("Isolados de"~italic("Colletotrichum lindemuthianum"))
sub <- list(
"Médias seguidas de mesma letra indicam diferença nula à 5%.",
font=1, cex=0.8)
p0 <- transform(p0, iso=reorder(iso, estimate))
p0 <- p0[order(p0$iso),]
segplot(iso~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)
})
