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

library(readxl)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ purrr::%||%()   masks base::%||%()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
# browseURL("https://cran.r-project.org/web/packages/epiphy/epiphy.pdf")
# browseURL("https://chgigot.github.io/epiphy/")
#
# browseURL("https://cran.r-project.org/web/packages/epifitter/epifitter.pdf")
# browseURL("https://github.com/AlvesKS/epifitter")
#
# browseURL("https://rdrr.io/cran/nlraa/")
# browseURL("https://www.openanalytics.eu/blog/2022/05/17/new-nonlinear-ls-solvers-in-r/")
# browseURL("https://www.statforbiology.com/2020/stat_nls_usefulfunctions/")

# To create rounding specifications for data display.
# browseURL("https://bcjaeger.github.io/table.glue/index.html")

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

tb <- read_excel("exvivo.xlsx", sheet = 1)
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 6 6 6 6 6 ...
##  $ sev : num [1:4480] 0 0 0 0 0 0 0 0 0 0 ...
##  $ time: num [1:4480] 12 12 12 12 12 24 24 24 24 24 ...
DataExplorer::plot_scatterplot(tb, by = "sev", ncol = 2, nrow = 5)

# Resumo numérico das variáveis.
# gtsummary::tbl_summary(tb)
gtsummary::tbl_summary(tb) |>
    gtsummary::as_kable()
Characteristic N = 4,480
exp
1 2,240 (50%)
2 2,240 (50%)
iso
CrAaPR-01 1,120 (25%)
CrAaPR-27 1,120 (25%)
CrAaPR-40 1,120 (25%)
CrAlPR-73 1,120 (25%)
rept
1 896 (20%)
2 896 (20%)
3 896 (20%)
4 896 (20%)
5 896 (20%)
tem
15 1,120 (25%)
21 1,120 (25%)
27 1,120 (25%)
33 1,120 (25%)
mol
6 1,120 (25%)
12 1,120 (25%)
24 1,120 (25%)
48 1,120 (25%)
sev 0.00 (0.00, 0.05)
time
12 640 (14%)
24 640 (14%)
48 640 (14%)
62 640 (14%)
84 640 (14%)
100 640 (14%)
124 640 (14%) 
# Número de repetições por condição experimental/avaliação
tb |>
    count(exp, iso, tem, mol, time)
## # A tibble: 896 × 6
##      exp iso         tem   mol  time     n
##    <dbl> <chr>     <dbl> <dbl> <dbl> <int>
##  1     1 CrAaPR-01    15     6    12     5
##  2     1 CrAaPR-01    15     6    24     5
##  3     1 CrAaPR-01    15     6    48     5
##  4     1 CrAaPR-01    15     6    62     5
##  5     1 CrAaPR-01    15     6    84     5
##  6     1 CrAaPR-01    15     6   100     5
##  7     1 CrAaPR-01    15     6   124     5
##  8     1 CrAaPR-01    15    12    12     5
##  9     1 CrAaPR-01    15    12    24     5
## 10     1 CrAaPR-01    15    12    48     5
## # ℹ 886 more rows
# Obervações no tempo.
tb |>
    count(time)
## # A tibble: 7 × 2
##    time     n
##   <dbl> <int>
## 1    12   640
## 2    24   640
## 3    48   640
## 4    62   640
## 5    84   640
## 6   100   640
## 7   124   640
# Converte variáveis qualitativas para fator.
tb <-
    tb |>
    mutate_at(c("exp", "iso"), factor)

ggplot(tb, aes(x = time, y = sev, color = iso,
               shape = exp, linetype = exp)) +
    facet_grid(tem ~ mol) +
    stat_summary(geom = "line", fun = "mean") +
    geom_point()

# Confere a média em cada combinação de fatores.
ggplot(tb, aes(x = time, y = sev)) +
    facet_grid(exp + tem ~ iso + mol) +
    stat_summary(geom = "line", fun = "mean") +
    geom_point()

# Destaca cada unidade experimental.
ggplot(tb, aes(x = time, y = sev)) +
    facet_grid(exp + tem ~ iso + mol) +
    geom_line(aes(group = rept)) +
    geom_point()

# ggsave("sev_time.png", width = 16, height = 9)

# NOTE: Muitas curvas com severidade igual a zero em todo o período.
# Isso impede que de ajustar qualquer modelo de regressão pois não há
# variação nos dados.

# QUESTION: Para que estimar uma curva para a severidade em função do
# tempo?

# Identifica a unidade experimental onde a avaliação da severidade é
# realizada ao longo do tempo.
tb <- tb |>
    mutate(ue = interaction(exp, iso, tem, mol, rept, drop = TRUE))

# Mostra uma inidade experimental aleatória.
tb |>
    filter(ue %in% sample(ue, size = 1)) |>
    ggplot(aes(x = time, y = sev)) +
    facet_grid(tem ~ mol) +
    stat_summary(geom = "line", fun = "mean") +
    geom_point()

#-----------------------------------------------------------------------
# Selecionar uma condição experimental.

# Mostra uma inidade experimental aleatória.
tb |>
    filter(tem == 27, mol == 24) |>
    ggplot(aes(x = time, y = sev)) +
    facet_grid(exp ~ iso) +
    geom_line(aes(group = rept), alpha = 0.15) +
    geom_point() +
    stat_summary(geom = "line", fun = "mean", color = "red") +
    labs(title = "Temperatura 27°C e Molhamento 24h")

# ggsave("sev_time.png", width = 16, height = 9)

#-----------------------------------------------------------------------
# Curva.

# Monomolecular
#
#                    y = K * (1 - B * exp(-r * t))
#
# K = maximum level of disease or asymptote of disease progress curve, y
# = disease at timeof observation, B = a parameter related to the level
# of initial disease or point of inflexion, r = rate of disease increase
# and t = time.
#
# Essa é uma reparametrização na qual B e substituido por y0
#
#                    y = 1 - (1 - y0) * exp(-r * t)
#
# em que y0 é o intercepto da curva e r é a taxa. A assintora superior é
# 1.

# Com y0 < 0 não tem sentido biológico.
with(list(y0 = -0.1, r = 0.1), {
    curve(1 - (1 - y0) * exp(-r * x),
          from = 0, to = 50,
          ylim = c(min(y0, 0), 1),
          col = "red",
          lwd = 2)
    abline(h = 0, col = "gray")
})

