AUCPD em função de temperatura e molhamento para experimento exvivo

1 Carregamento de pacotes

Código

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

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

# Área abaixo da curva.
auc <- function(time, resp) {
    i <- !is.na(resp)
    time <- time[i]
    resp <- resp[i]
    alt <- diff(time)
    bas <- resp[-length(resp)] + diff(resp)/2
    a <- sum(alt * bas)
    return(a)
}

# Contornos em gráficos 3D.
source("https://raw.githubusercontent.com/walmes/wzRfun/master/R/panel.3d.contour.R")

2 Importação e preparo dos dados

Código

#-----------------------------------------------------------------------
# Importação e preparo dos dados.

tb <- read_excel("exvivo.xlsx", na = "NA")
str(tb)

tibble [4,480 × 7] (S3: tbl_df/tbl/data.frame)
 $ exp : num [1:4480] 1 1 1 1 1 1 1 1 1 1 ...
 $ iso : chr [1:4480] "CrAaPR-01" "CrAaPR-01" "CrAaPR-01" "CrAaPR-01" ...
 $ rept: num [1:4480] 1 2 3 4 5 1 2 3 4 5 ...
 $ tem : num [1:4480] 15 15 15 15 15 15 15 15 15 15 ...
 $ mol : num [1:4480] 6 6 6 6 6 12 12 12 12 12 ...
 $ sev : num [1:4480] 0 0 0 0 0 0 0 0 0 0 ...
 $ time: num [1:4480] 12 12 12 12 12 12 12 12 12 12 ...

Código

# Contagem das combinações experimentais.
xtabs(~exp + iso + tem + mol, data = tb) |>
    ftable()

                  mol  6 12 24 48
exp iso       tem                
1   CrAaPR-01 15      35 35 35 35
              21      35 35 35 35
              27      35 35 35 35
              33      35 35 35 35
    CrAaPR-27 15      35 35 35 35
              21      35 35 35 35
              27      35 35 35 35
              33      35 35 35 35
    CrAaPR-40 15      35 35 35 35
              21      35 35 35 35
              27      35 35 35 35
              33      35 35 35 35
    CrAlPR-73 15      35 35 35 35
              21      35 35 35 35
              27      35 35 35 35
              33      35 35 35 35
2   CrAaPR-01 15      35 35 35 35
              21      35 35 35 35
              27      35 35 35 35
              33      35 35 35 35
    CrAaPR-27 15      35 35 35 35
              21      35 35 35 35
              27      35 35 35 35
              33      35 35 35 35
    CrAaPR-40 15      35 35 35 35
              21      35 35 35 35
              27      35 35 35 35
              33      35 35 35 35
    CrAlPR-73 15      35 35 35 35
              21      35 35 35 35
              27      35 35 35 35
              33      35 35 35 35

Código

#-----------------------------------------------------------------------
# Análise gráfica exploratória.

ggplot(data = tb,
       mapping = aes(x = time,
                     y = sev,
                     group = iso,
                     color = iso)) +
    geom_point() +
    facet_grid(facets = tem ~ mol) +
    stat_summary(fun = mean, geom = "line")

Código

gg1 <-
    ggplot(data = tb,
           mapping = aes(x = tem, y = sev, col = as.factor(mol))) +
    geom_point() +
    stat_summary(fun = mean, geom = "line") +
    facet_grid(facets = exp ~ iso) +
    labs(title = "A",
         y = "Severity (%)",
         x = "Temperature (°C)",
         color = "Wetness period (h)") +
    guides(color = guide_legend(nrow = 1, title.position = "left")) +
    theme(legend.position = "top")

gg2 <-
    ggplot(data = tb,
           mapping = aes(x = mol, y = sev, col = as.factor(tem))) +
    geom_point() +
    stat_summary(fun.y = mean, geom = "line") +
    facet_grid(facets = exp ~ iso) +
    labs(title = "B",
         y = "Severity (%)",
         x = "Wetness period (h)",
         color = "Temperature (°C)") +
    guides(color = guide_legend(nrow = 1, title.position = "left")) +
    theme(legend.position = "top")

