Nonlinear logistic regression of leaf area per plant as a function of accumulated growing degree-days (GDD), with ANOVA on model parameters and derivative-based critical-point analysis.

Overview

A three-parameter logistic model is fitted separately to each combination of season × cultivar × block. The estimated parameters are then subjected to a two-factor ANOVA in randomised complete blocks (FAT2DBC) to test the effects of sowing date and cultivar.

The logistic function is:

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

Parameter Interpretation
β₁ Upper asymptote (maximum leaf area, cm² plant⁻¹)
β₂ Inflection point (GDD at maximum growth rate)
β₃ Steepness coefficient (growth rate)

Required packages

Show code 
library(rio)        # Data import
library(tidyverse)  # Data wrangling and visualisation
library(metan)      # mean_by(), doo(), corr_plot()
library(broom)      # tidy() – extract NLS coefficients as a data frame
library(AgroR)      # FAT2DBC() – two-factor ANOVA in RCBD
library(patchwork)  # Combine ggplot objects
library(hydroGOF)   # gof() – goodness-of-fit statistics (R², RMSE, …)
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 area

Show code 
# Define the logistic formula for NLS fitting
formula_af <- af_planta_mean ~ b1 / (1 + exp((b3 * (b2 - gda))))

# Initial parameter guesses (required by nls())
#   b1 = 100  → expected maximum leaf area (~100 cm² plant⁻¹)
#   b2 = 400  → inflection point around 400 GDD
#   b3 = 0.01 → moderate steepness
start_af <- list(b1 = 100, b2 = 400, b3 = 0.01)

# Fit one model per season × cultivar × block combination.
# The date "04/11" is excluded because it falls outside the vegetative growth
# window (post-flowering measurement, distorts the sigmoid shape).
# metan::doo() is a tidy wrapper around dplyr::group_map() for NLS.
mod_af <-
  df_model |>
  filter(data != "04/11") |>
  group_by(epoca, cultivar, bloco) |>
  doo(~nls(formula_af, data = ., start = start_af))

Extract estimated parameters

Show code 
# broom::tidy() converts each NLS object into a data frame with columns:
#   term, estimate, std.error, statistic, p.value
# pivot_wider() reshapes so each parameter (b1, b2, b3) is its own column.
parameters_af <-
  mod_af |>
  mutate(data = map(data, ~.x |> tidy())) |>
  unnest(data) |>
  select(epoca, cultivar, bloco, term, estimate) |>
  pivot_wider(names_from  = term,
              values_from = estimate)

Goodness of fit

