Three-parameter logistic model fitted to plant height as a function of accumulated growing degree-days (GDD), with ANOVA on model parameters and all derivative-based critical points.

Overview

The same logistic framework applied to leaf area is used here for plant height (stem length, cm). Parameters β₁–β₃ are estimated independently for each season × cultivar × block combination and subjected to two-factor ANOVA in randomised complete blocks (FAT2DBC). Critical points derived from the first and second derivatives are also computed at the block level and analysed by ANOVA.

\[ y = \frac{\beta_1}{1 + e^{\;\beta_3(\beta_2 - x)}} \]

Parameter Interpretation
β₁ Upper asymptote (maximum plant height, cm)
β₂ Inflection point (GDD at maximum elongation rate)
β₃ Steepness coefficient (elongation rate)

Required packages

Show code 
library(rio)        # Data import
library(tidyverse)  # Data wrangling and visualisation
library(metan)      # mean_by(), doo()
library(broom)      # tidy() – NLS summary as data frame
library(AgroR)      # FAT2DBC() – two-factor ANOVA in RCBD
library(patchwork)  # Combine ggplot panels
library(hydroGOF)   # gof() – goodness-of-fit metrics
library(gt)         # Professional tables

Load data

Show code 
# Import the pre-computed plot-mean dataset.
# Each row = one season × cultivar × block × date mean (already averaged across plants).
# This dataset is used both for model fitting and for the ANOVA on critical points.
df_model <-
  import("../df_model_cresc_medias.xlsx")

Fit logistic models – plant height

Show code 
# Response variable: 'cp' (comprimento de planta = plant height, cm)
formula_ap <- cp_mean ~ b1 / (1 + exp((b3 * (b2 - gda))))

# Starting values:
#   b1 = 100  → expected maximum height (~100 cm)
#   b2 = 400  → inflection point around 400 GDD
#   b3 = 0.01 → moderate steepness
start_ap <- list(b1 = 100, b2 = 400, b3 = 0.01)

# Fit one model per season × cultivar × block combination.
# The date "04/11" falls outside the vegetative growth window (post-flowering)
# and is excluded to preserve the sigmoid shape of the curve.
mod_ap <-
  df_model |>
  filter(data != "04/11") |>
  group_by(epoca, cultivar, bloco) |>
  doo(~nls(formula_ap, data = ., start = start_ap))

Goodness of fit

Show code 
# Helper function: extract R² and AIC from a single NLS model object
get_gof <- function(model) {
  aic  <- AIC(model)
  bic  <- BIC(model)
  fit  <- model$m$fitted()
  res  <- model$m$resid()
  obs  <- fit + res
  
  g <- hydroGOF::gof(sim = fit, obs = obs, digits = 4)
  
  r2   <- g["R2", ]
  rmse <- g["RMSE", ]
  mae  <- g["MAE", ]
  d    <- g["d", ]
  
  data.frame(aic = aic, bic = bic, r2 = r2, rmse = rmse, mae = mae, d = d)
}

# Apply to all models
gof_ap <-
  mod_ap |>
  mutate(map_dfr(.x = data, .f = ~get_gof(.))) |>
  select(-data) |>
  mutate(variable = "plant_height")

# Display goodness of fit table with gt
gof_ap |> 
  gt() |> 
  tab_header(
    title = "Supplementary Table S4",
    subtitle = "Goodness-of-fit metrics for Plant Height logistic models (per block)"
  ) |> 
  fmt_number(columns = where(is.numeric), decimals = 4)
Supplementary Table S4
Goodness-of-fit metrics for Plant Height logistic models (per block)
epoca bloco cultivar aic bic r2 rmse mae d variable
E1 B1 D 66.3935 68.3331 0.9785 2.7569 2.3200 0.9946 plant_height
E1 B2 D 61.7407 63.6803 0.9889 2.2711 1.7191 0.9972 plant_height
E1 B3 D 72.5145 74.4541 0.9684 3.5579 2.8977 0.9921 plant_height
E1 B1 M 76.0310 77.9706 0.9713 4.1193 3.3420 0.9928 plant_height
E1 B2 M 83.5905 85.5302 0.9499 5.6444 4.5812 0.9876 plant_height
E1 B3 M 81.8164 83.7560 0.9551 5.2422 3.9575 0.9888 plant_height
E2 B1 D 60.7096 61.9200 0.9675 3.3755 2.5272 0.9918 plant_height
E2 B2 D 48.5308 49.7412 0.9880 1.8360 1.5194 0.9970 plant_height
E2 B3 D 74.0980 75.3084 0.9032 6.5927 5.4324 0.9757 plant_height
E2 B1 M 68.6819 69.8922 0.9411 5.0287 3.2778 0.9848 plant_height
E2 B2 M 59.2642 60.4745 0.9600 3.1401 2.7115 0.9900 plant_height
E2 B3 M 63.7990 65.0093 0.9472 3.9393 2.3528 0.9863 plant_height
Show code 

