Three-parameter logistic model fitted to cumulative leaf number per plant as a function of accumulated GDD, with ANOVA on parameters and all derivative-based critical points.

Overview

Leaf number (leaves plant⁻¹) follows a sigmoidal accumulation pattern driven by thermal time. The logistic model uses a higher asymptote (β₁ ≈ 150 leaves) and a shallower steepness (β₃ ≈ 0.005) relative to leaf area and plant height, reflecting slower saturation kinetics of leaf appearance.

Critical points from the first and second derivatives are computed at the block level and subjected to two-factor ANOVA, consistent with the approach used for leaf area and plant height.

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

Parameter Interpretation
β₁ Upper asymptote (maximum leaf count, leaves plant⁻¹)
β₂ Inflection point (GDD at maximum leaf appearance rate)
β₃ Steepness coefficient (leaf appearance 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).
df_model <-
  import("../df_model_cresc_medias.xlsx")

Fit logistic models – leaf number

Show code 
# Response variable: 'n' (number of leaves per plant)
formula_nf <- n_mean ~ b1 / (1 + exp((b3 * (b2 - gda))))

# Starting values:
#   b1 = 150  → higher asymptote than leaf area (more leaves)
#   b2 = 400  → inflection around 400 GDD
#   b3 = 0.005 → shallower sigmoid (slower saturation of leaf count)
start_nf <- list(b1 = 150, b2 = 400, b3 = 0.005)

# Increase maxiter: the default 50 iterations is often insufficient for
# flat curves with small b3 values.
mod_nf <-
  df_model |>
  filter(data != "04/11") |>
  group_by(epoca, cultivar, bloco) |>
  doo(~nls(formula_nf,
           data    = .,
           start   = start_nf,
           control = list(maxiter = 500)))

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

gof_nf <-
  mod_nf |>
  mutate(map_dfr(.x = data, .f = ~get_gof(.))) |>
  select(-data) |>
  mutate(variable = "leaf_number")

# Display goodness of fit table with gt
gof_nf |> 
  gt() |> 
  tab_header(
    title = "Supplementary Table S7",
    subtitle = "Goodness-of-fit metrics for Leaf Number logistic models (per block)"
  ) |> 
  fmt_number(columns = where(is.numeric), decimals = 4)
Supplementary Table S7
Goodness-of-fit metrics for Leaf Number logistic models (per block)
epoca bloco cultivar aic bic r2 rmse mae d variable
E1 B1 D 93.4669 95.4065 0.9478 8.5180 7.2684 0.9866 leaf_number
E1 B2 D 110.3749 112.3146 0.8012 17.2306 12.5903 0.9434 leaf_number
E1 B3 D 104.4329 106.3725 0.8700 13.4516 9.9129 0.9638 leaf_number
E1 B1 M 95.0313 96.9710 0.9374 9.0917 7.2052 0.9839 leaf_number
E1 B2 M 96.9291 98.8688 0.9479 9.8399 8.3223 0.9865 leaf_number
E1 B3 M 105.9092 107.8488 0.8256 14.3050 11.5220 0.9509 leaf_number
E2 B1 D 76.8690 78.0793 0.9293 7.5724 4.6942 0.9814 leaf_number
E2 B2 D 92.3730 93.5833 0.7167 16.4399 11.5470 0.9155 leaf_number
E2 B3 D 79.2999 80.5103 0.9137 8.5511 6.9266 0.9772 leaf_number
E2 B1 M 82.5137 83.7241 0.9058 10.0417 8.3373 0.9750 leaf_number
E2 B2 M 88.9461 90.1564 0.8313 13.8511 9.0146 0.9534 leaf_number
E2 B3 M 71.3668 72.5771 0.9656 5.7512 4.8502 0.9913 leaf_number
Show code 

saveRDS(gof_nf, "../results/gof_nf.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_nf <-
  mod_nf |>
  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 leaf appearance 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: acceleration / deceleration of leaf appearance rate
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 leaf appearance rate (inflection)
#   ypi  = dy(β₂)                          → maximum leaf appearance 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_nf <-
  parameters_nf |>
  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 = "leaf_number")

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

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

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

means_summary_nf <- 
  bind_rows(means_epoca_nf, means_cultivar_nf) |> 
  relocate(Factor, Level)



