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_milleflora.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 18 9 20 20 18 22 20 26 6 32 ...
## $ ivg : num 0.799 0.367 1.071 1.094 0.403 ...
##-----------------------------------------------------------------------------
## 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.14473 0.38158 22.5933 1.211e-08 ***
## temper 1 0.24249 0.24249 14.3578 0.0005115 ***
## embal 1 0.00055 0.00055 0.0324 0.8579980
## armaz:temper 2 0.28461 0.14230 8.4258 0.0009090 ***
## armaz:embal 2 0.08945 0.04472 2.6481 0.0834889 .
## temper:embal 1 0.00130 0.00130 0.0771 0.7827934
## armaz:temper:embal 2 0.06753 0.03377 1.9993 0.1490762
## Residuals 39 0.65867 0.01689
## ---
## 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.65867
## 2 45 0.81749 -6 -0.15883 1.5674 0.1826
anova(m1, test="F")
## Analysis of Variance Table
##
## Response: sqrt(ivg)
## Df Sum Sq Mean Sq F value Pr(>F)
## armaz 3 1.14473 0.38158 21.0044 1.173e-08 ***
## temper 1 0.24249 0.24249 13.3480 0.0006728 ***
## armaz:temper 2 0.28461 0.14230 7.8332 0.0012048 **
## Residuals 45 0.81749 0.01817
## ---
## 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.8951220 0.7593882 1.0308557
## 5 2 5 0.5929148 0.4969365 0.6888930
## 13 2 18 0.6236226 0.5276444 0.7196009
## 21 4 5 0.4657552 0.3697769 0.5617334
## 29 4 18 0.3519125 0.2559342 0.4478907
## 37 6 5 0.8864148 0.7904366 0.9823931
## 45 6 18 0.5430924 0.4471142 0.6390707
## 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 0.895 0.593 0.624 0.466 0.352 ...
## $ lwr : num 0.759 0.497 0.528 0.37 0.256 ...
## $ upr : num 1.031 0.689 0.72 0.562 0.448 ...
## $ pEstimate: num 0.801 0.352 0.389 0.217 0.124 ...
## $ plwr : num 0.5767 0.2469 0.2784 0.1367 0.0655 ...
## $ pupr : num 1.063 0.475 0.518 0.316 0.201 ...
## $ 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 0.8012433 0.57667045 1.0626635 0 meses d
## 5 2 5 0.3515479 0.24694588 0.4745736 2:5 bc
## 13 2 18 0.3889052 0.27840859 0.5178255 2:18 b
## 21 4 5 0.2169279 0.13673497 0.3155445 4:5 ac
## 29 4 18 0.1238424 0.06550232 0.2006061 4:18 a
## 37 6 5 0.7857312 0.62478994 0.9650962 6:5 d
## 45 6 18 0.2949494 0.19991109 0.4084114 6:18 bc
##----------------------------------------------------------------------
## 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.~milleflora)),
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] 3.066968
## 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 291.75
## armaz 3 91.983 48 199.77 9.9972 5.096e-05 ***
## temper 1 30.289 47 169.48 9.8757 0.003194 **
## embal 1 0.270 46 169.21 0.0881 0.768235
## armaz:temper 2 23.847 44 145.36 3.8877 0.028869 *
## armaz:embal 2 14.476 42 130.88 2.3600 0.107766
## temper:embal 1 0.075 41 130.81 0.0246 0.876135
## armaz:temper:embal 2 10.849 39 119.96 1.7687 0.183968
## ---
## 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 119.96
## 2 45 145.63 -6 -25.673 1.3951 0.2411
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 291.75
## armaz 3 91.983 48 199.77 9.7047 4.688e-05 ***
## temper 1 30.289 47 169.48 9.5868 0.003368 **
## armaz:temper 2 23.844 45 145.63 3.7735 0.030545 *
## ---
## 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.6856565 -1.2075749 -0.1637381
## 5 2 5 -0.5537276 -0.9155435 -0.1919118
## 13 2 18 -0.5537276 -0.9155435 -0.1919118
## 21 4 5 -1.3553321 -1.7868118 -0.9238525
## 29 4 18 -2.0163205 -2.5572547 -1.4753862
## 37 6 5 -0.4895482 -0.8484150 -0.1306815
## 45 6 18 -1.5505974 -2.0090292 -1.0921656
## 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.686 -0.554 -0.554 -1.355 -2.016 ...
## $ lwr : num -1.208 -0.916 -0.916 -1.787 -2.557 ...
## $ upr : num -0.164 -0.192 -0.192 -0.924 -1.475 ...
## $ pEstimate: num 33.5 36.5 36.5 20.5 11.8 ...
## $ plwr : num 23.01 28.59 28.59 14.35 7.19 ...
## $ pupr : num 45.9 45.2 45.2 28.4 18.6 ...
## $ 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 33.50 23.013042 45.91567 0 meses bc
## 5 2 5 36.50 28.586681 45.21688 2:5 c
## 13 2 18 36.50 28.586681 45.21688 2:18 c
## 21 4 5 20.50 14.346405 28.41736 4:5 ab
## 29 4 18 11.75 7.194062 18.61253 4:18 a
## 37 6 5 38.00 29.976546 46.73760 6:5 c
## 45 6 18 17.50 11.825817 25.12107 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.~milleflora)),
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))
})