# Mas com x0 > 0, a curva pode ser deslocada para a direita.
with(list(r = 0.1, x0 = 10), {
    # y0 é função de x0.
    y0 <- (1 - exp(-r * (0 - x0)))
    curve((1 - exp(-r * (x - x0))),
          from = 0, to = 50, n = 501,
          ylim = c(min(0, y0), 1),
          col = "red",
          lty = 2,
          lwd = 1)
    abline(h = 0, v = x0, col = "gray")
    abline(h = y0, col = "gray", lty = 2)
    # Dessa forna, f(x) é não negativa, o que é compatível com o
    # fenômeno.
    curve((1 - exp(-r * (x - x0))) * (x > x0),
          add = TRUE,
          col = "red",
          lwd = 2)
})

#-----------------------------------------------------------------------
# Prototipação do ajuste de modelo aos dados.

levels(tb$iso)
## [1] "CrAaPR-01" "CrAaPR-27" "CrAaPR-40" "CrAlPR-73"
tbi <- tb |>
    filter(tem == 27, mol == 24, exp == "2", iso == levels(iso)[4])

plot(tbi$time, tbi$sev)
with(list(r = 0.01, x0 = 50), {
    curve((1 - exp(-r * (x - x0))) * (x > x0),
          from = 0, to = max(tbi$time),
          ylim = c(0, 1),
          col = "red",
          add = TRUE,
          lwd = 2)
    abline(h = 0, v = x0, col = "gray")
})

n0 <- nls(sev ~ (1 - exp(-r * (time - x0))) * (time > x0),
          data = tbi,
          start = list(r = 0.01, x0 = 50))
summary(n0)
## 
## Formula: sev ~ (1 - exp(-r * (time - x0))) * (time > x0)
## 
## Parameters:
##     Estimate Std. Error t value Pr(>|t|)    
## r   0.013557   0.002897   4.680 4.72e-05 ***
## x0 45.629751   6.495558   7.025 4.90e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2142 on 33 degrees of freedom
## 
## Number of iterations to convergence: 4 
## Achieved convergence tolerance: 3.063e-08
plot(tbi$time, tbi$sev)
with(as.list(coef(n0)), {
    curve((1 - exp(-r * (x - x0))) * (x > x0),
          from = 0, to = max(tbi$time),
          ylim = c(0, 1),
          col = "red",
          add = TRUE,
          lwd = 2)
    abline(h = 0, v = x0, col = "gray")
})

broom::tidy(n0)
## # A tibble: 2 × 5
##   term  estimate std.error statistic      p.value
##   <chr>    <dbl>     <dbl>     <dbl>        <dbl>
## 1 r       0.0136   0.00290      4.68 0.0000472   
## 2 x0     45.6      6.50         7.02 0.0000000490
broom::glance(n0)
## # A tibble: 1 × 9
##   sigma isConv       finTol logLik   AIC    BIC deviance df.residual  nobs
##   <dbl> <lgl>         <dbl>  <dbl> <dbl>  <dbl>    <dbl>       <int> <int>
## 1 0.214 TRUE   0.0000000306   5.29 -4.58 0.0853     1.51          33    35
confint(n0, level = 0.95)
## Waiting for profiling to be done...
##            2.5%       97.5%
## r   0.007955212  0.02187156
## x0 23.725754619 59.70567769
confint.default(n0, level = 0.95)
##           2.5 %      97.5 %
## r   0.007879138  0.01923434
## x0 32.898690695 58.36081140
# Calculando R² como quadrado da correlação entre observado e predito.
cor(fitted(n0), fitted(n0) + residuals(n0))^2
## [1] 0.5830175
#-----------------------------------------------------------------------
# Apenas para aprendizado, ajustar usando o pacote `minpack.lm`.

n0 <-
    minpack.lm::nlsLM(
        formula = sev ~ (1 - exp(-r * (time - x0))) * (time > x0),
        data = tbi,
        start = list(r = 0.01, x0 = 50),
        trace = TRUE
    )
## It.    0, RSS =    1.79561, Par. =       0.01         50
## It.    1, RSS =    1.53759, Par. =  0.0121644     40.366
## It.    2, RSS =    1.51563, Par. =  0.0135581    46.3222
## It.    3, RSS =    1.51462, Par. =  0.0135568    45.6331
## It.    4, RSS =    1.51462, Par. =  0.0135567    45.6298
## It.    5, RSS =    1.51462, Par. =  0.0135567    45.6298
summary(n0)
## 
## Formula: sev ~ (1 - exp(-r * (time - x0))) * (time > x0)
## 
## Parameters:
##     Estimate Std. Error t value Pr(>|t|)    
## r   0.013557   0.002897   4.680 4.72e-05 ***
## x0 45.629750   6.495559   7.025 4.90e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2142 on 33 degrees of freedom
## 
## Number of iterations to convergence: 5 
## Achieved convergence tolerance: 1.49e-08
broom::tidy(n0)
## # A tibble: 2 × 5
##   term  estimate std.error statistic      p.value
##   <chr>    <dbl>     <dbl>     <dbl>        <dbl>
## 1 r       0.0136   0.00290      4.68 0.0000472   
## 2 x0     45.6      6.50         7.02 0.0000000490
broom::glance(n0)
## # A tibble: 1 × 9
##   sigma isConv       finTol logLik   AIC    BIC deviance df.residual  nobs
##   <dbl> <lgl>         <dbl>  <dbl> <dbl>  <dbl>    <dbl>       <int> <int>
## 1 0.214 TRUE   0.0000000149   5.29 -4.58 0.0853     1.51          33    35
confint(n0, level = 0.95)
## Waiting for profiling to be done...
##            2.5%       97.5%
## r   0.007955212  0.02187156
## x0 23.725754625 59.70567768
confint.default(n0, level = 0.95)
##           2.5 %      97.5 %
## r   0.007879138  0.01923434
## x0 32.898689173 58.36081106
#-----------------------------------------------------------------------
# Apenas para aprendizado, ajustar usando o pacote `gslnls`.