gridExtra::grid.arrange(gg1, gg2, ncol = 1)

Código

#-----------------------------------------------------------------------
# Cálculo da área abaixo da curva.

tb_auc <-
    tb |>
    group_by(exp, iso, tem, mol, rept) |>
    arrange(time) |>
    summarise(AUCPD = auc(time, sev)) |>
    ungroup()
str(tb_auc)

tibble [640 × 6] (S3: tbl_df/tbl/data.frame)
 $ exp  : num [1:640] 1 1 1 1 1 1 1 1 1 1 ...
 $ iso  : chr [1:640] "CrAaPR-01" "CrAaPR-01" "CrAaPR-01" "CrAaPR-01" ...
 $ tem  : num [1:640] 15 15 15 15 15 15 15 15 15 15 ...
 $ mol  : num [1:640] 6 6 6 6 6 12 12 12 12 12 ...
 $ rept : num [1:640] 1 2 3 4 5 1 2 3 4 5 ...
 $ AUCPD: num [1:640] 0 0 0 0 0 0 0 0 0 0 ...

Código

# Declara variáveis como fatores.
tb_auc <- tb_auc |>
    mutate_at(c("exp", "iso"), factor) |>
    rename(Exp = exp, Iso = iso) |>
    mutate(Mol = factor(mol),
           Tem = factor(tem))
str(tb_auc)

