#-----------------------------------------------------------------------
# Pacotes.
library(emmeans)
library(glmnet)
library(nlme)
library(tidyverse)
# Para não imprimir a matriz de correlação das estimativas.
# assignInNamespace("print.correlation",
# function(x, title) { return() },
# ns = "nlme")#-----------------------------------------------------------------------
# Importação e preparo dos dados.
# Importa da planilha.
tb <- readxl::read_xlsx("invivo_reorganized.xlsx",
sheet = 1,
na = c("na", ""),
guess_max = 5000)
str(tb)tibble [3,600 × 9] (S3: tbl_df/tbl/data.frame)
$ iso : chr [1:3600] "CrAaPR-40" "CrAaPR-40" "CrAaPR-40" "CrAaPR-40" ...
$ exp : num [1:3600] 1 1 1 1 1 1 1 1 1 1 ...
$ mol : num [1:3600] 36 36 36 36 36 24 24 24 24 24 ...
$ vaso : num [1:3600] 1 1 1 1 1 1 1 1 1 1 ...
$ ramo : num [1:3600] 1 1 1 1 1 1 1 1 1 1 ...
$ folha: num [1:3600] 1 2 3 4 5 1 2 3 4 5 ...
$ incid: num [1:3600] 0 0 0 0 0 0 0 0 0 0 ...
$ sev : num [1:3600] NA NA NA NA NA NA NA NA NA NA ...
$ time : num [1:3600] 12 12 12 12 12 12 12 12 12 12 ...xtabs(~iso + exp, data = tb) exp
iso 1 2
CrAaPR-01 450 450
CrAaPR-27 450 450
CrAaPR-40 450 450
CrAlPR-73 450 450xtabs(~iso + mol, data = tb) mol
iso 12 24 36
CrAaPR-01 300 300 300
CrAaPR-27 300 300 300
CrAaPR-40 300 300 300
CrAlPR-73 300 300 300tb |>
filter(time == 12) |>
xtabs(~exp + interaction(iso, mol, sep = "@") + vaso, data = _) |>
ftable() vaso 1 2 3 4 5 6 7 8 9 10
exp interaction(iso, mol, sep = "@")
1 CrAaPR-01@12 5 0 5 0 5 0 5 0 5 0
CrAaPR-27@12 5 0 5 0 5 0 5 0 5 0
CrAaPR-40@12 5 0 5 0 5 0 5 0 5 0
CrAlPR-73@12 0 5 0 5 0 5 0 5 0 5
CrAaPR-01@24 0 5 0 5 0 5 0 5 0 5
CrAaPR-27@24 0 5 0 5 0 5 0 5 0 5
CrAaPR-40@24 0 5 0 5 0 5 0 5 0 5
CrAlPR-73@24 5 0 5 0 5 0 5 0 5 0
CrAaPR-01@36 0 5 0 5 0 5 0 5 0 5
CrAaPR-27@36 0 5 0 5 0 5 0 5 0 5
CrAaPR-40@36 5 0 5 0 5 0 5 0 5 0
CrAlPR-73@36 5 0 5 0 5 0 5 0 5 0
2 CrAaPR-01@12 5 0 5 0 5 0 5 0 5 0
CrAaPR-27@12 5 0 5 0 5 0 5 0 5 0
CrAaPR-40@12 5 0 5 0 5 0 5 0 5 0
CrAlPR-73@12 0 5 0 5 0 5 0 5 0 5
CrAaPR-01@24 0 5 0 5 0 5 0 5 0 5
CrAaPR-27@24 0 5 0 5 0 5 0 5 0 5
CrAaPR-40@24 0 5 0 5 0 5 0 5 0 5
CrAlPR-73@24 5 0 5 0 5 0 5 0 5 0
CrAaPR-01@36 0 5 0 5 0 5 0 5 0 5
CrAaPR-27@36 0 5 0 5 0 5 0 5 0 5
CrAaPR-40@36 5 0 5 0 5 0 5 0 5 0
CrAlPR-73@36 5 0 5 0 5 0 5 0 5 0# Vaso é um bloco incompleto com 6 UE dos 12 tratamentos.
tb |>
filter(exp == 1, time == 12) |>
xtabs(~iso + mol + vaso, data = _) |>
ftable() vaso 1 2 3 4 5 6 7 8 9 10
iso mol
CrAaPR-01 12 5 0 5 0 5 0 5 0 5 0
24 0 5 0 5 0 5 0 5 0 5
36 0 5 0 5 0 5 0 5 0 5
CrAaPR-27 12 5 0 5 0 5 0 5 0 5 0
24 0 5 0 5 0 5 0 5 0 5
36 0 5 0 5 0 5 0 5 0 5
CrAaPR-40 12 5 0 5 0 5 0 5 0 5 0
24 0 5 0 5 0 5 0 5 0 5
36 5 0 5 0 5 0 5 0 5 0
CrAlPR-73 12 0 5 0 5 0 5 0 5 0 5
24 5 0 5 0 5 0 5 0 5 0
36 5 0 5 0 5 0 5 0 5 0# Das as folhas de um ramo são do mesmo tratamento (claro).
tb |>
filter(exp == 1, time == 12) |>
xtabs(~iso + mol + vaso + folha, data = _) |>
ftable() folha 1 2 3 4 5
iso mol vaso
CrAaPR-01 12 1 1 1 1 1 1
2 0 0 0 0 0
3 1 1 1 1 1
4 0 0 0 0 0
5 1 1 1 1 1
6 0 0 0 0 0
7 1 1 1 1 1
8 0 0 0 0 0
9 1 1 1 1 1
10 0 0 0 0 0
24 1 0 0 0 0 0
2 1 1 1 1 1
3 0 0 0 0 0
4 1 1 1 1 1
5 0 0 0 0 0
6 1 1 1 1 1
7 0 0 0 0 0
8 1 1 1 1 1
9 0 0 0 0 0
10 1 1 1 1 1
36 1 0 0 0 0 0
2 1 1 1 1 1
3 0 0 0 0 0
4 1 1 1 1 1
5 0 0 0 0 0
6 1 1 1 1 1
7 0 0 0 0 0
8 1 1 1 1 1
9 0 0 0 0 0
10 1 1 1 1 1
CrAaPR-27 12 1 1 1 1 1 1
2 0 0 0 0 0
3 1 1 1 1 1
4 0 0 0 0 0
5 1 1 1 1 1
6 0 0 0 0 0
7 1 1 1 1 1
8 0 0 0 0 0
9 1 1 1 1 1
10 0 0 0 0 0
24 1 0 0 0 0 0
2 1 1 1 1 1
3 0 0 0 0 0
4 1 1 1 1 1
5 0 0 0 0 0
6 1 1 1 1 1
7 0 0 0 0 0
8 1 1 1 1 1
9 0 0 0 0 0
10 1 1 1 1 1
36 1 0 0 0 0 0
2 1 1 1 1 1
3 0 0 0 0 0
4 1 1 1 1 1
5 0 0 0 0 0
6 1 1 1 1 1
7 0 0 0 0 0
8 1 1 1 1 1
9 0 0 0 0 0
10 1 1 1 1 1
CrAaPR-40 12 1 1 1 1 1 1
2 0 0 0 0 0
3 1 1 1 1 1
4 0 0 0 0 0
5 1 1 1 1 1
6 0 0 0 0 0
7 1 1 1 1 1
8 0 0 0 0 0
9 1 1 1 1 1
10 0 0 0 0 0
24 1 0 0 0 0 0
2 1 1 1 1 1
3 0 0 0 0 0
4 1 1 1 1 1
5 0 0 0 0 0
6 1 1 1 1 1
7 0 0 0 0 0
8 1 1 1 1 1
9 0 0 0 0 0
10 1 1 1 1 1
36 1 1 1 1 1 1
2 0 0 0 0 0
3 1 1 1 1 1
4 0 0 0 0 0
5 1 1 1 1 1
6 0 0 0 0 0
7 1 1 1 1 1
8 0 0 0 0 0
9 1 1 1 1 1
10 0 0 0 0 0
CrAlPR-73 12 1 0 0 0 0 0
2 1 1 1 1 1
3 0 0 0 0 0
4 1 1 1 1 1
5 0 0 0 0 0
6 1 1 1 1 1
7 0 0 0 0 0
8 1 1 1 1 1
9 0 0 0 0 0
10 1 1 1 1 1
24 1 1 1 1 1 1
2 0 0 0 0 0
3 1 1 1 1 1
4 0 0 0 0 0
5 1 1 1 1 1
6 0 0 0 0 0
7 1 1 1 1 1
8 0 0 0 0 0
9 1 1 1 1 1
10 0 0 0 0 0
36 1 1 1 1 1 1
2 0 0 0 0 0
3 1 1 1 1 1
4 0 0 0 0 0
5 1 1 1 1 1
6 0 0 0 0 0
7 1 1 1 1 1
8 0 0 0 0 0
9 1 1 1 1 1
10 0 0 0 0 0tb |>
filter(exp == 1, time == 12) |>
xtabs(~vaso + ramo, data = _) |>
ftable() ramo 1 2 3 4 5
vaso
1 30 0 0 0 0
2 30 0 0 0 0
3 0 30 0 0 0
4 0 30 0 0 0
5 0 0 30 0 0
6 0 0 30 0 0
7 0 0 0 30 0
8 0 0 0 30 0
9 0 0 0 0 30
10 0 0 0 0 30# NOTE: ramo tá representando algo equivocado aqui.#-----------------------------------------------------------------------
# Gráficos.
# Molhamento x incidência.
ggplot(data = tb,
mapping = aes(x = mol, y = incid, color = iso)) +
facet_grid(time ~ exp) +
# geom_point() +
geom_jitter(width = 1, height = 0.1) +
stat_summary(geom = "line", fun = "mean")
# Tempo x incidência.
ggplot(data = tb,
mapping = aes(x = time, y = incid, color = iso)) +
facet_grid(mol ~ exp) +
# geom_point() +
geom_jitter(width = 1, height = 0.1) +
stat_summary(geom = "line", fun = "mean")
# Efeito da posição da folha x molhamento.
ggplot(data = tb,
mapping = aes(x = time, y = incid)) +
facet_grid(mol ~ folha) +
# geom_point() +
geom_jitter(width = 1, height = 0.1) +
stat_summary(geom = "line", fun = "mean")
# Efeito da posição da folha x isolado.
ggplot(data = tb,
mapping = aes(x = time, y = incid)) +
facet_grid(iso ~ folha) +
# geom_point() +
geom_jitter(width = 1, height = 0.1) +
stat_summary(geom = "line", fun = "mean")
# Severidade.
ggplot(data = filter(tb, time == max(time)),
mapping = aes(x = mol, y = sev)) +
facet_grid(iso ~ exp) +
# geom_point() +
geom_jitter(width = 0.2, height = 0.1) +
stat_summary(geom = "line", fun = "mean")
#-----------------------------------------------------------------------
# Análise da severidade ------------------------------------------------
#-----------------------------------------------------------------------
# Criação da tabela para análise dos dados.
# Define os fatores.
tb <- tb |>
mutate(Iso = factor(iso),
Exp = factor(exp),
Mol = factor(mol),
Vaso = interaction(exp, vaso, sep = "_"),
Ramo = interaction(exp, vaso, iso, mol, sep = "_"),
Folha = factor(folha)) |>
droplevels()
str(tb)tibble [3,600 × 15] (S3: tbl_df/tbl/data.frame)
$ iso : chr [1:3600] "CrAaPR-40" "CrAaPR-40" "CrAaPR-40" "CrAaPR-40" ...
$ exp : num [1:3600] 1 1 1 1 1 1 1 1 1 1 ...
$ mol : num [1:3600] 36 36 36 36 36 24 24 24 24 24 ...
$ vaso : num [1:3600] 1 1 1 1 1 1 1 1 1 1 ...
$ ramo : num [1:3600] 1 1 1 1 1 1 1 1 1 1 ...
$ folha: num [1:3600] 1 2 3 4 5 1 2 3 4 5 ...
$ incid: num [1:3600] 0 0 0 0 0 0 0 0 0 0 ...
$ sev : num [1:3600] NA NA NA NA NA NA NA NA NA NA ...
$ time : num [1:3600] 12 12 12 12 12 12 12 12 12 12 ...
$ Iso : Factor w/ 4 levels "CrAaPR-01","CrAaPR-27",..: 3 3 3 3 3 4 4 4 4 4 ...
$ Exp : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
$ Mol : Factor w/ 3 levels "12","24","36": 3 3 3 3 3 2 2 2 2 2 ...
$ Vaso : Factor w/ 20 levels "1_1","2_1","1_2",..: 1 1 1 1 1 1 1 1 1 1 ...
$ Ramo : Factor w/ 120 levels "1_1_CrAaPR-01_12",..: 101 101 101 101 101 71 71 71 71 71 ...
$ Folha: Factor w/ 5 levels "1","2","3","4",..: 1 2 3 4 5 1 2 3 4 5 ...#-----------------------------------------------------------------------
# Qual a correlação entre as folhas?
tb |>
filter(time == max(time)) |>
mutate_at("sev", ~log(. + 1)) |>
lattice::xyplot(sev ~ folha | Ramo,
group = Iso, auto.key = TRUE,
# col = 1,
type = c("p", "l"), pch = 19,
data = _,
strip = FALSE,
as.table = TRUE)
tb_folhas <- tb |>
filter(time == max(time)) |>
filter(exp == 2) |>
select(Exp, Ramo, Folha, sev) |>
mutate_at("sev", ~log(. + 1)) |>
pivot_wider(names_from = Folha, values_from = sev) |>
select(`1`:last_col())
lattice::splom(tb_folhas, type = c("p", "r"),
pch = 19,
as.matrix = TRUE,
col = "black")
tb_folhas |>
cor(use = "pairwise.complete.obs") |>
round(4) |>
corrplot::corrplot(method = "ellipse")
# NOTE: As folhas de um mesmo ramo são bastante correlacionadas. Como
# não há hipótese para o efeito da folha, vamos trabalhar com valores
# agregados por ramo.
# Obtenção de um resumo da severidade por ramo.
tb_sev <- tb |>
filter(time == max(time)) |>
group_by(Exp, Iso, Mol, mol, Vaso, Ramo) |>
# summarize_at("sev", mean, na.rm = TRUE) |>
# summarize_at("sev", median, na.rm = TRUE) |>
summarize_at("sev", max, na.rm = TRUE) |>
ungroup()
str(tb_sev)tibble [120 × 7] (S3: tbl_df/tbl/data.frame)
$ Exp : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
$ Iso : Factor w/ 4 levels "CrAaPR-01","CrAaPR-27",..: 1 1 1 1 1 1 1 1 1 1 ...
$ Mol : Factor w/ 3 levels "12","24","36": 1 1 1 1 1 2 2 2 2 2 ...
$ mol : num [1:120] 12 12 12 12 12 24 24 24 24 24 ...
$ Vaso: Factor w/ 20 levels "1_1","2_1","1_2",..: 1 5 9 13 17 3 7 11 15 19 ...
$ Ramo: Factor w/ 120 levels "1_1_CrAaPR-01_12",..: 1 3 5 7 9 41 43 45 47 49 ...
$ sev : num [1:120] 2.5 0.15 0 3.5 6 3.5 3 0.3 8 8 ...tb_sev |>
select(Exp, Iso, Mol, sev) |>
summary() Exp Iso Mol sev
1:60 CrAaPR-01:30 12:40 Min. : 0.00
2:60 CrAaPR-27:30 24:40 1st Qu.: 3.00
CrAaPR-40:30 36:40 Median : 6.50
CrAlPR-73:30 Mean :14.63
3rd Qu.:25.00
Max. :55.00 # Confere o número de combinações.
xtabs(~Exp + Iso + Mol, data = tb_sev) |>
ftable() Mol 12 24 36
Exp Iso
1 CrAaPR-01 5 5 5
CrAaPR-27 5 5 5
CrAaPR-40 5 5 5
CrAlPR-73 5 5 5
2 CrAaPR-01 5 5 5
CrAaPR-27 5 5 5
CrAaPR-40 5 5 5
CrAlPR-73 5 5 5# NOTE: Vários testes foram feitos. A severidade máxima é a mais
# interessante.
# Define a escala da variável resposta.
# tb_sev$resp_sev <- tb_sev$sev
tb_sev$resp_sev <- log10(tb_sev$sev + 1)#-----------------------------------------------------------------------
# Experimento 1.
# Modelo de 1 estrato apenas para examinar os pressupostos.
m0 <- lm(resp_sev ~ Iso * Mol + Vaso,
data = filter(tb_sev, Exp == "1"))
# Exame dos pressupostos. Em caso de afastamento, tentar transformações
# para a resposta que reduzam violações.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Modelo de efeitos aleatórios.
m0 <- lme(resp_sev ~ Iso * Mol,
random = ~ 1 | Vaso,
data = filter(tb_sev, Exp == "1"),
method = "ML")
# Componentes de variância.
VarCorr(m0)Vaso = pdLogChol(1)
Variance StdDev
(Intercept) 0.004257444 0.06524909
Residual 0.034903253 0.18682412# Quadro de testes de Wald para os termos de efeito fixo.
anova(m0) numDF denDF F-value p-value
(Intercept) 1 39 693.3252 <.0001
Iso 3 39 79.6804 <.0001
Mol 2 39 27.2936 <.0001
Iso:Mol 6 39 0.3931 0.8789# Define o modelo com termos de efeito nulo abandonados.
m1 <- update(m0, . ~ . - Iso:Mol)
anova(m1) numDF denDF F-value p-value
(Intercept) 1 45 763.5690 <.0001
Iso 3 45 85.0680 <.0001
Mol 2 45 29.1325 <.0001# anova(m1, m0)
# Obtém as médias e aplica comparações múltiplas.
tb_means <-
emmeans(m1, specs = ~Iso) |>
multcomp::cld(Letters = letters, reverse = TRUE)
tb_means Iso emmean SE df lower.CL upper.CL .group
CrAaPR-40 1.456 0.0568 9 1.328 1.585 a
CrAlPR-73 1.226 0.0568 9 1.098 1.354 b
CrAaPR-01 0.613 0.0568 9 0.485 0.742 c
CrAaPR-27 0.442 0.0568 9 0.314 0.571 c
Results are averaged over the levels of: Mol
Degrees-of-freedom method: containment
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 4 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. # Gráfico com médias e IC.
tb_means |>
as.data.frame() |>
mutate_at(".group", trimws) |>
mutate_at(c("emmean", "lower.CL", "upper.CL"), ~10^(.) - 1) |>
ggplot(data = _,
mapping = aes(x = Iso, y = emmean)) +
geom_point() +
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
width = 0) +
geom_text(mapping = aes(label = sprintf("%0.1f %s", emmean, .group)),
nudge_x = 0.025, hjust = 0) +
labs(y = "Severidade",
x = "Isolado")
# Obtém as médias e aplica comparações múltiplas.
tb_means <-
emmeans(m1, specs = ~Mol) |>
multcomp::cld(Letters = letters, reverse = TRUE)
tb_means Mol emmean SE df lower.CL upper.CL .group
36 1.156 0.0501 9 1.043 1.270 a
24 0.979 0.0508 9 0.864 1.094 b
12 0.667 0.0508 9 0.552 0.782 c
Results are averaged over the levels of: Iso
Degrees-of-freedom method: containment
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 3 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. # Gráfico com médias e IC.
tb_means |>
as.data.frame() |>
mutate_at(".group", trimws) |>
mutate_at(c("emmean", "lower.CL", "upper.CL"), ~10^(.) - 1) |>
ggplot(data = _,
mapping = aes(x = Mol, y = emmean)) +
geom_point() +
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
width = 0) +
geom_text(mapping = aes(label = sprintf("%0.1f %s", emmean, .group)),
nudge_x = 0.025, hjust = 0) +
labs(y = "Severidade",
x = "Período de molhamento (h)")
#-----------------------------------------------------------------------
# Experimento 2.
# Modelo de 1 estrato apenas para examinar os pressupostos.
m0 <- lm(resp_sev ~ Iso * Mol + Vaso,
data = filter(tb_sev, Exp == "2"))
# Exame dos pressupostos. Em caso de afastamento, tentar transformações
# para a resposta que reduzam violações.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Modelo de efeitos aleatórios.
m0 <- lme(resp_sev ~ Iso * Mol,
random = ~ 1 | Vaso,
data = filter(tb_sev, Exp == "2"),
method = "ML")
# Componentes de variância.
VarCorr(m0)Vaso = pdLogChol(1)
Variance StdDev
(Intercept) 0.003122686 0.0558810
Residual 0.012845993 0.1133402# Quadro de testes de Wald para os termos de efeito fixo.
anova(m0) numDF denDF F-value p-value
(Intercept) 1 39 1248.1916 <.0001
Iso 3 39 224.6315 <.0001
Mol 2 39 136.8970 <.0001
Iso:Mol 6 39 1.1676 0.3433# Define o modelo com termos de efeito nulo abandonados.
m1 <- update(m0, . ~ . - Iso:Mol)
anova(m1) numDF denDF F-value p-value
(Intercept) 1 45 1277.7690 <.0001
Iso 3 45 218.7367 <.0001
Mol 2 45 133.4745 <.0001# anova(m1, m0)
# Obtém as médias e aplica comparações múltiplas.
tb_means <-
emmeans(m1, specs = ~Iso) |>
multcomp::cld(Letters = letters, reverse = TRUE)
tb_means Iso emmean SE df lower.CL upper.CL .group
CrAaPR-40 1.396 0.0386 9 1.309 1.483 a
CrAlPR-73 1.251 0.0386 9 1.163 1.338 b
CrAaPR-01 0.517 0.0386 9 0.430 0.604 c
CrAaPR-27 0.461 0.0386 9 0.374 0.548 c
Results are averaged over the levels of: Mol
Degrees-of-freedom method: containment
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 4 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. # Gráfico com médias e IC.
tb_means |>
as.data.frame() |>
mutate_at(".group", trimws) |>
mutate_at(c("emmean", "lower.CL", "upper.CL"), ~10^(.) - 1) |>
ggplot(data = _,
mapping = aes(x = Iso, y = emmean)) +
geom_point() +
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
width = 0) +
geom_text(mapping = aes(label = sprintf("%0.1f %s", emmean, .group)),
nudge_x = 0.025, hjust = 0) +
labs(y = "Severidade",
x = "Isolado")
# Obtém as médias e aplica comparações múltiplas.
tb_means <-
emmeans(m1, specs = ~Mol) |>
multcomp::cld(Letters = letters, reverse = TRUE)
tb_means Mol emmean SE df lower.CL upper.CL .group
36 1.157 0.0345 9 1.079 1.235 a
24 1.046 0.0352 9 0.966 1.126 b
12 0.516 0.0352 9 0.437 0.596 c
Results are averaged over the levels of: Iso
Degrees-of-freedom method: containment
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 3 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. # Gráfico com médias e IC.
tb_means |>
as.data.frame() |>
mutate_at(".group", trimws) |>
mutate_at(c("emmean", "lower.CL", "upper.CL"), ~10^(.) - 1) |>
ggplot(data = _,
mapping = aes(x = Mol, y = emmean)) +
geom_point() +
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
width = 0) +
geom_text(mapping = aes(label = sprintf("%0.1f %s", emmean, .group)),
nudge_x = 0.025, hjust = 0) +
labs(y = "Severidade",
x = "Período de molhamento (h)")
#-----------------------------------------------------------------------
# Análise da incidência ------------------------------------------------
#-----------------------------------------------------------------------
# Preparo dos dados.
tb |>
select(Exp, Iso, Mol, time, incid) |>
summary() Exp Iso Mol time incid
1:1800 CrAaPR-01:900 12:1200 Min. : 12 Min. :0.0000
2:1800 CrAaPR-27:900 24:1200 1st Qu.: 24 1st Qu.:0.0000
CrAaPR-40:900 36:1200 Median : 42 Median :0.0000
CrAlPR-73:900 Mean : 88 Mean :0.4019
3rd Qu.:120 3rd Qu.:1.0000
Max. :288 Max. :1.0000
NA's :10 # tb |>
# filter(exp == 2) |>
# lattice::xyplot(incid ~ sqrt(time) | Ramo,
# # group = Folha, auto.key = TRUE,
# col = 1,
# type = c("p", "a"), pch = 19,
# data = _,
# # strip = FALSE,
# as.table = TRUE)
# Não pode a incidência voltar para 0 depois que for 1.
tb |>
filter(Ramo == "2_6_CrAlPR-73_12") |>
select(folha, time, incid) |>
arrange(folha, time) |>
print(n = Inf)# A tibble: 30 × 3
folha time incid
<dbl> <dbl> <dbl>
1 1 12 0
2 1 24 1
3 1 36 1
4 1 48 1
5 1 120 1
6 1 288 1
7 2 12 0
8 2 24 1
9 2 36 1
10 2 48 1
11 2 120 1
12 2 288 1
13 3 12 0
14 3 24 1
15 3 36 1
16 3 48 1
17 3 120 1
18 3 288 1
19 4 12 0
20 4 24 1
21 4 36 1
22 4 48 1
23 4 120 1
24 4 288 0
25 5 12 0
26 5 24 1
27 5 36 1
28 5 48 1
29 5 120 1
30 5 288 0# Arrumar isso usando `cummax()`.
tb_incid <- tb |>
group_by(Ramo, folha) |>
arrange(time) |>
mutate(incid = cummax(incid)) |>
ungroup()
# Confere o resultado.
tb_incid |>
filter(exp == 2) |>
lattice::xyplot(incid ~ sqrt(time) | Ramo,
# group = Folha, auto.key = TRUE,
col = 1,
type = c("p", "a"), pch = 19,
data = _,
strip = FALSE,
as.table = TRUE)
# Contém NA na resposta?
tb_incid |>
count(incid)# A tibble: 3 × 2
incid n
<dbl> <int>
1 0 2107
2 1 1443
3 NA 50# Em quais ramos ocorrem?
tb_incid |>
filter(is.na(incid)) |>
count(Ramo)# A tibble: 2 × 2
Ramo n
<fct> <int>
1 1_2_CrAaPR-01_24 25
2 1_2_CrAaPR-27_24 25# # Como é o padrão?
# tb_incid |>
# filter(Ramo %in% c("1_2_CrAaPR-01_24", "1_2_CrAaPR-27_24")) |>
# select(Ramo, folha, time, incid) |>
# arrange(Ramo, folha, time) |>
# print(n = Inf)
# Abandona condições que só tiveram NA.
tb_incid <- tb_incid |>
filter(!Ramo %in% c("1_2_CrAaPR-01_24", "1_2_CrAaPR-27_24"))#-----------------------------------------------------------------------
# Determinar o tempo para incidencia de 50% usando modelo.
# Usar um modelo logístico pode não funcionar porque em muitos casos é
# linearmente separável. Então tem que contornar isso, por exemplo,
# usando regularização. Também é interessante usar uma escala que deixa
# mais regular o espaçamento no tempo. É possível também mudar a função
# de ligação. São possibilidades a serem exploradas.
##' plot(jitter(incid, 0.1) ~ jitter(log2(time), 0.2),
##' data = tb_incid)
##' plot(jitter(incid, 0.1) ~ jitter(sqrt(time), 0.2),
##' data = tb_incid)
# Define a escala da variável tempo para deixar o espaçamento mais
# regular.
tb_incid$f_time <- sqrt(tb_incid$time)
# library(glmnet)
# Cria a matriz de dados e a reposta usando todos os registros.
X <- model.matrix(~0 + Ramo + f_time, data = tb_incid)
y <- tb_incid$incid
table(y) |> prop.table()y
0 1
0.5923729 0.4076271 # Estimando com todos os dados (sem regularização).
coef(glm(incid ~ f_time, data = tb_incid, family = "binomial")) |>
setNames(c("b0", "b1")) |>
as.list() |>
with({
c(b0 = b0, b1 = b1, x = -b0/b1)
}) b0 b1 x
-2.0363590 0.2001821 10.1725351 ggplot(data = tb_incid,
mapping = aes(x = f_time, y = incid)) +
geom_jitter(height = 0.025, width = 0.1) +
geom_smooth(method = "glm", method.args = list(family = "binomial"))
# Determina o parâmetro lambda de regularização usando todos os dados.
cv_ridge <- glmnet::cv.glmnet(X, y, alpha = 0, family = "binomial")
cv_ridge$lambda.min[1] 0.02022836# Função para fazer o ajuste e retornar o período de incidencia.
get_incub_period <- function(tbl,
formula = y ~ x,
lambda = 0,
alpha = 1) {
X <- model.matrix(formula, data = tbl)
y <- tbl[[all.vars(formula)[1]]]
# Ajusta o modelo com regularização.
suppressWarnings({
m0 <- glmnet::glmnet(X, y,
family = "binomial",
alpha = alpha,
lambda = lambda)
})
# Extrai e renomeia os coeficientes.
m0_coef <- coef(m0) |>
unclass() |>
attr("x") |>
setNames(c("b0", "b1"))
# Valor em x que corresponde a \hat{p} = 0.5.
x <- unname(-m0_coef["b0"]/m0_coef["b1"])
return(append(m0_coef, c(x = x)))
}# Calcula o período de incidência para todos os ramos.
tb_incub <-
tb_incid |>
group_by(Ramo) |>
summarise(
prop_incid = mean(incid),
params = list({
params <- try(
silent = TRUE,
expr = {
get_incub_period(tbl = data.frame(time, incid),
formula = incid ~ f_time,
lambda = cv_ridge$lambda.min)
})
if (class(params) == "try-error") {
tibble(b0 = NA_real_, b1 = NA_real_,
x = NA_real_)
} else {
tibble(b0 = params["b0"], b1 = params["b1"],
x = params["x"])
}
})) |>
unnest(params) |>
mutate(incub_period = x^2)
tb_incub# A tibble: 118 × 6
Ramo prop_incid b0 b1 x incub_period
<fct> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1_1_CrAaPR-01_12 0.133 -4.47 0.255 17.5 308.
2 2_1_CrAaPR-01_12 0.133 -4.47 0.255 17.5 308.
3 1_3_CrAaPR-01_12 0.0667 -4.58 0.191 23.9 572.
4 2_3_CrAaPR-01_12 0.133 -4.47 0.255 17.5 308.
5 1_5_CrAaPR-01_12 0 NA NA NA NA
6 2_5_CrAaPR-01_12 0.133 -4.47 0.255 17.5 308.
7 1_7_CrAaPR-01_12 0.167 -5.03 0.331 15.2 231.
8 2_7_CrAaPR-01_12 0.0667 -7.24 0.384 18.8 354.
9 1_9_CrAaPR-01_12 0.3 -8.50 0.847 10.0 101.
10 2_9_CrAaPR-01_12 0.2 -4.57 0.321 14.2 202.
# ℹ 108 more rowstb_grid <- with(tb_incub,
data.frame(f_time = seq(0, max(x, na.rm = TRUE),
length.out = 51)))
str(tb_grid)'data.frame': 51 obs. of 1 variable:
$ f_time: num 0 0.479 0.957 1.436 1.914 ...tb_curves <- tb_incub |>
rowwise() |>
do({
data.frame(Ramo = .$Ramo,
tb_grid,
p = binomial()$linkinv(.$b0 + .$b1 * tb_grid$f_time))
}) |>
ungroup()
str(tb_curves)tibble [6,018 × 3] (S3: tbl_df/tbl/data.frame)
$ Ramo : Factor w/ 120 levels "1_1_CrAaPR-01_12",..: 1 1 1 1 1 1 1 1 1 1 ...
$ f_time: num [1:6018] 0 0.479 0.957 1.436 1.914 ...
$ p : num [1:6018] 0.0114 0.0128 0.0144 0.0163 0.0184 ...# Faz o gráfico com as curvas.
ggplot(data = tb_incid,
mapping = aes(x = f_time, y = incid)) +
facet_wrap(facets = ~Ramo) +
geom_jitter(width = 0.1, height = 0.05,
color = "gray",
pch = 19,
size = 0.75) +
geom_line(data = tb_curves,
mapping = aes(x = f_time, y = p),
inherit.aes = FALSE) +
geom_vline(data = tb_incub,
mapping = aes(xintercept = x),
linewidth = 0.25,
color = "red") +
geom_hline(yintercept = 0.5, color = "red", linewidth = 0.25) +
# theme_minimal() +
# theme_void() +
theme_bw() +
theme(strip.background = element_blank(),
strip.text.x = element_blank(),
# panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank())
# IMPORTANT: Aqui foi abandonado as situações em que não houve estimação
# de período de incubação, ou seja, onde estimou NA. Acontece que são
# situações onde o período de incubação é bem alto e não foi possível
# estimar a partir dos dados porque deveria se observar por mais tempo.
# Então na realidade o período de incubação ali é grande, porém
# desconhecido. Abandonar essas unidades experimentais leva a
# subestimação do período de incubação. Seria interessante empregar um
# modelo misto para resposta Gaussiana com cersura. Isso vai dificultar
# a análise, então iremos abandonar os NA. Para estudos futuros, ajustar
# o tempo de observação do experimento para evitar falhas na estimação
# do período de incubação.#-----------------------------------------------------------------------
# Análise exploratória.
# Estima as proporções de incidência em cada tempo.
tb_incid |>
group_by(Ramo, time) |>
summarise(prop_incid = mean(incid), .groups = "drop") |>
# ungroup() |>
pivot_wider(names_from = time, values_from = prop_incid) |>
left_join(tb_incub, by = "Ramo")# A tibble: 118 × 12
Ramo `12` `24` `36` `48` `120` `288` prop_incid b0 b1 x
<fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1_1_CrAaPR… 0 0 0 0 0.4 0.4 0.133 -4.47 0.255 17.5
2 2_1_CrAaPR… 0 0 0 0 0.4 0.4 0.133 -4.47 0.255 17.5
3 1_3_CrAaPR… 0 0 0 0 0.2 0.2 0.0667 -4.58 0.191 23.9
4 2_3_CrAaPR… 0 0 0 0 0.4 0.4 0.133 -4.47 0.255 17.5
5 1_5_CrAaPR… 0 0 0 0 0 0 0 NA NA NA
6 2_5_CrAaPR… 0 0 0 0 0.4 0.4 0.133 -4.47 0.255 17.5
7 1_7_CrAaPR… 0 0 0 0 0.4 0.6 0.167 -5.03 0.331 15.2
8 2_7_CrAaPR… 0 0 0 0 0 0.4 0.0667 -7.24 0.384 18.8
9 1_9_CrAaPR… 0 0 0 0 0.8 1 0.3 -8.50 0.847 10.0
10 2_9_CrAaPR… 0 0 0 0 0.6 0.6 0.2 -4.57 0.321 14.2
# ℹ 108 more rows
# ℹ 1 more variable: incub_period <dbl># NOTE: Foram poucas os ramos que tiveram incidência observada superior
# a 50%. Portanto, por estimativas amostrais, seria bem mais complicado.
# Junta os dados.
tb_incub <-
left_join(tb_incub,
distinct(tb_incid, Ramo, Vaso, Exp, Iso, Mol, mol),
by = "Ramo")
# Análise exploratória.
tb_incub |>
drop_na(incub_period) |>
ggplot(data = _,
mapping = aes(x = mol, y = incub_period)) +
facet_grid(facets = Exp ~ Iso) +
# geom_point() +
geom_jitter(width = 1, height = 1) +
stat_summary(geom = "line", fun = "mean")
# Define a escala em que a resposta será analisada.
tb_incub$resp_incub <- tb_incub$incub_period
# tb_incub$resp_incub <- log10(tb_incub$incub_period)#-----------------------------------------------------------------------
# Análise para o experimento 1.
# Modelo de 1 estrato apenas para examinar os pressupostos.
m0 <- lm(resp_incub ~ Iso * Mol + Vaso,
data = filter(tb_incub, Exp == "1"))
# Exame dos pressupostos. Em caso de afastamento, tentar transformações
# para a resposta que reduzam violações.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Modelo de efeitos aleatórios.
m0 <- lme(resp_incub ~ Iso * Mol,
random = ~ 1 | Vaso,
data = filter(tb_incub, Exp == "1") |> drop_na(incub_period),
method = "ML")
# Componentes de variância.
VarCorr(m0)Vaso = pdLogChol(1)
Variance StdDev
(Intercept) 1785.834 42.25913
Residual 11197.298 105.81728# Quadro de testes de Wald para os termos de efeito fixo.
anova(m0) numDF denDF F-value p-value
(Intercept) 1 30 42.07280 <.0001
Iso 3 30 13.58899 <.0001
Mol 2 30 0.68575 0.5114
Iso:Mol 6 30 0.82838 0.5573# Define o modelo com termos de efeito nulo abandonados.
m1 <- update(m0, . ~ Iso)
anova(m1) numDF denDF F-value p-value
(Intercept) 1 38 50.00417 <.0001
Iso 3 38 13.67980 <.0001# anova(m1, m0)
# Obtém as médias e aplica comparações múltiplas.
tb_means <-
emmeans(m1, specs = ~Iso) |>
multcomp::cld(Letters = letters, reverse = TRUE)
tb_means Iso emmean SE df lower.CL upper.CL .group
CrAaPR-27 305.2 45.4 9 202.4 408.0 a
CrAaPR-01 258.2 36.1 9 176.7 339.8 a
CrAlPR-73 97.8 33.7 9 21.6 174.1 b
CrAaPR-40 23.4 33.7 9 -52.9 99.6 b
Degrees-of-freedom method: containment
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 4 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. # Gráfico com médias e IC.
tb_means |>
as.data.frame() |>
mutate_at(".group", trimws) |>
ggplot(data = _,
mapping = aes(x = Iso, y = emmean)) +
geom_point() +
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
width = 0) +
geom_text(mapping = aes(label = sprintf("%0.1f %s", emmean, .group)),
nudge_x = 0.025, hjust = 0) +
# scale_y_log10() +
labs(y = "Período de incubação (dias)",
x = "Isolado")
#-----------------------------------------------------------------------
# Análise para o experimento 2.
# Modelo de 1 estrato apenas para examinar os pressupostos.
m0 <- lm(resp_incub ~ Iso * Mol + Vaso,
data = filter(tb_incub, Exp == "2"))
# Exame dos pressupostos. Em caso de afastamento, tentar transformações
# para a resposta que reduzam violações.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Modelo de efeitos aleatórios.
m0 <- lme(resp_incub ~ Iso * Mol,
random = ~ 1 | Vaso,
data = filter(tb_incub, Exp == "2") |> drop_na(incub_period),
method = "ML")
# Componentes de variância.
VarCorr(m0)Vaso = pdLogChol(1)
Variance StdDev
(Intercept) 6.417619e-04 0.02533302
Residual 4.128004e+03 64.24954116# Quadro de testes de Wald para os termos de efeito fixo.
anova(m0) numDF denDF F-value p-value
(Intercept) 1 39 286.72590 <.0001
Iso 3 39 51.41380 <.0001
Mol 2 39 1.81693 0.1760
Iso:Mol 6 39 1.08589 0.3878# Define o modelo com termos de efeito nulo abandonados.
m1 <- update(m0, . ~ Iso)
anova(m1) numDF denDF F-value p-value
(Intercept) 1 47 232.13330 <.0001
Iso 3 47 51.35395 <.0001# anova(m1, m0)
# Obtém as médias e aplica comparações múltiplas.
tb_means <-
emmeans(m1, specs = ~Iso) |>
multcomp::cld(Letters = letters, reverse = TRUE)
tb_means Iso emmean SE df lower.CL upper.CL .group
CrAaPR-27 274.2 19.1 9 231.0 317.5 a
CrAaPR-01 270.3 19.1 9 227.0 313.6 a
CrAlPR-73 63.7 19.1 9 20.4 106.9 b
CrAaPR-40 19.9 19.1 9 -23.4 63.2 b
Degrees-of-freedom method: containment
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 4 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. # Gráfico com médias e IC.
tb_means |>
as.data.frame() |>
mutate_at(".group", trimws) |>
ggplot(data = _,
mapping = aes(x = Iso, y = emmean)) +
geom_point() +
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
width = 0) +
geom_text(mapping = aes(label = sprintf("%0.1f %s", emmean, .group)),
nudge_x = 0.025, hjust = 0) +
# scale_y_log10() +
labs(y = "Período de incubação (dias)",
x = "Isolado")
#-----------------------------------------------------------------------