saveRDS(gof_ap, "../results/gof_ap.rds")

Extract estimated parameters and compute critical points

Show code 
# broom::tidy() converts each NLS object into a tidy data frame;
# pivot_wider() places b1, b2, b3 as separate columns.
parameters_ap <-
  mod_ap |>
  mutate(data = map(data, ~.x |> tidy())) |>
  unnest(data) |>
  select(epoca, cultivar, bloco, term, estimate) |>
  pivot_wider(names_from  = term,
              values_from = estimate)

Derivative functions

Show code 
# Symbolically verify the first derivative
D(expression(b1 / (1 + exp((b3 * (b2 - x))))), "x")
## b1 * (exp((b3 * (b2 - x))) * b3)/(1 + exp((b3 * (b2 - x))))^2

# First derivative: instantaneous elongation rate at a given GDD (x)
dy <- function(x, b1, b2, b3) {
  b1 * (exp((b3 * (b2 - x))) * b3) / (1 + exp((b3 * (b2 - x))))^2
}

# Symbolically verify the second derivative
D(expression(b1 * (exp((b3 * (b2 - x))) * b3) / (1 + exp((b3 * (b2 - x))))^2), "x")
## -(b1 * (exp((b3 * (b2 - x))) * b3 * b3)/(1 + exp((b3 * (b2 - 
##     x))))^2 - b1 * (exp((b3 * (b2 - x))) * b3) * (2 * (exp((b3 * 
##     (b2 - x))) * b3 * (1 + exp((b3 * (b2 - x))))))/((1 + exp((b3 * 
##     (b2 - x))))^2)^2)

# Second derivative: rate of change of elongation rate (acceleration / deceleration)
d2y <- function(x, b1, b2, b3) {
  -(b1 * (exp((b3 * (b2 - x))) * b3 * b3) / (1 + exp((b3 * (b2 - x))))^2 -
      b1 * (exp((b3 * (b2 - x))) * b3) *
      (2 * (exp((b3 * (b2 - x))) * b3 * (1 + exp((b3 * (b2 - x)))))) /
      ((1 + exp((b3 * (b2 - x))))^2)^2)
}

Critical points at block level

Show code 
# Compute ALL critical points at the individual block level.
# This preserves the replication structure required for ANOVA.
#
# Critical points derived from the logistic model:
#   xpi  = β₂                              → GDD at maximum elongation rate (inflection)
#   ypi  = dy(β₂)                          → maximum elongation rate value
#   xmap = β₂ + ln(2 + √3) / β₃           → GDD at maximum acceleration point (MAP)
#   xmdp = β₂ − ln(2 + √3) / β₃           → GDD at maximum deceleration point (MDP)
#   ymap = d²y(xmap)                        → second-derivative value at MAP
#   ymdp = d²y(xmdp)                        → second-derivative value at MDP
parameters_ap <-
  parameters_ap |>
  mutate(
    xpi  = b2,
    ypi  = dy(xpi, b1, b2, b3),
    xmap = b2 - (log(2 + sqrt(3)) / b3),
    xmdp = b2 + (log(2 + sqrt(3)) / b3),
    ymap = d2y(xmap, b1, b2, b3),
    ymdp = d2y(xmdp, b1, b2, b3)
  ) |>
  as.data.frame() |>
  mutate(variable = "plant_height")

# Block-level means for derivative plots (used in stat_function args)
parameters_ap_mean <-
  parameters_ap |>
  metan::mean_by(epoca, cultivar)