# GSL (Levenberg-Marquardt algorithm)
# library(gslnls)  ## v1.1.1
n0 <-
    gslnls::gsl_nls(
        fn = sev ~ (1 - exp(-r * (time - x0))) * (time > x0),
        data = tbi,
        start = list(r = 0.01, x0 = 50),
        algorithm = "lm",
        control = list(scale = "levenberg")
    )

summary(n0)
## 
## Formula: sev ~ (1 - exp(-r * (time - x0))) * (time > x0)
## 
## Parameters:
##     Estimate Std. Error t value Pr(>|t|)    
## r   0.013557   0.002897   4.680 4.72e-05 ***
## x0 45.629750   6.495559   7.025 4.90e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2142 on 33 degrees of freedom
## 
## Number of iterations to convergence: 16 
## Achieved convergence tolerance: -2.22e-16
broom::tidy(n0)
## # A tibble: 2 × 5
##   term  estimate std.error statistic      p.value
##   <chr>    <dbl>     <dbl>     <dbl>        <dbl>
## 1 r       0.0136   0.00290      4.68 0.0000472   
## 2 x0     45.6      6.50         7.02 0.0000000490
broom::glance(n0)
## # A tibble: 1 × 9
##   sigma isConv    finTol logLik   AIC    BIC deviance df.residual  nobs
##   <dbl> <lgl>      <dbl>  <dbl> <dbl>  <dbl>    <dbl>       <int> <int>
## 1 0.214 TRUE   -2.22e-16   5.29 -4.58 0.0853     1.51          33    35
confint(n0, level = 0.95)
##           2.5 %     97.5 %
## r   0.007663179  0.0194503
## x0 32.414436995 58.8450635
confint(n0, method = "profile", level = 0.95)
## Waiting for profiling to be done...
##            2.5%       97.5%
## r   0.007955212  0.02187156
## x0 23.725754623 59.70567769
confint(n0, method = "asymptotic", level = 0.95)
##           2.5 %     97.5 %
## r   0.007663179  0.0194503
## x0 32.414436995 58.8450635
confint.default(n0, level = 0.95)
##           2.5 %      97.5 %
## r   0.007879138  0.01923434
## x0 32.898689413 58.36081104
# Intervalos de confiança e predição.
predict(n0, interval = "confidence", level = 0.95)
##              fit         lwr       upr
##  [1,] 0.00000000  0.00000000 0.0000000
##  [2,] 0.00000000  0.00000000 0.0000000
##  [3,] 0.00000000  0.00000000 0.0000000
##  [4,] 0.00000000  0.00000000 0.0000000
##  [5,] 0.00000000  0.00000000 0.0000000
##  [6,] 0.00000000  0.00000000 0.0000000
##  [7,] 0.00000000  0.00000000 0.0000000
##  [8,] 0.00000000  0.00000000 0.0000000
##  [9,] 0.00000000  0.00000000 0.0000000
## [10,] 0.00000000  0.00000000 0.0000000
## [11,] 0.03162208 -0.13277606 0.1960202
## [12,] 0.03162208 -0.13277606 0.1960202
## [13,] 0.03162208 -0.13277606 0.1960202
## [14,] 0.03162208 -0.13277606 0.1960202
## [15,] 0.03162208 -0.13277606 0.1960202
## [16,] 0.19902632  0.09328288 0.3047698
## [17,] 0.19902632  0.09328288 0.3047698
## [18,] 0.19902632  0.09328288 0.3047698
## [19,] 0.19902632  0.09328288 0.3047698
## [20,] 0.19902632  0.09328288 0.3047698
## [21,] 0.40558374  0.30784437 0.5033231
## [22,] 0.40558374  0.30784437 0.5033231
## [23,] 0.40558374  0.30784437 0.5033231
## [24,] 0.40558374  0.30784437 0.5033231
## [25,] 0.40558374  0.30784437 0.5033231
## [26,] 0.52149242  0.40912674 0.6338581
## [27,] 0.52149242  0.40912674 0.6338581
## [28,] 0.52149242  0.40912674 0.6338581
## [29,] 0.52149242  0.40912674 0.6338581
## [30,] 0.52149242  0.40912674 0.6338581
## [31,] 0.65439022  0.52946935 0.7793111
## [32,] 0.65439022  0.52946935 0.7793111
## [33,] 0.65439022  0.52946935 0.7793111
## [34,] 0.65439022  0.52946935 0.7793111
## [35,] 0.65439022  0.52946935 0.7793111
predict(n0, interval = "prediction", level = 0.95)
##              fit         lwr       upr
##  [1,] 0.00000000 -0.43586889 0.4358689
##  [2,] 0.00000000 -0.43586889 0.4358689
##  [3,] 0.00000000 -0.43586889 0.4358689
##  [4,] 0.00000000 -0.43586889 0.4358689
##  [5,] 0.00000000 -0.43586889 0.4358689
##  [6,] 0.00000000 -0.43586889 0.4358689
##  [7,] 0.00000000 -0.43586889 0.4358689
##  [8,] 0.00000000 -0.43586889 0.4358689
##  [9,] 0.00000000 -0.43586889 0.4358689
## [10,] 0.00000000 -0.43586889 0.4358689
## [11,] 0.03162208 -0.43421956 0.4974637
## [12,] 0.03162208 -0.43421956 0.4974637
## [13,] 0.03162208 -0.43421956 0.4974637
## [14,] 0.03162208 -0.43421956 0.4974637
## [15,] 0.03162208 -0.43421956 0.4974637
## [16,] 0.19902632 -0.24948607 0.6475387
## [17,] 0.19902632 -0.24948607 0.6475387
## [18,] 0.19902632 -0.24948607 0.6475387
## [19,] 0.19902632 -0.24948607 0.6475387
## [20,] 0.19902632 -0.24948607 0.6475387
## [21,] 0.40558374 -0.04110930 0.8522768
## [22,] 0.40558374 -0.04110930 0.8522768
## [23,] 0.40558374 -0.04110930 0.8522768
## [24,] 0.40558374 -0.04110930 0.8522768
## [25,] 0.40558374 -0.04110930 0.8522768
## [26,] 0.52149242  0.07137273 0.9716121
## [27,] 0.52149242  0.07137273 0.9716121
## [28,] 0.52149242  0.07137273 0.9716121
## [29,] 0.52149242  0.07137273 0.9716121
## [30,] 0.52149242  0.07137273 0.9716121
## [31,] 0.65439022  0.20097329 1.1078072
## [32,] 0.65439022  0.20097329 1.1078072
## [33,] 0.65439022  0.20097329 1.1078072
## [34,] 0.65439022  0.20097329 1.1078072
## [35,] 0.65439022  0.20097329 1.1078072
# IC para funções de parâmetros usando método delta.
gslnls::confintd(n0, expr = c("log(r)", "sqrt(x0)"), level = 0.95)
##                fit       lwr       upr
## log(r)   -4.300872 -4.735604 -3.866139
## sqrt(x0)  6.754980  5.776789  7.733170
#-----------------------------------------------------------------------
# Função self-start.