tibble [640 × 8] (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 ...
 $ tem  : num [1:640] 15 15 15 15 15 15 15 15 15 15 ...
 $ mol  : num [1:640] 6 6 6 6 6 12 12 12 12 12 ...
 $ rept : num [1:640] 1 2 3 4 5 1 2 3 4 5 ...
 $ AUCPD: num [1:640] 0 0 0 0 0 0 0 0 0 0 ...
 $ Mol  : Factor w/ 4 levels "6","12","24",..: 1 1 1 1 1 2 2 2 2 2 ...
 $ Tem  : Factor w/ 4 levels "15","21","27",..: 1 1 1 1 1 1 1 1 1 1 ...

Código

# Frequência das combinações experimentais.
xtabs(~Exp + Iso + tem + mol, data = tb_auc) |>
    ftable()

                  mol 6 12 24 48
Exp Iso       tem               
1   CrAaPR-01 15      5  5  5  5
              21      5  5  5  5
              27      5  5  5  5
              33      5  5  5  5
    CrAaPR-27 15      5  5  5  5
              21      5  5  5  5
              27      5  5  5  5
              33      5  5  5  5
    CrAaPR-40 15      5  5  5  5
              21      5  5  5  5
              27      5  5  5  5
              33      5  5  5  5
    CrAlPR-73 15      5  5  5  5
              21      5  5  5  5
              27      5  5  5  5
              33      5  5  5  5
2   CrAaPR-01 15      5  5  5  5
              21      5  5  5  5
              27      5  5  5  5
              33      5  5  5  5
    CrAaPR-27 15      5  5  5  5
              21      5  5  5  5
              27      5  5  5  5
              33      5  5  5  5
    CrAaPR-40 15      5  5  5  5
              21      5  5  5  5
              27      5  5  5  5
              33      5  5  5  5
    CrAlPR-73 15      5  5  5  5
              21      5  5  5  5
              27      5  5  5  5
              33      5  5  5  5

Código

# Gráfico para mostrar a AUCPD em função de mol e tem.
ggplot(data = tb_auc,
       mapping = aes(x = mol, y = AUCPD, col = Tem)) +
    geom_point() +
    stat_summary(fun = mean, geom = "line") +
    facet_grid(facets = Exp ~ Iso) +
    labs(y = "Severity AUDPC",
        x = "Wetness period (h)",
         color = "Temperature (°C)") +
    guides(color = guide_legend(nrow = 1, title.position = "left")) +
    theme(legend.position = "top")

Código

# Gráfico para mostrar a AUCPD em função de mol e tem.
ggplot(data = tb_auc,
       mapping = aes(x = tem, y = AUCPD, col = Mol)) +
    geom_point() +
    stat_summary(fun = mean, geom = "line") +
    facet_grid(facets = Exp ~ Iso) +
    labs(y = "Severity AUDPC",
         color = "Wetness period (h)",
         x = "Temperature (°C)") +
    guides(color = guide_legend(nrow = 1, title.position = "left")) +
    theme(legend.position = "top")

Código

# IMPORTANT: Não existe um claro sinal de que a resposta como função da
# temperatura seja quadrática. O padrão desejável só é observado para
# dois isolados. Da mesma forma, não se observa padrão consistente com
# alguma função contínua para o efeito de molhamento. Dito isto,
# espera-se falta de ajuste ao declarar o modelo de segundo grau
# completo na temperatura e molhamento.

3 Análise do experimento 1

Código

#-----------------------------------------------------------------------
# Experimento 1 --------------------------------------------------------

tbi <- tb_auc |>
    filter(Exp == "1")

# Confere qual o intervalo de valores.
summary(tbi$AUCPD)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0     0.0   982.2  2248.3  4534.4  6978.0 

Código

# Soma uma constante para evitar valores negativos.
m0 <- lm((AUCPD + 500) ~ Tem * Iso * Mol, data = tbi)

# Exame dos pressupostos.
par(mfrow = c(2, 2))
plot(m0)

Código

layout(1)

# Transformação Box-Cox.
MASS::boxcox(m0)
abline(v = 0, col = "red")

Código

# Ajuste com os dados transformados.
m1 <- lm(log(AUCPD + 500) ~ Tem * Iso * Mol, data = tbi)

# Confere os pressupostos.
par(mfrow = c(2, 2))
plot(m1)

Código

layout(1)

#-----------------------------------------------------------------------
# Modelo que descreve uma superfície para temperatura e molhamento.

# Refina o modelo para descrever a superfície de resposta.
m2 <- lm(log(AUCPD + 500) ~
             Iso * (tem + I(tem^2) + mol + I(mol^0.5) + tem:mol),
         data = tbi)

# Teste da falta de ajuste.
anova(m2, m1)

Analysis of Variance Table

Model 1: log(AUCPD + 500) ~ Iso * (tem + I(tem^2) + mol + I(mol^0.5) + 
    tem:mol)
Model 2: log(AUCPD + 500) ~ Tem * Iso * Mol
  Res.Df     RSS Df Sum of Sq      F   Pr(>F)    
1    296 121.469                                 
2    256  90.257 40    31.212 2.2132 0.000118 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Código

# ATTENTION: Conforme antecipado, houve falta de ajuste. Isso quer dizer
# que o modelo considerado é uma alternativa muito grosseira para
# descrever o efeito do molhamento x temperatura. Portanto, o estudo dos
# fatores será feito considerando que temperatura e molhamento são
# fatores de níveis qualitativos.

#-----------------------------------------------------------------------
# Estudo com procedimento de comparação múltipla de médias.

anova(m1)

Analysis of Variance Table

Response: log(AUCPD + 500)
             Df  Sum Sq Mean Sq  F value    Pr(>F)    
Tem           3  81.595  27.198  77.1441 < 2.2e-16 ***
Iso           3 155.213  51.738 146.7461 < 2.2e-16 ***
Mol           3   5.260   1.753   4.9733  0.002266 ** 
Tem:Iso       9  19.333   2.148   6.0927 9.254e-08 ***
Tem:Mol       9   3.702   0.411   1.1667  0.316768    
Iso:Mol       9   2.244   0.249   0.7071  0.702243    
Tem:Iso:Mol  27  25.545   0.946   2.6836 3.145e-05 ***
Residuals   256  90.257   0.353                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Código

tb_means <- emmeans(m1, specs = ~Tem * Iso * Mol) |>
    as.data.frame() |>
    mutate(tem = as.numeric(as.character(Tem)),
           mol = as.numeric(as.character(Mol)))

ggplot(data = tb_means,
       mapping = aes(x = mol, y = emmean, group = 1)) +
    geom_line(color = "gray50", lty = 1, linewidth = 0.5) +
    geom_point() +
    geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                  width = 0) +
    facet_grid(facets = Iso ~ Tem) +
    labs(y = "Log of severity AUDPC",
         x = "Wetness period (h)")