# Display means summary using gt
means_summary_nf |> 
  gt(groupname_col = "Factor") |> 
  tab_header(
    title = "Supplementary Table S8", 
    subtitle = "Main Effect Means: Model parameters and critical points for Leaf Number"
  ) |> 
  fmt_number(columns = where(is.numeric), decimals = 5)
Supplementary Table S8
Main Effect Means: Model parameters and critical points for Leaf Number
Level b1 b2 b3 xpi ypi xmap xmdp ymap ymdp
Sowing Season
E1 120.03457 281.03311 0.00635 281.03311 0.18777 69.04974 493.01649 0.00046 −0.00046
E2 101.01244 172.69364 0.00894 172.69364 0.22367 20.65962 324.72765 0.00079 −0.00079
Cultivar
D 105.12950 218.11180 0.00786 218.11180 0.20241 44.37288 391.85072 0.00062 −0.00062
M 115.91751 235.61495 0.00743 235.61495 0.20904 45.33648 425.89342 0.00063 −0.00063
Show code 

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

ANOVA on model parameters and critical points

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

# 1. Run the ANOVA for all variables and store the full objects in a list
mod_anova_nf <- 
  vars_nf |> 
  set_names() |> 
  map(~with(parameters_nf, FAT2DBC(epoca, cultivar, bloco, get(.x),
                                   names.fat = c("Season", "Cultivar"))))
