#-----------------------------------------------------------------------
# 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á
# 2021-abr-15 · Curitiba/PR/Brazil
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
# Pacotes.
rm(list = objects())
library(readxl)
library(tidyverse)
library(emmeans)
#-----------------------------------------------------------------------
# Importação e preparo.
# Importação.
da <- read_excel(path = "factorial_audpc_pear_apple.xlsx")
str(da)
## tibble[,5] [153 × 5] (S3: tbl_df/tbl/data.frame)
## $ isol : chr [1:153] "mn16" "mn16" "mn16" "mn16" ...
## $ fruit: chr [1:153] "pear" "pear" "pear" "pear" ...
## $ temp : num [1:153] 5 5 5 5 5 5 10 10 10 10 ...
## $ rep : num [1:153] 1 2 3 4 5 6 1 2 3 4 ...
## $ AUDPC: num [1:153] 0 0 0 0 0 ...
# Conversão para fator.
da <- da %>%
mutate_at(c("isol", "fruit"), factor)
ggplot(data = da,
mapping = aes(x = temp, y = AUDPC, color = isol)) +
facet_wrap(facets = ~fruit, ncol = 1, scale = "free_y") +
geom_point()
Para maça, não houve crescimento para temperaturas diferentes de 25 °C. A única opção de análise é comparar os isolados na temperatura de 20 °C. Para a pera pode-se estudar o efeito de temperatura com ajuste de regressão.
# Aplica o filtro.
tb_apple <- da %>%
filter(fruit == "apple", temp == 25)
tb_apple
## # A tibble: 10 x 5
## isol fruit temp rep AUDPC
## <fct> <fct> <dbl> <dbl> <dbl>
## 1 mn16 apple 25 1 19.4
## 2 mn16 apple 25 2 28.7
## 3 mn16 apple 25 3 19.4
## 4 mn16 apple 25 4 53.1
## 5 mn16 apple 25 5 47.1
## 6 mn20 apple 25 1 0
## 7 mn20 apple 25 2 20.8
## 8 mn20 apple 25 3 5.10
## 9 mn20 apple 25 4 1.98
## 10 mn20 apple 25 5 0
# Análise com `sqrt()` da resposta é mais apropriada para atendimento
# dos pressupostos.
# m0 <- lm(AUDPC ~ isol, data = tb_apple)
m0 <- lm(sqrt(AUDPC) ~ isol, data = tb_apple)
# Análise de diagnóstico.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Médias na escala transformada.
emm <- emmeans(m0, specs = ~isol) %>%
multcomp::cld(Letters = letters)
emm
## isol emmean SE df lower.CL upper.CL .group
## mn20 1.65 0.736 8 -0.0523 3.34 a
## mn16 5.66 0.736 8 3.9655 7.36 b
##
## Results are given on the sqrt (not the response) scale.
## Confidence level used: 0.95
## Note: contrasts are still on the sqrt scale
## significance level used: alpha = 0.05
# Leva as estimativas para a escala natural.
emm_nat <- emm %>%
as.data.frame %>%
mutate_at(c("emmean", "lower.CL", "upper.CL"), ~ .^2) %>%
identity()
emm_nat
## isol emmean SE df lower.CL upper.CL .group
## 2 mn20 2.707939 0.736278 8 0.002733017 11.1786 a
## 1 mn16 32.073219 0.736278 8 15.724890485 54.1870 b
ggplot(data = emm_nat,
mapping = aes(x = isol,
y = emmean,
label = sprintf("%0.1f %s", emmean, .group),
ymin = lower.CL,
ymax = upper.CL)) +
geom_errorbar(width = 0.1) +
geom_point() +
geom_label(hjust = 0, nudge_x = 0.025) +
geom_point(data = tb_apple,
mapping = aes(x = isol, y = AUDPC),
position = position_nudge(x = -0.1),
color = "gray",
inherit.aes = FALSE) +
labs(x = "Isolates",
y = "Growth (mm)",
title = "Isolates colony growth at 25 °C for Apple",
caption = "Estimates follow by same letter are not different at 5% by LSD test.")
# Aplica o filtro.
tb_pear <- da %>%
filter(fruit == "pear")
tb_pear
## # A tibble: 83 x 5
## isol fruit temp rep AUDPC
## <fct> <fct> <dbl> <dbl> <dbl>
## 1 mn16 pear 5 1 0
## 2 mn16 pear 5 2 0
## 3 mn16 pear 5 3 0
## 4 mn16 pear 5 4 0
## 5 mn16 pear 5 5 0
## 6 mn16 pear 5 6 0
## 7 mn16 pear 10 1 30.0
## 8 mn16 pear 10 2 27.8
## 9 mn16 pear 10 3 17.0
## 10 mn16 pear 10 4 31.3
## # … with 73 more rows
ggplot(data = tb_pear,
mapping = aes(x = temp, y = AUDPC, color = isol)) +
geom_point() +
stat_summary(geom = "line", fun = "mean")
# ATTENTION: 1) remover temperaturas com valor 0 para a resposta porque
# isso prejudica a suposição de variância homogenea já que não há
# variância nesse caso ou 2) usar pesos na análise de regressão com
# pesos pequenos para os pontos com variância 0, pois isso evitar de
# eliminá-los. Aplicar alguma transformação na resposta para estabilizar
# a variância. Tem um valor bem baixo de crescimento na temperatura 25
# para o isolado mn16. Seria o caso de remover?
# Filtra para as temperaturas com crescimento e sem "outlier".
tb_pear <- da %>%
filter(fruit == "pear") %>%
filter(
# between(temp, left = 10, right = 25),
!(temp == 25 & AUDPC < 25)
)
ggplot(data = tb_pear,
mapping = aes(x = temp, y = sqrt(AUDPC), color = isol)) +
geom_point() +
stat_summary(geom = "line", fun = "mean")
# Peso 0.1 para pontos experimentais com variância 0. Essa escolha de
# peso é arbitrária e foi considerada para não eliminar as temperaturas
# da extremidade pois, apesar da variância 0, elas são úteis para dar
# forma de sino a curva de regressão ajustada. Então o peso reduz o
# efeito negativo da variância 0 sem precisar remover dados.
tb_pear <- tb_pear %>%
group_by(isol, temp) %>%
mutate(w = {
ifelse(var(AUDPC) == 0, 0.1, 1)
}) %>%
ungroup()
# Análise com `sqrt()` da resposta é mais apropriada para atendimento
# dos pressupostos.
# m0 <- lm(AUDPC ~ isol, data = tb_apple)
m0 <- lm(sqrt(AUDPC) ~ isol * poly(temp, degree = 3),
data = tb_pear, weights = w)
# Teste para a falta de ajuste do modelo de regressão.
m_full <- update(m0, formula = . ~ isol * factor(temp))
anova(m0, m_full)
## Analysis of Variance Table
##
## Model 1: sqrt(AUDPC) ~ isol * poly(temp, degree = 3)
## Model 2: sqrt(AUDPC) ~ isol + factor(temp) + isol:factor(temp)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 74 60.503
## 2 68 15.218 6 45.285 33.726 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Análise de diagnóstico.
par(mfrow = c(2, 2))
plot(m0)
# plot(m_full)
layout(1)
# Quadro de análise de variância.
anova(m0)
## Analysis of Variance Table
##
## Response: sqrt(AUDPC)
## Df Sum Sq Mean Sq F value Pr(>F)
## isol 1 0.82 0.817 0.9994 0.3207
## poly(temp, degree = 3) 3 494.07 164.691 201.4297 <2e-16 ***
## isol:poly(temp, degree = 3) 3 1.39 0.465 0.5683 0.6376
## Residuals 74 60.50 0.818
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m0)
##
## Call:
## lm(formula = sqrt(AUDPC) ~ isol * poly(temp, degree = 3), data = tb_pear,
## weights = w)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -2.3705 -0.6739 -0.2341 0.6066 2.3064
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.42947 0.23505 27.354 < 2e-16
## isolmn20 0.01978 0.33522 0.059 0.953
## poly(temp, degree = 3)1 4.06625 3.38218 1.202 0.233
## poly(temp, degree = 3)2 -31.72802 2.82678 -11.224 < 2e-16
## poly(temp, degree = 3)3 -11.93662 2.21592 -5.387 8.15e-07
## isolmn20:poly(temp, degree = 3)1 -1.97790 4.94179 -0.400 0.690
## isolmn20:poly(temp, degree = 3)2 -0.73609 4.07946 -0.180 0.857
## isolmn20:poly(temp, degree = 3)3 1.78982 3.14211 0.570 0.571
##
## (Intercept) ***
## isolmn20
## poly(temp, degree = 3)1
## poly(temp, degree = 3)2 ***
## poly(temp, degree = 3)3 ***
## isolmn20:poly(temp, degree = 3)1
## isolmn20:poly(temp, degree = 3)2
## isolmn20:poly(temp, degree = 3)3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9042 on 74 degrees of freedom
## Multiple R-squared: 0.8913, Adjusted R-squared: 0.8811
## F-statistic: 86.71 on 7 and 74 DF, p-value: < 2.2e-16
# NOTE: não há evidências sobre existência de interação
# isolado:temperatura. Dessa forma, os termos tem efeito aditivo.
# Modelo com termos aditivos.
m0 <- update(m0, formula = . ~ isol + poly(temp, degree = 3))
# Quadro de anova.
car::Anova(m0)
## Anova Table (Type II tests)
##
## Response: sqrt(AUDPC)
## Sum Sq Df F value Pr(>F)
## isol 0.15 1 0.1838 0.6693
## poly(temp, degree = 3) 494.07 3 204.8753 <2e-16 ***
## Residuals 61.90 77
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# NOTE: não há evidência sobre o efeito de isolado.
summary(m0)
##
## Call:
## lm(formula = sqrt(AUDPC) ~ isol + poly(temp, degree = 3), data = tb_pear,
## weights = w)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -2.2879 -0.6580 -0.2744 0.5651 2.5854
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.38672 0.20245 31.546 < 2e-16 ***
## isolmn20 0.09825 0.22916 0.429 0.669
## poly(temp, degree = 3)1 3.11200 2.44340 1.274 0.207
## poly(temp, degree = 3)2 -32.07800 2.01903 -15.888 < 2e-16 ***
## poly(temp, degree = 3)3 -10.95875 1.55529 -7.046 6.8e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8966 on 77 degrees of freedom
## Multiple R-squared: 0.8888, Adjusted R-squared: 0.8831
## F-statistic: 153.9 on 4 and 77 DF, p-value: < 2.2e-16
# range(tb_pear$temp)
grid <- with(tb_pear,
crossing(temp = seq(5, 35, length.out = 71),
isol = levels(isol)))
grid <- cbind(grid,
predict(m0, newdata = grid, interval = "confidence"))
grid <- grid %>%
mutate_at(c("fit", "lwr", "upr"), ~.^2)
head(grid)
## temp isol fit lwr upr
## 1 5.000000 mn16 4.329336 0.5749614 11.58142
## 2 5.000000 mn20 4.747865 0.7423046 12.22445
## 3 5.428571 mn16 5.045393 1.1032050 11.84774
## 4 5.428571 mn20 5.496442 1.3301004 12.50046
## 5 5.857143 mn16 5.886098 1.8154651 12.28407
## 6 5.857143 mn20 6.372505 2.1024175 12.95131
gg <-
ggplot(data = grid,
mapping = aes(x = temp,
y = fit,
color = isol,
fill = isol,
ymin = lwr,
ymax = upr)) +
geom_ribbon(alpha = 0.2,
show.legend = FALSE,
# color = NA,
size = 0.1) +
geom_line() +
geom_point(data = tb_pear,
mapping = aes(x = temp, y = AUDPC, color = isol),
inherit.aes = FALSE) +
labs(x = "Temperature",
y = "Growth (mm)",
color = "Isolates",
title = "Colony growth as a function of temperature for Pear")
gg
# Determinar o ponto de máximo da função de forma numérica.
fun_opt <- function(temp, isol = "mn16") {
-1 * predict(m0, newdata = list(temp = temp, isol = isol))
}
# Temperatura para o máximo crescimento. A temperatura de máximo será a
# mesma porque o efeito de isolado e temperatura é aditivo.
opt <- optimise(f = fun_opt,
lower = 20, upper = 30,
isol = "mn16",
maximum = FALSE)
opt
## $minimum
## [1] 23.6539
##
## $objective
## 1
## -11.50009
# opt <- nlm(f = fun_opt, p = 25, isol = "mn16", hessian = FALSE)
# opt <- opt[c("estimate", "minimum")]
# opt
gg +
geom_vline(xintercept = opt[[1]], color = "red") +
geom_label(data = data.frame(x = opt[[1]], y = 5),
mapping = aes(x = x,
y = y,
label = sprintf("%0.2f", x)),
inherit.aes = FALSE)