Código

ggplot(data = tb_means,
       mapping = aes(x = tem, y = emmean, group = 1)) +
    geom_line(color = "gray50", lty = 1, linewidth = 0.5) +
    geom_point() +
    geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                  width = 0.1) +
    facet_grid(facets = Iso ~ Mol) +
    labs(y = "Log of severity AUDPC",
         x = "Temperature (°C)")

4 Análise do experimento 2

Código

#-----------------------------------------------------------------------
# Experimento 2 --------------------------------------------------------

tbi <- tb_auc |>
    filter(Exp == "2")

# Confere qual o intervalo de valores.
summary(tbi$AUCPD)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0     0.0   295.5  1695.2  3513.8  8213.0 

Código

# Soma uma constante para evitar valores negativos.
m0 <- lm((AUCPD + 50) ~ Tem * Iso * Mol, data = tbi)

# Exame dos pressupostos.
par(mfrow = c(2, 2))
plot(m0)

Código

layout(1)

# Transformação Box-Cox.
MASS::boxcox(m0)
abline(v = 0, col = "red")

Código

# Ajuste com os dados transformados.
m1 <- lm(log(AUCPD + 50) ~ Tem * Iso * Mol, data = tbi)

# Confere os pressupostos.
par(mfrow = c(2, 2))
plot(m1)

Código

layout(1)

#-----------------------------------------------------------------------
# Modelo que descreve uma superfície para temperatura e molhamento.

# Refina o modelo para descrever a superfície de resposta.
m2 <- lm(log(AUCPD + 50) ~
             Iso * (tem + I(tem^2) + mol + I(mol^0.5) + tem:mol),
         data = tbi)

# Teste da falta de ajuste.
anova(m2, m1)

Analysis of Variance Table

Model 1: log(AUCPD + 50) ~ Iso * (tem + I(tem^2) + mol + I(mol^0.5) + 
    tem:mol)
Model 2: log(AUCPD + 50) ~ Tem * Iso * Mol
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1    296 778.26                                  
2    256 567.32 40    210.93 2.3796 2.687e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Código

# ATTENTION: Conforme antecipado, houve falta de ajuste. Isso quer dizer
# que o modelo considerado é uma alternativa muito grosseira para
# descrever o efeito do molhamento x temperatura. Portanto, o estudo dos
# fatores será feito considerando que temperatura e molhamento são
# fatores de níveis qualitativos.

#-----------------------------------------------------------------------
# Estudo com procedimento de comparação múltipla de médias.

anova(m1)

Analysis of Variance Table

Response: log(AUCPD + 50)
             Df Sum Sq Mean Sq F value    Pr(>F)    
Tem           3 217.50  72.500 32.7149 < 2.2e-16 ***
Iso           3   5.16   1.719  0.7758 0.5084294    
Mol           3  51.45  17.151  7.7394 5.737e-05 ***
Tem:Iso       9  94.60  10.511  4.7431 7.562e-06 ***
Tem:Mol       9  75.84   8.427  3.8027 0.0001596 ***
Iso:Mol       9  18.95   2.105  0.9499 0.4823343    
Tem:Iso:Mol  27 104.63   3.875  1.7486 0.0147600 *  
Residuals   256 567.32   2.216                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Código

tb_means <- emmeans(m1, specs = ~Tem * Iso * Mol) |>
    as.data.frame() |>
    mutate(tem = as.numeric(as.character(Tem)),
           mol = as.numeric(as.character(Mol)))