## 
## -----------------------------------------------------------------
## Normality of errors
## -----------------------------------------------------------------
##                          Method Statistic   p.value
##  Shapiro-Wilk normality test(W) 0.9256118 0.3359078
## 
## -----------------------------------------------------------------
## Homogeneity of Variances
## -----------------------------------------------------------------
##                               Method Statistic   p.value
##  Bartlett test(Bartlett's K-squared)   3.14787 0.3693836
## 
## -----------------------------------------------------------------
## Independence from errors
## -----------------------------------------------------------------
##                  Method Statistic   p.value
##  Durbin-Watson test(DW)  3.169048 0.9189987
## 
## -----------------------------------------------------------------
## Additional Information
## -----------------------------------------------------------------
## 
## CV (%) =  11.97
## Mean =  110.5235
## Median =  108.0218
## Possible outliers =  No discrepant point
## 
## -----------------------------------------------------------------
## Analysis of Variance
## -----------------------------------------------------------------
##                   Df      Sum Sq     Mean.Sq    F value      Pr(F)
## Season             1 1085.524288 1085.524288 6.19968509 0.04716755
## Cultivar           1  349.143354  349.143354 1.99404000 0.20762177
## Block              2  100.698804   50.349402 0.28755730 0.75987782
## Season × Cultivar  1    7.895831    7.895831 0.04509495 0.83886001
## Residuals          6 1050.560735  175.093456                      
## 
## -----------------------------------------------------------------
## No significant interaction
## -----------------------------------------------------------------
## 
## -----------------------------------------------------------------
## Season
## -----------------------------------------------------------------
## Multiple Comparison Test: Tukey HSD 
##        resp groups
## E1 120.0346      a
## E2 101.0124      b
## NULL
## 
## -----------------------------------------------------------------
## Isolated factors 2 not significant
## -----------------------------------------------------------------
##       Mean
## D 105.1295
## M 115.9175
## 
## -----------------------------------------------------------------
## Normality of errors
## -----------------------------------------------------------------
##                          Method Statistic   p.value
##  Shapiro-Wilk normality test(W) 0.9479183 0.6067783
## 
## -----------------------------------------------------------------
## Homogeneity of Variances
## -----------------------------------------------------------------
##                               Method Statistic   p.value
##  Bartlett test(Bartlett's K-squared)  2.672334 0.4449493
## 
## -----------------------------------------------------------------
## Independence from errors
## -----------------------------------------------------------------
##                  Method Statistic   p.value
##  Durbin-Watson test(DW)  2.989822 0.8376624
## 
## -----------------------------------------------------------------
## Additional Information
## -----------------------------------------------------------------
## 
## CV (%) =  19.16
## Mean =  226.8634
## Median =  212.1108
## Possible outliers =  No discrepant point
## 
## -----------------------------------------------------------------
## Analysis of Variance
## -----------------------------------------------------------------
##                   Df     Sum Sq    Mean.Sq     F value       Pr(F)
## Season             1 35212.3268 35212.3268 18.64514586 0.004993589
## Cultivar           1   919.0813   919.0813  0.48665924 0.511534775
## Block              2  2108.3689  1054.1845  0.55819722 0.599340623
## Season × Cultivar  1   165.0161   165.0161  0.08737704 0.777501528
## Residuals          6 11331.3118  1888.5520                        
## 
## -----------------------------------------------------------------
## No significant interaction
## -----------------------------------------------------------------
## 
## -----------------------------------------------------------------
## Season
## -----------------------------------------------------------------
## Multiple Comparison Test: Tukey HSD 
##        resp groups
## E1 281.0331      a
## E2 172.6936      b
## NULL
## 
## -----------------------------------------------------------------
## Isolated factors 2 not significant
## -----------------------------------------------------------------
##       Mean
## D 218.1118
## M 235.6150
## 
## -----------------------------------------------------------------
## Normality of errors
## -----------------------------------------------------------------
##                          Method Statistic   p.value
##  Shapiro-Wilk normality test(W) 0.9793568 0.9809638
## 
## -----------------------------------------------------------------
## Homogeneity of Variances
## -----------------------------------------------------------------
##                               Method Statistic  p.value
##  Bartlett test(Bartlett's K-squared)  5.808614 0.121302
## 
## -----------------------------------------------------------------
## Independence from errors
## -----------------------------------------------------------------
##                  Method Statistic   p.value
##  Durbin-Watson test(DW)  3.036684 0.8610195
## 
## -----------------------------------------------------------------
## Additional Information
## -----------------------------------------------------------------
## 
## CV (%) =  18.88
## Mean =  0.0076
## Median =  0.0073
## Possible outliers =  No discrepant point
## 
## -----------------------------------------------------------------
## Analysis of Variance
## -----------------------------------------------------------------
##                   Df       Sum Sq      Mean.Sq   F value      Pr(F)
## Season             1 2.008052e-05 2.008052e-05 9.6291730 0.02103336
## Cultivar           1 5.392787e-07 5.392787e-07 0.2585993 0.62923574
## Block              2 7.826049e-06 3.913025e-06 1.8764053 0.23284365
## Season × Cultivar  1 3.685094e-07 3.685094e-07 0.1767106 0.68885473
## Residuals          6 1.251230e-05 2.085384e-06                     
## 
## -----------------------------------------------------------------
## No significant interaction
## -----------------------------------------------------------------
## 
## -----------------------------------------------------------------
## Season
## -----------------------------------------------------------------
## Multiple Comparison Test: Tukey HSD 
##           resp groups
## E2 0.008940489      a
## E1 0.006353308      b
## NULL
## 
## -----------------------------------------------------------------
## Isolated factors 2 not significant
## -----------------------------------------------------------------
##          Mean
## D 0.007858889
## M 0.007434908
## 
## -----------------------------------------------------------------
## Normality of errors
## -----------------------------------------------------------------
##                          Method Statistic   p.value
##  Shapiro-Wilk normality test(W) 0.9479183 0.6067783
## 
## -----------------------------------------------------------------
## Homogeneity of Variances
## -----------------------------------------------------------------
##                               Method Statistic   p.value
##  Bartlett test(Bartlett's K-squared)  2.672334 0.4449493
## 
## -----------------------------------------------------------------
## Independence from errors
## -----------------------------------------------------------------
##                  Method Statistic   p.value
##  Durbin-Watson test(DW)  2.989822 0.8376624
## 
## -----------------------------------------------------------------
## Additional Information
## -----------------------------------------------------------------
## 
## CV (%) =  19.16
## Mean =  226.8634
## Median =  212.1108
## Possible outliers =  No discrepant point
## 
## -----------------------------------------------------------------
## Analysis of Variance
## -----------------------------------------------------------------
##                   Df     Sum Sq    Mean.Sq     F value       Pr(F)
## Season             1 35212.3268 35212.3268 18.64514586 0.004993589
## Cultivar           1   919.0813   919.0813  0.48665924 0.511534775
## Block              2  2108.3689  1054.1845  0.55819722 0.599340623
## Season × Cultivar  1   165.0161   165.0161  0.08737704 0.777501528
## Residuals          6 11331.3118  1888.5520                        
## 
## -----------------------------------------------------------------
## No significant interaction
## -----------------------------------------------------------------
## 
## -----------------------------------------------------------------
## Season
## -----------------------------------------------------------------
## Multiple Comparison Test: Tukey HSD 
##        resp groups
## E1 281.0331      a
## E2 172.6936      b
## NULL
## 
## -----------------------------------------------------------------
## Isolated factors 2 not significant
## -----------------------------------------------------------------
##       Mean
## D 218.1118
## M 235.6150
## 
## -----------------------------------------------------------------
## Normality of errors
## -----------------------------------------------------------------
##                          Method Statistic   p.value
##  Shapiro-Wilk normality test(W) 0.9873152 0.9987483
## 
## -----------------------------------------------------------------
## Homogeneity of Variances
## -----------------------------------------------------------------
##                               Method Statistic   p.value
##  Bartlett test(Bartlett's K-squared)  5.143787 0.1615661
## 
## -----------------------------------------------------------------
## Independence from errors
## -----------------------------------------------------------------
##                  Method Statistic   p.value
##  Durbin-Watson test(DW)  3.006183 0.8459505
## 
## -----------------------------------------------------------------
## Additional Information
## -----------------------------------------------------------------
## 
## CV (%) =  10.57
## Mean =  0.2057
## Median =  0.1988
## Possible outliers =  No discrepant point
## 
## -----------------------------------------------------------------
## Analysis of Variance
## -----------------------------------------------------------------
##                   Df       Sum Sq      Mean.Sq   F value      Pr(F)
## Season             1 0.0038652014 0.0038652014 8.1694080 0.02886902
## Cultivar           1 0.0001318426 0.0001318426 0.2786597 0.61652476
## Block              2 0.0034542905 0.0017271453 3.6504577 0.09179291
## Season × Cultivar  1 0.0009344047 0.0009344047 1.9749380 0.20953082
## Residuals          6 0.0028387869 0.0004731312                     
## 
## -----------------------------------------------------------------
## No significant interaction
## -----------------------------------------------------------------
## 
## -----------------------------------------------------------------
## Season
## -----------------------------------------------------------------
## Multiple Comparison Test: Tukey HSD 
##         resp groups
## E2 0.2236683      a
## E1 0.1877740      b
## NULL
## 
## -----------------------------------------------------------------
## Isolated factors 2 not significant
## -----------------------------------------------------------------
##        Mean
## D 0.2024065
## M 0.2090358
## 
## -----------------------------------------------------------------
## Normality of errors
## -----------------------------------------------------------------
##                          Method Statistic   p.value
##  Shapiro-Wilk normality test(W) 0.9614341 0.8040929
## 
## -----------------------------------------------------------------
## Homogeneity of Variances
## -----------------------------------------------------------------
##                               Method Statistic   p.value
##  Bartlett test(Bartlett's K-squared)  1.409281 0.7033606
## 
## -----------------------------------------------------------------
## Independence from errors
## -----------------------------------------------------------------
##                  Method Statistic   p.value
##  Durbin-Watson test(DW)  2.519047 0.5616332
## 
## -----------------------------------------------------------------
## Additional Information
## -----------------------------------------------------------------
## 
## CV (%) =  27.82
## Mean =  44.8547
## Median =  39.6692
## Possible outliers =  No discrepant point
## 
## -----------------------------------------------------------------
## Analysis of Variance
## -----------------------------------------------------------------
##                   Df      Sum Sq     Mean.Sq     F value      Pr(F)
## Season             1 7024.809063 7024.809063 45.09861723 0.00052990
## Cultivar           1    2.785608    2.785608  0.01788334 0.89798988
## Block              2 2264.116490 1132.058245  7.26770806 0.02494269
## Season × Cultivar  1   31.290353   31.290353  0.20088114 0.66974070
## Residuals          6  934.593053  155.765509                       
## 
## -----------------------------------------------------------------
## No significant interaction
## -----------------------------------------------------------------
## 
## -----------------------------------------------------------------
## Season
## -----------------------------------------------------------------
## Multiple Comparison Test: Tukey HSD 
##        resp groups
## E1 69.04974      a
## E2 20.65962      b
## NULL
## 
## -----------------------------------------------------------------
## Isolated factors 2 not significant
## -----------------------------------------------------------------
##       Mean
## D 44.37288
## M 45.33648
## 
## -----------------------------------------------------------------
## Normality of errors
## -----------------------------------------------------------------
##                          Method Statistic   p.value
##  Shapiro-Wilk normality test(W) 0.9398903 0.4966522
## 
## -----------------------------------------------------------------
## Homogeneity of Variances
## -----------------------------------------------------------------
##                               Method Statistic   p.value
##  Bartlett test(Bartlett's K-squared)  2.810452 0.4217822
## 
## -----------------------------------------------------------------
## Independence from errors
## -----------------------------------------------------------------
##                  Method Statistic   p.value
##  Durbin-Watson test(DW)  3.026438 0.8560152
## 
## -----------------------------------------------------------------
## Additional Information
## -----------------------------------------------------------------
## 
## CV (%) =  19.15
## Mean =  408.8721
## Median =  402.7028
## Possible outliers =  No discrepant point
## 
## -----------------------------------------------------------------
## Analysis of Variance
## -----------------------------------------------------------------
##                   Df     Sum Sq    Mean.Sq    F value       Pr(F)
## Season             1 84963.4022 84963.4022 13.8621127 0.009815345
## Cultivar           1  3476.7169  3476.7169  0.5672400 0.479856076
## Block              2  5604.8426  2802.4213  0.4572260 0.653401965
## Season × Cultivar  1   978.7822   978.7822  0.1596922 0.703277019
## Residuals          6 36775.0878  6129.1813                       
## 
## -----------------------------------------------------------------
## No significant interaction
## -----------------------------------------------------------------
## 
## -----------------------------------------------------------------
## Season
## -----------------------------------------------------------------
## Multiple Comparison Test: Tukey HSD 
##        resp groups
## E1 493.0165      a
## E2 324.7277      b
## NULL
## 
## -----------------------------------------------------------------
## Isolated factors 2 not significant
## -----------------------------------------------------------------
##       Mean
## D 391.8507
## M 425.8934

