Análise de área abaixo da curva

Jhulia Gelain & Walmes Marques Zeviani

2021-04-15

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á
#                                       2021-abr-15 · Curitiba/PR/Brazil
#-----------------------------------------------------------------------

#-----------------------------------------------------------------------
# Pacotes.

rm(list = objects())

library(readxl)
library(tidyverse)
library(emmeans)

2 Importação e preparo

#-----------------------------------------------------------------------
# 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.

3 Maça

# 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.")

4 Pera

# 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)