Show code 
# Helper function: extract R² and AIC from a single NLS object
get_gof <- function(model) {
  aic  <- AIC(model)
  bic  <- BIC(model)
  fit  <- model$m$fitted()   # fitted values
  res  <- model$m$resid()    # residuals
  obs  <- fit + res           # observed = fitted + residual
  
  # Compute GOF metrics using hydroGOF::gof
  g <- hydroGOF::gof(sim = fit, obs = obs, digits = 4)
  
  # Extract specific metrics
  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 the helper to every model in the nested data frame
gof_af <-
  mod_af |>
  mutate(map_dfr(.x = data, .f = ~get_gof(.))) |>
  select(-data) |>
  mutate(variable = "leaf_area")

# Display goodness of fit table with gt
gof_af |> 
  gt() |> 
  tab_header(
    title = "Supplementary Table S1",
    subtitle = "Goodness-of-fit metrics for Leaf Area logistic models (per block)"
  ) |> 
  fmt_number(columns = where(is.numeric), decimals = 4)
Supplementary Table S1
Goodness-of-fit metrics for Leaf Area logistic models (per block)
epoca bloco cultivar aic bic r2 rmse mae d variable
E1 B1 D 92.1544 94.0940 0.9512 8.0647 5.8542 0.9875 leaf_area
E1 B2 D 96.2183 98.1579 0.8796 9.5527 7.6902 0.9681 leaf_area
E1 B3 D 91.6287 93.5684 0.9259 7.8900 6.3099 0.9810 leaf_area
E1 B1 M 98.0237 99.9633 0.9402 10.2990 8.1658 0.9849 leaf_area
E1 B2 M 99.4571 101.3968 0.9518 10.9329 8.6882 0.9883 leaf_area
E1 B3 M 100.9898 102.9294 0.8989 11.6538 9.1617 0.9732 leaf_area
E2 B1 D 58.6348 59.8451 0.9750 3.0429 2.7395 0.9937 leaf_area
E2 B2 D 70.9637 72.1740 0.9185 5.6364 4.5106 0.9790 leaf_area
E2 B3 D 65.4323 66.6426 0.9478 4.2745 3.2005 0.9864 leaf_area
E2 B1 M 79.8645 81.0749 0.8642 8.7959 7.1313 0.9625 leaf_area
E2 B2 M 78.0279 79.2382 0.9025 8.0242 6.2223 0.9732 leaf_area
E2 B3 M 74.8398 76.0502 0.9147 6.8418 5.1728 0.9780 leaf_area
Show code 

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

Fitted logistic curves

Show code 
# Logistic formula used inside geom_smooth (x/y notation for ggplot internals)
formula_xy <- y ~ b1 / (1 + exp((b3 * (b2 - x))))

ggplot(
  df_model |> filter(!data %in% c("04/11")),
  aes(gda, af_planta_mean, color = cultivar)
) +
  # Fitted NLS curves via geom_smooth (one per 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)
  ) +
  # One panel per sowing season
  facet_wrap(~epoca, ncol = 1) +
  # Treatment means overlaid as points
  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     = expression(Leaf ~ area ~ (cm^2 ~ plant^{-1})),
    color = "Cultivar"
  ) +
  scale_color_manual(values = c("gold", "brown")) +
  # Add the excluded 04/11 point as a transparent reference
  stat_summary(
    data          = df_model |> filter(data == "04/11"),
    aes(x = gda, y = af_planta_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 growth curves fitted to leaf area per plant (cm² plant⁻¹) as a function of accumulated GDD, faceted by sowing season. Points = treatment means. The semi-transparent point at the far right (04/11) was excluded from fitting.
Show code 

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

Derivative analysis

First derivative – growth rate

The first derivative of the logistic function with respect to GDD gives the instantaneous growth rate of leaf area:

\[ \frac{dy}{dx} = \frac{\beta_1 \cdot \beta_3 \cdot e^{\beta_3(\beta_2 - x)}} {\left[1 + e^{\beta_3(\beta_2 - x)}\right]^2} \]

This reaches its maximum at the inflection point x = β₂.

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

# Numerical first-derivative function
dy <- function(x, b1, b2, b3) {
  b1 * (exp((b3 * (b2 - x))) * b3) / (1 + exp((b3 * (b2 - x))))^2
}

# Average parameters over blocks (mean_by collapses replications)
# and compute the inflection-point coordinates (xpi, ypi)
parameters_af <-
  parameters_af |>
  mutate(
    xpi = b2,                        # x at inflection = β₂
    ypi = dy(xpi, b1, b2, b3)        # maximum growth rate value
  ) |>
  as.data.frame() |>
  mutate(variable = "leaf_area")

Second derivative – acceleration / deceleration

The second derivative identifies the critical acceleration point (MAP: maximum acceleration point) and the maximum deceleration point (MDP):

\[ x_{MAP} = \beta_2 + \frac{\ln(2 + \sqrt{3})}{\beta_3}, \qquad x_{MDP} = \beta_2 - \frac{\ln(2 + \sqrt{3})}{\beta_3} \]

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

# Numerical second-derivative function
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)
}

# Compute MAP and MDP x-coordinates and their second-derivative values
parameters_af <-
  parameters_af |>
  mutate(
    xmap = b2 - (log(2 + sqrt(3)) / b3),   # GDD at maximum acceleration point
    xmdp = b2 + (log(2 + sqrt(3)) / b3),   # GDD at maximum deceleration point
    ymap = d2y(xmap, b1, b2, b3),
    ymdp = d2y(xmdp, b1, b2, b3)
  )

parameters_af_mean <- 
  parameters_af |> 
  metan::mean_by(epoca, cultivar) 

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

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

means_summary_af <- 
  bind_rows(means_epoca_af, means_cultivar_af) |> 
  relocate(Factor, Level)


