--- title: "Modeling soybean canopy cover" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Modeling soybean canopy cover} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} 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. ```{r setup, warning=FALSE, message=FALSE} library(flexFitR) library(dplyr) library(ggpubr) library(ggplot2) ``` ```{r} 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. ```{r} ex <- explorer(dt_soybean_22, x = time_since_sowing, y = Canopy_cover, id = plot.UID) names(ex) ``` ```{r, fig.width= 8, fig.height=4, fig.alt="plot evolution"} 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} ```{r} 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. ```{r, fig.width= 8, fig.height=4, fig.alt="plot fn"} 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. ```{r, warning=FALSE, message=FALSE} 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. ```{r, fig.width= 8, fig.height=5, fig.alt="plot fit"} plot(mod_1, id = plots_ids[1:4]) ``` ```{r} 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. ```{r} coef(mod_1, id = plots_ids[1]) ``` ```{r} confint(mod_1, id = plots_ids[1]) ``` ```{r} vcov(mod_1, id = plots_ids[1])$FPSB0160001 |> round(digits = 3) ``` ```{r} 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: ```{r, fig.width= 8, fig.height=5, fig.alt="plot coef"} 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: ```{r, fig.width= 8, fig.height=4, fig.alt="plot curves"} 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: ```{r, fig.width= 8, fig.height=4, fig.alt="plot fit and derivative"} 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: ```{r} 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) ) ``` ```{r} 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)`: ```{r} 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: ```{r} 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. ```{r, eval = FALSE} 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.