# xx <- tbi$time
# yy <- tbi$sev
# plot(log(1/(1 - yy * 0.97)) ~ xx)
# coef(lm(log(1/(1 - yy * 0.97)) ~ 0 + xx))

model.init <- function(mCall, data, LHS, ...){
    xy <- data.frame(sortedXyData(mCall[["time"]], LHS, data))
    xx <- xy[["x"]]
    yy <- xy[["y"]]
    my <- tapply(yy, xx, mean, na.rm = TRUE)
    mx <- as.numeric(names(my))
    # Primeiro x com y > 0.
    x0 <- min(xx[yy > 0])
    if (max(yy) == 1) {
        k <- 0.97
        log_y <- log(1/(1 - yy * k))
    } else {
        log_y <- log(1/(1 - yy))
    }
    # Taxa obtida com modelo linearizado.
    xx_pos <- xx[yy > 0]
    log_y_pos <- log_y[yy > 0]
    r <- coef(lm(log_y_pos ~ 0 + xx_pos))
    value <- c(x0, r)
    names(value) <- mCall[c("x0", "r")]
    return(value)
}

# model <- expression((1 - (exp(-r * (time - x0)))))
# pars <- c("x0", "r")
# sapply(pars, D, expr = model)
modelfun <- function(time, x0, r){
    value <- (1 - (exp(-r * (time - x0)))) * (time > x0)
    attr(value, "gradient") <-
        cbind(
            x0 = (-(exp(-r * (time - x0)) * r)) * (time > x0),
            r = (exp(-r * (time - x0)) * (time - x0)) * (time > x0)
        )
    return(value)
}
# attr(modelfun, "initial") <- model.init
str(modelfun)
## function (time, x0, r)  
##  - attr(*, "srcref")= 'srcref' int [1:8] 273 13 281 1 13 1 273 281
##   ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' <environment: 0x55cba6bd5008>
# Testa o gradiente.
with(unique(tbi[, c("time")]), {
    cbind(time = time,
          attr(modelfun(time, x0 = 50, r = 0.01),
               which = "gradient"))
})
##      time           x0        r
## [1,]   12  0.000000000  0.00000
## [2,]   24  0.000000000  0.00000
## [3,]   48  0.000000000  0.00000
## [4,]   62 -0.008869204 10.64305
## [5,]   84 -0.007117703 24.20019
## [6,]  100 -0.006065307 30.32653
## [7,]  124 -0.004771139 35.30643
# Saí o gradiente atributo gradiente junto.
SSmodel <- selfStart(
    model = modelfun,
    parameters = c("x0", "r"),
    initial = model.init,
    template = function(time, x0, r) NULL
)
# SSmodel

# As pré-estimativas avaliadas.
getInitial(sev ~ SSmodel(time, x0, r),
           data = tbi)
##          x0           r 
## 24.00000000  0.00674073
# NOTE: Depois de testes vi que essa self-start não é tão eficiente.
# Preciso dedicar mais tempo para fazer tratamentos de exceções e deixar
# ela mais "esperta".

#-----------------------------------------------------------------------
# Ajustes por unidade experimental.

tbi <-
    tb |>
    filter(tem == 27, mol == 24) |>
    droplevels() |>
    mutate(cond = interaction(exp, iso, tem, mol, drop = TRUE))
str(tbi)
## tibble [280 × 9] (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 ...
##  $ rept: num [1:280] 1 2 3 4 5 1 2 3 4 5 ...
##  $ tem : num [1:280] 27 27 27 27 27 27 27 27 27 27 ...
##  $ mol : num [1:280] 24 24 24 24 24 24 24 24 24 24 ...
##  $ sev : num [1:280] 0 0 0 0 0 0 0 0 0 0 ...
##  $ time: num [1:280] 12 12 12 12 12 24 24 24 24 24 ...
##  $ ue  : Factor w/ 40 levels "1.CrAaPR-01.27.24.1",..: 1 9 17 25 33 1 9 17 25 33 ...
##  $ cond: Factor w/ 8 levels "1.CrAaPR-01.27.24",..: 1 1 1 1 1 1 1 1 1 1 ...
n0_ue <-
    nlme::nlsList(sev ~ (1 - exp(-r * (time - x0))) * (time > x0) | ue,
                  start = list(r = 0.02, x0 = 55),
                  data = tbi)
## Warning: 8 errors caught in nls(model, data = data, control = controlvals, start = start).  The error messages and their frequencies are
## 
##                                                singular gradient 
##                                                                4 
## step factor 0.000488281 reduced below 'minFactor' of 0.000976562 
##                                                                4
# n0_ue <-
#     nlme::nlsList(sev ~ SSmodel(time, x0, r) | ue,
#                   # start = list(x0 = 55, r = 0.02),
#                   data = tbi)