ggplot(data = tb_means,
       mapping = aes(x = mol, y = emmean, group = 1)) +
    geom_line(color = "gray50", lty = 1, linewidth = 0.5) +
    geom_point() +
    geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                  width = 0) +
    facet_grid(facets = Iso ~ Tem) +
    labs(y = "Log of severity AUDPC",
         x = "Wetness period (h)")

Código

ggplot(data = tb_means,
       mapping = aes(x = tem, y = emmean, group = 1)) +
    geom_line(color = "gray50", lty = 1, linewidth = 0.5) +
    geom_point() +
    geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                  width = 0.1) +
    facet_grid(facets = Iso ~ Mol) +
    labs(y = "Log of severity AUDPC",
         x = "Temperature (°C)")

5 Ajuste de modelos de spline e polinomial

Código

#-----------------------------------------------------------------------
# Fazer ajuste usando splines por isolado x experimento para mostrar a
# forma da relação temperatura x molhamento.

library(rms)

# Definir 3 knots para cada variável. Usar o ponto médio entre níveis já
# que são 4 níveis por fator.
get_knots <- function(x) {
    u <- x |>
        unique() |>
        sort()
    head(u, n = -1) + diff(u)/2
}

# Knots.
tem_knots <- get_knots(tb_auc$tem)
mol_knots <- get_knots(tb_auc$mol)

seq_expand <- function(x, f = 0.05, length.out = 31) {
    er <- extendrange(x, f = f)
    seq(er[1], er[2], length.out = length.out)
}
seq_expand(0:10, f = 0.05, length.out = 31)

 [1] -0.5000000 -0.1333333  0.2333333  0.6000000  0.9666667  1.3333333
 [7]  1.7000000  2.0666667  2.4333333  2.8000000  3.1666667  3.5333333
[13]  3.9000000  4.2666667  4.6333333  5.0000000  5.3666667  5.7333333
[19]  6.1000000  6.4666667  6.8333333  7.2000000  7.5666667  7.9333333
[25]  8.3000000  8.6666667  9.0333333  9.4000000  9.7666667 10.1333333
[31] 10.5000000

Código

# Para obter os valores preditos.
pred <- with(tb_auc, {
    expand.grid(tem = seq_expand(tem, f = 0.05, length.out = 31),
                mol = seq_expand(mol, f = 0.05, length.out = 31),
                KEEP.OUT.ATTRS = FALSE)
})

# rcs(): rcs is a linear tail-restricted cubic spline function.
##' help(rcs, help_type = "html")

# FIXME TODO: incluir p-valor contra o modelo nulo. Incluir p-valor para a
# hipótese de LOF.

tb_auc$y <- log(tb_auc$AUCPD + 500)

# Ajustes para todas as combinações.
tb_fits <-
    tb_auc |>
    group_by(Exp, Iso) |>
    do(
        model_null = {
            # Modelo nulo.
            m0_null <- lm(y ~ 1, data = .)
            m0_null
        },
        model_sat = {
            # Modelo saturado.
            m0_sat <- lm(y ~ Tem * Mol,
                         data = .)
            m0_sat
        },
        model_2nd = {
            # Modelo de segunda ordem completo.
            m0_2nd <- lm(y ~
                             # tem + I(tem^2) + mol + I(mol^2) + tem:mol,
                             tem + I(tem^2) + mol + sqrt(mol) + tem:mol,
                             # poly(tem, 2) * poly(mol, 2),
                         data = .)
            m0_2nd
        },
        model_rcs = {
            # Modelo de splines.
            m0_rcs <- lm(y ~
                         rcs(tem, tem_knots) * rcs(mol, mol_knots),
                     data = .)
            m0_rcs
        }
    )

# Testes da falta de ajuste para o modelo de segunda ordem.
tb_fits |>
    mutate(
        signif_2nd = {
            anova(model_2nd, model_null)[["Pr(>F)"]][2]
        },
        lof_2nd_pvalue = {
            anova(model_2nd, model_sat)[["Pr(>F)"]][2]
        },
        signif_rcs = {
            anova(model_rcs, model_null)[["Pr(>F)"]][2]
        },
        lof_rcs_pvalue = {
            anova(model_rcs, model_sat)[["Pr(>F)"]][2]
        }
    ) |>
    select(Exp, Iso, starts_with("signif"), starts_with("lof"))

