Grasiela Bruzamarello Tognon
Walmes Marques Zeviani
Francine Lorena Cuquel
##----------------------------------------------------------------------
## Definições da sessão.
## Lista de pacotes a serem usados na sessão.
pkg <- c("lattice", "latticeExtra", "gridExtra",
"doBy", "plyr", "multcomp", "wzRfun")
sapply(pkg, require, character.only=TRUE)
## lattice latticeExtra gridExtra doBy plyr multcomp
## TRUE TRUE TRUE TRUE TRUE TRUE
## wzRfun
## TRUE
sessionInfo()
## R version 3.2.0 (2015-04-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.2 LTS
##
## 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 grid stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] wzRfun_0.4 multcomp_1.4-0 TH.data_1.0-6 mvtnorm_1.0-2
## [5] plyr_1.8.3 doBy_4.5-13 survival_2.38-2 gridExtra_0.9.1
## [9] latticeExtra_0.6-26 RColorBrewer_1.1-2 lattice_0.20-31 knitr_1.10.5
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.11.6 codetools_0.2-11 zoo_1.7-12 MASS_7.3-41 formatR_1.2
## [6] magrittr_1.5 evaluate_0.7 stringi_0.4-1 Matrix_1.2-1 sandwich_2.3-3
## [11] splines_3.2.0 tools_3.2.0 stringr_1.0.0
citation()
##
## To cite R in publications use:
##
## R Core Team (2015). R: A language and environment for statistical computing. R
## Foundation for Statistical Computing, Vienna, Austria. URL
## http://www.R-project.org/.
##
## A BibTeX entry for LaTeX users is
##
## @Manual{,
## title = {R: A Language and Environment for Statistical Computing},
## author = {{R Core Team}},
## organization = {R Foundation for Statistical Computing},
## address = {Vienna, Austria},
## year = {2015},
## url = {http://www.R-project.org/},
## }
##
## We have invested a lot of time and effort in creating R, please cite it when
## using it for data analysis. See also 'citation("pkgname")' for citing R
## packages.
##----------------------------------------------------------------------
## Definições gráficas trellis.
## display.brewer.all()
## mycol <- brewer.pal(6, "Set1")
mycol <- c(brewer.pal(6, "Dark2"), brewer.pal(6, "Set1"))
## names(trellis.par.get())
## dput(trellis.par.get())
ps <- list(box.rectangle=list(col=1, fill=c("gray70")),
box.umbrella=list(col=1, lty=1),
dot.symbol=list(col=1),
dot.line=list(col=1, lty=3),
plot.symbol=list(col=1, cex=0.7),
plot.line=list(col=1),
plot.polygon=list(col="gray80"),
superpose.line=list(col=mycol),
superpose.symbol=list(col=mycol),
superpose.polygon=list(col=mycol),
strip.background=list(col=c("gray90","gray50"))
)
trellis.par.set(ps)
## show.settings()
## show.settings(ps)
##-----------------------------------------------------------------------------
## Ler.
rm(list=ls())
## list.files(pattern="*.txt")
da <- read.table("ivg_tridentata.txt", header=TRUE, sep="\t")
da$armaz <- factor(da$armaz)
da$temper <- factor(da$temper, levels=c("NULA","5","18"))
str(da)
## 'data.frame': 52 obs. of 6 variables:
## $ armaz : Factor w/ 4 levels "0","2","4","6": 1 1 1 1 2 2 2 2 2 2 ...
## $ temper : Factor w/ 3 levels "NULA","5","18": 1 1 1 1 2 2 2 2 2 2 ...
## $ embal : Factor w/ 3 levels "NULA","Papel Kraft",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ rept : int 1 2 3 4 1 2 3 4 1 2 ...
## $ ngerm50: int 31 30 29 32 31 28 31 24 21 28 ...
## $ ivg : num 2.86 2.52 2.67 2.85 2.27 ...
##-----------------------------------------------------------------------------
## Ver.
ftable(xtabs(~temper+embal+armaz, da))
## armaz 0 2 4 6
## temper embal
## NULA NULA 4 0 0 0
## Papel Kraft 0 0 0 0
## Polietileno 0 0 0 0
## 5 NULA 0 0 0 0
## Papel Kraft 0 4 4 4
## Polietileno 0 4 4 4
## 18 NULA 0 0 0 0
## Papel Kraft 0 4 4 4
## Polietileno 0 4 4 4
##----------------------------------------------------------------------
## Fatorial incompleto.
xyplot(ngerm50/50~armaz|temper, groups=embal, data=da,
type=c("p","a"), jitter.x=TRUE)
xyplot(ivg~armaz|temper, groups=embal, data=da,
type=c("p","a"), jitter.x=TRUE)
##-----------------------------------------------------------------------------
m0 <- lm(ivg~armaz*temper*embal, data=da)
## par(mfrow=c(2,2)); plot(m0); layout(1)
## MASS::boxcox(m0)
m0 <- update(m0, sqrt(.)~.)
## Quadro de testes F baseados na deviance para termos do modelo.
anova(m0, test="F")
## Analysis of Variance Table
##
## Response: sqrt(ivg)
## Df Sum Sq Mean Sq F value Pr(>F)
## armaz 3 1.02467 0.34156 17.7566 2e-07 ***
## temper 1 0.05003 0.05003 2.6011 0.114851
## embal 1 0.10713 0.10713 5.5692 0.023378 *
## armaz:temper 2 0.22847 0.11424 5.9388 0.005604 **
## armaz:embal 2 0.07998 0.03999 2.0790 0.138697
## temper:embal 1 0.02079 0.02079 1.0808 0.304929
## armaz:temper:embal 2 0.01889 0.00945 0.4911 0.615706
## Residuals 39 0.75019 0.01924
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Modelo com abandono de termos não significativos.
m1 <- update(m0, .~armaz*temper)
anova(m0, m1, test="F")
## Analysis of Variance Table
##
## Model 1: sqrt(ivg) ~ armaz + temper + embal + armaz:temper + armaz:embal +
## temper:embal + armaz:temper:embal
## Model 2: sqrt(ivg) ~ armaz + temper + armaz:temper
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 39 0.75019
## 2 45 0.97697 -6 -0.22679 1.965 0.09445 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(m1, test="F")
## Analysis of Variance Table
##
## Response: sqrt(ivg)
## Df Sum Sq Mean Sq F value Pr(>F)
## armaz 3 1.02467 0.34156 15.7323 3.89e-07 ***
## temper 1 0.05003 0.05003 2.3046 0.135986
## armaz:temper 2 0.22847 0.11424 5.2618 0.008842 **
## Residuals 45 0.97697 0.02171
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Não existe efeito de níveis de embalagem. Existe interação entre
## temperatura e período de armazenamento.
##-----------------------------------------------------------------------------
## Ajustar modelo a partir de uma matriz de delineamento com colunas
## apenas para os efeitos estimáveis.
X <- model.matrix(m1)
X <- X[,!is.na(coef(m1))]
ncol(X)
## [1] 7
m2 <- update(m1, .~0+X)
## anova(m2, m1)
##----------------------------------------------------------------------
## Gerar matriz para estimar as médias ajustadas.
grid <- unique(subset(da, select=c(armaz, temper)))
grid
## armaz temper
## 1 0 NULA
## 5 2 5
## 13 2 18
## 21 4 5
## 29 4 18
## 37 6 5
## 45 6 18
Xm <- model.matrix(m1)[rownames(grid),]
Xm <- Xm[,!is.na(coef(m1))]
ci <- confint(glht(m2, linfct=Xm), calpha=univariate_calpha())
grid <- cbind(grid, ci$confint)
grid
## armaz temper Estimate lwr upr
## 1 0 NULA 1.649982 1.501598 1.798366
## 5 2 5 1.350851 1.245928 1.455775
## 13 2 18 1.380471 1.275548 1.485395
## 21 4 5 1.607449 1.502526 1.712372
## 29 4 18 1.643785 1.538861 1.748708
## 37 6 5 1.832501 1.727577 1.937424
## 45 6 18 1.572829 1.467906 1.677753
## Médias amostrais, só para conferir.
## aggregate(sqrt(ivg)~temper+armaz, data=da, FUN=mean)
## segplot(armaz~lwr+upr, centers=Estimate, data=grid, draw=FALSE,
## groups=temper, f=0.1, panel=panel.segplotBy,
## ylab=expression(Período~de~armazenamento~(meses)),
## xlab="logit da Germinação")
## Passa a inversa nos extremos do IC.
ci <- sapply(grid[,3:5], function(x) x^2)
colnames(ci) <- paste0("p", colnames(ci))
grid <- cbind(grid, ci)
key <- list(type="o", divide=1,
title=expression(Temperatura~(degree*C)), cex.title=1.1,
text=list(levels(grid$temper)[-1]),
lines=list(pch=2:3, lty=1))
## segplot(armaz~plwr+pupr, centers=pEstimate, data=grid, draw=FALSE,
## groups=temper, pch=grid$temper, f=0.1, panel=panel.segplotBy,
## ylab=expression(Período~de~armazenamento~(meses)),
## ## xlab="IVG",
## key=key, par.settings=ps)
##-----------------------------------------------------------------------------
## Comparações múltiplas, contraste entre os sete parâmetros.
lev <- apply(grid[,1:2], 1, paste, collapse=":")
## lev <- with(grid, paste(armaz, "meses :", temper, " °C"))
lev[grep("NULA", lev)] <- "0 meses"
grid$trat <- factor(lev)
str(grid)
## 'data.frame': 7 obs. of 9 variables:
## $ armaz : Factor w/ 4 levels "0","2","4","6": 1 2 2 3 3 4 4
## $ temper : Factor w/ 3 levels "NULA","5","18": 1 2 3 2 3 2 3
## $ Estimate : num 1.65 1.35 1.38 1.61 1.64 ...
## $ lwr : num 1.5 1.25 1.28 1.5 1.54 ...
## $ upr : num 1.8 1.46 1.49 1.71 1.75 ...
## $ pEstimate: num 2.72 1.82 1.91 2.58 2.7 ...
## $ plwr : num 2.25 1.55 1.63 2.26 2.37 ...
## $ pupr : num 3.23 2.12 2.21 2.93 3.06 ...
## $ trat : Factor w/ 7 levels "0 meses","2:18",..: 1 3 2 5 4 7 6
Xc <- apc(Xm, lev=lev)
g0 <- summary(glht(m2, linfct=Xc), test=adjusted(type="fdr"))
source("/home/walmes/Dropbox/Public/mycld.R")
## g0
g0$focus <- "Tratamentos"
g0$Type <- "Tukey"
grid$cld <- my.cld(g0)$mcletters$Letters
## Tabelas com as estimativas de % germinação.
grid[,-c(3:5)]
## armaz temper pEstimate plwr pupr trat cld
## 1 0 NULA 2.722439 2.254795 3.234119 0 meses bc
## 5 2 5 1.824799 1.552337 2.119280 2:5 a
## 13 2 18 1.905701 1.627023 2.206397 2:18 a
## 21 4 5 2.583892 2.257583 2.932218 4:5 c
## 29 4 18 2.702028 2.368094 3.057979 4:18 c
## 37 6 5 3.358059 2.984524 3.753612 6:5 b
## 45 6 18 2.473793 2.154749 2.814854 6:18 c
##----------------------------------------------------------------------
## Gráfico.
## segplot(trat~plwr+pupr, centers=pEstimate, data=grid, draw=FALSE,
## xlab=expression(IVG~de~italic(B.~milleflora)~"(%)"),
## ylab=expression(Armazenamento~(meses)*" : "*Temperatura~(degree*C)),
## txt=with(grid, sprintf("%0.1f %s", pEstimate, cld)),
## panel=function(z, centers, txt, ...){
## panel.segplot(z=z, centers=centers, ...)
## panel.text(x=centers, y=z, labels=txt, pos=3)
## })
segplot(armaz~plwr+pupr, centers=pEstimate, data=grid, draw=FALSE,
groups=temper, pch=grid$temper, f=0.1,
txt=with(grid, sprintf("%0.1f %s", pEstimate, cld)),
ylab=expression(Período~de~armazenamento~(meses)),
xlab=expression(IVG~de~italic(B.~tridentata)),
key=key,
panel=function(x, y, z, centers, subscripts,
groups, f, txt, ...){
d <- 2*((as.numeric(groups)-1)/(nlevels(groups)-1))-1
z <- as.numeric(z)+f*d
panel.segplot(x, y, z, centers=centers,
subscripts=subscripts, ...)
panel.text(x=centers, y=z, labels=txt, pos=c(3,1))
})
##-----------------------------------------------------------------------------
m0 <- glm(cbind(ngerm50, 50-ngerm50)~armaz*temper*embal,
data=da, family="quasibinomial")
## par(mfrow=c(2,2)); plot(m0); layout(1)
summary(m0)$dispersion
## [1] 1.573133
## Quadro de testes F baseados na deviance para termos do modelo.
anova(m0, test="F")
## Analysis of Deviance Table
##
## Model: quasibinomial, link: logit
##
## Response: cbind(ngerm50, 50 - ngerm50)
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 51 96.732
## armaz 3 0.9844 48 95.748 0.2086 0.88984
## temper 1 0.5597 47 95.188 0.3558 0.55430
## embal 1 7.5308 46 87.657 4.7871 0.03473 *
## armaz:temper 2 15.9279 44 71.730 5.0625 0.01110 *
## armaz:embal 2 5.8306 42 65.899 1.8532 0.17027
## temper:embal 1 1.2140 41 64.685 0.7717 0.38508
## armaz:temper:embal 2 2.8063 39 61.879 0.8919 0.41806
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Modelo com abandono de termos não significativos.
m1 <- update(m0, .~armaz*temper)
anova(m0, m1, test="F")
## Analysis of Deviance Table
##
## Model 1: cbind(ngerm50, 50 - ngerm50) ~ armaz * temper * embal
## Model 2: cbind(ngerm50, 50 - ngerm50) ~ armaz + temper + armaz:temper
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 39 61.879
## 2 45 79.311 -6 -17.432 1.8468 0.1151
anova(m1, test="F")
## Analysis of Deviance Table
##
## Model: quasibinomial, link: logit
##
## Response: cbind(ngerm50, 50 - ngerm50)
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 51 96.732
## armaz 3 0.9844 48 95.748 0.1903 0.90247
## temper 1 0.5597 47 95.188 0.3246 0.57171
## armaz:temper 2 15.8776 45 79.311 4.6037 0.01517 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Não existe efeito de níveis de embalagem. Existe interação entre
## temperatura e período de armazenamento.
##-----------------------------------------------------------------------------
## Ajustar modelo a partir de uma matriz de delineamento com colunas
## apenas para os efeitos estimáveis.
X <- model.matrix(m1)
X <- X[,!is.na(coef(m1))]
ncol(X)
## [1] 7
m2 <- update(m1, .~0+X)
## anova(m2, m1)
##----------------------------------------------------------------------
## Gerar matriz para estimar as médias ajustadas.
grid <- unique(subset(da, select=c(armaz, temper)))
grid
## armaz temper
## 1 0 NULA
## 5 2 5
## 13 2 18
## 21 4 5
## 29 4 18
## 37 6 5
## 45 6 18
Xm <- model.matrix(m1)[rownames(grid),]
Xm <- Xm[,!is.na(coef(m1))]
ci <- confint(glht(m2, linfct=Xm), calpha=univariate_calpha())
grid <- cbind(grid, ci$confint)
grid
## armaz temper Estimate lwr upr
## 1 0 NULA 0.4473122 0.074181368 0.8204431
## 5 2 5 0.2411621 -0.018090488 0.5004146
## 13 2 18 0.5001731 0.234703198 0.7656429
## 21 4 5 0.2614797 0.001897728 0.5210617
## 29 4 18 0.5971325 0.328196276 0.8660688
## 37 6 5 0.5429569 0.276034857 0.8098789
## 45 6 18 0.1402293 -0.117782723 0.3982414
## segplot(armaz~lwr+upr, centers=Estimate, data=grid, draw=FALSE,
## groups=temper, f=0.1, panel=panel.segplotBy,
## ylab=expression(Período~de~armazenamento~(meses)),
## xlab="logit da Germinação")
## Passa a inversa nos extremos do IC.
ci <- sapply(grid[,3:5], m2$family$linkinv)
colnames(ci) <- paste0("p", colnames(ci))
grid <- cbind(grid, 100*ci)
key <- list(type="o", divide=1,
title=expression(Temperatura~(degree*C)), cex.title=1.1,
text=list(levels(grid$temper)[-1]),
lines=list(pch=2:3, lty=1))
## segplot(armaz~plwr+pupr, centers=pEstimate, data=grid, draw=FALSE,
## groups=temper, pch=grid$temper, f=0.1, panel=panel.segplotBy,
## ylab=expression(Período~de~armazenamento~(meses)),
## xlab="Germinação (%)",
## key=key, par.settings=ps)
##-----------------------------------------------------------------------------
## Comparações múltiplas, contraste entre os sete parâmetros.
lev <- apply(grid[,1:2], 1, paste, collapse=":")
## lev <- with(grid, paste(armaz, "meses :", temper, " °C"))
lev[grep("NULA", lev)] <- "0 meses"
grid$trat <- factor(lev)
str(grid)
## 'data.frame': 7 obs. of 9 variables:
## $ armaz : Factor w/ 4 levels "0","2","4","6": 1 2 2 3 3 4 4
## $ temper : Factor w/ 3 levels "NULA","5","18": 1 2 3 2 3 2 3
## $ Estimate : num 0.447 0.241 0.5 0.261 0.597 ...
## $ lwr : num 0.0742 -0.0181 0.2347 0.0019 0.3282 ...
## $ upr : num 0.82 0.5 0.766 0.521 0.866 ...
## $ pEstimate: num 61 56 62.3 56.5 64.5 ...
## $ plwr : num 51.9 49.5 55.8 50 58.1 ...
## $ pupr : num 69.4 62.3 68.3 62.7 70.4 ...
## $ trat : Factor w/ 7 levels "0 meses","2:18",..: 1 3 2 5 4 7 6
Xc <- apc(Xm, lev=lev)
g0 <- summary(glht(m2, linfct=Xc), test=adjusted(type="fdr"))
source("/home/walmes/Dropbox/Public/mycld.R")
## g0
g0$focus <- "Tratamentos"
g0$Type <- "Tukey"
grid$cld <- my.cld(g0)$mcletters$Letters
## Tabelas com as estimativas de % germinação.
grid[,-c(3:5)]
## armaz temper pEstimate plwr pupr trat cld
## 1 0 NULA 61.00 51.85368 69.43304 0 meses a
## 5 2 5 56.00 49.54775 62.25568 2:5 a
## 13 2 18 62.25 55.84079 68.25776 2:18 a
## 21 4 5 56.50 50.04744 62.73960 4:5 a
## 29 4 18 64.50 58.13204 70.39270 4:18 a
## 37 6 5 63.25 56.85738 69.20837 6:5 a
## 45 6 18 53.50 47.05883 59.82651 6:18 a
##----------------------------------------------------------------------
## Gráfico.
## segplot(trat~plwr+pupr, centers=pEstimate, data=grid, draw=FALSE,
## xlab=expression(Germinação~de~italic(B.~milleflora)~"(%)"),
## ylab=expression(Armazenamento~(meses)*" : "*Temperatura~(degree*C)),
## txt=with(grid, sprintf("%0.1f %s", pEstimate, cld)),
## panel=function(z, centers, txt, ...){
## panel.segplot(z=z, centers=centers, ...)
## panel.text(x=centers, y=z, labels=txt, pos=3)
## })
segplot(armaz~plwr+pupr, centers=pEstimate, data=grid, draw=FALSE,
groups=temper, pch=grid$temper, f=0.1,
txt=with(grid, sprintf("%0.1f %s", pEstimate, cld)),
ylab=expression(Período~de~armazenamento~(meses)),
xlab=expression(Germinação~de~italic(B.~tridentata)),
key=key,
panel=function(x, y, z, centers, subscripts,
groups, f, txt, ...){
d <- 2*((as.numeric(groups)-1)/(nlevels(groups)-1))-1
z <- as.numeric(z)+f*d
panel.segplot(x, y, z, centers=centers,
subscripts=subscripts, ...)
panel.text(x=centers, y=z, labels=txt, pos=c(3,1))
})
