Material suplementar do artigo
Germination and Vigor of caiaué seeds: American oil palm tree
Download the dataset: caiuae.txt
Download the R script: caiuae.R
Definições da sessão
#-----------------------------------------------------------------------
# Required packages and session definition.
library(lattice)
library(latticeExtra)
library(doBy)
library(multcomp)
library(plyr)
library(wzRfun)
trellis.par.set(canonical.theme(color = FALSE))Importação dos dados para análise
#-----------------------------------------------------------------------
# Importação.
da <- read.table(file = "caiuae.txt",
header = TRUE,
sep = "\t",
encoding = "latin1")
str(da)'data.frame': 45 obs. of 14 variables:
$ temp : int 55 55 55 55 55 55 55 55 55 55 ...
$ umid : Factor w/ 5 levels "18-19","19-20",..: 1 1 1 2 2 2 3 3 3 4 ...
$ rept : int 1 2 3 1 2 3 1 2 3 1 ...
$ X2012.11.22: int 6 5 0 21 12 22 51 84 78 23 ...
$ X2012.11.30: int 30 19 12 71 58 54 86 125 120 78 ...
$ X2012.12.07: int 30 40 21 81 41 39 63 63 48 42 ...
$ X2012.12.14: int 29 41 50 51 60 64 50 54 69 30 ...
$ X2012.12.21: int 12 21 34 17 23 42 22 18 16 15 ...
$ germini : int 36 24 12 92 70 76 137 209 198 101 ...
$ germfin : int 107 126 117 241 194 221 272 344 331 188 ...
$ ngermfin : int 203 224 263 159 175 167 183 96 90 205 ...
$ pgerm : num 34.52 36 30.79 6.03 5.26 ...
$ pcont : num 11.61 6.86 3.16 23 18.97 ...
$ ivg : num 5.81 6.04 4.61 14.24 10.61 ...#-----------------------------------------------------------------------
# Propriedades.
# Todos casos completos (sem dados perdidos).
all(complete.cases(da))[1] TRUE# Frequência.
xtabs(~temp + umid, data = da) umid
temp 18-19 19-20 20-21 21-22 22-23
55 3 3 3 3 3
75 3 3 3 3 3
100 3 3 3 3 3Visualização e análise exploratória dos dados
#-----------------------------------------------------------------------
# Germinação.
leg <- list("Heat treatment period (days)",
"Moisture content (%)")
xyplot(germfin/(germfin + ngermfin) ~ temp,
groups = umid, data = da,
type = c("p", "a"),
auto.key = list(columns = 2, title = leg[[2]],
cex.title = 1.1),
xlab = leg[[1]],
ylab = "Germination")
xyplot(germfin/(germfin + ngermfin) ~ umid,
groups = temp, data = da,
type = c("p", "a"),
auto.key = list(columns = 3, title = leg[[1]],
cex.title = 1.1),
xlab = leg[[2]],
ylab = "Germination")
#-----------------------------------------------------------------------
# Primeira germinação (aos 15 dias).
xyplot(pcont ~ temp, groups = umid, data = da,
type = c("p", "a"),
auto.key = list(columns = 2, title = leg[[2]], cex.title = 1.1),
xlab = leg[[1]],
ylab = "Germination in the first count (%)")
xyplot(pcont ~ umid, groups = temp, data = da,
type = c("p", "a"),
auto.key = list(columns = 3, title = leg[[1]], cex.title = 1.1),
xlab = leg[[2]], ylab = "Germination in the first count (%)")
#-----------------------------------------------------------------------
# IVG.
xyplot(ivg ~ umid, groups = temp, data = da,
type = c("p", "a"),
auto.key = list(columns = 3, title = leg[[1]], cex.title = 1.1),
xlab = leg[[2]], ylab = "Germination speed index")
xyplot(ivg ~ temp, groups = umid, data = da,
type = c("p", "a"),
auto.key = list(columns = 3, title = leg[[1]], cex.title = 1.1),
xlab = leg[[2]], ylab = "Germination speed index")
Germinação final (36 dias)
Para análise da germinação final (15+3 * 7 = 36 dias após estÃmulo) considerou distribuição binomial sendo a germinação de uma semente o sucesso. A função de ligação foi a logit. Para essa análise é necessário o par de variáveis: número de sementes germinadas e número de sementes não germinadas (\(y\), \(n-y\)). Havendo alguma dispersão extra à quela prevista pela distribuição binomial, declarou-se o modelo quasi-binomial para correção dos erros padrões dos efeitos estimados, bem como o uso da estatÃstica no quadro de análise de deviance para os termos do modelo estatÃstico (perÃodo, umidade e periodo*umidade). Tal análise foi feita por meio da função glm() do programa R para computação estatÃstica.
Representou-se o efeito de intervalos de umidade como uma variável categórica, ao invés de contÃnua visto que 1) a umidade não é precisamente fixada sob um valor mas oscila dentro do intervalo e 2) os aparelhos permitem regular apenas nestes 5 nÃveis, então qualquer interpolação númericamente estimada não se pode aplicar na prática, o que faz a decisão por um desses nÃveis algo mais simples (Wanderlei, rever esse argumento).
#-----------------------------------------------------------------------
# Modelo.
da$Temp <- as.factor(da$temp)
m0 <- glm(cbind(germfin, ngermfin) ~ Temp * umid,
data = da, family = quasibinomial)
deviance(m0)/df.residual(m0)[1] 14.44564par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Testes hipóteses sequenciais pela estatÃstica F de Wald.
anova(m0, test = "F")Analysis of Deviance Table
Model: quasibinomial, link: logit
Response: cbind(germfin, ngermfin)
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev F Pr(>F)
NULL 44 5162.4
Temp 2 143.34 42 5019.1 5.0332 0.01304 *
umid 4 3008.32 38 2010.7 52.8182 3.711e-13 ***
Temp:umid 8 1577.38 30 433.4 13.8473 3.552e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Testes hipóteses marginais pela estatÃstica F da razão de
# verossimilhanças.
drop1(m0, test = "F", scope = . ~ .)Single term deletions
Model:
cbind(germfin, ngermfin) ~ Temp * umid
Df Deviance F value Pr(>F)
<none> 433.37
Temp 2 941.16 17.5760 8.869e-06 ***
umid 4 872.17 7.5941 0.0002375 ***
Temp:umid 8 2010.75 13.6492 4.173e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1O modelo quasi-binomial foi ajustado devido à variação extra binomial, o que pode ser justificado pela não uniformidade das sementes quanto ao tamanho e demais caracterÃsticas, o que faz com que apresentem um potencial germinativo mais variável que o previsto pela binomial que assume probabilidade de sucesso constante em todas as unidades experimentais.
O quadro acima indica que o efeito da interação foi significativo pela estatÃstica F. Trata-se de um quadro de deviance, uma vez o modelo quasi-binomial foi considerado para a germinação de sementes.
O estudo da interação (ou de efeitos principais, quando o caso) foi feito por comparações múltiplas baseadas em contrastes dois à dois (contrastes de Tukey, all pairwise comparisons). Métodos que se baseiam em diferença mÃnima sigficativa entre médias, como Tukey, não podem ser aplicados pois os contrastes de médias estimadas (probabilidade de germinação) nesse modelo tem precisão diferente, mas o número de repetições sendo iguais. Sendo assim, aplicou-se o teste t para cada contraste entre médias e corrigiu-se o p-valor pelo método fdr (false discovery rate). O estudo da interação preconizou o efeito da faixa de umidade e cada perÃodo de permanência.
Citar esse livro sobre comparações múltiplas: http://www.ievbras.ru/ecostat/Kiril/R/Biblio/R_eng/Bretz%20Multiple%20Comparisons.pdf, https://books.google.com.br/books?id=U8Xc9zujgcsC&source=gbs_navlinks_s;
#-----------------------------------------------------------------------
# Comparações múltiplas de umidade dentro de perÃodo.
L <- LSmatrix(m0, effect = c("Temp", "umid"))
grid <- equallevels(attr(L, "grid"), da)
Lu <- by(data = L, INDICES = grid$Temp, FUN = as.matrix)
Lu <- lapply(Lu, FUN = "rownames<-", levels(grid$umid))
Lu <- ldply(lapply(Lu, apmc, model = m0, test = "fdr", focus = "umid"),
.id = "Temp")
Lu[, c("fit", "lwr", "upr")] <-
apply(Lu[, c("fit", "lwr", "upr")],
MARGIN = 2,
FUN = m0$family$linkinv)
segplot(umid ~ lwr + upr | Temp, centers = fit,
data = Lu, horizontal = FALSE, draw = FALSE,
as.table = TRUE, layout = c(NA, 1),
txt = Lu$cld,
ylab = "Germination",
xlab = leg[[2]],
panel = function(z, centers, subscripts, txt, ...) {
panel.segplot(z = z, centers = centers,
subscripts = subscripts, ...)
panel.text(y = centers[subscripts],
x = as.integer(z[subscripts]), pos = 2,
labels = txt[subscripts])
})
#-----------------------------------------------------------------------
# Comparações múltiplas dentro de umidade.
Lt <- by(data = L, INDICES = grid$umid, FUN = as.matrix)
Lt <- lapply(Lt, FUN = "rownames<-", levels(grid$Temp))
Lt <- ldply(lapply(Lt, apmc, model = m0, test = "fdr", focus = "Temp"),
.id = "umid")
Lt[, c("fit", "lwr", "upr")] <-
apply(Lt[, c("fit", "lwr", "upr")],
MARGIN = 2, FUN = m0$family$linkinv)
Lt$Temp <- factor(Lt$Temp, levels = levels(da$Temp))
segplot(Temp ~ lwr + upr | umid, centers = fit,
data = Lt, horizontal = FALSE,
draw = FALSE, layout = c(NA, 1), as.table = TRUE,
txt = Lt$cld,
ylab = "Germination", xlab = leg[[1]],
panel = function(z, centers, subscripts, txt, ...) {
panel.segplot(z = z, centers = centers,
subscripts = subscripts, ...)
panel.text(y = centers[subscripts],
x = as.integer(z[subscripts]), pos = 4,
labels = txt[subscripts])
})
Nestes gráficos tem-se as estimativas pontuais de germinação (pontos preenchidos) junto ao intervalo de confiança de 95% para tais estimativas. Os intervalos são assimétricos devido emprego da função inversa da função de ligação logit à s estimativas na escala do preditor linear para determinar as probabilidades estimadas. Letras diferentes ao lado das estimativas indicam contrastes não nulos (à 5%) entre parâmetros para um nÃvel fixo do fator indicado nas tarjas.
Germinação inicial (aos 15 dias)
Assim como para a germinação final, para a germinação inicial considerou-se distribuição quasi binomial também, pelos mesmos motivos.
#-----------------------------------------------------------------------
# Modelo.
# aggregate(germini~Temp+umid, data=da, FUN=sum)
# subset(da, Temp=="100" & umid=="18-19")
da$germini[31] <- 1
m0 <- glm(cbind(germini, ngermfin) ~ Temp * umid,
data = da,
family = quasibinomial)
deviance(m0)/df.residual(m0)[1] 11.53034par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Testes hipóteses sequenciais pela estatÃstica F de Wald.
anova(m0, test = "F")Analysis of Deviance Table
Model: quasibinomial, link: logit
Response: cbind(germini, ngermfin)
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev F Pr(>F)
NULL 44 4907.6
Temp 2 292.6 42 4615.1 13.329 7.207e-05 ***
umid 4 3252.6 38 1362.5 74.096 4.129e-15 ***
Temp:umid 8 1016.6 30 345.9 11.579 2.514e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Testes hipóteses marginais pela estatÃstica F da razão de
# verossimilhanças.
drop1(m0, test = "F", scope = . ~ .)Single term deletions
Model:
cbind(germini, ngermfin) ~ Temp * umid
Df Deviance F value Pr(>F)
<none> 345.91
Temp 2 464.02 5.1217 0.0122 *
umid 4 893.64 11.8759 6.689e-06 ***
Temp:umid 8 1362.49 11.0207 4.241e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#-----------------------------------------------------------------------
# Comparações múltiplas de umidade dentro de perÃodo.
L <- LSmatrix(m0, effect = c("Temp", "umid"))
grid <- equallevels(attr(L, "grid"), da)
Lu <- by(data = L, INDICES = grid$Temp, FUN = as.matrix)
Lu <- lapply(Lu, FUN = "rownames<-", levels(grid$umid))
Lu <- ldply(lapply(Lu, apmc, model = m0, test = "fdr", focus = "umid"),
.id = "Temp")
Lu[, c("fit", "lwr", "upr")] <-
apply(Lu[, c("fit", "lwr", "upr")],
MARGIN = 2, FUN = m0$family$linkinv)
segplot(umid ~ lwr + upr | Temp,
centers = fit, data = Lu,
horizontal = FALSE, draw = FALSE,
layout = c(NA, 1), as.table = TRUE,
txt = Lu$cld,
ylab = "Germination in the first count",
xlab = leg[[2]],
panel = function(z, centers, subscripts, txt, ...) {
panel.segplot(z = z, centers = centers,
subscripts = subscripts, ...)
panel.text(y = centers[subscripts],
x = as.integer(z[subscripts]), pos = 2,
labels = txt[subscripts])
})
#-----------------------------------------------------------------------
# Comparações múltiplas dentro de umidade.
Lt <- by(data = L, INDICES = grid$umid, FUN = as.matrix)
Lt <- lapply(Lt, FUN = "rownames<-", levels(grid$Temp))
Lt <- ldply(lapply(Lt, apmc, model = m0, test = "fdr", focus = "Temp"),
.id = "umid")
Lt[, c("fit", "lwr", "upr")] <-
apply(Lt[, c("fit", "lwr", "upr")],
MARGIN = 2, FUN = m0$family$linkinv)
Lt$Temp <- factor(Lt$Temp, levels = levels(da$Temp))
segplot(Temp ~ lwr + upr | umid,
centers = fit, data = Lt, horizontal = FALSE,
draw = FALSE, layout = c(NA, 1), as.table = TRUE,
txt = Lt$cld,
ylab = "Germination in the first count",
xlab = leg[[1]],
panel = function(z, centers, subscripts, txt, ...) {
panel.segplot(z = z, centers = centers,
subscripts = subscripts, ...)
panel.text(y = centers[subscripts],
x = as.integer(z[subscripts]), pos = 4,
labels = txt[subscripts])
})
A germinação inicial se aplicam as mesmas interpretações para os gráficos. Vale salientar que a cela perÃodo de 100 dias e umidade 18-19 não apresentou germinação em nenhuma repetição. Reescreveu-se umas das repetições como tendo 1 semente germinada para evitar problemas de estimação, visto que a probabilidade de sucesso tem qus ser maior que zero (\(p > 0\)).
IVG - Ãndice de velocidade de germinação
O indice de velocidade de gerimanação, embora seja uma função de variáveis discretas, assumiu-se distribuição normal (mesmo porque trata-se de uma soma de termos). No entanto, se a análise gráfica dos resÃduos indicar fuga dos pressupostos, uma transformação da famÃlia Box-Cox será aplicada para fazer com que os pressupostos, na escala transformada, sejam atendidos. Todas as inferências serão feitos na escala na qual os pressupostos são atendidos e para facilidade de exposição, os valores médios serão transformados pra escala original da variável.
#-----------------------------------------------------------------------
# Model.
m0 <- lm(ivg + 0.5 ~ Temp * umid, data = da)
par(mfrow = c(2, 2))
plot(m0)
layout(1)
MASS::boxcox(m0)
abline(v = 0.5, col = "red")
m0 <- lm(sqrt(ivg + 0.5) ~ Temp * umid, data = da)
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Testes hipóteses sequenciais pela estatÃstica F de Wald.
anova(m0)Analysis of Variance Table
Response: sqrt(ivg + 0.5)
Df Sum Sq Mean Sq F value Pr(>F)
Temp 2 5.777 2.8883 23.338 7.706e-07 ***
umid 4 60.934 15.2336 123.093 < 2.2e-16 ***
Temp:umid 8 14.661 1.8327 14.809 1.662e-08 ***
Residuals 30 3.713 0.1238
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Para atender os pressupostos, analisou-se a raÃz quadrada do IVG mais 0.5, para tornar a variável positiva. Nessa escala, tem-se interação entre perÃodo e umidade. Para comparações múltiplas, considerou-se o ts etet correção do p-valor pelo método fdr.
#-----------------------------------------------------------------------
# Comparações múltiplas
L <- LSmatrix(m0, effect = c("Temp", "umid"))
grid <- equallevels(attr(L, "grid"), da)
Lu <- by(data = L, INDICES = grid$Temp, FUN = as.matrix)
Lu <- lapply(Lu, FUN = "rownames<-", levels(grid$umid))
Lu <- ldply(lapply(Lu, apmc, model = m0, test = "fdr", focus = "umid"),
.id = "Temp")
Lu[, c("fit", "lwr", "upr")] <-
apply(Lu[, c("fit", "lwr", "upr")],
MARGIN = 2, FUN = function(x) x^2 - 0.5)
segplot(umid ~ lwr + upr | Temp, centers = fit, data = Lu,
horizontal = FALSE,
draw = FALSE, layout = c(NA, 1),
txt = Lu$cld, as.table = TRUE,
ylab = "Germination speed index",
xlab = leg[[2]],
panel = function(z, centers, subscripts, txt, ...) {
panel.segplot(z = z, centers = centers,
subscripts = subscripts, ...)
panel.text(y = centers[subscripts],
x = as.integer(z[subscripts]), pos = 2,
labels = txt[subscripts])
})
Lt <- by(data = L, INDICES = grid$umid, FUN = as.matrix)
Lt <- lapply(Lt, FUN = "rownames<-", levels(grid$Temp))
Lt <- ldply(lapply(Lt, apmc, model = m0, test = "fdr", focus = "Temp"),
.id = "umid")
Lt[, c("fit", "lwr", "upr")] <-
apply(Lt[, c("fit", "lwr", "upr")],
MARGIN = 2, FUN = function(x) x^2 - 0.5)
Lt$Temp <- factor(Lt$Temp, levels = levels(da$Temp))
segplot(Temp ~ lwr + upr | umid,
centers = fit, data = Lt, horizontal = FALSE,
draw = FALSE, layout = c(NA, 1), as.table = TRUE,
txt = Lt$cld,
ylab = "Germination speed index",
xlab = leg[[1]],
panel = function(z, centers, subscripts, txt, ...) {
panel.segplot(z = z, centers = centers,
subscripts = subscripts, ...)
panel.text(y = centers[subscripts],
x = as.integer(z[subscripts]), pos = 4,
labels = txt[subscripts])
})
Detalhes da sessão
options(width=90)
devtools::session_info()
Sys.time() setting value
version R version 3.3.0 (2016-05-03)
system x86_64, linux-gnu
ui X11
language en_US
collate en_US.UTF-8
tz <NA>
date 2016-06-23
package * version date source
codetools 0.2-14 2015-07-15 CRAN (R 3.2.1)
devtools 1.11.1 2016-04-21 CRAN (R 3.3.0)
digest 0.6.9 2016-01-08 CRAN (R 3.3.0)
doBy * 4.5-15 2016-03-31 CRAN (R 3.3.0)
evaluate 0.9 2016-04-29 CRAN (R 3.3.0)
formatR 1.4 2016-05-09 CRAN (R 3.3.0)
htmltools 0.3.5 2016-03-21 CRAN (R 3.3.0)
knitr * 1.13 2016-05-09 CRAN (R 3.3.0)
lattice * 0.20-33 2015-07-14 CRAN (R 3.2.1)
latticeExtra * 0.6-28 2016-02-09 CRAN (R 3.3.0)
magrittr 1.5 2014-11-22 CRAN (R 3.3.0)
MASS * 7.3-45 2015-11-10 CRAN (R 3.2.5)
Matrix 1.2-6 2016-05-02 CRAN (R 3.3.0)
memoise 1.0.0 2016-01-29 CRAN (R 3.3.0)
multcomp * 1.4-5 2016-05-04 CRAN (R 3.3.0)
mvtnorm * 1.0-5 2016-02-02 CRAN (R 3.3.0)
plyr * 1.8.4 2016-06-08 CRAN (R 3.3.0)
RColorBrewer * 1.1-2 2014-12-07 CRAN (R 3.3.0)
Rcpp 0.12.3 2016-01-10 CRAN (R 3.2.3)
rmarkdown * 0.9.6 2016-05-01 CRAN (R 3.3.0)
rpanel 1.1-3 2014-02-28 CRAN (R 3.3.0)
sandwich 2.3-4 2015-09-24 CRAN (R 3.3.0)
stringi 1.1.1 2016-05-27 CRAN (R 3.3.0)
stringr 1.0.0 2015-04-30 CRAN (R 3.3.0)
survival * 2.39-4 2016-05-11 CRAN (R 3.3.0)
TH.data * 1.0-7 2016-01-28 CRAN (R 3.3.0)
withr 1.0.1 2016-02-04 CRAN (R 3.3.0)
wzRfun * 0.61 2016-06-20 local (walmes/wzRfun)
yaml 2.1.13 2014-06-12 CRAN (R 3.3.0)
zoo 1.7-13 2016-05-03 CRAN (R 3.3.0)
[1] "2016-06-23 15:47:21 BRT"