summary(n0_ue)
## Call:
##   Model: sev ~ (1 - exp(-r * (time - x0))) * (time > x0) | ue 
##    Data: tbi 
## 
## Coefficients:
##    r 
##                        Estimate  Std. Error   t value     Pr(>|t|)
## 1.CrAaPR-01.27.24.1 0.152936035 0.140265783 1.0903303 1.363726e-03
## 2.CrAaPR-01.27.24.1 0.152936035 0.140265783 1.0903303 1.363726e-03
## 1.CrAaPR-27.27.24.1          NA          NA        NA           NA
## 2.CrAaPR-27.27.24.1 0.122442825 0.028324826 4.3228094 2.940508e-02
## 1.CrAaPR-40.27.24.1 0.152936035 0.140265783 1.0903303 1.363726e-03
## 2.CrAaPR-40.27.24.1          NA          NA        NA           NA
## 1.CrAlPR-73.27.24.1 0.152936038 0.140265797 1.0903302 1.402076e-03
## 2.CrAlPR-73.27.24.1 0.020899975 0.002897012 7.2143223 5.304805e-03
## 1.CrAaPR-01.27.24.2 0.104029731 0.030542441 3.4060713 8.635194e-07
## 2.CrAaPR-01.27.24.2 0.015828678 0.005442375 2.9084140 4.421903e-05
## 1.CrAaPR-27.27.24.2          NA          NA        NA           NA
## 2.CrAaPR-27.27.24.2 0.026223977 0.004724773 5.5503146 8.852566e-03
## 1.CrAaPR-40.27.24.2 0.114340041 0.025135780 4.5488957 1.274414e-06
## 2.CrAaPR-40.27.24.2          NA          NA        NA           NA
## 1.CrAlPR-73.27.24.2 0.241806322 0.220746004 1.0954052 3.246980e-04
## 2.CrAlPR-73.27.24.2 0.019473337 0.003382701 5.7567411 3.746884e-02
## 1.CrAaPR-01.27.24.3 0.238786699 0.220626488 1.0823120 7.898440e-03
## 2.CrAaPR-01.27.24.3 0.031618992 0.004269411 7.4059377 7.049363e-02
## 1.CrAaPR-27.27.24.3          NA          NA        NA           NA
## 2.CrAaPR-27.27.24.3 0.017624294 0.003079018 5.7239990 1.233996e-02
## 1.CrAaPR-40.27.24.3          NA          NA        NA           NA
## 2.CrAaPR-40.27.24.3 0.051085642 0.008488375 6.0183062 5.384525e-03
## 1.CrAlPR-73.27.24.3 0.160896893 0.050238646 3.2026518 2.044013e-05
## 2.CrAlPR-73.27.24.3 0.008230541 0.001813034 4.5396518 6.912598e-03
## 1.CrAaPR-01.27.24.4 0.126341104 0.130519052 0.9679897 2.069618e-03
## 2.CrAaPR-01.27.24.4 0.114340041 0.025135780 4.5488957 1.274414e-06
## 1.CrAaPR-27.27.24.4 0.233776634 0.220408148 1.0606533 3.974868e-04
## 2.CrAaPR-27.27.24.4 0.032956049 0.004913402 6.7073791 7.365514e-03
## 1.CrAaPR-40.27.24.4 0.114340041 0.025135780 4.5488957 1.274414e-06
## 2.CrAaPR-40.27.24.4 0.089065001 0.022179226 4.0156948 3.449326e-03
## 1.CrAlPR-73.27.24.4 0.267300948 0.252349232 1.0592501 4.593439e-03
## 2.CrAlPR-73.27.24.4 0.001460673 0.001656226 0.8819283 2.274943e-03
## 1.CrAaPR-01.27.24.5 0.243867714 0.220822855 1.1043590 3.126276e-04
## 2.CrAaPR-01.27.24.5 0.114340041 0.025135780 4.5488957 1.274414e-06
## 1.CrAaPR-27.27.24.5          NA          NA        NA           NA
## 2.CrAaPR-27.27.24.5 0.091004538 0.022419352 4.0591957 3.367861e-03
## 1.CrAaPR-40.27.24.5          NA          NA        NA           NA
## 2.CrAaPR-40.27.24.5 0.188991312 0.057565469 3.2830674 6.626691e-02
## 1.CrAlPR-73.27.24.5 0.160896893 0.050238646 3.2026518 2.044013e-05
## 2.CrAlPR-73.27.24.5 0.060706988 0.015994271 3.7955458 9.499104e-05
##    x0 
##                     Estimate Std. Error    t value     Pr(>|t|)
## 1.CrAaPR-01.27.24.1 59.28266  2.6133904  22.684196 4.471753e-10
## 2.CrAaPR-01.27.24.1 59.28266  2.6133904  22.684196 4.471753e-10
## 1.CrAaPR-27.27.24.1       NA         NA         NA           NA
## 2.CrAaPR-27.27.24.1 47.97243  0.5841565  82.122561 3.041825e-08
## 1.CrAaPR-40.27.24.1 59.28266  2.6133904  22.684196 4.471753e-10
## 2.CrAaPR-40.27.24.1       NA         NA         NA           NA
## 1.CrAlPR-73.27.24.1 59.28266  2.6133905  22.684196 4.610664e-10
## 2.CrAlPR-73.27.24.1 47.13088  3.1361109  15.028447 1.878725e-04
## 1.CrAaPR-01.27.24.2 61.97049  0.6869774  90.207470 6.709870e-14
## 2.CrAaPR-01.27.24.2 97.74921  5.1352291  19.035024 3.906232e-09
## 1.CrAaPR-27.27.24.2       NA         NA         NA           NA
## 2.CrAaPR-27.27.24.2 53.36052  4.3287621  12.326969 2.504811e-04
## 1.CrAaPR-40.27.24.2 47.99796  0.6223463  77.124205 9.229316e-13
## 2.CrAaPR-40.27.24.2       NA         NA         NA           NA
## 1.CrAlPR-73.27.24.2 83.85265  0.3358530 249.670697 6.047381e-16
## 2.CrAlPR-73.27.24.2 55.95992  4.6082262  12.143483 1.943645e-03
## 1.CrAaPR-01.27.24.3 83.65080  0.4621590 181.000043 1.015799e-13
## 2.CrAaPR-01.27.24.3 21.99398  2.4778576   8.876208 4.046655e-02
## 1.CrAaPR-27.27.24.3       NA         NA         NA           NA
## 2.CrAaPR-27.27.24.3 57.11650  4.7235967  12.091739 4.717468e-04
## 1.CrAaPR-40.27.24.3       NA         NA         NA           NA
## 2.CrAaPR-40.27.24.3 60.91400  1.5329330  39.736895 6.590929e-07
## 1.CrAlPR-73.27.24.3 47.77796  0.4706508 101.514663 6.677260e-13
## 2.CrAlPR-73.27.24.3 37.15748  9.4743483   3.921904 1.242669e-02
## 1.CrAaPR-01.27.24.4 54.54476  7.9239054   6.883571 1.516146e-07
## 2.CrAaPR-01.27.24.4 47.99796  0.6223463  77.124205 9.229316e-13
## 1.CrAaPR-27.27.24.4 83.30480  0.7543102 110.438379 3.764804e-14
## 2.CrAaPR-27.27.24.4 60.52787  2.3327732  25.946744 1.356317e-05
## 1.CrAaPR-40.27.24.4 47.99796  0.6223463  77.124205 9.229316e-13
## 2.CrAaPR-40.27.24.4 61.57250  0.8513217  72.325765 2.615327e-09
## 1.CrAlPR-73.27.24.4 47.39199  0.6603615  71.766743 4.853509e-12
## 2.CrAlPR-73.27.24.4 57.90667 45.4259548   1.274749 4.205155e-04
## 1.CrAaPR-01.27.24.5 83.98767  0.2930996 286.549905 3.038926e-16
## 2.CrAaPR-01.27.24.5 47.99796  0.6223463  77.124205 9.229316e-13
## 1.CrAaPR-27.27.24.5       NA         NA         NA           NA
## 2.CrAaPR-27.27.24.5 61.97588  0.7842564  79.025025 1.724507e-09
## 1.CrAaPR-40.27.24.5       NA         NA         NA           NA
## 2.CrAaPR-40.27.24.5 11.99870  0.3766389  31.857311 3.071365e-06
## 1.CrAlPR-73.27.24.5 47.77796  0.4706508 101.514663 6.677260e-13
## 2.CrAlPR-73.27.24.5 54.98117  2.8930070  19.004851 3.266359e-08
## 
## Residual standard error: 0.07115649 on 160 degrees of freedom
colMeans(coef(n0_ue), na.rm = TRUE)
##          r         x0 
##  0.1132631 56.6156026
# Ajustes que falharam.
coef(n0_ue)[[1]] |>
    is.na() |>
    sum()
