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.
## list.files(pattern="*.txt")
da <- read.table("germini.txt", header=TRUE, sep="\t")
da$temper <- factor(da$temper, levels=c("20", "25", "30", "20/30"))
str(da)
## 'data.frame': 32 obs. of 5 variables:
## $ especie: Factor w/ 2 levels "B. milleflora",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ temper : Factor w/ 4 levels "20","25","30",..: 4 4 4 4 1 1 1 1 2 2 ...
## $ rept : int 1 2 3 4 1 2 3 4 1 2 ...
## $ ngerm50: int 18 26 18 24 24 33 31 30 34 33 ...
## $ ivg : num 1.21 2.29 1.68 2.19 1.75 ...
##----------------------------------------------------------------------
## Ver.
## xtabs(~temper+especie, da)
xtabs(~especie+temper, da)
## temper
## especie 20 25 30 20/30
## B. milleflora 4 4 4 4
## B. tridentata 4 4 4 4
xyplot(ngerm50/50~temper|especie, data=da, type=c("p","a"))
xyplot(ivg~temper|especie, data=da, type=c("p","a"))
##----------------------------------------------------------------------
## B. milleflora.
m0 <- glm(cbind(ngerm50, 50-ngerm50)~temper,
data=subset(da, especie=="B. milleflora"),
family="quasibinomial")
## par(mfrow=c(2,2)); plot(m0); layout(1)
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 15 191.696
## temper 3 171.7 12 19.995 38.147 2.054e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m0)
##
## Call:
## glm(formula = cbind(ngerm50, 50 - ngerm50) ~ temper, family = "quasibinomial",
## data = subset(da, especie == "B. milleflora"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.44935 -0.70799 -0.00012 0.53617 1.72562
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.293 1.228 -4.311 0.00101 **
## temper25 -17.337 4311.541 -0.004 0.99686
## temper30 1.817 1.329 1.368 0.19652
## temper20/30 4.608 1.242 3.711 0.00298 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 1.500353)
##
## Null deviance: 191.696 on 15 degrees of freedom
## Residual deviance: 19.995 on 12 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 18
##----------------------------------------------------------------------
## Intervalos de confiança para a germianação.
m1 <- update(m0, .~0+temper)
ci <- confint(m1)
## Waiting for profiling to be done...
colnames(ci) <- c("lwr", "upr")
ci <- cbind(fit=coef(m1), ci)
grid <- equallevels(cbind(data.frame(temper=levels(da$temper)), ci), da)
grid[,c("fit","lwr","upr")] <-
sapply(grid[,c("fit","lwr","upr")], m0$family$linkinv)
## colwise(.fun=round, digits=4, .cols=is.numeric)(grid)
##----------------------------------------------------------------------
## Comparações de germinação.
## **NÃO** É UM TESTE DE TUKEY. São contrastes 2 a 2 com p-valor do
## teste corrigido por bonferroni.
g0 <- summary(glht(m0, linfct=mcp(temper="Tukey")),
test=adjusted(type="bonferroni"))
g0
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glm(formula = cbind(ngerm50, 50 - ngerm50) ~ temper, family = "quasibinomial",
## data = subset(da, especie == "B. milleflora"))
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 25 - 20 == 0 -17.3374 4311.5406 -0.004 1.00000
## 30 - 20 == 0 1.8172 1.3288 1.368 1.00000
## 20/30 - 20 == 0 4.6076 1.2416 3.711 0.00124 **
## 30 - 25 == 0 19.1546 4311.5404 0.004 1.00000
## 20/30 - 25 == 0 21.9451 4311.5404 0.005 1.00000
## 20/30 - 30 == 0 2.7904 0.5399 5.169 1.41e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)
grid$cld <- cld(g0)$mcletters$Letters
grid
## temper fit lwr upr cld
## temper20 20 5.000000e-03 1.054221e-04 2.770487e-02 b
## temper25 25 1.484581e-10 2.220446e-16 2.220446e-16 ab
## temper30 30 3.000000e-02 9.407174e-03 6.830164e-02 b
## temper20/30 20/30 3.350000e-01 2.585579e-01 4.177213e-01 a
##----------------------------------------------------------------------
## Gráfico.
segplot(temper~lwr+upr, centers=fit, data=grid, draw=FALSE,
xlab=expression(Germinação~de~italic(B.~milleflora)),
ylab=expression(Temperatura~(degree*C)),
txt=with(grid, sprintf("%0.3f %s", fit, cld)),
panel=function(z, centers, txt, ...){
panel.segplot(z=z, centers=centers, ...)
panel.text(x=centers, y=z, labels=txt, pos=3)
})
##----------------------------------------------------------------------
## B. tridentata.
m0 <- glm(cbind(ngerm50, 50-ngerm50)~temper,
data=subset(da, especie=="B. tridentata"),
family="quasibinomial")
## par(mfrow=c(2,2)); plot(m0); layout(1)
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 15 32.646
## temper 3 17.024 12 15.622 4.3481 0.02721 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m0)
##
## Call:
## glm(formula = cbind(ngerm50, 50 - ngerm50) ~ temper, family = "quasibinomial",
## data = subset(da, especie == "B. tridentata"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1336 -0.5769 0.1446 0.7878 1.3127
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.640e-01 1.642e-01 2.216 0.0468 *
## temper25 -1.460e-12 2.323e-01 0.000 1.0000
## temper30 8.335e-02 2.332e-01 0.357 0.7270
## temper20/30 -6.458e-01 2.315e-01 -2.790 0.0164 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 1.305073)
##
## Null deviance: 32.646 on 15 degrees of freedom
## Residual deviance: 15.622 on 12 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 3
##----------------------------------------------------------------------
## Intervalos de confiança para a germianação.
m1 <- update(m0, .~0+temper)
ci <- confint(m1)
## Waiting for profiling to be done...
colnames(ci) <- c("lwr", "upr")
ci <- cbind(fit=coef(m1), ci)
grid <- equallevels(cbind(data.frame(temper=levels(da$temper)), ci), da)
grid[,c("fit","lwr","upr")] <-
sapply(grid[,c("fit","lwr","upr")], m0$family$linkinv)
## colwise(.fun=round, digits=4, .cols=is.numeric)(grid)
##----------------------------------------------------------------------
## Comparações de germinação.
## **NÃO** É UM TESTE DE TUKEY. São contrastes 2 a 2 com p-valor do
## teste corrigido por bonferroni.
g0 <- summary(glht(m0, linfct=mcp(temper="Tukey")),
test=adjusted(type="bonferroni"))
g0
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glm(formula = cbind(ngerm50, 50 - ngerm50) ~ temper, family = "quasibinomial",
## data = subset(da, especie == "B. tridentata"))
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 25 - 20 == 0 -1.460e-12 2.323e-01 0.000 1.0000
## 30 - 20 == 0 8.335e-02 2.332e-01 0.357 1.0000
## 20/30 - 20 == 0 -6.458e-01 2.315e-01 -2.790 0.0317 *
## 30 - 25 == 0 8.335e-02 2.332e-01 0.357 1.0000
## 20/30 - 25 == 0 -6.458e-01 2.315e-01 -2.790 0.0317 *
## 20/30 - 30 == 0 -7.292e-01 2.325e-01 -3.136 0.0103 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)
grid$cld <- cld(g0)$mcletters$Letters
grid
## temper fit lwr upr cld
## temper20 20 0.59 0.5111164 0.6659010 b
## temper25 25 0.59 0.5111164 0.6659010 b
## temper30 30 0.61 0.5314239 0.6849307 b
## temper20/30 20/30 0.43 0.3532603 0.5090596 a
##----------------------------------------------------------------------
## Gráfico.
segplot(temper~lwr+upr, centers=fit, data=grid, draw=FALSE,
xlab=expression(Germinação~de~italic(B.~tridentata)),
ylab=expression(Temperatura~(degree*C)),
txt=with(grid, sprintf("%0.3f %s", fit, cld)),
panel=function(z, centers, txt, ...){
panel.segplot(z=z, centers=centers, ...)
panel.text(x=centers, y=z, labels=txt, pos=3)
})
##----------------------------------------------------------------------
## B. milleflora.
## Análise da raíz quadrada melhor atende os pressupostos.
m0 <- lm(sqrt(ivg)~temper, data=subset(da, especie=="B. milleflora"))
## par(mfrow=c(2,2)); plot(m0); layout(1)
anova(m0)
## Analysis of Variance Table
##
## Response: sqrt(ivg)
## Df Sum Sq Mean Sq F value Pr(>F)
## temper 3 2.01977 0.67326 20.155 5.601e-05 ***
## Residuals 12 0.40085 0.03340
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m0)
##
## Call:
## lm(formula = sqrt(ivg) ~ temper, data = subset(da, especie ==
## "B. milleflora"))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.28923 -0.07215 0.00000 0.07659 0.35417
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.07215 0.09138 0.790 0.445
## temper25 -0.07215 0.12924 -0.558 0.587
## temper30 0.13270 0.12924 1.027 0.325
## temper20/30 0.82297 0.12924 6.368 3.57e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1828 on 12 degrees of freedom
## Multiple R-squared: 0.8344, Adjusted R-squared: 0.793
## F-statistic: 20.16 on 3 and 12 DF, p-value: 5.601e-05
##----------------------------------------------------------------------
m1 <- update(m0, .~0+temper)
ci <- confint(m1)
colnames(ci) <- c("lwr", "upr")
ci <- cbind(fit=coef(m1), ci)
grid <- equallevels(cbind(data.frame(temper=levels(da$temper)), ci), da)
grid[,c("fit","lwr","upr")] <-
sapply(grid[,c("fit","lwr","upr")],
function(x) x^2)
## colwise(.fun=round, digits=4, .cols=is.numeric)(grid)
##----------------------------------------------------------------------
## Comparações de germinação.
## **NÃO** É UM TESTE DE TUKEY. São contrastes 2 a 2 com p-valor do
## teste corrigido por bonferroni.
g0 <- summary(glht(m0, linfct=mcp(temper="Tukey")),
test=adjusted(type="bonferroni"))
g0
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = sqrt(ivg) ~ temper, data = subset(da, especie ==
## "B. milleflora"))
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 25 - 20 == 0 -0.07215 0.12924 -0.558 1.000000
## 30 - 20 == 0 0.13270 0.12924 1.027 1.000000
## 20/30 - 20 == 0 0.82297 0.12924 6.368 0.000214 ***
## 30 - 25 == 0 0.20485 0.12924 1.585 0.833592
## 20/30 - 25 == 0 0.89512 0.12924 6.926 9.55e-05 ***
## 20/30 - 30 == 0 0.69027 0.12924 5.341 0.001057 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)
grid$cld <- cld(g0)$mcletters$Letters
grid
## temper fit lwr upr cld
## temper20 20 0.00520625 1.611722e-02 0.07358324 a
## temper25 25 0.00000000 3.964398e-02 0.03964398 a
## temper30 30 0.04196366 3.297483e-05 0.16318230 a
## temper20/30 20/30 0.80124334 4.844355e-01 1.19733914 b
##----------------------------------------------------------------------
## Gráfico.
segplot(temper~lwr+upr, centers=fit, data=grid, draw=FALSE,
xlab=expression(IVG~de~italic(B.~milleflora)),
ylab=expression(Temperatura~(degree*C)),
txt=with(grid, sprintf("%0.3f %s", fit, cld)),
panel=function(z, centers, txt, ...){
panel.segplot(z=z, centers=centers, ...)
panel.text(x=centers, y=z, labels=txt, pos=3)
})
##----------------------------------------------------------------------
## B. tridentata.
## Análise da raíz quadrada melhor atende os pressupostos.
m0 <- lm(sqrt(ivg)~temper, data=subset(da, especie=="B. tridentata"))
## par(mfrow=c(2,2)); plot(m0); layout(1)
anova(m0)
## Analysis of Variance Table
##
## Response: sqrt(ivg)
## Df Sum Sq Mean Sq F value Pr(>F)
## temper 3 0.25722 0.085739 4.3235 0.02767 *
## Residuals 12 0.23797 0.019831
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m0)
##
## Call:
## lm(formula = sqrt(ivg) ~ temper, data = subset(da, especie ==
## "B. tridentata"))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.24540 -0.06382 0.03313 0.07554 0.19912
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.46994 0.07041 20.877 8.45e-11 ***
## temper25 0.17306 0.09958 1.738 0.1078
## temper30 0.18004 0.09958 1.808 0.0957 .
## temper20/30 -0.12313 0.09958 -1.237 0.2399
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1408 on 12 degrees of freedom
## Multiple R-squared: 0.5194, Adjusted R-squared: 0.3993
## F-statistic: 4.323 on 3 and 12 DF, p-value: 0.02767
##----------------------------------------------------------------------
m1 <- update(m0, .~0+temper)
ci <- confint(m1)
colnames(ci) <- c("lwr", "upr")
ci <- cbind(fit=coef(m1), ci)
grid <- equallevels(cbind(data.frame(temper=levels(da$temper)), ci), da)
grid[,c("fit","lwr","upr")] <-
sapply(grid[,c("fit","lwr","upr")],
function(x) x^2)
## colwise(.fun=round, digits=4, .cols=is.numeric)(grid)
##----------------------------------------------------------------------
## Comparações de germinação.
## **NÃO** É UM TESTE DE TUKEY. São contrastes 2 a 2 com p-valor do
## teste corrigido por bonferroni.
g0 <- summary(glht(m0, linfct=mcp(temper="Tukey")),
test=adjusted(type="bonferroni"))
g0
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = sqrt(ivg) ~ temper, data = subset(da, especie ==
## "B. tridentata"))
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 25 - 20 == 0 0.173057 0.099576 1.738 0.6467
## 30 - 20 == 0 0.180043 0.099576 1.808 0.5742
## 20/30 - 20 == 0 -0.123131 0.099576 -1.237 1.0000
## 30 - 25 == 0 0.006986 0.099576 0.070 1.0000
## 20/30 - 25 == 0 -0.296188 0.099576 -2.974 0.0696 .
## 20/30 - 30 == 0 -0.303174 0.099576 -3.045 0.0611 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)
grid$cld <- cld(g0)$mcletters$Letters
grid
## temper fit lwr upr cld
## temper20 20 2.160719 1.733240 2.635269 a
## temper25 25 2.699434 2.218857 3.227082 a
## temper30 30 2.722439 2.239719 3.252231 a
## temper20/30 20/30 1.813890 1.424191 2.250660 a
##----------------------------------------------------------------------
## Gráfico.
segplot(temper~lwr+upr, centers=fit, data=grid, draw=FALSE,
xlab=expression(IVG~de~italic(B.~milleflora)),
ylab=expression(Temperatura~(degree*C)),
txt=with(grid, sprintf("%0.3f %s", fit, cld)),
panel=function(z, centers, txt, ...){
panel.segplot(z=z, centers=centers, ...)
panel.text(x=centers, y=z, labels=txt, pos=3)
})