# 2. Extract and combine tidy ANOVA tables for the summary
anova_nf <- 
  mod_anova_nf |> 
  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_nf |> 
  filter(term != "Residuals", Parameter != "xpi") |> 
  gt() |> 
  tab_header(
    title = "Supplementary Table S9",
    subtitle = "ANOVA Summary: Model parameters and critical points for Leaf Number"
  ) |> 
  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 S9
ANOVA Summary: Model parameters and critical points for Leaf Number
Parameter term GL SQ QM Fcal p.value
b1 Season 1.00000 1,085.52429 1,085.52429 6.19969 0.04717
b1 Cultivar 1.00000 349.14335 349.14335 1.99404 0.20762
b1 bloco 2.00000 100.69880 50.34940 0.28756 0.75988
b1 Season:Cultivar 1.00000 7.89583 7.89583 0.04509 0.83886
b2 Season 1.00000 35,212.32683 35,212.32683 18.64515 0.00499
b2 Cultivar 1.00000 919.08126 919.08126 0.48666 0.51153
b2 bloco 2.00000 2,108.36892 1,054.18446 0.55820 0.59934
b2 Season:Cultivar 1.00000 165.01608 165.01608 0.08738 0.77750
b3 Season 1.00000 0.00002 0.00002 9.62917 0.02103
b3 Cultivar 1.00000 0.00000 0.00000 0.25860 0.62924
b3 bloco 2.00000 0.00001 0.00000 1.87641 0.23284
b3 Season:Cultivar 1.00000 0.00000 0.00000 0.17671 0.68885
ypi Season 1.00000 0.00387 0.00387 8.16941 0.02887
ypi Cultivar 1.00000 0.00013 0.00013 0.27866 0.61652
ypi bloco 2.00000 0.00345 0.00173 3.65046 0.09179
ypi Season:Cultivar 1.00000 0.00093 0.00093 1.97494 0.20953
xmap Season 1.00000 7,024.80906 7,024.80906 45.09862 0.00053
xmap Cultivar 1.00000 2.78561 2.78561 0.01788 0.89799
xmap bloco 2.00000 2,264.11649 1,132.05824 7.26771 0.02494
xmap Season:Cultivar 1.00000 31.29035 31.29035 0.20088 0.66974
xmdp Season 1.00000 84,963.40218 84,963.40218 13.86211 0.00982
xmdp Cultivar 1.00000 3,476.71695 3,476.71695 0.56724 0.47986
xmdp bloco 2.00000 5,604.84264 2,802.42132 0.45723 0.65340
xmdp Season:Cultivar 1.00000 978.78218 978.78218 0.15969 0.70328
Show code 

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