## [1] 8
# Condições experimentais para a predição.
x_seq <- c(seq(0, max(tbi$time), length.out = 51),
           coef(n0_ue)[["x0"]]) |>
    sort()
tb_pred <- tbi |>
    distinct(exp, iso, tem, mol, cond, ue)
tb_pred <- cross_join(tb_pred,
                      data.frame(time = x_seq))
str(tb_pred)
## tibble [3,320 × 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 : num [1:3320] 27 27 27 27 27 27 27 27 27 27 ...
##  $ mol : num [1:3320] 24 24 24 24 24 24 24 24 24 24 ...
##  $ cond: Factor w/ 8 levels "1.CrAaPR-01.27.24",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ ue  : Factor w/ 40 levels "1.CrAaPR-01.27.24.1",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ time: num [1:3320] 0 2.48 4.96 7.44 9.92 ...
# Predição.
tb_pred$fit <- predict(n0_ue, newdata = tb_pred)

tbi |>
    ggplot(aes(x = time, y = sev)) +
    facet_wrap(~ue) +
    # facet_grid(exp ~ iso) +
    geom_line(aes(group = rept), alpha = 0.10) +
    stat_summary(geom = "line", fun = "mean",
                 color = "darkgreen",
                 alpha = 0.40) +
    geom_point() +
    labs(title = "Temperatura 27°C e Molhamento 24h") +
    geom_line(data = tb_pred, aes(y = fit, group = ue),
              color = "red",
              linewidth = 1)
## Warning: Removed 664 rows containing missing values or values outside the scale range
## (`geom_line()`).

#-----------------------------------------------------------------------
# Ajustes por condição experimental.

# Ajuste com valores individuais de cada repetição.
n0_cond <-
    tbi |>
    nlme::nlsList(sev ~ (1 - exp(-r * (time - x0))) * (time > x0) | cond,
                  data = _,
                  start = list(r = 0.01, x0 = 50))
summary(n0_cond)
## Call:
##   Model: sev ~ (1 - exp(-r * (time - x0))) * (time > x0) | cond 
##    Data: tbi 
## 
## Coefficients:
##    r 
##                     Estimate  Std. Error  t value     Pr(>|t|)
## 1.CrAaPR-01.27.24 0.04606052 0.011782116 3.909359 2.304052e-04
## 2.CrAaPR-01.27.24 0.03972980 0.008665458 4.584847 1.034662e-03
## 1.CrAaPR-27.27.24 0.01180054 0.004978534 2.370283 3.915898e-02
## 2.CrAaPR-27.27.24 0.03004052 0.007820522 3.841242 2.142797e-05
## 1.CrAaPR-40.27.24 0.10781814 0.031760021 3.394776 1.997721e-09
## 2.CrAaPR-40.27.24 0.05014405 0.013225029 3.791602 5.129499e-03
## 1.CrAlPR-73.27.24 0.05640563 0.013498387 4.178694 1.486783e-04
## 2.CrAlPR-73.27.24 0.01355674 0.002989674 4.534521 4.722470e-05
##    x0 
##                   Estimate Std. Error   t value     Pr(>|t|)
## 1.CrAaPR-01.27.24 57.71086  3.1612759 18.255559 1.517437e-19
## 2.CrAaPR-01.27.24 45.90430  2.8162358 16.299877 2.429773e-14
## 1.CrAaPR-27.27.24 77.61032 10.2381487  7.580503 7.656407e-08
## 2.CrAaPR-27.27.24 52.77070  5.8156301  9.073944 2.818170e-13
## 1.CrAaPR-40.27.24 47.99845  0.9168719 52.350226 7.566490e-46
## 2.CrAaPR-40.27.24 44.03573  2.9348010 15.004672 1.894637e-13
## 1.CrAlPR-73.27.24 46.93523  1.9410264 24.180624 6.537040e-23
## 2.CrAlPR-73.27.24 45.62975  6.7038392  6.806510 4.896645e-08
## 
## Residual standard error: 0.2211068 on 264 degrees of freedom
colMeans(coef(n0_cond), na.rm = TRUE)
##           r          x0 
##  0.04444449 52.32441664
# Ajuste com as médias das repetições em cada tempo.
n0_condm <-
    tbi |>
    group_by(cond, time) |>
    summarise(across("sev", mean)) |>
    ungroup() |>
    nlme::nlsList(sev ~ (1 - exp(-r * (time - x0))) * (time > x0) | cond,
                  data = _,
                  start = list(r = 0.01, x0 = 50))