# Factor-level means (Sowing Season and Cultivar)
# Since interaction is often non-significant, we present main effects
means_epoca_ap <- 
  parameters_ap |> 
  metan::mean_by(epoca) |> 
  mutate(Factor = "Sowing Season") |> 
  rename(Level = epoca)

means_cultivar_ap <- 
  parameters_ap |> 
  metan::mean_by(cultivar) |> 
  mutate(Factor = "Cultivar") |> 
  rename(Level = cultivar)

means_summary_ap <- 
  bind_rows(means_epoca_ap, means_cultivar_ap) |> 
  relocate(Factor, Level)



# Display means summary using gt
means_summary_ap |> 
  gt(groupname_col = "Factor") |> 
  tab_header(
    title = "Supplementary Table S5", 
    subtitle = "Main Effect Means: Model parameters and critical points for Plant Height"
  ) |> 
  fmt_number(columns = where(is.numeric), decimals = 5)
Supplementary Table S5
Main Effect Means: Model parameters and critical points for Plant Height
Level b1 b2 b3 xpi ypi xmap xmdp ymap ymdp
Sowing Season
E1 92.66452 685.49906 0.00456 685.49906 0.10362 386.76668 984.23144 0.00018 −0.00018
E2 106.82804 587.90582 0.00507 587.90582 0.12165 307.24299 868.56865 0.00023 −0.00023
Cultivar
D 116.20156 722.50312 0.00458 722.50312 0.12090 406.22671 1,038.77953 0.00021 −0.00021
M 83.29100 550.90176 0.00506 550.90176 0.10437 287.78296 814.02056 0.00020 −0.00020
Show code 

saveRDS(means_summary_ap, "../results/means_summary_ap.rds")

ANOVA on model parameters and critical points

Show code 
# List of variables to analyze
vars_ap <- c("b1", "b2", "b3", "xpi", "ypi", "xmap", "xmdp")

# 1. Run the ANOVA for all variables and store the full objects in a list
# This allows access to means tests and other FAT2DBC outputs later.
mod_anova_ap <- 
  vars_ap |> 
  set_names() |> 
  map(~with(parameters_ap, FAT2DBC(epoca, cultivar, bloco, get(.x))))
