Diferenciação de isolados em relação às respostas para as condições de desenvolvimento
Glória Soriano Vidal, Henrique da Silva Silveira Duarte e Walmes Marques Zeviani
2020-08-19
Importação e análise exploratória
Os seguintes pacotes são utilizados para a realização das análises presentes nesse relatório. Exceto o wzRfun, todos os pacotes são oficiais e estão disponíveis no CRAN. O wzRfun está no github e pode ser instalado com devtools::install_github("walmes/wzRfun").
#-----------------------------------------------------------------------
# Pacotes.
library(lattice)
library(latticeExtra)
library(splines)
library(mgcv)
library(car)
library(wzRfun)
library(doBy)
library(emmeans)
library(tidyverse)
devtools::session_info()## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 4.0.2 (2020-06-22)
## os Ubuntu 20.04.1 LTS
## system x86_64, linux-gnu
## ui X11
## language en_US
## collate en_US.UTF-8
## ctype pt_BR.UTF-8
## tz America/Sao_Paulo
## date 2020-08-19
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date lib source
## abind 1.4-5 2016-07-21 [3] CRAN (R 4.0.2)
## assertthat 0.2.1 2019-03-21 [3] CRAN (R 4.0.1)
## backports 1.1.8 2020-06-17 [3] CRAN (R 4.0.1)
## blob 1.2.1 2020-01-20 [3] CRAN (R 4.0.1)
## broom 0.7.0 2020-07-09 [3] CRAN (R 4.0.2)
## callr 3.4.3 2020-03-28 [3] CRAN (R 4.0.1)
## car * 3.0-8 2020-05-21 [3] CRAN (R 4.0.2)
## carData * 3.0-4 2020-05-22 [3] CRAN (R 4.0.2)
## cellranger 1.1.0 2016-07-27 [3] CRAN (R 4.0.1)
## cli 2.0.2 2020-02-28 [3] CRAN (R 4.0.1)
## coda 0.19-3 2019-07-05 [3] CRAN (R 4.0.2)
## codetools 0.2-16 2018-12-24 [4] CRAN (R 4.0.0)
## colorspace 1.4-1 2019-03-18 [3] CRAN (R 4.0.1)
## crayon 1.3.4 2017-09-16 [3] CRAN (R 4.0.1)
## curl 4.3 2019-12-02 [3] CRAN (R 4.0.1)
## data.table 1.12.8 2019-12-09 [3] CRAN (R 4.0.2)
## DBI 1.1.0 2019-12-15 [3] CRAN (R 4.0.1)
## dbplyr 1.4.4 2020-05-27 [3] CRAN (R 4.0.1)
## Deriv 4.0 2019-12-10 [3] CRAN (R 4.0.2)
## desc 1.2.0 2018-05-01 [3] CRAN (R 4.0.1)
## devtools 2.3.0 2020-04-10 [3] CRAN (R 4.0.2)
## digest 0.6.25 2020-02-23 [3] CRAN (R 4.0.1)
## doBy * 4.6.7 2020-07-08 [3] CRAN (R 4.0.2)
## dplyr * 1.0.0 2020-05-29 [3] CRAN (R 4.0.1)
## ellipsis 0.3.1 2020-05-15 [3] CRAN (R 4.0.1)
## emmeans * 1.4.8 2020-06-26 [3] CRAN (R 4.0.2)
## estimability 1.3 2018-02-11 [3] CRAN (R 4.0.2)
## evaluate 0.14 2019-05-28 [3] CRAN (R 4.0.1)
## fansi 0.4.1 2020-01-08 [3] CRAN (R 4.0.1)
## forcats * 0.5.0 2020-03-01 [3] CRAN (R 4.0.1)
## foreign 0.8-79 2020-04-26 [4] CRAN (R 4.0.0)
## fs 1.4.1 2020-04-04 [3] CRAN (R 4.0.1)
## generics 0.0.2 2018-11-29 [3] CRAN (R 4.0.1)
## ggplot2 * 3.3.2 2020-06-19 [3] CRAN (R 4.0.1)
## glue 1.4.1 2020-05-13 [3] CRAN (R 4.0.1)
## gtable 0.3.0 2019-03-25 [3] CRAN (R 4.0.1)
## haven 2.3.1 2020-06-01 [3] CRAN (R 4.0.1)
## hms 0.5.3 2020-01-08 [3] CRAN (R 4.0.1)
## htmltools 0.5.0 2020-06-16 [3] CRAN (R 4.0.1)
## httr 1.4.2 2020-07-20 [3] CRAN (R 4.0.2)
## jpeg 0.1-8.1 2019-10-24 [3] CRAN (R 4.0.2)
## jsonlite 1.7.0 2020-06-25 [3] CRAN (R 4.0.2)
## knitr * 1.29 2020-06-23 [3] CRAN (R 4.0.2)
## lattice * 0.20-41 2020-04-02 [4] CRAN (R 4.0.0)
## latticeExtra * 0.6-29 2019-12-19 [3] CRAN (R 4.0.2)
## lifecycle 0.2.0 2020-03-06 [3] CRAN (R 4.0.1)
## lubridate 1.7.9 2020-06-08 [3] CRAN (R 4.0.1)
## magrittr 1.5 2014-11-22 [3] CRAN (R 4.0.1)
## MASS 7.3-51.6 2020-04-26 [4] CRAN (R 4.0.0)
## Matrix 1.2-18 2019-11-27 [4] CRAN (R 4.0.0)
## memoise 1.1.0 2017-04-21 [3] CRAN (R 4.0.2)
## mgcv * 1.8-31 2019-11-09 [4] CRAN (R 4.0.0)
## modelr 0.1.8 2020-05-19 [3] CRAN (R 4.0.1)
## multcomp 1.4-13 2020-04-08 [3] CRAN (R 4.0.2)
## munsell 0.5.0 2018-06-12 [3] CRAN (R 4.0.1)
## mvtnorm 1.1-1 2020-06-09 [3] CRAN (R 4.0.2)
## nlme * 3.1-148 2020-05-24 [3] CRAN (R 4.0.2)
## openxlsx 4.1.5 2020-05-06 [3] CRAN (R 4.0.2)
## pillar 1.4.6 2020-07-10 [3] CRAN (R 4.0.2)
## pkgbuild 1.1.0 2020-07-13 [3] CRAN (R 4.0.2)
## pkgconfig 2.0.3 2019-09-22 [3] CRAN (R 4.0.1)
## pkgload 1.1.0 2020-05-29 [3] CRAN (R 4.0.1)
## png 0.1-7 2013-12-03 [3] CRAN (R 4.0.2)
## prettyunits 1.1.1 2020-01-24 [3] CRAN (R 4.0.1)
## processx 3.4.3 2020-07-05 [3] CRAN (R 4.0.2)
## ps 1.3.3 2020-05-08 [3] CRAN (R 4.0.1)
## purrr * 0.3.4 2020-04-17 [3] CRAN (R 4.0.1)
## R6 2.4.1 2019-11-12 [3] CRAN (R 4.0.1)
## RColorBrewer 1.1-2 2014-12-07 [3] CRAN (R 4.0.2)
## Rcpp 1.0.5 2020-07-06 [3] CRAN (R 4.0.2)
## readr * 1.3.1 2018-12-21 [3] CRAN (R 4.0.1)
## readxl 1.3.1 2019-03-13 [3] CRAN (R 4.0.1)
## remotes 2.1.1 2020-02-15 [3] CRAN (R 4.0.2)
## reprex 0.3.0 2019-05-16 [3] CRAN (R 4.0.1)
## rio 0.5.16 2018-11-26 [3] CRAN (R 4.0.2)
## rlang 0.4.7 2020-07-09 [3] CRAN (R 4.0.2)
## rmarkdown * 2.3 2020-06-18 [3] CRAN (R 4.0.1)
## rprojroot 1.3-2 2018-01-03 [3] CRAN (R 4.0.1)
## rstudioapi 0.11 2020-02-07 [3] CRAN (R 4.0.1)
## rvest 0.3.5 2019-11-08 [3] CRAN (R 4.0.1)
## sandwich 2.5-1 2019-04-06 [3] CRAN (R 4.0.2)
## scales 1.1.1 2020-05-11 [3] CRAN (R 4.0.1)
## sessioninfo 1.1.1 2018-11-05 [3] CRAN (R 4.0.2)
## stringi 1.4.6 2020-02-17 [3] CRAN (R 4.0.1)
## stringr * 1.4.0 2019-02-10 [3] CRAN (R 4.0.1)
## survival 3.1-12 2020-04-10 [4] CRAN (R 4.0.0)
## testthat 2.3.2 2020-03-02 [3] CRAN (R 4.0.1)
## TH.data 1.0-10 2019-01-21 [3] CRAN (R 4.0.2)
## tibble * 3.0.3 2020-07-10 [3] CRAN (R 4.0.2)
## tidyr * 1.1.0 2020-05-20 [3] CRAN (R 4.0.1)
## tidyselect 1.1.0 2020-05-11 [3] CRAN (R 4.0.1)
## tidyverse * 1.3.0 2019-11-21 [3] CRAN (R 4.0.1)
## usethis 1.6.1 2020-04-29 [3] CRAN (R 4.0.2)
## vctrs 0.3.2 2020-07-15 [3] CRAN (R 4.0.2)
## withr 2.2.0 2020-04-20 [3] CRAN (R 4.0.1)
## wzRfun * 0.91 2020-07-21 [3] Github (walmes/wzRfun@a90b379)
## xfun 0.15 2020-06-21 [3] CRAN (R 4.0.1)
## xml2 1.3.2 2020-04-23 [3] CRAN (R 4.0.1)
## xtable 1.8-4 2019-04-21 [3] CRAN (R 4.0.2)
## yaml 2.2.1 2020-02-01 [3] CRAN (R 4.0.1)
## zip 2.0.4 2019-09-01 [3] CRAN (R 4.0.2)
## zoo 1.8-8 2020-05-02 [3] CRAN (R 4.0.2)
##
## [1] /home/walmes/R/x86_64-pc-linux-gnu-library/4.0
## [2] /usr/local/lib/R/site-library
## [3] /usr/lib/R/site-library
## [4] /usr/lib/R/libraryOs dados dos 3 ensaios (germinação, latência e pústulas) foram arquivos em documentos de texto pleno com valores separados por delimitador (tabulação ou vírgula). O código abaixo importa para o R as tabelas de dados.
#-----------------------------------------------------------------------
# Leitura dos conjuntos de dados.
# Para não mostrar a mensagem de parse das colunas ao importar.
options(readr.num_columns = 0)
ger <- read_tsv("germinacao.txt", progress = FALSE)
ger$rept <- rep(1:3, times = nrow(ger)/3)
# dput(substr(tolower(names(ger)), 1, 3))
# str(ger)
ger <- ger %>%
gather(key = "molha_rept",
value = "germ",
-c(experimento, origem, temperatura, rept)) %>%
separate(col = molha_rept, into = c("molha", "leaf")) %>%
# replace_na(replace = list(leaf = 1)) %>%
mutate(molha = as.integer(molha),
origem = factor(origem),
experimento = factor(experimento)) %>%
group_by(experimento, origem, temperatura, molha, rept) %>%
summarise_at("germ", "mean") %>%
ungroup()
names(ger) <- c("exp", "est", "tem", "mol", "rept", "ger")
str(ger)## tibble [648 × 6] (S3: tbl_df/tbl/data.frame)
## $ exp : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
## $ est : Factor w/ 3 levels "PR","SC","SP": 1 1 1 1 1 1 1 1 1 1 ...
## $ tem : num [1:648] 10 10 10 10 10 10 10 10 10 10 ...
## $ mol : int [1:648] 3 3 3 6 6 6 9 9 9 12 ...
## $ rept: int [1:648] 1 2 3 1 2 3 1 2 3 1 ...
## $ ger : num [1:648] 15 15.7 15.7 25.3 27 ...lat <- read_csv("latencia.csv")
# dput(substr(tolower(names(lat)), 1, 3))
names(lat) <- c("est", "tem", "rep", "pl", "fi")
lat <- lat %>%
mutate(fi = as.numeric(sub(",", ".", fi)),
est = factor(est))
str(lat)## tibble [108 × 5] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ est: Factor w/ 3 levels "PR","SC","SP": 3 3 3 3 3 3 3 3 3 3 ...
## $ tem: num [1:108] 10 10 10 10 13 13 13 13 18 18 ...
## $ rep: num [1:108] 1 2 3 4 1 2 3 4 1 2 ...
## $ pl : num [1:108] 0 0 0 0 15 17 15 16 15 12 ...
## $ fi : num [1:108] 0 0 0 0 0.5 0.75 0.5 0.75 0.75 1.5 ...
## - attr(*, "spec")=
## .. cols(
## .. ESTADO = col_character(),
## .. TEMPERATURA = col_double(),
## .. REP = col_double(),
## .. PL = col_double(),
## .. FI = col_character()
## .. )pus <- read_csv("pustulas.csv")
# dput(substr(tolower(names(pus)), 1, 3))
names(pus) <- c("est", "rep", "fol", "dia", "npu")
pus <- pus %>%
mutate(est = factor(est))
str(pus)## tibble [312 × 5] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ est: Factor w/ 3 levels "PR","SC","SP": 3 3 3 3 3 3 3 3 3 3 ...
## $ rep: num [1:312] 1 1 2 2 3 3 4 4 1 1 ...
## $ fol: num [1:312] 1 2 1 2 1 2 1 2 1 2 ...
## $ dia: num [1:312] 0 0 0 0 0 0 0 0 1 1 ...
## $ npu: num [1:312] 0 0 0 0 0 0 0 0 18 2 ...
## - attr(*, "spec")=
## .. cols(
## .. ESTADO = col_character(),
## .. REP = col_double(),
## .. FOLHA = col_double(),
## .. DIA = col_double(),
## .. NPUSTULAS = col_double()
## .. )A seguir são feitos gráficos para análise exploratória dos dados. Neles é possível verificar a existência do efeito dos fatores experimentais avaliados (temperatura e período de molhamento) sob a germinação dos isolados, o período de latência e a frequência de infecção.
#-----------------------------------------------------------------------
# Análise exploratória para inspeção.
gg1 <-
ggplot(data = ger,
mapping = aes(x = tem, y = ger, col = as.factor(mol))) +
geom_point() +
stat_summary(fun.y = mean, geom = "line") +
facet_grid(facets = exp ~ est) +
ylab("Germination (%)") +
xlab("Temperature (°C)") +
# labs(title = "A: Temperature effect on germination") +
labs(title = "A") +
scale_color_discrete("Wetness period (h)") +
theme(legend.position = "top") +
guides(color = guide_legend(nrow = 1, title.position = "left"))
# gg1
gg2 <-
ggplot(data = ger,
mapping = aes(x = mol, y = ger, col = as.factor(tem))) +
geom_point() +
stat_summary(fun.y = mean, geom = "line") +
facet_grid(facets = exp ~ est) +
ylab("Germination (%)") +
xlab("Wetness period (h)") +
# labs(title = "B: Wetness period effect on germination") +
labs(title = "B") +
scale_color_discrete("Temperature (°C)") +
theme(legend.position = "top") +
guides(color = guide_legend(nrow = 1, title.position = "left"))
# gg2
ggplot(data = lat,
mapping = aes(x = tem, y = pl, col = est)) +
geom_point() +
stat_summary(fun.y = mean, geom = "line") +
geom_jitter(position = position_jitter(height = 0)) +
ylab("Período de latência (dias)") +
xlab("Temperatura (°C)") +
scale_color_discrete("Origem") +
theme(legend.position = "top") +
guides(color = guide_legend(nrow = 1, title.position = "left"))
ggplot(data = lat,
mapping = aes(x = tem, y = fi, col = est)) +
geom_point() +
scale_y_log10() +
stat_summary(fun.y = mean, geom = "line") +
ylab("Frequência de infecção") +
xlab("Temperatura (°C)") +
scale_color_discrete("Origem") +
theme(legend.position = "top") +
guides(color = guide_legend(nrow = 1, title.position = "left"))
# ATTENTION: Não fazer a análise, esperar próximo experimento.
ggplot(data = pus,
mapping = aes(x = dia, y = npu, col = est)) +
geom_point() +
scale_y_log10() +
geom_smooth(se = FALSE) +
ylab("Número de pústulas") +
xlab("Tempo após o período de latência (dias)") +
scale_color_discrete("Origem") +
theme(legend.position = "top") +
guides(color = guide_legend(nrow = 1, title.position = "left"))
gridExtra::grid.arrange(gg1, gg2, ncol = 1)
Verificou-se pelos gráficos que algumas variáveis resposta possivelmente terão que ser transformadas para os pressupostos dos modelos sejam atendidos. Da mesma forma, os níveis mais mais extremos de temepratura não produziram valores não nulos em algumas respostas e portanto poderão ser removidos da análise.
Germinação
Para análise de germinação será considerado o modelo linear com resposta gaussiana (embora pudesse ser resposta binomial). Os efeitos de temperatura e período de molhamento, que são fatores métricos, sob a germinação média serão acomodados utilizando polinômios. Baseado nos gráficos de análise exploratória, polinômios de grau inferior a 3 devem ser suficientes para acomodar tais efeitos.
#-----------------------------------------------------------------------
# Análise da germinação.
# NOTE: para não precisar usar um grau alto de polinômio, optou-se por
# remover da análise temperaturas inferiores à 15 °C porque, caso fossem
# mantidas seria necessário usar um modelo complexo para acomodar o
# comportamento da germinação na proximidade da temperatura 15. Pela
# análise exploratória viu-se que o máximo está distante de 15, sendo
# portanto essa remoção não prejudicial para determinação da temperatura
# de ótimo.
unique(ger$tem)## [1] 10 15 20 25 30 35ger <- subset(ger, tem >= 15)
# Modelo linear com interações e termos não lineares.
m0 <- lm(log(ger) ~ exp + est * (tem + I(tem^2) + mol + I(mol^0.5) + tem:mol),
data = ger)
# Quadro de análise de variância.
anova(m0)## Analysis of Variance Table
##
## Response: log(ger)
## Df Sum Sq Mean Sq F value Pr(>F)
## exp 1 0.787 0.787 28.6107 1.327e-07 ***
## est 2 2.376 1.188 43.1677 < 2.2e-16 ***
## tem 1 7.029 7.029 255.4577 < 2.2e-16 ***
## I(tem^2) 1 42.733 42.733 1553.0828 < 2.2e-16 ***
## mol 1 22.581 22.581 820.6961 < 2.2e-16 ***
## I(mol^0.5) 1 12.328 12.328 448.0431 < 2.2e-16 ***
## tem:mol 1 0.313 0.313 11.3619 0.0008053 ***
## est:tem 2 0.001 0.000 0.0122 0.9878774
## est:I(tem^2) 2 0.590 0.295 10.7291 2.715e-05 ***
## est:mol 2 0.054 0.027 0.9753 0.3777812
## est:I(mol^0.5) 2 0.018 0.009 0.3230 0.7241071
## est:tem:mol 2 0.011 0.005 0.1984 0.8201089
## Residuals 521 14.335 0.028
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# summary(m0)
# Análise de diagnóstico do ajuste pelos residuos.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Ajuste do modelo reduzido eliminando termos não relevantes.
m0 <- lm(log(ger) ~ exp + est * (tem + I(tem^2)) + mol + I(mol^0.5) +
tem:mol,
data = ger)
# Quadro de ANOVA.
anova(m0)## Analysis of Variance Table
##
## Response: log(ger)
## Df Sum Sq Mean Sq F value Pr(>F)
## exp 1 0.787 0.787 28.7748 1.219e-07 ***
## est 2 2.376 1.188 43.4154 < 2.2e-16 ***
## tem 1 7.029 7.029 256.9235 < 2.2e-16 ***
## I(tem^2) 1 42.733 42.733 1561.9943 < 2.2e-16 ***
## mol 1 22.581 22.581 825.4052 < 2.2e-16 ***
## I(mol^0.5) 1 12.328 12.328 450.6139 < 2.2e-16 ***
## est:tem 2 0.001 0.000 0.0123 0.9878082
## est:I(tem^2) 2 0.590 0.295 10.7906 2.553e-05 ***
## tem:mol 1 0.313 0.313 11.4271 0.0007775 ***
## Residuals 527 14.418 0.027
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Estimativas dos coeficientes.
summary(m0)##
## Call:
## lm(formula = log(ger) ~ exp + est * (tem + I(tem^2)) + mol +
## I(mol^0.5) + tem:mol, data = ger)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.49687 -0.10271 -0.00352 0.11877 0.36502
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.211e-01 1.835e-01 -0.660 0.509548
## exp2 7.636e-02 1.424e-02 5.364 1.22e-07 ***
## estSC -1.245e+00 2.481e-01 -5.021 7.06e-07 ***
## estSP -4.188e-01 2.481e-01 -1.688 0.091959 .
## tem 2.743e-01 1.488e-02 18.433 < 2e-16 ***
## I(tem^2) -5.879e-03 2.947e-04 -19.949 < 2e-16 ***
## mol -4.897e-02 3.169e-03 -15.452 < 2e-16 ***
## I(mol^0.5) 5.132e-01 2.418e-02 21.228 < 2e-16 ***
## estSC:tem 9.472e-02 2.098e-02 4.514 7.86e-06 ***
## estSP:tem 3.144e-02 2.098e-02 1.498 0.134711
## estSC:I(tem^2) -1.901e-03 4.168e-04 -4.562 6.30e-06 ***
## estSP:I(tem^2) -6.349e-04 4.168e-04 -1.523 0.128291
## tem:mol 2.215e-04 6.552e-05 3.380 0.000777 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1654 on 527 degrees of freedom
## Multiple R-squared: 0.8602, Adjusted R-squared: 0.8571
## F-statistic: 270.3 on 12 and 527 DF, p-value: < 2.2e-16# Análise de diagnóstico do ajuste pelos residuos.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Ainda existe alguma falta de ajuste mas que não é sistemática ou
# parsimoniamente corrigível.
residualPlots(m0)
## Test stat Pr(>|Test stat|)
## exp
## est
## tem 0.0535 0.9573171
## I(tem^2) 20.5554 < 2.2e-16 ***
## mol 3.9028 0.0001075 ***
## I(mol^0.5) 0.0918 0.9268960
## Tukey test -4.6426 3.44e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#-----------------------------------------------------------------------
# Predição.
# Valores únicos que foram avaliados dos fatores experimentais.
sort(unique(ger$tem))## [1] 15 20 25 30 35sort(unique(ger$mol))## [1] 3 6 9 12 24 48# Grid de valores para a predição e então construção da superfície.
# grid <- with(ger,
# expand.grid(est = unique(est),
# tem = seq(min(tem), max(tem), length.out = 30),
# mol = seq(min(mol), max(mol), length.out = 30),
# KEEP.OUT.ATTRS = FALSE))
# grid$ger <- predict(m0, newdata = grid) # Modelo linear.
# grid$ger <- predict(m1, newdata = grid) # GAM.
# grid$ger <- predict(m2, newdata = grid) # Splines.
# str(grid)
# Grid de valores para a predição e então construção da superfície.
grid <- emmeans(m0, specs = ~est + tem + mol,
at = with(ger,
list(tem = seq(min(tem),
max(tem),
length.out = 30),
mol = seq(min(mol),
max(mol),
length.out = 30))))
grid <- as.data.frame(grid)
str(grid)## 'data.frame': 2700 obs. of 8 variables:
## $ est : Factor w/ 3 levels "PR","SC","SP": 1 2 3 1 2 3 1 2 3 1 ...
## $ tem : num 15 15 15 15.7 15.7 ...
## $ mol : num 3 3 3 3 3 3 3 3 3 3 ...
## $ emmean : num 3.46 3.21 3.37 3.53 3.3 ...
## $ SE : num 0.0304 0.0304 0.0304 0.0278 0.0278 ...
## $ df : num 527 527 527 527 527 527 527 527 527 527 ...
## $ lower.CL: num 3.4 3.15 3.31 3.47 3.24 ...
## $ upper.CL: num 3.52 3.27 3.43 3.58 3.35 ...
## - attr(*, "estName")= chr "emmean"
## - attr(*, "clNames")= chr [1:2] "lower.CL" "upper.CL"
## - attr(*, "pri.vars")= chr [1:3] "est" "tem" "mol"
## - attr(*, "adjust")= chr "none"
## - attr(*, "side")= num 0
## - attr(*, "delta")= num 0
## - attr(*, "type")= chr "link"
## - attr(*, "mesg")= chr [1:3] "Results are averaged over the levels of: exp" "Results are given on the log (not the response) scale." "Confidence level used: 0.95"colr <- RColorBrewer::brewer.pal(11, "Spectral")
colr <- colorRampPalette(colr, space = "rgb")
# Gráfico de níveis com linhas de contornos de nível.
levelplot(exp(emmean) ~ mol + tem | est,
data = grid,
xlab = "Período de molhamento (h)",
ylab = "Temperatura (°C)",
as.table = TRUE,
contour = TRUE,
labels = TRUE,
at = seq(from = 0, to = 110, by = 5),
col.regions = colr(101))
# Gráfico 3D da superfície.
# x11(width = 10, height = 4)
# dev.off()
wf <-
wireframe(exp(emmean) ~ mol + tem | est,
layout = c(NA, 1),
data = grid,
# col = "white",
xlab = list("Wetness period (h)", rot = 30),
ylab = list("Temperature (°C)", rot = -40),
zlab = list("Germination (%)", rot = 90),
zlim = c(0, 105),
as.table = TRUE,
panel.3d.wireframe = panel.3d.contour,
col.contour = "black",
type = c("on", "bottom")[2],
drape = TRUE,
col.regions = colr(101),
# zoom = 0.8,
scales = list(arrows = FALSE,
z = list(distance = 1.25)))
wf
#-----------------------------------------------------------------------
# Testando hipótestes.
# Termos presentes no modelo considerado.
formula(m0)## log(ger) ~ exp + est * (tem + I(tem^2)) + mol + I(mol^0.5) +
## tem:mol# Teste entre modelos aninhados para avaliar o efeito de localidade dos
# isolados, ou seja, se as curvas de germinação independem da
# localidade.
anova(m0,
update(m0, # Modelo com os termos de interação.
. ~ . -est:tem - est:I(tem^2))) # Modelo sem ele.## Analysis of Variance Table
##
## Model 1: log(ger) ~ exp + est * (tem + I(tem^2)) + mol + I(mol^0.5) +
## tem:mol
## Model 2: log(ger) ~ exp + est + tem + I(tem^2) + mol + I(mol^0.5) + tem:mol
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 527 14.418
## 2 531 15.009 -4 -0.59109 5.4014 0.0002871 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# NOTE: conforme o teste acima, os modelos tem curvas que diferem à 5%.
anova(m0) # O termo `est:I(tem^2)` é significativo.## Analysis of Variance Table
##
## Response: log(ger)
## Df Sum Sq Mean Sq F value Pr(>F)
## exp 1 0.787 0.787 28.7748 1.219e-07 ***
## est 2 2.376 1.188 43.4154 < 2.2e-16 ***
## tem 1 7.029 7.029 256.9235 < 2.2e-16 ***
## I(tem^2) 1 42.733 42.733 1561.9943 < 2.2e-16 ***
## mol 1 22.581 22.581 825.4052 < 2.2e-16 ***
## I(mol^0.5) 1 12.328 12.328 450.6139 < 2.2e-16 ***
## est:tem 2 0.001 0.000 0.0123 0.9878082
## est:I(tem^2) 2 0.590 0.295 10.7906 2.553e-05 ***
## tem:mol 1 0.313 0.313 11.4271 0.0007775 ***
## Residuals 527 14.418 0.027
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m0)##
## Call:
## lm(formula = log(ger) ~ exp + est * (tem + I(tem^2)) + mol +
## I(mol^0.5) + tem:mol, data = ger)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.49687 -0.10271 -0.00352 0.11877 0.36502
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.211e-01 1.835e-01 -0.660 0.509548
## exp2 7.636e-02 1.424e-02 5.364 1.22e-07 ***
## estSC -1.245e+00 2.481e-01 -5.021 7.06e-07 ***
## estSP -4.188e-01 2.481e-01 -1.688 0.091959 .
## tem 2.743e-01 1.488e-02 18.433 < 2e-16 ***
## I(tem^2) -5.879e-03 2.947e-04 -19.949 < 2e-16 ***
## mol -4.897e-02 3.169e-03 -15.452 < 2e-16 ***
## I(mol^0.5) 5.132e-01 2.418e-02 21.228 < 2e-16 ***
## estSC:tem 9.472e-02 2.098e-02 4.514 7.86e-06 ***
## estSP:tem 3.144e-02 2.098e-02 1.498 0.134711
## estSC:I(tem^2) -1.901e-03 4.168e-04 -4.562 6.30e-06 ***
## estSP:I(tem^2) -6.349e-04 4.168e-04 -1.523 0.128291
## tem:mol 2.215e-04 6.552e-05 3.380 0.000777 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1654 on 527 degrees of freedom
## Multiple R-squared: 0.8602, Adjusted R-squared: 0.8571
## F-statistic: 270.3 on 12 and 527 DF, p-value: < 2.2e-16#-----------------------------------------------------------------------
# Estimação intervalar do ponto de máximo.
# A curvatura ao redor do ponto de máximo não a mesma para todos as
# localidades pois não houve interação `est:I(tem^2)` (termo quadrático
# comum), portanto, o decaimento da germinação com o afastamento da
# tempetura ótima é igual para todas as localidades.
# Ponto de máximo de uma curva quadrática y = c + b * x + a * x^2 é dado
# por x = -b/a. Sendo assim, o ponto de máximo para cada localidade é
# calculado pelos código abaixo.
coef(m0)## (Intercept) exp2 estSC estSP tem
## -0.1211401989 0.0763625934 -1.2454100138 -0.4187815554 0.2742774681
## I(tem^2) mol I(mol^0.5) estSC:tem estSP:tem
## -0.0058791299 -0.0489666400 0.5132082073 0.0947178280 0.0314358378
## estSC:I(tem^2) estSP:I(tem^2) tem:mol
## -0.0019014829 -0.0006348602 0.0002214988# Extraí os coeficientes que são usados para determinar o coeficientes
# da curva com relação a temperatura (`tem`) para cada localidade
# (`est`). Está fixando o molhamento em 48 horas.
lsm <- LE_matrix(m0, effect = "est", at = list(tem = 1, mol = 48))
lsm## (Intercept) exp2 estSC estSP tem I(tem^2) mol I(mol^0.5) estSC:tem
## [1,] 1 0.5 0 0 1 1 48 6.928203 0
## [2,] 1 0.5 1 0 1 1 48 6.928203 1
## [3,] 1 0.5 0 1 1 1 48 6.928203 0
## estSP:tem estSC:I(tem^2) estSP:I(tem^2) tem:mol
## [1,] 0 0 0 48
## [2,] 0 1 0 48
## [3,] 1 0 1 48# Cria o vetor com os coeficientes da equação de segundo (a, b, c) grau
# para cada localidade.
# Multiplica estimativas pelo valor das variáveis.
coefs <- sweep(lsm, 2, rbind(coef(m0)), FUN = "*")
# Agrupa para somar coeficientes de mesma ordem.
coeflist <- list(
tapply(coefs[1, as.logical(lsm[1, ])], c(1, 1, 2, 3, 1, 1, 2), sum),
tapply(coefs[2, as.logical(lsm[2, ])], c(1, 1, 1, 2, 3, 1, 1, 2, 3, 2), sum),
tapply(coefs[3, as.logical(lsm[3, ])], c(1, 1, 1, 2, 3, 1, 1, 2, 3, 2), sum))
names(coeflist) <- levels(ger$est)
coeflist <- lapply(coeflist, FUN = setNames, nm = c("c", "b", "a"))
coeflist## $PR
## c b a
## 1.12225314 0.28490941 -0.00587913
##
## $SC
## c b a
## -0.123156874 0.379627240 -0.007780613
##
## $SP
## c b a
## 0.70347158 0.31634525 -0.00651399# Obtém o ponto crítico.
temmax <- sapply(coeflist,
FUN = function(x) {
x <- -0.5 * x["b"]/x["a"]
names(x) <- NULL
return(x)
})
temmax## PR SC SP
## 24.23058 24.39572 24.28199# Exibe os valores.
levelplot(emmean ~ mol + tem | est,
data = grid,
xlab = "Wet períod (h)",
ylab = "Temperature (°C)",
as.table = TRUE,
contour = TRUE,
col.regions = heat.colors(101)) +
latticeExtra::layer({
panel.abline(h = temmax[which.packet()], lty = 2)
})
#-----------------------------------------------------------------------
# Calculando erro padrão para as temperaturas no ponto de máximo (método
# delta).
# Nomes para os parâmetros para usar nas contas no método delta.
# dput(names(coef(m0)))
pn <- c("(Intercept)",
"exp2",
"estSC",
"estSP",
"tem",
"tem.squared",
"mol",
"mol.root",
"estSC.tem",
"estSP.tem",
"estSC.tem2",
"estSP.tem2",
"tem.mol")
# Expressões que serão avaliadas para obter a estimativa e erro padrão.
g <- list("-0.5 * (tem + 48 * tem.mol)/tem.squared",
"-0.5 * (tem + estSC.tem + 48 * tem.mol)/(tem.squared + estSC.tem2)",
"-0.5 * (tem + estSP.tem + 48 * tem.mol)/(tem.squared + estSP.tem2)")
names(g) <- levels(ger$est)
# Avalia as expressões e organiza os resultados em tabela.
pars <- lapply(g,
FUN = deltaMethod, # Aplica o método delta.
object = m0,
parameterNames = pn) %>%
do.call(what = rbind) %>%
as.data.frame() %>%
setNames(nm = c("estimate", "se", "lwr", "upr")) %>%
mutate(est = factor(rownames(.)))
pars## estimate se lwr upr est
## 1 24.23058 0.2309048 23.77801 24.68314 PR
## 2 24.39572 0.1735395 24.05558 24.73585 SC
## 3 24.28199 0.2080248 23.87427 24.68971 SP# Gráfico de segmentos para exibir os intervalos de confiança.
segplot(factor(est) ~ lwr + upr,
data = pars,
draw = FALSE,
centers = estimate,
xlab = "Temperature (°C)",
ylab = "Origin") +
latticeExtra::layer({
panel.text(x = centers,
y = z,
labels = sprintf("%0.2f", centers),
pos = 3)
})
# Adiciona a germinação para o molhamento = 48h.
X <- lapply(1:nlevels(ger$est),
FUN = function(level) {
X <- LE_matrix(m0,
at = list(tem = pars$estimate[level],
mol = 48,
est = pars$est[level]))
return(X)
})
X <- do.call(rbind, X)
pars$ger48 <- exp(X %*% coef(m0))
pars## estimate se lwr upr est ger48
## 1 24.23058 0.2309048 23.77801 24.68314 PR 96.93230
## 2 24.39572 0.1735395 24.05558 24.73585 SC 90.69319
## 3 24.28199 0.2080248 23.87427 24.68971 SP 94.08675#-----------------------------------------------------------------------
# Gráficos com os perfis para a temperatura e molhamento.
grid2 <- list(tem = with(ger,
expand.grid(est = unique(est),
tem = seq(min(tem),
max(tem),
length.out = 60),
mol = 48,
exp = "1",
KEEP.OUT.ATTRS = FALSE)),
mol = with(ger,
expand.grid(est = unique(est),
tem = 0,
mol = seq(min(mol),
max(mol),
length.out = 60),
exp = "1",
KEEP.OUT.ATTRS = FALSE)))
grid2$mol$tem <- pars$estimate[match(grid2$mol$est, table = pars$est)]
# Predição.
grid2 <- lapply(grid2,
FUN = function(x) {
x$ger <- exp(predict(m0, newdata = x))
return(x)
})
xyplot(ger ~ tem,
groups = est,
data = grid2$tem,
type = "l",
auto.key = list(corner = c(0.05, 0.95),
columns = 1,
lines = TRUE,
points = FALSE,
title = "Origem",
cex.title = 1),
ylab = expression("Germinação" ~ ("%") ~
"com 48 horas de molhamento"),
# sub = list("48 horas de molhamento", font = 1, cex = 0.8),
xlab = expression("Temperatura" ~ (degree * C))) +
latticeExtra::layer({
panel.segments(pars$estimate,
0,
pars$estimate,
pars$ger48,
lty = 2)
})
xyplot(ger ~ mol,
groups = est,
data = grid2$mol,
type = "l",
auto.key = list(corner = c(0.05, 0.95),
columns = 1,
lines = TRUE,
points = FALSE,
title = "Origem",
cex.title = 1),
ylab = expression("Germinação" ~ ("%") ~
"na temperatura ótima"),
# sub = list("Temperatura ótima", font = 1, cex = 0.8),
xlab = "Molhamento (horas)")
Período de latência
A análise do período de latência e frequência de infecção considerou polinômio de 3 grau pois tem a flexibilidade necessária para acomodar o efeito da temperatura para ambas respostas. O ponto crítico na curva estimada será encontrado fazendo-se minimização numérica (embora seja possível por expressões analíticas, o resultado será o mesmo).
#-----------------------------------------------------------------------
# Período de latência.
# Para um ajuste adequado do modelo, serão removidas as temperaturas
# onde não foi determinado período de latência (valor registrado foi 0).
# Para um melhor ajuste do modelo, será usando apenas valores de
# temperatura superior a 15 pois pela análise exploratória o ponto de
# mínimo está distante da borda esquerda e não haverá prejuízo com a
# remoção.
lat <- subset(lat, pl > 0 & tem > 15)
# # Usando base-splines.
# m0 <- lm(pl ~ est * bs(tem),
# data = lat)
# Usar polinômio de grau 3 dá igual aos bs() sem knots internos.
m0 <- lm(pl ~ est * poly(tem, degree = 3),
data = lat)
# Quadro de ANOVA e coeficientes estimados.
anova(m0)## Analysis of Variance Table
##
## Response: pl
## Df Sum Sq Mean Sq F value Pr(>F)
## est 2 4.933 2.4667 1.6186 0.208803
## poly(tem, degree = 3) 3 80.846 26.9485 17.6834 7.209e-08 ***
## est:poly(tem, degree = 3) 6 29.805 4.9675 3.2597 0.009063 **
## Residuals 48 73.149 1.5239
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m0)##
## Call:
## lm(formula = pl ~ est * poly(tem, degree = 3), data = lat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.63040 -0.81572 -0.06264 0.81506 2.72749
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.4000 0.2760 41.299 < 2e-16 ***
## estSC -0.3000 0.3904 -0.768 0.44596
## estSP 0.4000 0.3904 1.025 0.31066
## poly(tem, degree = 3)1 -3.6500 2.1382 -1.707 0.09427 .
## poly(tem, degree = 3)2 5.5557 2.1382 2.598 0.01241 *
## poly(tem, degree = 3)3 -3.1473 2.1382 -1.472 0.14756
## estSC:poly(tem, degree = 3)1 8.3929 3.0238 2.776 0.00783 **
## estSP:poly(tem, degree = 3)1 -1.0273 3.0238 -0.340 0.73555
## estSC:poly(tem, degree = 3)2 7.7157 3.0238 2.552 0.01396 *
## estSP:poly(tem, degree = 3)2 1.5194 3.0238 0.502 0.61762
## estSC:poly(tem, degree = 3)3 2.2059 3.0238 0.729 0.46924
## estSP:poly(tem, degree = 3)3 0.6149 3.0238 0.203 0.83973
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.234 on 48 degrees of freedom
## Multiple R-squared: 0.6124, Adjusted R-squared: 0.5236
## F-statistic: 6.895 on 11 and 48 DF, p-value: 8.225e-07# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Grid de valores para a predição.
grid <- with(lat,
expand.grid(est = unique(est),
tem = seq(min(tem), max(tem), length.out = 50),
KEEP.OUT.ATTRS = FALSE))
grid <- cbind(grid,
as.data.frame(predict(m0,
newdata = grid,
interval = "confidence")))
str(grid)## 'data.frame': 150 obs. of 5 variables:
## $ est: Factor w/ 3 levels "PR","SC","SP": 3 1 2 3 1 2 3 1 2 3 ...
## $ tem: num 18 18 18 18.2 18.2 ...
## $ fit: num 13.9 13.2 12.3 13.6 12.9 ...
## $ lwr: num 12.7 12 11.1 12.5 11.8 ...
## $ upr: num 15.2 14.4 13.5 14.7 13.9 ...# Função para obter o ponto de mínimo com avaliação a predição.
predfun <- function(x, est) {
predict(m0, newdata = list(est = est, tem = x))
}
# Determinar o número mínimo de dias para cada origem. Fazer minimização
# numérica.
opts <- sapply(levels(lat$est),
simplify = FALSE,
FUN = function(est) {
# Minimização numérica em uma dimensão.
opt <- optimise(f = predfun,
interval = c(18, 28),
est = est)
setNames(unlist(opt, use.names = FALSE),
nm = c("tem", "pl"))
}) %>%
do.call(what = rbind) %>%
as.data.frame() %>%
mutate(est = rownames(.))
opts## tem pl est
## 1 21.90253 10.331821 PR
## 2 22.14304 8.941689 SC
## 3 22.43209 10.596390 SPp1 <-
xyplot(pl ~ tem,
groups = est,
jitter.x = TRUE,
xlab = "Temperature (°C)",
ylab = "Latent period (days)",
auto.key = list(columns = 3),
data = lat,
subset = pl > 0 & tem > 15) +
latticeExtra::as.layer({
xyplot(fit ~ tem,
type = "l",
groups = est,
cty = "bands",
ly = grid$lwr,
uy = grid$upr,
alpha = 0.45,
prepanel = prepanel.cbH,
panel = panel.superpose,
panel.groups = panel.cbH,
data = grid)
}, under = TRUE) +
latticeExtra::layer({
panel.segments(opts$tem,
0,
opts$tem,
opts$pl,
lty = 2)
})
p1
# Teste de hipótese para verificar se o número de dias difere.
aux <- predict(m0, newdata = opts, interval = "confidence")
opts <- opts %>% cbind(aux) %>% mutate(est = factor(est))
# Criar a matriz X para aplicar o teste formal.
cfs <- attr(m0$model[[3]], "coefs")
X <- model.matrix(~est * poly(tem, degree = 3, coefs = cfs),
data = opts)
rownames(X) <- levels(lat$est)
# Avalia os contrastes par a par.
opts <- apmc(X, model = m0, focus = "est") %>%
select(est, cld) %>%
inner_join(opts)
p2 <-
segplot(factor(est) ~ lwr + upr,
data = opts,
draw = FALSE,
cld = opts$cld,
centers = fit,
xlab = "Shortest latent period (days)",
ylab = "Origin") +
latticeExtra::layer({
panel.text(x = centers,
y = z,
labels = sprintf("%0.2f %s", centers, cld),
pos = 3)
})
p2
gridExtra::grid.arrange(p1, p2, ncol = 1, heights = c(2, 1.5))
Frequência de infecção
A análise segue as mesmas etapas da análise para o período de latência.
#-----------------------------------------------------------------------
# Frequência de infecção.
# ATTENTION: remover a temperatura 15 para ter um melhor ajuste do
# modelo proposto. Pela análise exploratória tem-se que a tempeturatura
# de ótimo está na região [15, 30]. Então, removendo-se valores de
# tempetura inferior a 15 beneficia o ajuste do modelo local para
# determinação do ponto de máximo.
lat <- subset(lat, fi > 0 & tem > 15)
# m0 <- lm(log10(fi) ~ est * bs(tem),
# data = lat)
m0 <- lm(log10(fi) ~ est * poly(tem, degree = 3),
data = lat)
anova(m0)## Analysis of Variance Table
##
## Response: log10(fi)
## Df Sum Sq Mean Sq F value Pr(>F)
## est 2 0.5567 0.2783 4.8230 0.01234 *
## poly(tem, degree = 3) 3 22.3027 7.4342 128.8184 < 2.2e-16 ***
## est:poly(tem, degree = 3) 6 2.1556 0.3593 6.2252 6.908e-05 ***
## Residuals 48 2.7701 0.0577
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m0)##
## Call:
## lm(formula = log10(fi) ~ est * poly(tem, degree = 3), data = lat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.41924 -0.15064 -0.02524 0.15069 0.40989
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.07358 0.05372 19.986 < 2e-16 ***
## estSC 0.19396 0.07597 2.553 0.01391 *
## estSP -0.01935 0.07597 -0.255 0.80000
## poly(tem, degree = 3)1 1.41715 0.41609 3.406 0.00134 **
## poly(tem, degree = 3)2 -3.29446 0.41609 -7.918 2.93e-10 ***
## poly(tem, degree = 3)3 1.15522 0.41609 2.776 0.00782 **
## estSC:poly(tem, degree = 3)1 -1.99904 0.58844 -3.397 0.00138 **
## estSP:poly(tem, degree = 3)1 0.19411 0.58844 0.330 0.74293
## estSC:poly(tem, degree = 3)2 -2.62417 0.58844 -4.460 4.94e-05 ***
## estSP:poly(tem, degree = 3)2 -1.14729 0.58844 -1.950 0.05707 .
## estSC:poly(tem, degree = 3)3 -0.29212 0.58844 -0.496 0.62186
## estSP:poly(tem, degree = 3)3 -0.29531 0.58844 -0.502 0.61806
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2402 on 48 degrees of freedom
## Multiple R-squared: 0.9003, Adjusted R-squared: 0.8775
## F-statistic: 39.4 on 11 and 48 DF, p-value: < 2.2e-16# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Grid de valores para a predição.
grid <- with(lat,
expand.grid(est = unique(est),
tem = seq(min(tem), max(tem), length.out = 50),
KEEP.OUT.ATTRS = FALSE))
grid <- cbind(grid,
as.data.frame(predict(m0,
newdata = grid,
interval = "confidence")))
str(grid)## 'data.frame': 150 obs. of 5 variables:
## $ est: Factor w/ 3 levels "PR","SC","SP": 3 1 2 3 1 2 3 1 2 3 ...
## $ tem: num 18 18 18 18.2 18.2 ...
## $ fit: num 0.0257 0.2135 0.4061 0.1936 0.3685 ...
## $ lwr: num -0.213 -0.0252 0.1674 -0.0162 0.1588 ...
## $ upr: num 0.264 0.452 0.645 0.403 0.578 ...# Função para obter o ponto de mínimo com avaliação a predição.
predfun <- function(x, est) {
-predict(m0, newdata = list(est = est, tem = x))
}
# Determinar o número mínimo de dias. Fazer minimização numérica.
opts <- sapply(levels(lat$est),
simplify = FALSE,
FUN = function(est) {
# Minimização numérica em uma dimensão.
opt <- optimise(f = predfun,
interval = c(18, 28),
est = est)
setNames(unlist(opt, use.names = FALSE),
nm = c("tem", "pi"))
}) %>%
do.call(what = rbind) %>%
as.data.frame() %>%
mutate(est = rownames(.))
opts## tem pi est
## 1 22.19972 -1.645343 PR
## 2 22.22259 -2.241017 SC
## 3 22.64792 -1.773904 SPp1 <-
xyplot(#log10(fi) ~ tem,
fi ~ tem,
scale = list(y = list(log = 10)),
yscale.components = yscale.components.log10.3,
groups = est,
xlab = "Temperature (°C)",
# ylab = expression(log[10] ~ "of infection frequency"),
ylab = expression("Number of pustules" ~ cm^{-2} ~ (log[10])),
jitter.x = TRUE,
auto.key = list(columns = 3),
data = lat,
ylim = c(NA, 290),
subset = fi > 0 & tem > 15) +
latticeExtra::as.layer({
xyplot(fit ~ tem,
type = "l",
groups = est,
cty = "bands",
ly = grid$lwr,
uy = grid$upr,
alpha = 0.45,
prepanel = prepanel.cbH,
panel = panel.superpose,
panel.groups = panel.cbH,
data = grid)
}, under = TRUE) +
latticeExtra::layer({
panel.segments(opts$tem,
-1,
opts$tem,
-opts$pi,
lty = 2)
})
p1
# p1 <-
# xyplot(fi ~ tem,
# groups = est,
# xlab = "Temperature (°C)",
# ylab = expression(log[10] ~ "of infection frequency"),
# jitter.x = TRUE,
# auto.key = list(columns = 3),
# data = lat,
# ylim = c(NA, 300),
# subset = fi > 0 & tem > 15) +
# latticeExtra::as.layer({
# xyplot(10^fit ~ tem,
# type = "l",
# groups = est,
# cty = "bands",
# ly = 10^grid$lwr,
# uy = 10^grid$upr,
# alpha = 0.45,
# prepanel = prepanel.cbH,
# panel = panel.superpose,
# panel.groups = panel.cbH,
# data = grid)
# }, under = TRUE) +
# latticeExtra::layer({
# panel.segments(opts$tem,
# -1,
# opts$tem,
# -opts$pi,
# lty = 2)
# })
# p1
# Teste de hipótese para verificar se o número de dias difere.
aux <- predict(m0, newdata = opts, interval = "confidence")
opts <- opts %>% cbind(aux) %>% mutate(est = factor(est))
# Criar a matriz X para aplicar o teste formal.
cfs <- attr(m0$model[[3]], "coefs")
X <- model.matrix(~est * poly(tem, degree = 3, coefs = cfs),
data = opts)
rownames(X) <- levels(lat$est)
# Avalia os contrastes par a par.
opts <- apmc(X, model = m0, focus = "est") %>%
select(est, cld) %>%
inner_join(opts)
opts## est cld tem pi fit lwr upr
## 1 PR b 22.19972 -1.645343 1.645343 1.459979 1.830708
## 2 SC a 22.22259 -2.241017 2.241017 2.056329 2.425705
## 3 SP b 22.64792 -1.773904 1.773904 1.599927 1.947880p2 <-
segplot(factor(est) ~ lwr + upr,
# data = opts,
data = mutate_at(opts, c("fit", "lwr", "upr"), ~10^.),
# scale = list(x = list(log = 10)),
# xscale.components = xscale.components.log10.3,
draw = FALSE,
cld = opts$cld,
centers = fit,
# xlab = expression(log[10] ~ "of infection frequency"),
xlab = expression("Number of pustules" ~ cm^{-2}),
ylab = "Origin") +
latticeExtra::layer({
panel.text(x = centers,
y = z,
labels = sprintf("%0.2f %s", centers, cld),
pos = 3)
})
p2
gridExtra::grid.arrange(p1, p2, ncol = 1, heights = c(2, 1.5))