## `summarise()` has grouped output by 'cond'. You can override using the
## `.groups` argument.
summary(n0_condm)
## Call:
##   Model: sev ~ (1 - exp(-r * (time - x0))) * (time > x0) | cond 
##    Data: ungroup(summarise(group_by(tbi, cond, time), across("sev", mean))) 
## 
## Coefficients:
##    r 
##                     Estimate  Std. Error  t value     Pr(>|t|)
## 1.CrAaPR-01.27.24 0.04606052 0.008023731 5.740537 3.841292e-03
## 2.CrAaPR-01.27.24 0.03973002 0.005901293 6.732427 2.921148e-03
## 1.CrAaPR-27.27.24 0.01180054 0.003390428 3.480545 4.323747e-03
## 2.CrAaPR-27.27.24 0.03004052 0.005325849 5.640514 6.848684e-03
## 1.CrAaPR-40.27.24 0.10781814 0.021628870 4.984918 1.040934e-06
## 2.CrAaPR-40.27.24 0.05014498 0.009006562 5.567605 8.683786e-03
## 1.CrAlPR-73.27.24 0.05640597 0.009192594 6.136024 3.191565e-04
## 2.CrAlPR-73.27.24 0.01355674 0.002035996 6.658530 6.032308e-04
##    x0 
##                   Estimate Std. Error   t value     Pr(>|t|)
## 1.CrAaPR-01.27.24 57.71086  2.1528583 26.806622 2.484269e-06
## 2.CrAaPR-01.27.24 45.90434  1.9178682 23.935084 7.013542e-06
## 1.CrAaPR-27.27.24 77.61032  6.9722747 11.131276 1.848883e-05
## 2.CrAaPR-27.27.24 52.77070  3.9604982 13.324259 1.379152e-04
## 1.CrAaPR-40.27.24 47.99845  0.6243983 76.871528 1.209691e-12
## 2.CrAaPR-40.27.24 44.03587  1.9985777 22.033603 1.479934e-05
## 1.CrAlPR-73.27.24 46.93525  1.3218467 35.507337 5.606377e-08
## 2.CrAlPR-73.27.24 45.62975  4.5653769  9.994739 8.712354e-05
## 
## Residual standard error: 0.06733953 on 40 degrees of freedom
colMeans(coef(n0_condm), na.rm = TRUE)
##           r          x0 
##  0.04444468 52.32444268
# Comparando as estimativas pontuais.
cbind(coef(n0_cond)[["x0"]], coef(n0_condm)[["x0"]])
##          [,1]     [,2]
## [1,] 57.71086 57.71086
## [2,] 45.90430 45.90434
## [3,] 77.61032 77.61032
## [4,] 52.77070 52.77070
## [5,] 47.99845 47.99845
## [6,] 44.03573 44.03587
## [7,] 46.93523 46.93525
## [8,] 45.62975 45.62975
cbind(coef(n0_cond)[["r"]], coef(n0_condm)[["r"]])
##            [,1]       [,2]
## [1,] 0.04606052 0.04606052
## [2,] 0.03972980 0.03973002
## [3,] 0.01180054 0.01180054
## [4,] 0.03004052 0.03004052
## [5,] 0.10781814 0.10781814
## [6,] 0.05014405 0.05014498
## [7,] 0.05640563 0.05640597
## [8,] 0.01355674 0.01355674
# NOTE: As estimativas são iguais porque o número de observações em cada
# tempo é o mesmo.

# Define o modelo para a predição.
# n0_final <- n0_cond
n0_final <- n0_condm

# plot(residuals(n0_final) ~ fitted(n0_final))

# Condições experimentais para a predição.
x_seq <- c(seq(0, max(tbi$time), length.out = 51),
           coef(n0_final)[["x0"]]) |>
    sort()
tb_pred <- tbi |>
    distinct(exp, iso, tem, mol, cond)
tb_pred <- cross_join(tb_pred,
                      data.frame(time = x_seq))
str(tb_pred)
## tibble [472 × 6] (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:472] 27 27 27 27 27 27 27 27 27 27 ...
##  $ mol : num [1:472] 24 24 24 24 24 24 24 24 24 24 ...
##  $ cond: Factor w/ 8 levels "1.CrAaPR-01.27.24",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ time: num [1:472] 0 2.48 4.96 7.44 9.92 ...
# Predição.
tb_pred$fit <- predict(n0_final, newdata = tb_pred)

tbi |>
    ggplot(aes(x = time, y = sev)) +
    facet_grid(exp ~ iso) +
    geom_line(aes(group = rept), alpha = 0.10) +
    stat_summary(geom = "line", fun = "mean",
                 color = "darkgreen",
                 alpha = 0.40) +
    geom_point() +
    labs(title = "Temperatura 27°C e Molhamento 24h") +
    geom_line(data = tb_pred, aes(y = fit, group = cond),
              color = "red",
              linewidth = 1)

# ggsave("sev_time.png", width = 16, height = 9)

#-----------------------------------------------------------------------
# Fitted-Model-Based Annotations.
# browseURL("https://docs.r4photobiology.info/ggpmisc/articles/model-based-annotations.html")

library(ggpmisc)
## Loading required package: ggpp
## Registered S3 methods overwritten by 'ggpp':
##   method                  from   
##   heightDetails.titleGrob ggplot2
##   widthDetails.titleGrob  ggplot2
## 
## Attaching package: 'ggpp'
## 
## The following object is masked from 'package:ggplot2':
## 
##     annotate
# Com esse pacote dá para adicionar a equação do modelo no gráfico além
# de outras medidas de ajuste. Ele refaz o ajuste para cada painel,
# portanto ele não aproveita os ajustes já feitos anteriormente.

# micmen.formula <- y ~ SSmicmen(x, Vm, K)
# ggplot(Puromycin, aes(conc, rate, colour = state)) +
#   facet_wrap(~state) +
#   geom_point() +
#   geom_smooth(method = "nls",
#               formula = micmen.formula,
#               se = FALSE) +
#   stat_fit_tb(method = "nls",
#               method.args = list(formula = micmen.formula),
#               tb.type = "fit.coefs",
#               label.x = 0.9,
#               label.y = c(0.75, 0.2)) +
#   theme(legend.position = "none") +
#   labs(x = "C", y = "V")

