Modeling soybean canopy cover"

In flexFitR: Flexible Non-Linear Least Square Model Fitting

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Estimating soybean canopy development phases

We use canopy-cover data from the 2022 soybean season collected with ETH Zurich's Field Phenotyping Platform (Keller et al., 2026). The dataset contains 78 plots measured on 30 dates, with canopy cover expressed as a proportion from 0 to 1. The seasonal trajectories include canopy expansion, a period of maximum cover, and canopy decline. We use a piecewise function to estimate four transition times: t1, the onset of canopy expansion; t2, the start of the maximum-cover plateau; t3, the onset of canopy decline; and t4, the end of canopy decline. The parameters k and n represent maximum and terminal canopy cover, respectively.

library(flexFitR)
library(dplyr)
library(ggpubr)
library(ggplot2)

data(dt_soybean_22)
head(dt_soybean_22)

1. Exploring data

We start with explorer(), which summarizes the series and lets us look at the temporal evolution of every plot before model fitting.

ex <- explorer(dt_soybean_22, x = time_since_sowing, y = Canopy_cover, id = plot.UID)
names(ex)

plot(ex, type = "evolution", add_avg = TRUE)

The curve starts flat, rises steeply, plateaus near 0.87, then declines through senescence but stops above zero. The curve does not return completely to zero, indicating that some green canopy cover remains at the end of the observed period. This remaining cover is represented by the parameter n.

2. Regression function

The function takes time (t) first and the parameters after:

\begin{equation} f(t; t_1, t_2, t_3, t_4, k, n) = \begin{cases} 0 & \text{if } t < t_1 \ \dfrac{k}{t_2 - t_1} \cdot (t - t_1) & \text{if } t_1 \leq t \leq t_2 \ k & \text{if } t_2 < t \leq t_3 \ n + (k - n) \cdot \dfrac{t_4 - t}{t_4 - t_3} & \text{if } t_3 < t \leq t_4 \ n & \text{if } t > t_4 \end{cases} \end{equation}

fn_piecewise <- function(t, t1, t2, t3, t4, k, n) {
  ifelse(
    test = t < t1, yes = 0,
    no = ifelse(
      test = t <= t2, yes = k / (t2 - t1) * (t - t1),
      no = ifelse(
        test = t <= t3, yes = k,
        no = ifelse(
          test = t <= t4, yes = n + (k - n) * (t4 - t) / (t4 - t3),
          no = n
        )
      )
    )
  )
}

Before fitting anything, plot_fn() lets us draw the function at our proposed initial values.

initial_vals <- c(t1 = 25, t2 = 62, t3 = 100, t4 = 120, k = 1, n = 0.05)

plot_fn(
  fn = "fn_piecewise",
  params = initial_vals,
  interval = c(0, 151),
  color = "black",
  base_size = 15
)

That is a good starting guess: it reproduces the four breakpoints seen in the evolution plot.

3. Fitting models

We fit 10 plots first. Scaling to the full trial will be shown at the end.

plots_ids <- unique(dt_soybean_22$plot.UID)[1:10]

mod_1 <- dt_soybean_22 |>
  modeler(
    x = time_since_sowing,
    y = Canopy_cover,
    grp = plot.UID,
    fn = "fn_piecewise",
    parameters = initial_vals,
    subset = plots_ids,
    method = c("BFGS", "subplex")
  )
print(mod_1)

Passing more than one optimizer to method makes modeler() try each and keep the best solution per plot. Use list_methods() to see the full set of available optimizers.

plot(mod_1, id = plots_ids[1:4])

knitr::kable(mutate_if(mod_1$param, is.numeric, round, 2))

3.1. Extracting model coefficients and uncertainty measures

coef(), confint(), and vcov() return the parameter estimates, confidence intervals, and variance-covariance matrices, respectively.

coef(mod_1, id = plots_ids[1])

confint(mod_1, id = plots_ids[1])

vcov(mod_1, id = plots_ids[1])$FPSB0160001 |> round(digits = 3)

knitr::kable(mutate_if(metrics(mod_1), is.numeric, round, 2))

4. Plotting options

type = 2 shows the coefficients with their confidence intervals. Restricting parm to the four time parameters keeps them on a common scale:

mod_1 |>
  plot(type = 2, id = plots_ids, parm = c("t1", "t2", "t3", "t4"), label_size = 10) +
  theme(axis.text.x = element_text(angle = 65, hjust = 1))

type = 3 overlays every fitted curve, which is a good way to spot a plot that behaved differently from the rest:

plot(mod_1, type = 3, id = plots_ids)

type = 4 adds confidence (blue) and prediction (red) intervals, and type = 5 plots the first derivative. For this function, the canopy expansion and senescence rates appear as two flat steps:

a <- plot(mod_1, type = 4, id = plots_ids[1], color = "black")
b <- plot(mod_1, type = 5, id = plots_ids[1], color = "black")
ggarrange(a, b)

5. Deriving canopy development traits

The fitted parameters are already stage estimates, but some other interesting quantities we usually compare are differences between them. predict.modeler() accepts a formula involving the fitted parameters and propagates their uncertainty:

durations <- rbind(
  predict(mod_1, formula = ~ t2 - t1, id = plots_ids),
  predict(mod_1, formula = ~ t3 - t2, id = plots_ids),
  predict(mod_1, formula = ~ t4 - t3, id = plots_ids)
)

durations |>
  mutate_if(is.numeric, round, 2) |>
  filter(uid %in% "FPSB0160001") |>
  select(-fn_name) |>
  knitr::kable()

These three read as the duration of canopy expansion, the length of the full-canopy plateau, and the duration of senescence.

We can also get rates rather than durations — the slope of the expansion phase is k / (t2 - t1):

predict(mod_1, formula = ~ k / (t2 - t1), id = plots_ids[1:2]) |>
  mutate_if(is.numeric, round, 3) |>
  knitr::kable()

Integrating the fitted curve provides the area under the canopy-cover curve, expressed in canopy-cover days:

predict(mod_1, x = c(0, 151), type = "auc", id = plots_ids[1:3]) |>
  mutate_if(is.numeric, round, 2) |>
  knitr::kable()

6. Modeling all plots using parallel processing

Finally, the same call scales to all 78 plots by adding the options argument.

mod <- dt_soybean_22 |>
  modeler(
    x = time_since_sowing,
    y = Canopy_cover,
    grp = plot.UID,
    keep = c(location, Year),
    fn = "fn_piecewise",
    parameters = initial_vals,
    method = c("BFGS", "subplex"),
    options = list(progress = TRUE, parallel = TRUE, workers = 5)
  )

7. Conclusion

Using the publicly available soybean ground cover dataset from Keller et al. (2026), we demonstrated how flexFitR can transform time-series observations into interpretable growth parameters. We thank the authors for making the dataset publicly available.

References

Keller, B., Kirchgessner, N., Oppliger, C., Kronenberg, L., Roth, L., Zumsteg, O., Corrado, S., Liebisch, F., Aasen, H., Storni, N., Tschurr, F., Zellweger, H., Betrix, C. A., Barendregt, C., Hund, A., & Walter, A. (2026). FIP 1.0 soybean data: Insights on soybean growth from eight years of high-throughput image field phenotyping. Scientific Data, 13(1), 476. https://doi.org/10.1038/s41597-026-06663-z

flexFitR documentation built on Aug. 22, 2026, 1:09 a.m.