## 
## -----------------------------------------------------------------
## Normality of errors
## -----------------------------------------------------------------
##                          Method Statistic   p.value
##  Shapiro-Wilk normality test(W) 0.9713285 0.9242106
## 
## -----------------------------------------------------------------
## Homogeneity of Variances
## -----------------------------------------------------------------
##                               Method Statistic   p.value
##  Bartlett test(Bartlett's K-squared)  2.498853 0.4754984
## 
## -----------------------------------------------------------------
## Independence from errors
## -----------------------------------------------------------------
##                  Method Statistic    p.value
##  Durbin-Watson test(DW)  1.522037 0.04147825
## 
## -----------------------------------------------------------------
## Additional Information
## -----------------------------------------------------------------
## 
## CV (%) =  34.18
## Mean =  99.7463
## Median =  90.1007
## Possible outliers =  No discrepant point
## 
## -----------------------------------------------------------------
## Analysis of Variance
## -----------------------------------------------------------------
##           Df    Sum Sq   Mean.Sq   F value     Pr(F)
## F1         1  601.8165  601.8165 0.5176079 0.4989154
## F2         1 3249.3144 3249.3144 2.7946572 0.1456110
## Block      2 5952.2038 2976.1019 2.5596737 0.1571143
## F1 × F2    1 2274.2512 2274.2512 1.9560288 0.2114456
## Residuals  6 6976.1279 1162.6880
## -----------------------------------------------------------------
## No significant interaction
## -----------------------------------------------------------------
## 
## -----------------------------------------------------------------
## Isolated factors not significant
## -----------------------------------------------------------------
##            D       M
## E1  95.35313 89.9759
## E2 137.04998 76.6061
## 
## -----------------------------------------------------------------
## Normality of errors
## -----------------------------------------------------------------
##                          Method Statistic   p.value
##  Shapiro-Wilk normality test(W) 0.9291302 0.3709843
## 
## -----------------------------------------------------------------
## Homogeneity of Variances
## -----------------------------------------------------------------
##                               Method Statistic   p.value
##  Bartlett test(Bartlett's K-squared)   2.84019 0.4169267
## 
## -----------------------------------------------------------------
## Independence from errors
## -----------------------------------------------------------------
##                  Method Statistic    p.value
##  Durbin-Watson test(DW)  1.548727 0.04802587
## 
## -----------------------------------------------------------------
## Additional Information
## -----------------------------------------------------------------
## 
## CV (%) =  22.07
## Mean =  636.7024
## Median =  592.9882
## Possible outliers =  No discrepant point
## 
## -----------------------------------------------------------------
## Analysis of Variance
## -----------------------------------------------------------------
##           Df     Sum Sq   Mean.Sq   F value      Pr(F)
## F1         1  28573.321 28573.321 1.4475777 0.27423233
## F2         1  88341.080 88341.080 4.4755236 0.07876703
## Block      2 125466.557 62733.278 3.1781847 0.11449357
## F1 × F2    1   9600.018  9600.018 0.4863548 0.51166195
## Residuals  6 118432.284 19738.714
## -----------------------------------------------------------------
## No significant interaction
## -----------------------------------------------------------------
## 
## -----------------------------------------------------------------
## Isolated factors not significant
## -----------------------------------------------------------------
##           D        M
## E1 743.0154 627.9827
## E2 701.9908 473.8208
## 
## -----------------------------------------------------------------
## Normality of errors
## -----------------------------------------------------------------
##                          Method Statistic   p.value
##  Shapiro-Wilk normality test(W) 0.9487436 0.6186954
## 
## -----------------------------------------------------------------
## Homogeneity of Variances
## -----------------------------------------------------------------
##                               Method Statistic    p.value
##  Bartlett test(Bartlett's K-squared)  7.243356 0.06453268
## 
## -----------------------------------------------------------------
## Independence from errors
## -----------------------------------------------------------------
##                  Method Statistic   p.value
##  Durbin-Watson test(DW)  2.071861 0.2845967
## 
## -----------------------------------------------------------------
## Additional Information
## -----------------------------------------------------------------
## 
## CV (%) =  28.33
## Mean =  0.0048
## Median =  0.0047
## Possible outliers =  No discrepant point
## 
## -----------------------------------------------------------------
## Analysis of Variance
## -----------------------------------------------------------------
##           Df       Sum Sq      Mean.Sq   F value     Pr(F)
## F1         1 7.815089e-07 7.815089e-07 0.4191803 0.5413109
## F2         1 6.933254e-07 6.933254e-07 0.3718810 0.5643458
## Block      2 5.706489e-06 2.853244e-06 1.5304032 0.2903710
## F1 × F2    1 1.372507e-06 1.372507e-06 0.7361757 0.4238321
## Residuals  6 1.118625e-05 1.864374e-06                    
## 
## -----------------------------------------------------------------
## No significant interaction
## -----------------------------------------------------------------
## 
## -----------------------------------------------------------------
## Isolated factors not significant
## -----------------------------------------------------------------
##              D           M
## E1 0.003985501 0.005142627
## E2 0.005172285 0.004976633
## 
## -----------------------------------------------------------------
## Normality of errors
## -----------------------------------------------------------------
##                          Method Statistic   p.value
##  Shapiro-Wilk normality test(W) 0.9291302 0.3709843
## 
## -----------------------------------------------------------------
## Homogeneity of Variances
## -----------------------------------------------------------------
##                               Method Statistic   p.value
##  Bartlett test(Bartlett's K-squared)   2.84019 0.4169267
## 
## -----------------------------------------------------------------
## Independence from errors
## -----------------------------------------------------------------
##                  Method Statistic    p.value
##  Durbin-Watson test(DW)  1.548727 0.04802587
## 
## -----------------------------------------------------------------
## Additional Information
## -----------------------------------------------------------------
## 
## CV (%) =  22.07
## Mean =  636.7024
## Median =  592.9882
## Possible outliers =  No discrepant point
## 
## -----------------------------------------------------------------
## Analysis of Variance
## -----------------------------------------------------------------
##           Df     Sum Sq   Mean.Sq   F value      Pr(F)
## F1         1  28573.321 28573.321 1.4475777 0.27423233
## F2         1  88341.080 88341.080 4.4755236 0.07876703
## Block      2 125466.557 62733.278 3.1781847 0.11449357
## F1 × F2    1   9600.018  9600.018 0.4863548 0.51166195
## Residuals  6 118432.284 19738.714
## -----------------------------------------------------------------
## No significant interaction
## -----------------------------------------------------------------
## 
## -----------------------------------------------------------------
## Isolated factors not significant
## -----------------------------------------------------------------
##           D        M
## E1 743.0154 627.9827
## E2 701.9908 473.8208
## 
## -----------------------------------------------------------------
## Normality of errors
## -----------------------------------------------------------------
##                          Method Statistic   p.value
##  Shapiro-Wilk normality test(W)  0.970736 0.9183601
## 
## -----------------------------------------------------------------
## Homogeneity of Variances
## -----------------------------------------------------------------
##                               Method Statistic   p.value
##  Bartlett test(Bartlett's K-squared)  4.453592 0.2164673
## 
## -----------------------------------------------------------------
## Independence from errors
## -----------------------------------------------------------------
##                  Method Statistic   p.value
##  Durbin-Watson test(DW)  2.491596 0.5443447
## 
## -----------------------------------------------------------------
## Additional Information
## -----------------------------------------------------------------
## 
## CV (%) =  18.55
## Mean =  0.1126
## Median =  0.1106
## Possible outliers =  No discrepant point
## 
## -----------------------------------------------------------------
## Analysis of Variance
## -----------------------------------------------------------------
##           Df       Sum Sq      Mean.Sq   F value      Pr(F)
## F1         1 0.0009746050 0.0009746050  2.233622 0.18565893
## F2         1 0.0008200382 0.0008200382  1.879382 0.21947169
## Block      2 0.0010670424 0.0005335212  1.222736 0.35857660
## F1 × F2    1 0.0044215256 0.0044215256 10.133354 0.01899663
## Residuals  6 0.0026180034 0.0004363339                     
## 
## -----------------------------------------------------------------
## 
## Significant interaction: analyzing the interaction
## 
## -----------------------------------------------------------------
## 
## -----------------------------------------------------------------
## Analyzing  F1  inside of each level of  F2
## -----------------------------------------------------------------
##              Df    Sum Sq   Mean Sq F value  Pr(>F)  
## Block         2 0.0010670 0.0005335  1.2227 0.35858  
## F2            1 0.0008200 0.0008200  1.8794 0.21947  
## F1 × F2 + F1  2 0.0053961 0.0026981  6.1835 0.03486 *
##    F1:F2 D    1 0.0047739 0.0047739 10.9410 0.01625 *
##    F1:F2 M    1 0.0006222 0.0006222  1.4260 0.27749  
## Residuals     6 0.0026180 0.0004363                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -----------------------------------------------------------------
## Analyzing  F2  inside of the level of  F1
## -----------------------------------------------------------------
## 
##              Df    Sum Sq   Mean Sq F value  Pr(>F)  
## Block         2 0.0010670 0.0005335  1.2227 0.35858  
## F1            1 0.0009746 0.0009746  2.2336 0.18566  
## F1 × F2 + F2  2 0.0052416 0.0026208  6.0064 0.03696 *
##    F2:F1 E1   1 0.0007166 0.0007166  1.6424 0.24730  
##    F2:F1 E2   1 0.0045249 0.0045249 10.3704 0.01813 *
## Residuals     6 0.0026180 0.0004363                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## -----------------------------------------------------------------
## Final table
## -----------------------------------------------------------------
##            D         M
## E1 0.0927 bA 0.1146 aA
## E2 0.1491 aA 0.0942 aB
## 
## -----------------------------------------------------------------
## Normality of errors
## -----------------------------------------------------------------
##                          Method Statistic   p.value
##  Shapiro-Wilk normality test(W) 0.9450673 0.5663506
## 
## -----------------------------------------------------------------
## Homogeneity of Variances
## -----------------------------------------------------------------
##                               Method Statistic   p.value
##  Bartlett test(Bartlett's K-squared)  2.022598 0.5677298
## 
## -----------------------------------------------------------------
## Independence from errors
## -----------------------------------------------------------------
##                  Method Statistic    p.value
##  Durbin-Watson test(DW)  1.438158 0.02480983
## 
## -----------------------------------------------------------------
## Additional Information
## -----------------------------------------------------------------
## 
## CV (%) =  24.96
## Mean =  347.0048
## Median =  352.688
## Possible outliers =  No discrepant point
## 
## -----------------------------------------------------------------
## Analysis of Variance
## -----------------------------------------------------------------
##           Df   Sum Sq   Mean.Sq  F value      Pr(F)
## F1         1 18972.05 18972.052 2.528423 0.16291476
## F2         1 42086.77 42086.766 5.608943 0.05564813
## Block      2 51134.36 25567.181 3.407362 0.10264219
## F1 × F2    1 19967.73 19967.733 2.661119 0.15394946
## Residuals  6 45021.07  7503.511
## -----------------------------------------------------------------
## No significant interaction
## -----------------------------------------------------------------
## 
## -----------------------------------------------------------------
## Isolated factors not significant
## -----------------------------------------------------------------
##           D        M
## E1 405.1967 368.3367
## E2 407.2567 207.2292
## 
## -----------------------------------------------------------------
## Normality of errors
## -----------------------------------------------------------------
##                          Method Statistic   p.value
##  Shapiro-Wilk normality test(W) 0.9371275 0.4617801
## 
## -----------------------------------------------------------------
## Homogeneity of Variances
## -----------------------------------------------------------------
##                               Method Statistic   p.value
##  Bartlett test(Bartlett's K-squared)  3.316516 0.3453506
## 
## -----------------------------------------------------------------
## Independence from errors
## -----------------------------------------------------------------
##                  Method Statistic   p.value
##  Durbin-Watson test(DW)  1.652061 0.0790389
## 
## -----------------------------------------------------------------
## Additional Information
## -----------------------------------------------------------------
## 
## CV (%) =  21.45
## Mean =  926.4
## Median =  857.6216
## Possible outliers =  No discrepant point
## 
## -----------------------------------------------------------------
## Analysis of Variance
## -----------------------------------------------------------------
##           Df     Sum Sq    Mean.Sq    F value      Pr(F)
## F1         1  40133.643  40133.643 1.01665553 0.35223181
## F2         1 151549.781 151549.781 3.83902160 0.09778429
## Block      2 232680.404 116340.202 2.94710125 0.12836534
## F1 × F2    1   2986.855   2986.855 0.07566228 0.79249092
## Residuals  6 236856.882  39476.147                      
## 
## -----------------------------------------------------------------
## No significant interaction
## -----------------------------------------------------------------
## 
## -----------------------------------------------------------------
## Isolated factors not significant
## -----------------------------------------------------------------
##            D        M
## E1 1080.8342 887.6287
## E2  996.7248 740.4125