# A tibble: 8 × 6
# Rowwise: 
  Exp   Iso       signif_2nd signif_rcs lof_2nd_pvalue lof_rcs_pvalue
  <fct> <fct>          <dbl>      <dbl>          <dbl>          <dbl>
1 1     CrAaPR-01   7.36e- 6   9.85e- 6      0.00142        0.00234  
2 1     CrAaPR-27   1.69e- 4   5.01e- 4      0.0210         0.0165   
3 1     CrAaPR-40   2.99e-18   3.85e-19      0.0235         0.512    
4 1     CrAlPR-73   6.84e-11   9.43e-10      0.468          0.480    
5 2     CrAaPR-01   5.42e- 4   1.12e- 3      0.0278         0.0269   
6 2     CrAaPR-27   1.07e- 4   6.68e- 4      0.0000558      0.0000145
7 2     CrAaPR-40   8.89e- 6   3.60e- 5      0.135          0.132    
8 2     CrAlPR-73   4.64e- 5   1.82e- 4      0.496          0.534    

Código

# Obter o R2 para os modelos.
tb_fits |>
    mutate(
        R2_sat = {
            summary(model_sat)$r.squared
        },
        R2_2nd = {
            summary(model_2nd)$r.squared
        },
        R2_rcs = {
            summary(model_rcs)$r.squared
        }
    ) |>
    select(Exp, Iso, starts_with("R2"))

# A tibble: 8 × 5
# Rowwise: 
  Exp   Iso       R2_sat R2_2nd R2_rcs
  <fct> <fct>      <dbl>  <dbl>  <dbl>
1 1     CrAaPR-01  0.568  0.342  0.397
2 1     CrAaPR-27  0.470  0.278  0.314
3 1     CrAaPR-40  0.782  0.704  0.760
4 1     CrAlPR-73  0.589  0.526  0.547
5 2     CrAaPR-01  0.444  0.252  0.295
6 2     CrAaPR-27  0.588  0.288  0.307
7 2     CrAaPR-40  0.469  0.338  0.371
8 2     CrAlPR-73  0.395  0.305  0.337

Código

# Tabela com valores preditos no grid para fazer os gráficos.
tb_pred <-
    tb_fits |>
    mutate(pred = {
        list(cbind(pred, fit = predict(model_rcs, newdata = pred)))
        # list(cbind(pred, fit = predict(model_2nd, newdata = pred)))
    }) |>
    ungroup() |>
    select(Exp, Iso, pred) |>
    unnest(pred)
str(tb_pred)