# Mostra as estimativas dos parâmetros em uma tabela dentro do gráfico.
ggplot(tbi,
       aes(x = time, y = sev)) +
    # facet_wrap(~cond) +
    facet_grid(exp ~ iso) +
    geom_smooth(
        method = "nls",
        formula = y ~ (1 - (1 * exp(-r * (x - x0)))) * (x > x0),
        method.args = list(start = list(r = 0.01, x0 = 50)),
        se = FALSE
    ) +
    ggpmisc::stat_fit_tb(
        method = "nls",
        method.args = list(
            formula = y ~ (1 - (1 * exp(-r * (x - x0)))) * (x > x0),
            start = list(r = 0.01, x0 = 50)
        ),
        tb.type = "fit.coefs",
        label.x = 0.05,
        label.y = 0.95
    ) +
    geom_point()

#-----------------------------------------------------------------------
# Criando uma tabela com resumo dos ajustes.

tb_r2 <-
    with(tbi, {
    data.frame(cond = cond, obs = sev, fit = predict(n0_cond, tbi)) |>
        group_by(cond) |>
        summarise(r2 = cor(obs, fit)^2)
    })
tb_r2
## # A tibble: 8 × 2
##   cond                 r2
##   <fct>             <dbl>
## 1 1.CrAaPR-01.27.24 0.807
## 2 2.CrAaPR-01.27.24 0.645
## 3 1.CrAaPR-27.27.24 0.286
## 4 2.CrAaPR-27.27.24 0.810
## 5 1.CrAaPR-40.27.24 0.964
## 6 2.CrAaPR-40.27.24 0.652
## 7 1.CrAlPR-73.27.24 0.804
## 8 2.CrAlPR-73.27.24 0.583
tb_coef <-
    n0_cond |>
    coef() |>
    as.data.frame() |>
    rownames_to_column(var = "cond") |>
    left_join(distinct(tbi, cond, exp, iso), by = "cond") |>
    left_join(tb_r2, by = "cond")
tb_coef
##                cond          r       x0 exp       iso        r2
## 1 1.CrAaPR-01.27.24 0.04606052 57.71086   1 CrAaPR-01 0.8070173
## 2 2.CrAaPR-01.27.24 0.03972980 45.90430   2 CrAaPR-01 0.6450326
## 3 1.CrAaPR-27.27.24 0.01180054 77.61032   1 CrAaPR-27 0.2857190
## 4 2.CrAaPR-27.27.24 0.03004052 52.77070   2 CrAaPR-27 0.8100359
## 5 1.CrAaPR-40.27.24 0.10781814 47.99845   1 CrAaPR-40 0.9644831
## 6 2.CrAaPR-40.27.24 0.05014405 44.03573   2 CrAaPR-40 0.6523180
## 7 1.CrAlPR-73.27.24 0.05640563 46.93523   1 CrAlPR-73 0.8038452
## 8 2.CrAlPR-73.27.24 0.01355674 45.62975   2 CrAlPR-73 0.5830175
# Criando a equação.
#   (1 - exp(-r * (time - x0))) * (time > x0)
tb_coef <-
    tb_coef |>
    mutate(
        eq = sprintf(
            "y == (1 - exp(-%0.3f * (x - %0.1f))) * (x > %0.1f) ~~ (R^2 == %0.1f)",
            r, x0, x0, 100 * r2),
        x = 5,
        y = 0.92)
tb_coef
##                cond          r       x0 exp       iso        r2
## 1 1.CrAaPR-01.27.24 0.04606052 57.71086   1 CrAaPR-01 0.8070173
## 2 2.CrAaPR-01.27.24 0.03972980 45.90430   2 CrAaPR-01 0.6450326
## 3 1.CrAaPR-27.27.24 0.01180054 77.61032   1 CrAaPR-27 0.2857190
## 4 2.CrAaPR-27.27.24 0.03004052 52.77070   2 CrAaPR-27 0.8100359
## 5 1.CrAaPR-40.27.24 0.10781814 47.99845   1 CrAaPR-40 0.9644831
## 6 2.CrAaPR-40.27.24 0.05014405 44.03573   2 CrAaPR-40 0.6523180
## 7 1.CrAlPR-73.27.24 0.05640563 46.93523   1 CrAlPR-73 0.8038452
## 8 2.CrAlPR-73.27.24 0.01355674 45.62975   2 CrAlPR-73 0.5830175
##                                                                  eq x    y
## 1 y == (1 - exp(-0.046 * (x - 57.7))) * (x > 57.7) ~~ (R^2 == 80.7) 5 0.92
## 2 y == (1 - exp(-0.040 * (x - 45.9))) * (x > 45.9) ~~ (R^2 == 64.5) 5 0.92
## 3 y == (1 - exp(-0.012 * (x - 77.6))) * (x > 77.6) ~~ (R^2 == 28.6) 5 0.92
## 4 y == (1 - exp(-0.030 * (x - 52.8))) * (x > 52.8) ~~ (R^2 == 81.0) 5 0.92
## 5 y == (1 - exp(-0.108 * (x - 48.0))) * (x > 48.0) ~~ (R^2 == 96.4) 5 0.92
## 6 y == (1 - exp(-0.050 * (x - 44.0))) * (x > 44.0) ~~ (R^2 == 65.2) 5 0.92
## 7 y == (1 - exp(-0.056 * (x - 46.9))) * (x > 46.9) ~~ (R^2 == 80.4) 5 0.92
## 8 y == (1 - exp(-0.014 * (x - 45.6))) * (x > 45.6) ~~ (R^2 == 58.3) 5 0.92
tbi |>
    ggplot(aes(x = time, y = sev)) +
    facet_grid(exp ~ iso) +
    geom_line(aes(group = rept), alpha = 0.10) +
    stat_summary(geom = "line", fun = "mean",
                 color = "darkgreen",
                 alpha = 0.40) +
    geom_point() +
    labs(title = "Temperatura 27°C e Molhamento 24h") +
    geom_line(data = tb_pred, aes(y = fit, group = cond),
              color = "red",
              linewidth = 1) +
    geom_label(data = tb_coef,
              mapping = aes(x = x, y = y, label = eq), parse = TRUE,
              hjust = 0, size = 2.75)

#-----------------------------------------------------------------------