---
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.
<https://doi.org/10.1038/s41597-026-06663-z>

<br>