Fitted logistic curves

Show code 
formula_xy <- y ~ b1 / (1 + exp((b3 * (b2 - x))))

ggplot(
  df_model |> filter(!data %in% c("04/11")),
  aes(gda, n_mean, color = cultivar)
) +
  geom_smooth(
    method      = "nls",
    method.args = list(
      formula = formula_xy,
      start   = c(b1 = 100, 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     = "Number of leaves per plant",
    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 = n_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 cumulative leaf number per plant 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_nf.jpg", width = 6, height = 9)

Derivative plots

Show code 

# --- First derivative plot ---
plot_pi_nf <-
  ggplot() +
  stat_function(fun = dy, aes(color = "D", linetype = "E1"), n = 500, linewidth = 1,
                xlim = c(0, 1000),
                args = c(b1 = parameters_nf_mean[[1, 3]],
                         b2 = parameters_nf_mean[[1, 4]],
                         b3 = parameters_nf_mean[[1, 5]])) +
  stat_function(fun = dy, aes(color = "M", linetype = "E1"), linewidth = 1,
                xlim = c(0, 1000),
                args = c(b1 = parameters_nf_mean[[2, 3]],
                         b2 = parameters_nf_mean[[2, 4]],
                         b3 = parameters_nf_mean[[2, 5]])) +
  stat_function(fun = dy, aes(color = "D", linetype = "E2"), linewidth = 1,
                xlim = c(0, 1000),
                args = c(b1 = parameters_nf_mean[[3, 3]],
                         b2 = parameters_nf_mean[[3, 4]],
                         b3 = parameters_nf_mean[[3, 5]])) +
  stat_function(fun = dy, aes(color = "M", linetype = "E2"), linewidth = 1,
                xlim = c(0, 1000),
                args = c(b1 = parameters_nf_mean[[4, 3]],
                         b2 = parameters_nf_mean[[4, 4]],
                         b3 = parameters_nf_mean[[4, 5]])) +
  # Mark the inflection point for each treatment
  geom_point(aes(xpi, ypi, shape = epoca, color = cultivar),
             data = parameters_nf_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(Number ~ of ~ leaves ~ plant^{-1} ~ degree ~ 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_nf <-
  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_nf_mean[[1, 3]],
                         b2 = parameters_nf_mean[[1, 4]],
                         b3 = parameters_nf_mean[[1, 5]])) +
  stat_function(fun = d2y, aes(color = "M", linetype = "E1"), n = 500, linewidth = 1,
                xlim = c(0, 1000),
                args = c(b1 = parameters_nf_mean[[2, 3]],
                         b2 = parameters_nf_mean[[2, 4]],
                         b3 = parameters_nf_mean[[2, 5]])) +
  stat_function(fun = d2y, aes(color = "D", linetype = "E2"), n = 500, linewidth = 1,
                xlim = c(0, 1000),
                args = c(b1 = parameters_nf_mean[[3, 3]],
                         b2 = parameters_nf_mean[[3, 4]],
                         b3 = parameters_nf_mean[[3, 5]])) +
  stat_function(fun = d2y, aes(color = "M", linetype = "E2"), n = 500, linewidth = 1,
                xlim = c(0, 1000),
                args = c(b1 = parameters_nf_mean[[4, 3]],
                         b2 = parameters_nf_mean[[4, 4]],
                         b3 = parameters_nf_mean[[4, 5]])) +
  # MAP (triangle up) and MDP (triangle down) markers
  geom_point(aes(xmap, ymap, fill = cultivar), data = parameters_nf_mean,
             size = 3, shape = 24, show.legend = FALSE) +
  geom_point(aes(xmdp, ymdp, fill = cultivar), data = parameters_nf_mean,
             size = 3, shape = 25, show.legend = FALSE) +
  theme_bw(base_size = 20) +
  labs(
    x        = "Accumulated Degree Days",
    y        = expression(Number ~ of ~ leaves ~ plant^{-1} ~ degree ~ 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_nf / plot_sd_nf +
  plot_layout(guides = "collect") &
  theme(legend.position = "bottom")

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

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

Correlation between model parameters

Show code 
# Correlation between parameters and critical points

  parameters_nf |> 
  select(all_of(vars_nf)) |> 
  metan::corr_coef() |> 
  plot()