# 2. Extract and combine tidy ANOVA tables for the summary
anova_ap <- 
  mod_anova_ap |> 
  map_dfr(~{
    .x$anova |> 
      broom::tidy()
  }, .id = "Parameter") |> 
  mutate(term = str_replace(term, "Fator1", "Season"),
         term = str_replace(term, "Fator2", "Cultivar"))

# Display the summary table with gt
anova_ap |> 
  filter(term != "Residuals", Parameter != "xpi") |> 
  gt() |> 
  tab_header(
    title = "Supplementary Table S6",
    subtitle = "ANOVA Summary: Model parameters and critical points for Plant Height"
  ) |> 
  fmt_number(columns = where(is.numeric), decimals = 5) |> 
  tab_style(
    style = cell_text(color = "red", weight = "bold"),
    locations = cells_body(
      rows = p.value < 0.05
    )
  )
Supplementary Table S6
ANOVA Summary: Model parameters and critical points for Plant Height
Parameter term GL SQ QM Fcal p.value
b1 Season 1.00000 601.81645 601.81645 0.51761 0.49892
b1 Cultivar 1.00000 3,249.31436 3,249.31436 2.79466 0.14561
b1 bloco 2.00000 5,952.20377 2,976.10189 2.55967 0.15711
b1 Season:Cultivar 1.00000 2,274.25119 2,274.25119 1.95603 0.21145
b2 Season 1.00000 28,573.32138 28,573.32138 1.44758 0.27423
b2 Cultivar 1.00000 88,341.07965 88,341.07965 4.47552 0.07877
b2 bloco 2.00000 125,466.55664 62,733.27832 3.17818 0.11449
b2 Season:Cultivar 1.00000 9,600.01759 9,600.01759 0.48635 0.51166
b3 Season 1.00000 0.00000 0.00000 0.41918 0.54131
b3 Cultivar 1.00000 0.00000 0.00000 0.37188 0.56435
b3 bloco 2.00000 0.00001 0.00000 1.53040 0.29037
b3 Season:Cultivar 1.00000 0.00000 0.00000 0.73618 0.42383
ypi Season 1.00000 0.00097 0.00097 2.23362 0.18566
ypi Cultivar 1.00000 0.00082 0.00082 1.87938 0.21947
ypi bloco 2.00000 0.00107 0.00053 1.22274 0.35858
ypi Season:Cultivar 1.00000 0.00442 0.00442 10.13335 0.01900
xmap Season 1.00000 18,972.05156 18,972.05156 2.52842 0.16291
xmap Cultivar 1.00000 42,086.76633 42,086.76633 5.60894 0.05565
xmap bloco 2.00000 51,134.36127 25,567.18064 3.40736 0.10264
xmap Season:Cultivar 1.00000 19,967.73288 19,967.73288 2.66112 0.15395
xmdp Season 1.00000 40,133.64309 40,133.64309 1.01666 0.35223
xmdp Cultivar 1.00000 151,549.78106 151,549.78106 3.83902 0.09778
xmdp bloco 2.00000 232,680.40415 116,340.20207 2.94710 0.12837
xmdp Season:Cultivar 1.00000 2,986.85534 2,986.85534 0.07566 0.79249
Show code 

