1 Definições da sessĂ£o
#-----------------------------------------------------------------------
# Prof. Dr. Walmes M. Zeviani
# leg.ufpr.br/~walmes · github.com/walmes
# walmes@ufpr.br · @walmeszeviani
# Laboratory of Statistics and Geoinformation (LEG)
# Department of Statistics · Federal University of ParanĂ¡
# 2020-jul-05 · Curitiba/PR/Brazil
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
# Pacotes.
rm(list = ls())
library(emmeans)
library(tidyverse)
library(readxl)
as_letters <- function(x) {
xs <- strsplit(trimws(as.character(x)), "")
xs <- lapply(xs, as.integer)
m <- max(unlist(xs))
sapply(xs, function(i) paste(rev(letters[m - i + 1]), collapse = ""))
}2 AnĂ¡lise para insetos
2.1 ImportaĂ§Ă£o e preparo
#-----------------------------------------------------------------------
# ImportaĂ§Ă£o da planilha.
# Importa.
tb <- gdata::read.xls("DADOS 2 Experimento.xlsx", sheet = 1, skip = 1, dec = ",") %>%
as_tibble()
# tb <- read_excel("DADOS 2 Experimento.xlsx", sheet = 1, skip = 1)
attr(tb, "spec") <- NULL
str(tb)## tibble [327 Ă— 7] (S3: tbl_df/tbl/data.frame)
## $ ensaio : int [1:327] 1 1 1 1 1 1 1 1 1 1 ...
## $ tratamento: chr [1:327] "NAC" "NAC" "NAC" "NAC" ...
## $ exposicao : int [1:327] 24 24 24 24 24 24 24 24 24 24 ...
## $ inseto : int [1:327] 1 2 3 4 5 6 7 8 9 10 ...
## $ teste : int [1:327] 1 1 1 1 1 1 2 2 2 2 ...
## $ qpcr : num [1:327] 27.1 30 31.1 34.4 33.5 ...
## $ morte : int [1:327] 0 0 0 0 0 1 0 0 1 0 ...# Declara os fatores.
tb <- tb %>%
mutate(expos = exposicao,
tratamento = trimws(tratamento)) %>%
mutate_at(c("tratamento", "exposicao", "ensaio"), factor)
# Tabela de frequĂªncia das condições experimentais (insetos).
# NĂºmero de insetos.
tb %>%
xtabs(formula = ~tratamento + exposicao + ensaio) %>%
ftable()## ensaio 1 2 3 4
## tratamento exposicao
## CONTROLE 24 9 24 27 15
## 72 6 6 15 15
## 120 3 6 15 9
## NAC 24 15 21 18 21
## 72 18 15 6 15
## 120 15 3 15 15cell_counts <- tb %>%
count(tratamento, exposicao, ensaio, teste) %>%
print(n = Inf)## # A tibble: 25 x 5
## tratamento exposicao ensaio teste n
## <fct> <fct> <fct> <int> <int>
## 1 CONTROLE 24 1 3 9
## 2 CONTROLE 24 2 9 24
## 3 CONTROLE 24 3 15 27
## 4 CONTROLE 24 4 21 15
## 5 CONTROLE 72 1 4 6
## 6 CONTROLE 72 2 11 6
## 7 CONTROLE 72 3 19 15
## 8 CONTROLE 72 4 23 15
## 9 CONTROLE 120 1 7 3
## 10 CONTROLE 120 2 12 6
## 11 CONTROLE 120 3 18 15
## 12 CONTROLE 120 4 25 9
## 13 NAC 24 1 1 6
## 14 NAC 24 1 2 9
## 15 NAC 24 2 8 21
## 16 NAC 24 3 14 18
## 17 NAC 24 4 20 21
## 18 NAC 72 1 5 18
## 19 NAC 72 2 10 15
## 20 NAC 72 3 16 6
## 21 NAC 72 4 22 15
## 22 NAC 120 1 6 15
## 23 NAC 120 2 13 3
## 24 NAC 120 3 17 15
## 25 NAC 120 4 24 15# OBS: com o ensaio, mudou-se a planta teste. O nĂºmero de insetos Ă©
# variĂ¡vel ao final. A anĂ¡lise irĂ¡ SUPOR que o nĂºmero de insetos nĂ£o
# depende das condições experimentais.
# NĂºmero de plantas teste.
cell_counts %>%
# xtabs(formula = n ~ tratamento + exposicao + ensaio) %>%
xtabs(formula = ~tratamento + exposicao + ensaio) %>%
ftable()## ensaio 1 2 3 4
## tratamento exposicao
## CONTROLE 24 1 1 1 1
## 72 1 1 1 1
## 120 1 1 1 1
## NAC 24 2 1 1 1
## 72 1 1 1 1
## 120 1 1 1 1# OBS: HĂ¡ praticamente uma planta por combinaĂ§Ă£o experimental em cada
# ensaio. Assim, o efeito de planta teste dentro de ensaio confunde-se
# com o efeito de interaĂ§Ă£o entre os fatores. NĂ£o serĂ¡ possĂvel declarar
# esse efeito. Se houvessem mais plantas teste com combinaĂ§Ă£o
# experimental, aĂ seria interessante.
ggplot(data = cell_counts,
mapping = aes(y = n,
x = exposicao,
color = tratamento,
group = tratamento)) +
facet_wrap(facets = ~ensaio) +
geom_point() +
geom_line()
ggplot(data = cell_counts,
mapping = aes(y = n,
x = exposicao,
color = tratamento,
group = tratamento)) +
# facet_wrap(facets = ~ensaio) +
geom_point() +
stat_summary(geom = "line", fun = "mean")
# OBS: O `n` cai com a exposiĂ§Ă£o, o que tem justificativa biolĂ³gica mas
# isso mostra que a frequĂªncia depende dos fatores experimentais. Pode
# haver vĂes.2.2 AnĂ¡lise exploratĂ³ria
#-----------------------------------------------------------------------
# AnĂ¡lise exploratĂ³ria das variĂ¡veis resposta.
# qPCR.
ggplot(data = tb,
mapping = aes(x = expos,
y = qpcr,
color = tratamento)) +
facet_wrap(facets = ~ensaio) +
# geom_point() +
geom_jitter(width = 3, height = 0) +
stat_summary(geom = "line", fun = "mean") +
# scale_y_log10() +
theme()
# Mortalidade.
ggplot(data = tb,
mapping = aes(x = expos,
y = morte,
color = tratamento)) +
facet_wrap(facets = ~ensaio) +
# geom_point() +
geom_jitter(width = 3, height = 0.1) +
stat_summary(geom = "line", fun = "mean") +
theme()
2.3 AnĂ¡lise das variĂ¡veis resposta
2.3.1 AnĂ¡lise para qPCR
#-----------------------------------------------------------------------
# AnĂ¡lise para qPCR.
# Dicotomiza.
tb$qpcr <- ifelse(tb$qpcr <= 30, 1, 0)
tb_means <- tb %>%
group_by(ensaio, tratamento, exposicao, expos) %>%
do({
bind_cols(
summarise_at(., c("qpcr", "morte"), mean),
summarise_at(., c("qpcr"), c(n = length)))
}) %>%
ungroup()
tb_means## # A tibble: 24 x 7
## ensaio tratamento exposicao expos qpcr morte n
## <fct> <fct> <fct> <int> <dbl> <dbl> <int>
## 1 1 CONTROLE 24 24 0 0.556 9
## 2 1 CONTROLE 72 72 0.167 0.167 6
## 3 1 CONTROLE 120 120 0.333 0.333 3
## 4 1 NAC 24 24 0.133 0.333 15
## 5 1 NAC 72 72 0.111 0.389 18
## 6 1 NAC 120 120 0.0667 0.6 15
## 7 2 CONTROLE 24 24 0 0.0417 24
## 8 2 CONTROLE 72 72 0.333 0.167 6
## 9 2 CONTROLE 120 120 0.167 0 6
## 10 2 NAC 24 24 0.238 0.190 21
## # … with 14 more rows# Modelo saturado.
m0 <- glm(qpcr ~ ensaio + tratamento * exposicao,
data = tb_means,
weights = n,
family = quasibinomial)
anova(m0, test = "F")## Analysis of Deviance Table
##
## Model: quasibinomial, link: logit
##
## Response: qpcr
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 23 172.032
## ensaio 3 108.137 20 63.895 11.7342 0.0003232
## tratamento 1 0.043 19 63.852 0.0140 0.9073788
## exposicao 2 8.350 17 55.501 1.3592 0.2867544
## tratamento:exposicao 2 5.078 15 50.424 0.8265 0.4565597
##
## NULL
## ensaio ***
## tratamento
## exposicao
## tratamento:exposicao
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Usa termo linear para efeito de exposiĂ§Ă£o.
m1 <- update(m0, . ~ ensaio + tratamento * expos)
anova(m0, m1, test = "F")## Analysis of Deviance Table
##
## Model 1: qpcr ~ ensaio + tratamento * exposicao
## Model 2: qpcr ~ ensaio + tratamento + expos + tratamento:expos
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 15 50.424
## 2 17 51.364 -2 -0.93987 0.153 0.8595anova(m1, test = "F")## Analysis of Deviance Table
##
## Model: quasibinomial, link: logit
##
## Response: qpcr
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 23 172.032
## ensaio 3 108.137 20 63.895 13.1053 0.0001111 ***
## tratamento 1 0.043 19 63.852 0.0156 0.9019515
## expos 1 8.045 18 55.807 2.9249 0.1054094
## tratamento:expos 1 4.443 17 51.364 1.6155 0.2208373
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Remove a interaĂ§Ă£o.
m2 <- update(m0, . ~ ensaio + tratamento + expos)
anova(m1, m2, test = "F")## Analysis of Deviance Table
##
## Model 1: qpcr ~ ensaio + tratamento + expos + tratamento:expos
## Model 2: qpcr ~ ensaio + tratamento + expos
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 17 51.364
## 2 18 55.807 -1 -4.4434 1.6155 0.2208anova(m2, test = "F")## Analysis of Deviance Table
##
## Model: quasibinomial, link: logit
##
## Response: qpcr
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 23 172.032
## ensaio 3 108.137 20 63.895 12.6379 0.0001103 ***
## tratamento 1 0.043 19 63.852 0.0151 0.9036271
## expos 1 8.045 18 55.807 2.8206 0.1103367
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# ObteĂ§Ă£o de valores preditos.
expos_seq <- with(tb_means, seq(min(expos),
max(expos),
length.out = 31))
# # Valores preditos na média dos ensaios.
# pred <- emmeans(m1,
# specs = ~expos + tratamento,
# at = list(expos = expos_seq)) %>%
# as.data.frame()
# # attributes(pred) <- NULL
# str(pred)
#
# # Curvas com bandas de confiança.
# ggplot(data = pred,
# mapping = aes(x = expos, y = emmean, color = tratamento)) +
# geom_ribbon(mapping = aes(ymin = lower.CL,
# ymax = upper.CL,
# fill = tratamento),
# alpha = 0.3,
# show.legend = FALSE) +
# geom_line() +
# labs(x = "PerĂodo de exposiĂ§Ă£o (h)",
# y = "qPCR",
# color = "Tratamento")
# Valores preditos na média dos ensaios.
pred <- emmeans(m2,
specs = ~expos + tratamento,
at = list(expos = expos_seq),
type = "response") %>%
as.data.frame()
# attributes(pred) <- NULL
# str(pred)
ggplot(data = pred,
mapping = aes(x = expos, y = prob, color = tratamento)) +
geom_ribbon(mapping = aes(ymin = asymp.LCL,
ymax = asymp.UCL,
fill = tratamento),
alpha = 0.3,
show.legend = FALSE) +
geom_line() +
scale_x_continuous(breaks = c(24, 72, 120)) +
labs(x = "ExposiĂ§Ă£o",
y = "Resultado do qPCR",
color = "Tratamento")
2.3.2 AnĂ¡lise para mortalidade
#-----------------------------------------------------------------------
# AnĂ¡lise para mortalidade.
# Modelo saturado.
m0 <- glm(morte ~ ensaio + tratamento * exposicao,
data = tb_means,
weights = n,
family = quasibinomial)
anova(m0, test = "F")## Analysis of Deviance Table
##
## Model: quasibinomial, link: logit
##
## Response: morte
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 23 54.930
## ensaio 3 18.3150 20 36.615 4.7274 0.01627 *
## tratamento 1 4.4612 19 32.154 3.4545 0.08281 .
## exposicao 2 4.5328 17 27.621 1.7550 0.20660
## tratamento:exposicao 2 5.9677 15 21.653 2.3105 0.13344
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Usa termo linear para efeito de exposiĂ§Ă£o.
m1 <- update(m0, . ~ ensaio + tratamento * expos)
anova(m0, m1, test = "F")## Analysis of Deviance Table
##
## Model 1: morte ~ ensaio + tratamento * exposicao
## Model 2: morte ~ ensaio + tratamento + expos + tratamento:expos
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 15 21.653
## 2 17 23.961 -2 -2.308 0.8936 0.4299anova(m1, test = "F")## Analysis of Deviance Table
##
## Model: quasibinomial, link: logit
##
## Response: morte
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 23 54.930
## ensaio 3 18.3150 20 36.615 4.8002 0.01337 *
## tratamento 1 4.4612 19 32.154 3.5077 0.07838 .
## expos 1 4.5237 18 27.630 3.5568 0.07651 .
## tratamento:expos 1 3.6687 17 23.961 2.8846 0.10766
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Remove a interaĂ§Ă£o.
m2 <- update(m0, . ~ ensaio + tratamento + expos)
anova(m1, m2, test = "F")## Analysis of Deviance Table
##
## Model 1: morte ~ ensaio + tratamento + expos + tratamento:expos
## Model 2: morte ~ ensaio + tratamento + expos
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 17 23.961
## 2 18 27.630 -1 -3.6687 2.8846 0.1077anova(m2, test = "F")## Analysis of Deviance Table
##
## Model: quasibinomial, link: logit
##
## Response: morte
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 23 54.930
## ensaio 3 18.3150 20 36.615 4.3859 0.01748 *
## tratamento 1 4.4612 19 32.154 3.2050 0.09025 .
## expos 1 4.5237 18 27.630 3.2499 0.08820 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# NOTE: Apenas efeitos significativos para o nĂvel de 10%.
# Valores preditos na média dos ensaios.
pred <- emmeans(m2,
specs = ~expos + tratamento,
at = list(expos = expos_seq),
type = "response") %>%
as.data.frame()
# attributes(pred) <- NULL
str(pred)## 'data.frame': 62 obs. of 7 variables:
## $ expos : num 24 27.2 30.4 33.6 36.8 40 43.2 46.4 49.6 52.8 ...
## $ tratamento: Factor w/ 2 levels "CONTROLE","NAC": 1 1 1 1 1 1 1 1 1 1 ...
## $ prob : num 0.183 0.186 0.189 0.193 0.196 ...
## $ SE : num 0.0429 0.0426 0.0423 0.0421 0.0419 ...
## $ df : num Inf Inf Inf Inf Inf ...
## $ asymp.LCL : num 0.113 0.116 0.12 0.123 0.126 ...
## $ asymp.UCL : num 0.282 0.284 0.286 0.288 0.291 ...
## - attr(*, "estName")= chr "prob"
## - attr(*, "clNames")= chr [1:2] "asymp.LCL" "asymp.UCL"
## - attr(*, "pri.vars")= chr [1:2] "expos" "tratamento"
## - attr(*, "adjust")= chr "none"
## - attr(*, "side")= num 0
## - attr(*, "delta")= num 0
## - attr(*, "type")= chr "response"
## - attr(*, "mesg")= chr [1:3] "Results are averaged over the levels of: ensaio" "Confidence level used: 0.95" "Intervals are back-transformed from the logit scale"ggplot(data = pred,
mapping = aes(x = expos, y = prob, color = tratamento)) +
geom_ribbon(mapping = aes(ymin = asymp.LCL,
ymax = asymp.UCL,
fill = tratamento),
alpha = 0.3,
show.legend = FALSE) +
geom_line() +
scale_x_continuous(breaks = c(24, 72, 120)) +
labs(x = "ExposiĂ§Ă£o",
y = "Mortalidade",
color = "Tratamento")
3 AnĂ¡lise para plantas teste
3.1 ImportaĂ§Ă£o e preparo
tb <- read_excel("DADOS 2 Experimento.xlsx", sheet = 2, skip = 1)
attr(tb, "spec") <- NULL
str(tb)## tibble [327 Ă— 6] (S3: tbl_df/tbl/data.frame)
## $ ensaio : num [1:327] 1 1 1 1 1 1 1 1 1 1 ...
## $ tratamento: chr [1:327] "NAC" "NAC" "NAC" "NAC" ...
## $ exposicao : num [1:327] 24 24 24 24 24 24 24 24 24 24 ...
## $ inseto : num [1:327] 1 2 3 4 5 6 7 8 9 10 ...
## $ teste : num [1:327] 1 1 1 2 2 2 3 3 3 4 ...
## $ qpcr : num [1:327] 0 0 0 1 1 1 0 0 0 0 ...# Declara os fatores.
tb <- tb %>%
mutate(expos = exposicao,
tratamento = trimws(tratamento)) %>%
mutate_at(c("tratamento", "exposicao", "ensaio"), factor)
# DANGER: por alguma razĂ£o os dados estĂ£o ao nĂvel de inseto mas
# duplicados. Precisa criar a variĂ¡vel ensaio.
# xtabs(~fonte, data = tb)
tb <- tb %>%
group_by(ensaio, tratamento, exposicao, teste) %>%
slice(1) %>%
ungroup()
# Tabela de frequĂªncia das condições experimentais (nĂºmero de plantas teste).
tb %>%
xtabs(formula = ~tratamento + exposicao + ensaio) %>%
ftable()## ensaio 1 2 3 4
## tratamento exposicao
## CONTROLE 24 3 8 9 5
## 72 2 2 5 5
## 120 1 2 5 3
## NAC 24 5 7 6 7
## 72 6 5 2 5
## 120 5 1 5 5ggplot(data = tb,
mapping = aes(x = teste, y = qpcr, color = tratamento)) +
facet_wrap(facets = ~ensaio, scale = "free_x") +
geom_point()
3.2 AnĂ¡lise
#-----------------------------------------------------------------------
# Ajuste do modelos aos dados de PCR.
m0 <- glm(qpcr ~ ensaio + tratamento * exposicao,
data = tb,
family = quasibinomial)
anova(m0, test = "F")## Analysis of Deviance Table
##
## Model: quasibinomial, link: logit
##
## Response: qpcr
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 108 75.582
## ensaio 3 17.8398 105 57.742 8.1827 6.34e-05
## tratamento 1 0.0566 104 57.686 0.0779 0.7808
## exposicao 2 0.9857 102 56.700 0.6782 0.5099
## tratamento:exposicao 2 2.3436 100 54.356 1.6124 0.2045
##
## NULL
## ensaio ***
## tratamento
## exposicao
## tratamento:exposicao
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Usa termo linear para efeito de exposiĂ§Ă£o.
m1 <- update(m0, . ~ ensaio + tratamento * expos)
anova(m0, m1, test = "F")## Analysis of Deviance Table
##
## Model 1: qpcr ~ ensaio + tratamento * exposicao
## Model 2: qpcr ~ ensaio + tratamento + expos + tratamento:expos
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 100 54.356
## 2 102 55.942 -2 -1.5852 1.0906 0.34anova(m1, test = "F")## Analysis of Deviance Table
##
## Model: quasibinomial, link: logit
##
## Response: qpcr
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 108 75.582
## ensaio 3 17.8398 105 57.742 9.1850 1.956e-05 ***
## tratamento 1 0.0566 104 57.686 0.0874 0.7681
## expos 1 0.2144 103 57.471 0.3312 0.5662
## tratamento:expos 1 1.5297 102 55.942 2.3627 0.1274
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Remove a interaĂ§Ă£o.
m2 <- update(m0, . ~ ensaio + tratamento + expos)
anova(m1, m2, test = "F")## Analysis of Deviance Table
##
## Model 1: qpcr ~ ensaio + tratamento + expos + tratamento:expos
## Model 2: qpcr ~ ensaio + tratamento + expos
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 102 55.942
## 2 103 57.471 -1 -1.5297 2.3627 0.1274anova(m2, test = "F")## Analysis of Deviance Table
##
## Model: quasibinomial, link: logit
##
## Response: qpcr
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 108 75.582
## ensaio 3 17.8398 105 57.742 7.4223 0.0001503 ***
## tratamento 1 0.0566 104 57.686 0.0706 0.7909322
## expos 1 0.2144 103 57.471 0.2676 0.6060297
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m2)##
## Call:
## glm(formula = qpcr ~ ensaio + tratamento + expos, family = quasibinomial,
## data = tb)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.04438 -0.40724 -0.33859 -0.00008 2.40317
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.362e-01 7.851e-01 -0.556 0.57966
## ensaio2 -1.925e+00 7.948e-01 -2.422 0.01718 *
## ensaio3 -2.105e+00 7.895e-01 -2.667 0.00889 **
## ensaio4 -1.899e+01 1.753e+03 -0.011 0.99138
## tratamentoNAC 2.159e-01 6.554e-01 0.329 0.74250
## expos -4.204e-03 8.201e-03 -0.513 0.60929
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 0.8011769)
##
## Null deviance: 75.582 on 108 degrees of freedom
## Residual deviance: 57.471 on 103 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 18# Modelo nulo.
m_null <- update(m0, . ~ ensaio)
summary(m_null)##
## Call:
## glm(formula = qpcr ~ ensaio, family = quasibinomial, data = tb)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.95077 -0.40837 -0.35927 -0.00008 2.35482
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.5596 0.3844 -1.456 0.14846
## ensaio2 -1.8827 0.7461 -2.523 0.01312 *
## ensaio3 -2.1484 0.7410 -2.899 0.00455 **
## ensaio4 -19.0065 1703.0556 -0.011 0.99112
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 0.752381)
##
## Null deviance: 75.582 on 108 degrees of freedom
## Residual deviance: 57.742 on 105 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 18# Modelo nulo contra modelo com termos experimentais.
anova(m_null, m2, test = "F")## Analysis of Deviance Table
##
## Model 1: qpcr ~ ensaio
## Model 2: qpcr ~ ensaio + tratamento + expos
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 105 57.742
## 2 103 57.471 2 0.27102 0.1691 0.8446# Estimativa de probabilidade de transmissĂ£o.
emmeans(m_null, specs = ~1, type = "response")## 1 prob SE df asymp.LCL asymp.UCL
## overall 0.0018 0.764 Inf 2.22e-16 1
##
## Results are averaged over the levels of: ensaio
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale# NOTE: Apenas efeitos significativos para o nĂvel de 10%.
# Valores preditos na média dos ensaios.
pred <- emmeans(m2,
specs = ~expos + tratamento,
at = list(expos = expos_seq),
type = "response") %>%
as.data.frame()
# attributes(pred) <- NULL
str(pred)## 'data.frame': 62 obs. of 7 variables:
## $ expos : num 24 27.2 30.4 33.6 36.8 40 43.2 46.4 49.6 52.8 ...
## $ tratamento: Factor w/ 2 levels "CONTROLE","NAC": 1 1 1 1 1 1 1 1 1 1 ...
## $ prob : num 0.00185 0.00182 0.0018 0.00177 0.00175 ...
## $ SE : num 0.808 0.797 0.787 0.776 0.766 ...
## $ df : num Inf Inf Inf Inf Inf ...
## $ asymp.LCL : num 2.22e-16 2.22e-16 2.22e-16 2.22e-16 2.22e-16 ...
## $ asymp.UCL : num 1 1 1 1 1 1 1 1 1 1 ...
## - attr(*, "estName")= chr "prob"
## - attr(*, "clNames")= chr [1:2] "asymp.LCL" "asymp.UCL"
## - attr(*, "pri.vars")= chr [1:2] "expos" "tratamento"
## - attr(*, "adjust")= chr "none"
## - attr(*, "side")= num 0
## - attr(*, "delta")= num 0
## - attr(*, "type")= chr "response"
## - attr(*, "mesg")= chr [1:3] "Results are averaged over the levels of: ensaio" "Confidence level used: 0.95" "Intervals are back-transformed from the logit scale"ggplot(data = pred,
mapping = aes(x = expos, y = prob, color = tratamento)) +
# geom_ribbon(mapping = aes(ymin = asymp.LCL,
# ymax = asymp.UCL,
# fill = tratamento),
# alpha = 0.3,
# show.legend = FALSE) +
geom_line() +
scale_x_continuous(breaks = c(24, 72, 120)) +
labs(x = "ExposiĂ§Ă£o",
y = "Resultado do qPCR",
color = "Tratamento")
#