tibble [7,688 × 5] (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 ...
 $ tem: num [1:7688] 14.1 14.8 15.4 16.1 16.7 ...
 $ mol: num [1:7688] 3.9 3.9 3.9 3.9 3.9 3.9 3.9 3.9 3.9 3.9 ...
 $ fit: num [1:7688] 6.24 6.26 6.29 6.31 6.33 ...

Código

# Escala de cores.
colr <- RColorBrewer::brewer.pal(11, "Spectral")
colr <- colorRampPalette(colr, space = "rgb")

tb_pred$fit |> summary()

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  5.696   6.813   7.285   7.381   7.898   9.059 

Código

# Superfície.
lattice::wireframe(fit ~ tem + mol | Iso + Exp,
                   data = tb_pred,
                   xlab = list("Temperature (°C)", rot = 30),
                   ylab = list("Wetness period (h)", rot = -40),
                   zlab = list("AUCPD", rot = 90),
                   # zlim = c(-3000, 7000),
                   # zoom = 1,
                   as.table = TRUE,
                   panel.3d.wireframe = panel.3d.contour,
                   col.contour = "black",
                   type = c("on", "bottom")[2],
                   drape = TRUE,
                   col.regions = colr(101),
                   scales = list(arrows = FALSE,
                                 z = list(distance = 1.25)))

Código

# Gráfico de níveis com linhas de contornos de nível.
levelplot(fit ~ tem + mol | Iso + Exp,
          data = tb_pred,
          ylab = "Período de molhamento (h)",
          xlab = "Temperatura (°C)",
          as.table = TRUE,
          contour = TRUE,
          # labels = TRUE,
          # at = seq(from = 0, to = 110, by = 5),
          col.regions = colr(101))

Código

# Fazer o gráfico de níveis com ggplot.
ggplot(data = tb_pred,
       mapping = aes(x = tem, y = mol, z = fit, group = 1)) +
    geom_tile(mapping = aes(fill = fit)) +
    geom_contour(color = "black", bins = 20) +
    facet_grid(facets = Exp ~ Iso) +
    # scale_fill_gradientn(colours = colr(101)) +
    scale_fill_distiller(palette = "Spectral") +
    labs(y = "Wetness period (h)",
         x = "Temperature (°C)",
         fill = "AUCPD") +
    theme_bw()

Código

# NOTE: Os resultados realmente mostram que nem todas ascurvas não podem
# ser simplificadas por modelos de segunda ordem ou splines.

6 GAM com bivariate thin-plate spline

Código

#-----------------------------------------------------------------------
# Usando GAM com bivariate thin-plate spline.
# https://stackoverflow.com/questions/52279218/rough-thin-plate-spline-fitting-thin-plate-spline-interpolation-in-r-with-mgcv

library(mgcv)

# Serão necessários para voltar para escala original ao final.
scaled_data <-
    tb_auc |>
    select(tem, mol) |>
    distinct() |>
    reframe(across(everything(), range))
scaled_data

# A tibble: 2 × 2
    tem   mol
  <dbl> <dbl>
1    15     6
2    33    48

Código

# Cria cópias com escalas padronizadas.
for (v in names(scaled_data)) {
    tb_auc[[paste0(v, "_sc")]] <-
        (tb_auc[[v]] - scaled_data[[v]][1])/diff(scaled_data[[v]])
}
str(tb_auc)

tibble [640 × 11] (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 ...
 $ tem   : num [1:640] 15 15 15 15 15 15 15 15 15 15 ...
 $ mol   : num [1:640] 6 6 6 6 6 12 12 12 12 12 ...
 $ rept  : num [1:640] 1 2 3 4 5 1 2 3 4 5 ...
 $ AUCPD : num [1:640] 0 0 0 0 0 0 0 0 0 0 ...
 $ Mol   : Factor w/ 4 levels "6","12","24",..: 1 1 1 1 1 2 2 2 2 2 ...
 $ Tem   : Factor w/ 4 levels "15","21","27",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ y     : num [1:640] 6.21 6.21 6.21 6.21 6.21 ...
 $ tem_sc: num [1:640] 0 0 0 0 0 0 0 0 0 0 ...
 $ mol_sc: num [1:640] 0 0 0 0 0 ...

Código

# Confere se estão na mesma escala.
tb_auc |>
    distinct(tem_sc, mol_sc) |>
    plot(tem_sc ~ mol_sc, data = _, asp = 1)

Código

# Para obter os valores preditos.
pred <- with(tb_auc, {
    expand.grid(tem_sc = seq_expand(c(0, 1), f = 0.05, length.out = 31),
                mol_sc = seq_expand(c(0, 1), f = 0.05, length.out = 31),
                KEEP.OUT.ATTRS = FALSE)
})
str(pred)

'data.frame':   961 obs. of  2 variables:
 $ tem_sc: num  -0.05 -0.0133 0.0233 0.06 0.0967 ...
 $ mol_sc: num  -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 ...

Código

# Volta para a escala original.
for (v in names(scaled_data)) {
    pred[[v]] <-
        pred[[paste0(v, "_sc")]] *
        diff(scaled_data[[v]]) +
        scaled_data[[v]][1]
}
str(pred)

'data.frame':   961 obs. of  4 variables:
 $ tem_sc: num  -0.05 -0.0133 0.0233 0.06 0.0967 ...
 $ mol_sc: num  -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 ...
 $ tem   : num  14.1 14.8 15.4 16.1 16.7 ...
 $ mol   : num  3.9 3.9 3.9 3.9 3.9 3.9 3.9 3.9 3.9 3.9 ...

Código

# Número de combinações únicas.
k <- nrow(distinct(tb_auc, tem_sc, mol_sc))
k

[1] 16

Código

# Define a variável resposta.
tb_auc$y <- log(tb_auc$AUCPD + 500)

# Ajustes para todas as combinações.
tb_fits <-
    tb_auc |>
    group_by(Exp, Iso) |>
    do(
        model_tps = {
            # Com parâmetro de penalidade.
            # gam(y ~ s(tem_sc, mol_sc, k = k, sp = 1), data = .)
            gam(y ~ s(tem_sc, mol_sc, k = k), data = .)
        }
    )

# Tabela com valores preditos no grid para fazer os gráficos.
tb_pred <-
    tb_fits |>
    mutate(pred = {
        # list(cbind(pred, fit = predict(model_rcs, newdata = pred)))
        list(cbind(pred, fit = predict(model_tps, newdata = pred)))
    }) |>
    ungroup() |>
    select(Exp, Iso, pred) |>
    unnest(pred)
str(tb_pred)

tibble [7,688 × 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 ...
 $ tem_sc: num [1:7688] -0.05 -0.0133 0.0233 0.06 0.0967 ...
 $ mol_sc: num [1:7688] -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 ...
 $ tem   : num [1:7688] 14.1 14.8 15.4 16.1 16.7 ...
 $ mol   : num [1:7688] 3.9 3.9 3.9 3.9 3.9 3.9 3.9 3.9 3.9 3.9 ...
 $ fit   : num [1:7688] 6.16 6.19 6.21 6.24 6.27 ...

Código

# Escala de cores.
colr <- RColorBrewer::brewer.pal(11, "Spectral")
colr <- colorRampPalette(colr, space = "rgb")

tb_pred$fit |> summary()

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  5.999   6.779   7.280   7.383   7.940   8.874 

Código

# Superfície.
lattice::wireframe(fit ~ tem + mol | Iso + Exp,
                   data = tb_pred,
                   xlab = list("Temperature (°C)", rot = 30),
                   ylab = list("Wetness period (h)", rot = -40),
                   zlab = list("AUCPD", rot = 90),
                   # zlim = c(-3000, 7000),
                   # zoom = 1,
                   as.table = TRUE,
                   panel.3d.wireframe = panel.3d.contour,
                   col.contour = "black",
                   type = c("on", "bottom")[2],
                   drape = TRUE,
                   col.regions = colr(101),
                   scales = list(arrows = FALSE,
                                 z = list(distance = 1.25)))

Código

# Gráfico de níveis com linhas de contornos de nível.
levelplot(fit ~ tem + mol | Iso + Exp,
          data = tb_pred,
          ylab = "Período de molhamento (h)",
          xlab = "Temperatura (°C)",
          as.table = TRUE,
          contour = TRUE,
          # labels = TRUE,
          # at = seq(from = 0, to = 110, by = 5),
          col.regions = colr(101))

Código

# Fazer o gráfico de níveis com ggplot.
ggplot(data = tb_pred,
       mapping = aes(x = tem, y = mol, z = fit, group = 1)) +
    geom_tile(mapping = aes(fill = fit)) +
    geom_contour(color = "black", bins = 20) +
    facet_grid(facets = Exp ~ Iso) +
    # scale_fill_gradientn(colours = colr(101)) +
    scale_fill_distiller(palette = "Spectral") +
    labs(y = "Wetness period (h)",
         x = "Temperature (°C)",
         fill = "AUCPD") +
    theme_bw()