saveRDS(anova_ap, "../results/anova_ap.rds")

Fitted logistic curves

Show code 
# Logistic formula in x/y notation (required by geom_smooth internals)
formula_xy <- y ~ b1 / (1 + exp((b3 * (b2 - x))))

ggplot(
  df_model |> filter(!data %in% c("04/11")),
  aes(gda, cp_mean, color = cultivar)
) +
  geom_smooth(
    method      = "nls",
    method.args = list(
      formula = formula_xy,
      start   = c(b1 = 60, b2 = 400, b3 = 0.01)
    ),
    se  = FALSE,
    aes(color = cultivar)
  ) +
  facet_wrap(~epoca, ncol = 1) +
  stat_summary(
    fun      = mean,
    geom     = "point",
    aes(color = cultivar),
    size     = 2.5,
    position = position_dodge(width = 0.8)
  ) +
  scale_y_continuous(breaks = seq(0, 150, by = 25)) +
  labs(
    x     = "Accumulated Degree Days",
    y     = "Plant height (cm)",
    color = "Cultivar"
  ) +
  scale_color_manual(values = c("gold", "brown")) +
  # Excluded point shown transparently for reference
  stat_summary(
    data        = df_model |> filter(data == "04/11"),
    aes(x = gda, y = cp_mean),
    fun         = mean,
    alpha       = 0.4,
    shape       = 16,
    show.legend = FALSE
  ) +
  theme_bw(base_size = 18) +
  theme(
    panel.grid.major  = element_blank(),
    legend.background = element_rect(fill = "transparent"),
    legend.position   = "bottom"
  )