# Display means summary using gt
means_summary_af |> 
  gt(groupname_col = "Factor") |> 
  tab_header(
    title = "Supplementary Table S2", 
    subtitle = "Main Effect Means: Model parameters and critical points for Leaf Area"
  ) |> 
  fmt_number(columns = where(is.numeric), decimals = 5)
Supplementary Table S2
Main Effect Means: Model parameters and critical points for Leaf Area
Level b1 b2 b3 xpi ypi xmap xmdp ymap ymdp
Sowing Season
E1 96.86681 467.21673 0.00899 467.21673 0.21646 318.26514 616.16831 0.00076 −0.00076
E2 65.60805 327.41369 0.00866 327.41369 0.13720 161.08138 493.74601 0.00048 −0.00048
Cultivar
D 69.47417 379.84506 0.00893 379.84506 0.15161 224.16396 535.52616 0.00054 −0.00054
M 93.00068 414.78536 0.00872 414.78536 0.20204 255.18256 574.38815 0.00070 −0.00070
Show code 

saveRDS(means_summary_af, "../results/means_summary_af.rds")
Show code 

# --- First derivative plot ---
plot_pi <-
  ggplot() +
  # Each combination: cultivar (colour) × season (linetype)
  stat_function(fun = dy, aes(color = "D", linetype = "E1"), n = 500, linewidth = 1,
                xlim = c(0, 1000),
                args = c(b1 = parameters_af_mean[[1, 3]],
                         b2 = parameters_af_mean[[1, 4]],
                         b3 = parameters_af_mean[[1, 5]])) +
  stat_function(fun = dy, aes(color = "M", linetype = "E1"), linewidth = 1,
                xlim = c(0, 1000),
                args = c(b1 = parameters_af_mean[[2, 3]],
                         b2 = parameters_af_mean[[2, 4]],
                         b3 = parameters_af_mean[[2, 5]])) +
  stat_function(fun = dy, aes(color = "D", linetype = "E2"), linewidth = 1,
                xlim = c(0, 1000),
                args = c(b1 = parameters_af_mean[[3, 3]],
                         b2 = parameters_af_mean[[3, 4]],
                         b3 = parameters_af_mean[[3, 5]])) +
  stat_function(fun = dy, aes(color = "M", linetype = "E2"), linewidth = 1,
                xlim = c(0, 1000),
                args = c(b1 = parameters_af_mean[[4, 3]],
                         b2 = parameters_af_mean[[4, 4]],
                         b3 = parameters_af_mean[[4, 5]])) +
  # Mark the inflection point for each treatment
  geom_point(aes(xpi, ypi, shape = epoca, color = cultivar),
             data = parameters_af_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(Leaf ~ area ~ (cm^2 ~ 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 <-
  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_af_mean[[1, 3]],
                         b2 = parameters_af_mean[[1, 4]],
                         b3 = parameters_af_mean[[1, 5]])) +
  stat_function(fun = d2y, aes(color = "M", linetype = "E1"), n = 500, linewidth = 1,
                xlim = c(0, 1000),
                args = c(b1 = parameters_af_mean[[2, 3]],
                         b2 = parameters_af_mean[[2, 4]],
                         b3 = parameters_af_mean[[2, 5]])) +
  stat_function(fun = d2y, aes(color = "D", linetype = "E2"), n = 500, linewidth = 1,
                xlim = c(0, 1000),
                args = c(b1 = parameters_af_mean[[3, 3]],
                         b2 = parameters_af_mean[[3, 4]],
                         b3 = parameters_af_mean[[3, 5]])) +
  stat_function(fun = d2y, aes(color = "M", linetype = "E2"), n = 500, linewidth = 1,
                xlim = c(0, 1000),
                args = c(b1 = parameters_af_mean[[4, 3]],
                         b2 = parameters_af_mean[[4, 4]],
                         b3 = parameters_af_mean[[4, 5]])) +
  # MAP (triangle up) and MDP (triangle down) markers
  geom_point(aes(xmap, ymap, fill = cultivar), data = parameters_af_mean,
             size = 3, shape = 24, show.legend = FALSE) +
  geom_point(aes(xmdp, ymdp, fill = cultivar), data = parameters_af_mean,
             size = 3, shape = 25, show.legend = FALSE) +
  theme_bw(base_size = 20) +
  labs(
    x       = "Accumulated Degree Days",
    y       = expression(cm^2 ~ 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"))

# Combine both plots with a shared legend
plot_pi / plot_sd +
  plot_layout(guides = "collect") &
  theme(legend.position = "bottom")

First (left) and second (right) derivatives of the fitted logistic curves for leaf area. Points mark the inflection, acceleration, and deceleration thresholds for each cultivar × season combination.
Show code 

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

ANOVA on model parameters and critical points

Show code 
# List of variables to analyze
vars_af <- 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_af <- 
  vars_af |> 
  set_names() |> 
  map(~with(parameters_af, FAT2DBC(epoca, cultivar, bloco, get(.x),
                                   names.fat = c("Season", "Cultivar"))))
## 
## -----------------------------------------------------------------
## Normality of errors
## -----------------------------------------------------------------
##                          Method Statistic   p.value
##  Shapiro-Wilk normality test(W) 0.9099109 0.2127872
## 
## -----------------------------------------------------------------
## Homogeneity of Variances
## -----------------------------------------------------------------
##                               Method Statistic   p.value
##  Bartlett test(Bartlett's K-squared)  2.288671 0.5146948
## 
## -----------------------------------------------------------------
## Independence from errors
## -----------------------------------------------------------------
##                  Method Statistic   p.value
##  Durbin-Watson test(DW)  2.611349 0.6193389
## 
## -----------------------------------------------------------------
## Additional Information
## -----------------------------------------------------------------
## 
## CV (%) =  16.53
## Mean =  81.2374
## Median =  70.2282
## Possible outliers =  No discrepant point
## 
## -----------------------------------------------------------------
## Analysis of Variance
## -----------------------------------------------------------------
##                   Df    Sum Sq   Mean.Sq   F value       Pr(F)
## Season             1 2931.3288 2931.3288 16.263905 0.006858053
## Cultivar           1 1660.4898 1660.4898  9.212903 0.022941726
## Block              2  578.7362  289.3681  1.605502 0.276396627
## Season × Cultivar  1  241.7477  241.7477  1.341290 0.290818627
## Residuals          6 1081.4114  180.2352                      
## 
## -----------------------------------------------------------------
## No significant interaction
## -----------------------------------------------------------------
## 
## -----------------------------------------------------------------
## Season
## -----------------------------------------------------------------
## Multiple Comparison Test: Tukey HSD 
##        resp groups
## E1 96.86681      a
## E2 65.60805      b
## NULL
## 
## -----------------------------------------------------------------
## Cultivar
## -----------------------------------------------------------------
## Multiple Comparison Test: Tukey HSD 
##       resp groups
## M 93.00068      a
## D 69.47417      b
## NULL
## 
## -----------------------------------------------------------------
## Normality of errors
## -----------------------------------------------------------------
##                          Method Statistic   p.value
##  Shapiro-Wilk normality test(W) 0.8410314 0.0284992
## 
## -----------------------------------------------------------------
## Homogeneity of Variances
## -----------------------------------------------------------------
##                               Method Statistic   p.value
##  Bartlett test(Bartlett's K-squared)  1.168717 0.7605167
## 
## -----------------------------------------------------------------
## Independence from errors
## -----------------------------------------------------------------
##                  Method Statistic   p.value
##  Durbin-Watson test(DW)  1.874465 0.1755968
## 
## -----------------------------------------------------------------
## Additional Information
## -----------------------------------------------------------------
## 
## CV (%) =  14.38
## Mean =  397.3152
## Median =  390.9044
## Possible outliers =  No discrepant point
## 
## -----------------------------------------------------------------
## Analysis of Variance
## -----------------------------------------------------------------
##                   Df    Sum Sq   Mean.Sq   F value       Pr(F)
## Season             1 58634.667 58634.667 17.955465 0.005455459
## Cultivar           1  3662.472  3662.472  1.121545 0.330346269
## Block              2 20225.434 10112.717  3.096778 0.119141390
## Season × Cultivar  1  4259.832  4259.832  1.304472 0.296919138
## Residuals          6 19593.366  3265.561
## -----------------------------------------------------------------
## No significant interaction
## -----------------------------------------------------------------
## 
## -----------------------------------------------------------------
## Season
## -----------------------------------------------------------------
## Multiple Comparison Test: Tukey HSD 
##        resp groups
## E1 467.2167      a
## E2 327.4137      b
## NULL
## 
## -----------------------------------------------------------------
## Isolated factors 2 not significant
## -----------------------------------------------------------------
##       Mean
## D 379.8451
## M 414.7854
## 
## -----------------------------------------------------------------
## Normality of errors
## -----------------------------------------------------------------
##                          Method Statistic   p.value
##  Shapiro-Wilk normality test(W) 0.9654659 0.8580232
## 
## -----------------------------------------------------------------
## Homogeneity of Variances
## -----------------------------------------------------------------
##                               Method Statistic   p.value
##  Bartlett test(Bartlett's K-squared)  1.427575 0.6990839
## 
## -----------------------------------------------------------------
## Independence from errors
## -----------------------------------------------------------------
##                  Method Statistic   p.value
##  Durbin-Watson test(DW)  1.782741 0.1311485
## 
## -----------------------------------------------------------------
## Additional Information
## -----------------------------------------------------------------
## 
## CV (%) =  20.61
## Mean =  0.0088
## Median =  0.0092
## Possible outliers =  No discrepant point
## 
## -----------------------------------------------------------------
## Analysis of Variance
## -----------------------------------------------------------------
##                   Df       Sum Sq      Mean.Sq     F value      Pr(F)
## Season             1 3.292602e-07 3.292602e-07 0.099510351 0.76310173
## Cultivar           1 1.342966e-07 1.342966e-07 0.040587681 0.84699241
## Block              2 2.570451e-05 1.285226e-05 3.884261213 0.08275452
## Season × Cultivar  1 1.214815e-08 1.214815e-08 0.003671465 0.95365153
## Residuals          6 1.985282e-05 3.308803e-06                       
## 
## -----------------------------------------------------------------
## No significant interaction
## -----------------------------------------------------------------
## 
## -----------------------------------------------------------------
## Isolated factors not significant
## -----------------------------------------------------------------
##              D           M
## E1 0.009130135 0.008854921
## E2 0.008735210 0.008587266
## 
## -----------------------------------------------------------------
## Normality of errors
## -----------------------------------------------------------------
##                          Method Statistic   p.value
##  Shapiro-Wilk normality test(W) 0.8410314 0.0284992
## 
## -----------------------------------------------------------------
## Homogeneity of Variances
## -----------------------------------------------------------------
##                               Method Statistic   p.value
##  Bartlett test(Bartlett's K-squared)  1.168717 0.7605167
## 
## -----------------------------------------------------------------
## Independence from errors
## -----------------------------------------------------------------
##                  Method Statistic   p.value
##  Durbin-Watson test(DW)  1.874465 0.1755968
## 
## -----------------------------------------------------------------
## Additional Information
## -----------------------------------------------------------------
## 
## CV (%) =  14.38
## Mean =  397.3152
## Median =  390.9044
## Possible outliers =  No discrepant point
## 
## -----------------------------------------------------------------
## Analysis of Variance
## -----------------------------------------------------------------
##                   Df    Sum Sq   Mean.Sq   F value       Pr(F)
## Season             1 58634.667 58634.667 17.955465 0.005455459
## Cultivar           1  3662.472  3662.472  1.121545 0.330346269
## Block              2 20225.434 10112.717  3.096778 0.119141390
## Season × Cultivar  1  4259.832  4259.832  1.304472 0.296919138
## Residuals          6 19593.366  3265.561
## -----------------------------------------------------------------
## No significant interaction
## -----------------------------------------------------------------
## 
## -----------------------------------------------------------------
## Season
## -----------------------------------------------------------------
## Multiple Comparison Test: Tukey HSD 
##        resp groups
## E1 467.2167      a
## E2 327.4137      b
## NULL
## 
## -----------------------------------------------------------------
## Isolated factors 2 not significant
## -----------------------------------------------------------------
##       Mean
## D 379.8451
## M 414.7854
## 
## -----------------------------------------------------------------
## Normality of errors
## -----------------------------------------------------------------
##                          Method Statistic   p.value
##  Shapiro-Wilk normality test(W) 0.9808685 0.9868669
## 
## -----------------------------------------------------------------
## Homogeneity of Variances
## -----------------------------------------------------------------
##                               Method Statistic   p.value
##  Bartlett test(Bartlett's K-squared)  4.499618 0.2123244
## 
## -----------------------------------------------------------------
## Independence from errors
## -----------------------------------------------------------------
##                  Method Statistic   p.value
##  Durbin-Watson test(DW)  2.727901 0.6905783
## 
## -----------------------------------------------------------------
## Additional Information
## -----------------------------------------------------------------
## 
## CV (%) =  21.27
## Mean =  0.1768
## Median =  0.1701
## Possible outliers =  No discrepant point
## 
## -----------------------------------------------------------------
## Analysis of Variance
## -----------------------------------------------------------------
##                   Df      Sum Sq     Mean.Sq   F value      Pr(F)
## Season             1 0.018847429 0.018847429 13.329511 0.01069407
## Cultivar           1 0.007627941 0.007627941  5.394727 0.05923325
## Block              2 0.005837123 0.002918561  2.064101 0.20790108
## Season × Cultivar  1 0.001734116 0.001734116  1.226423 0.31051411
## Residuals          6 0.008483775 0.001413963                     
## 
## -----------------------------------------------------------------
## No significant interaction
## -----------------------------------------------------------------
## 
## -----------------------------------------------------------------
## Season
## -----------------------------------------------------------------
## Multiple Comparison Test: Tukey HSD 
##         resp groups
## E1 0.2164574      a
## E2 0.1371954      b
## NULL
## 
## -----------------------------------------------------------------
## Isolated factors 2 not significant
## -----------------------------------------------------------------
##        Mean
## D 0.1516141
## M 0.2020387
## 
## -----------------------------------------------------------------
## Normality of errors
## -----------------------------------------------------------------
##                          Method Statistic   p.value
##  Shapiro-Wilk normality test(W) 0.9419661 0.5239558
## 
## -----------------------------------------------------------------
## Homogeneity of Variances
## -----------------------------------------------------------------
##                               Method Statistic   p.value
##  Bartlett test(Bartlett's K-squared)  1.034669 0.7928639
## 
## -----------------------------------------------------------------
## Independence from errors
## -----------------------------------------------------------------
##                  Method Statistic   p.value
##  Durbin-Watson test(DW)  2.333798 0.4447802
## 
## -----------------------------------------------------------------
## Additional Information
## -----------------------------------------------------------------
## 
## CV (%) =  14.8
## Mean =  239.6733
## Median =  215.6151
## Possible outliers =  No discrepant point
## 
## -----------------------------------------------------------------
## Analysis of Variance
## -----------------------------------------------------------------
##                   Df    Sum Sq   Mean.Sq   F value        Pr(F)
## Season             1 74120.208 74120.208 58.919373 0.0002557639
## Cultivar           1  2886.460  2886.460  2.294495 0.1806092329
## Block              2  4413.819  2206.910  1.754309 0.2512474277
## Season × Cultivar  1  4475.296  4475.296  3.557487 0.1082312616
## Residuals          6  7547.963  1257.994                       
## 
## -----------------------------------------------------------------
## No significant interaction
## -----------------------------------------------------------------
## 
## -----------------------------------------------------------------
## Season
## -----------------------------------------------------------------
## Multiple Comparison Test: Tukey HSD 
##        resp groups
## E1 318.2651      a
## E2 161.0814      b
## NULL
## 
## -----------------------------------------------------------------
## Isolated factors 2 not significant
## -----------------------------------------------------------------
##       Mean
## D 224.1640
## M 255.1826
## 
## -----------------------------------------------------------------
## Normality of errors
## -----------------------------------------------------------------
##                          Method Statistic    p.value
##  Shapiro-Wilk normality test(W) 0.8811267 0.09060961
## 
## -----------------------------------------------------------------
## Homogeneity of Variances
## -----------------------------------------------------------------
##                               Method Statistic   p.value
##  Bartlett test(Bartlett's K-squared)  1.510064 0.6799495
## 
## -----------------------------------------------------------------
## Independence from errors
## -----------------------------------------------------------------
##                  Method Statistic    p.value
##  Durbin-Watson test(DW)  1.670053 0.08535883
## 
## -----------------------------------------------------------------
## Additional Information
## -----------------------------------------------------------------
## 
## CV (%) =  16.34
## Mean =  554.9572
## Median =  560.4246
## Possible outliers =  No discrepant point
## 
## -----------------------------------------------------------------
## Analysis of Variance
## -----------------------------------------------------------------
##                   Df    Sum Sq   Mean.Sq   F value      Pr(F)
## Season             1 44961.664 44961.664 5.4685954 0.05796251
## Cultivar           1  4530.762  4530.762 0.5510674 0.48591633
## Block              2 53508.458 26754.229 3.2540622 0.11037664
## Season × Cultivar  1  4049.683  4049.683 0.4925547 0.50908425
## Residuals          6 49330.764  8221.794                     
## 
## -----------------------------------------------------------------
## No significant interaction
## -----------------------------------------------------------------
## 
## -----------------------------------------------------------------
## Isolated factors not significant
## -----------------------------------------------------------------
##           D        M
## E1 578.3669 653.9698
## E2 492.6855 494.8065

# 2. Extract and combine tidy ANOVA tables for the summary
anova_af <- 
  mod_anova_af |> 
  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_af |> 
  filter(term != "Residuals", Parameter != "xpi") |> 
  gt() |> 
  tab_header(
    title = "Supplementary Table S3",
    subtitle = "ANOVA Summary: Model parameters and critical points for Leaf Area"
  ) |> 
  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 S3
ANOVA Summary: Model parameters and critical points for Leaf Area
Parameter term GL SQ QM Fcal p.value
b1 Season 1.00000 2,931.32875 2,931.32875 16.26390 0.00686
b1 Cultivar 1.00000 1,660.48982 1,660.48982 9.21290 0.02294
b1 bloco 2.00000 578.73620 289.36810 1.60550 0.27640
b1 Season:Cultivar 1.00000 241.74771 241.74771 1.34129 0.29082
b2 Season 1.00000 58,634.66721 58,634.66721 17.95547 0.00546
b2 Cultivar 1.00000 3,662.47227 3,662.47227 1.12154 0.33035
b2 bloco 2.00000 20,225.43353 10,112.71677 3.09678 0.11914
b2 Season:Cultivar 1.00000 4,259.83208 4,259.83208 1.30447 0.29692
b3 Season 1.00000 0.00000 0.00000 0.09951 0.76310
b3 Cultivar 1.00000 0.00000 0.00000 0.04059 0.84699
b3 bloco 2.00000 0.00003 0.00001 3.88426 0.08275
b3 Season:Cultivar 1.00000 0.00000 0.00000 0.00367 0.95365
ypi Season 1.00000 0.01885 0.01885 13.32951 0.01069
ypi Cultivar 1.00000 0.00763 0.00763 5.39473 0.05923
ypi bloco 2.00000 0.00584 0.00292 2.06410 0.20790
ypi Season:Cultivar 1.00000 0.00173 0.00173 1.22642 0.31051
xmap Season 1.00000 74,120.20836 74,120.20836 58.91937 0.00026
xmap Cultivar 1.00000 2,886.46026 2,886.46026 2.29449 0.18061
xmap bloco 2.00000 4,413.81913 2,206.90956 1.75431 0.25125
xmap Season:Cultivar 1.00000 4,475.29629 4,475.29629 3.55749 0.10823
xmdp Season 1.00000 44,961.66445 44,961.66445 5.46860 0.05796
xmdp Cultivar 1.00000 4,530.76244 4,530.76244 0.55107 0.48592
xmdp bloco 2.00000 53,508.45801 26,754.22900 3.25406 0.11038
xmdp Season:Cultivar 1.00000 4,049.68342 4,049.68342 0.49255 0.50908
Show code 

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

Correlation between model parameters

Show code 
# Correlation between parameters and critical points

  parameters_af |> 
  select(all_of(vars_af)) |> 
  metan::corr_coef() |> 
  plot()