inst/doc/modeling_soybean.R

In flexFitR: Flexible Non-Linear Least Square Model Fitting

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

## ----setup, warning=FALSE, message=FALSE--------------------------------------
library(flexFitR)
library(dplyr)
library(ggpubr)
library(ggplot2)

## -----------------------------------------------------------------------------
data(dt_soybean_22)
head(dt_soybean_22)

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

## ----fig.width= 8, fig.height=4, fig.alt="plot evolution"---------------------
plot(ex, type = "evolution", add_avg = TRUE)

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

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

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

## ----fig.width= 8, fig.height=5, fig.alt="plot fit"---------------------------
plot(mod_1, id = plots_ids[1:4])

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

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

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

## ----fig.width= 8, fig.height=4, fig.alt="plot curves"------------------------
plot(mod_1, type = 3, id = plots_ids)

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

## -----------------------------------------------------------------------------
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()

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

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

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