Logistic curves fitted to plant height (cm) as a function of accumulated GDD, faceted by sowing season. Points = treatment means. The semi-transparent point (04/11) was excluded from fitting.
Show code 

ggsave("../figs/mod_logistico_ap.jpg", width = 6, height = 9)

Derivative plots

Show code 

# --- First derivative plot ---
plot_pi_ap <-
  ggplot() +
  stat_function(fun = dy, aes(color = "D", linetype = "E1"), n = 500, linewidth = 1,
                xlim = c(0, 1000),
                args = c(b1 = parameters_ap_mean[[1, 3]],
                         b2 = parameters_ap_mean[[1, 4]],
                         b3 = parameters_ap_mean[[1, 5]])) +
  stat_function(fun = dy, aes(color = "M", linetype = "E1"), linewidth = 1,
                xlim = c(0, 1000),
                args = c(b1 = parameters_ap_mean[[2, 3]],
                         b2 = parameters_ap_mean[[2, 4]],
                         b3 = parameters_ap_mean[[2, 5]])) +
  stat_function(fun = dy, aes(color = "D", linetype = "E2"), linewidth = 1,
                xlim = c(0, 1000),
                args = c(b1 = parameters_ap_mean[[3, 3]],
                         b2 = parameters_ap_mean[[3, 4]],
                         b3 = parameters_ap_mean[[3, 5]])) +
  stat_function(fun = dy, aes(color = "M", linetype = "E2"), linewidth = 1,
                xlim = c(0, 1000),
                args = c(b1 = parameters_ap_mean[[4, 3]],
                         b2 = parameters_ap_mean[[4, 4]],
                         b3 = parameters_ap_mean[[4, 5]])) +
  # Mark the inflection point for each treatment
  geom_point(aes(xpi, ypi, shape = epoca, color = cultivar),
             data = parameters_ap_mean, size = 3, show.legend = FALSE) +
  theme_bw(base_size = 20) +
  scale_x_continuous(breaks = seq(0, 1200, by = 200)) +
  labs(
    x        = "Accumulated Degree Days",
    y        = expression(Plant ~ height ~ (cm ~ plant^{-1} ~ degree ~ C ~ day^{-1})),
    color    = "Cultivar",
    linetype = "Season"
  ) +
  theme(legend.position = "bottom", panel.grid.minor = element_blank()) +
  scale_color_manual(values = c("gold", "brown"))

# --- Second derivative plot ---
plot_sd_ap <-
  ggplot() +
  geom_hline(yintercept = 0) +   # zero line: above = acceleration, below = deceleration
  stat_function(fun = d2y, aes(color = "D", linetype = "E1"), n = 500, linewidth = 1,
                xlim = c(0, 1000),
                args = c(b1 = parameters_ap_mean[[1, 3]],
                         b2 = parameters_ap_mean[[1, 4]],
                         b3 = parameters_ap_mean[[1, 5]])) +
  stat_function(fun = d2y, aes(color = "M", linetype = "E1"), n = 500, linewidth = 1,
                xlim = c(0, 1000),
                args = c(b1 = parameters_ap_mean[[2, 3]],
                         b2 = parameters_ap_mean[[2, 4]],
                         b3 = parameters_ap_mean[[2, 5]])) +
  stat_function(fun = d2y, aes(color = "D", linetype = "E2"), n = 500, linewidth = 1,
                xlim = c(0, 1000),
                args = c(b1 = parameters_ap_mean[[3, 3]],
                         b2 = parameters_ap_mean[[3, 4]],
                         b3 = parameters_ap_mean[[3, 5]])) +
  stat_function(fun = d2y, aes(color = "M", linetype = "E2"), n = 500, linewidth = 1,
                xlim = c(0, 1000),
                args = c(b1 = parameters_ap_mean[[4, 3]],
                         b2 = parameters_ap_mean[[4, 4]],
                         b3 = parameters_ap_mean[[4, 5]])) +
  # MAP (triangle up) and MDP (triangle down) markers
  geom_point(aes(xmap, ymap, fill = cultivar), data = parameters_ap_mean,
             size = 3, shape = 24, show.legend = FALSE) +
  geom_point(aes(xmdp, ymdp, fill = cultivar), data = parameters_ap_mean,
             size = 3, shape = 25, show.legend = FALSE) +
  theme_bw(base_size = 20) +
  labs(
    x        = "Accumulated Degree Days",
    y        = expression(Plant ~ height ~ (cm ~ plant^{-1} ~ degree ~ C ~ day^{-2})),
    color    = "Cultivar",
    linetype = "Season"
  ) +
  theme(legend.position = "bottom", panel.grid.minor = element_blank()) +
  scale_color_manual(values = c("gold", "brown")) +
  scale_fill_manual(values  = c("gold", "brown"))

plot_pi_ap / plot_sd_ap +
  plot_layout(guides = "collect") &
  theme(legend.position = "bottom")

First (left) and second (right) derivatives of the fitted logistic curves for plant height. Points mark the inflection (●), maximum acceleration (▲), and maximum deceleration (▽) thresholds per cultivar × season combination.
Show code 

ggsave("../figs/deriv_ap.jpg", width = 8, height = 12)

Correlation between model parameters

Show code 
# Correlation between parameters and critical points

  parameters_ap |> 
  select(all_of(vars_ap)) |> 
  metan::corr_coef() |> 
  plot()