add generic function to model curves

DESCRIPTION

@@ -1,37 +1,38 @@
-Package: exploreHTP
-Type: Package
-Title: Explore High-Throughput Phenotypic (HTP) Data
-Version: 0.1.0
-Authors@R: 
-    c(person(given = "Johan",
-             family = "Aparicio",
-             role = c("aut", "cre"),
-             email = "[email protected]"),
-      person(given = "Jeffrey",
-             family = "Endelman",
-             role = "ctb",
-             email = "[email protected]"),
-      person(given = "University of Wisconsin Madison",
-             role = "cph"))
-Maintainer: Johan Aparicio <[email protected]> 
-Description: Utilities designed to explore high-throughput phenotypic (HTP) data
-    for field trials in plant breeding.
-License: MIT + file LICENSE
-Encoding: UTF-8
-LazyData: true
-Depends: 
-    R (>= 2.10)
-RoxygenNote: 7.3.1
-Author: Johan Aparicio [cre, aut], Jeffrey Endelman [ctb],
- University of Wisconsin Madison [cph]
-Imports: 
-    agriutilities,
-    dplyr,
-    ggplot2,
-    optimx,
-    subplex,
-    tibble,
-    tidyr
-URL: https://apariciojohan.github.io/exploreHTP/,
-     https://github.com/AparicioJohan/exploreHTP
-BugReports: https://github.com/AparicioJohan/exploreHTP/issues
+Package: exploreHTP
+Type: Package
+Title: Explore High-Throughput Phenotypic (HTP) Data
+Version: 0.1.0
+Authors@R: 
+    c(person(given = "Johan",
+             family = "Aparicio",
+             role = c("aut", "cre"),
+             email = "[email protected]"),
+      person(given = "Jeffrey",
+             family = "Endelman",
+             role = "ctb",
+             email = "[email protected]"),
+      person(given = "University of Wisconsin Madison",
+             role = "cph"))
+Maintainer: Johan Aparicio <[email protected]> 
+Description: Utilities designed to explore high-throughput phenotypic (HTP) data
+    for field trials in plant breeding.
+License: MIT + file LICENSE
+Encoding: UTF-8
+LazyData: true
+Depends: 
+    R (>= 2.10)
+RoxygenNote: 7.3.1
+Author: Johan Aparicio [cre, aut], Jeffrey Endelman [ctb],
+ University of Wisconsin Madison [cph]
+Imports: 
+    agriutilities,
+    dplyr,
+    ggplot2,
+    optimx,
+    rlang,
+    subplex,
+    tibble,
+    tidyr
+URL: https://apariciojohan.github.io/exploreHTP/,
+     https://github.com/AparicioJohan/exploreHTP
+BugReports: https://github.com/AparicioJohan/exploreHTP/issues

NAMESPACE

@@ -3,6 +3,7 @@
 S3method(plot,canopy_HTP)
 S3method(plot,height_HTP)
 S3method(plot,maturity_HTP)
+S3method(plot,modeler_HTP)
 S3method(plot,read_HTP)
 export(canopy_HTP)
 export(fn_exp1_exp)
@@ -15,6 +16,7 @@ export(fn_lin_pl_lin3)
 export(fn_piwise)
 export(height_HTP)
 export(maturity_HTP)
+export(modeler_HTP)
 export(read_HTP)
 export(sse_exp1_exp)
 export(sse_exp1_lin)

R/04_maturity.R

@@ -3,9 +3,9 @@
 #' @param results An object of class \code{read_HTP}, containing the results of the \code{read_HTP()} function.
 #' @param canopy An object of class \code{canopy_HTP}, containing the results of the \code{canopy_HTP()} function.
 #' @param index A string specifying the trait to be modeled. Default is \code{"GLI"}.
-#' @param check_negative Logical. Convert negative values to zero. \code{TRUE} by default.
 #' @param plot_id An optional vector of plot IDs to filter the data. Default is \code{NULL}, meaning all plots are used.
 #' @param add_zero \code{TRUE} or \code{FALSE}. Add zero to the time series.\code{TRUE} by default.
+#' @param check_negative Logical. Convert negative values to zero. \code{TRUE} by default.
 #' @param method A vector of the methods to be used, each as a character string.
 #' See optimx package. c("subplex", "pracmanm", "anms") by default.
 #' @param return_method \code{TRUE} or \code{FALSE}. To return the method selected for the
@@ -69,9 +69,9 @@
 maturity_HTP <- function(results,
                          canopy,
                          index = "GLI",
-                         check_negative = TRUE,
-                         add_zero = TRUE,
                          plot_id = NULL,
+                         add_zero = TRUE,
+                         check_negative = TRUE,
                          method = c("subplex", "pracmanm", "anms"),
                          return_method = FALSE,
                          parameters = c(t1 = 38.7, t2 = 62, t3 = 90, k = 0.32, beta = -0.01),

R/05_growth_curves.R

@@ -738,24 +738,26 @@ sse_lin_pl_lin3 <- function(params, t, y) {
 }
 
 
+#' @examples
+#' params <- c(34.9, 61.8, 100)
+#' t <- c(0, 29, 36, 42, 56, 76, 92, 100, 108)
+#' y <- c(0, 0, 4.379, 26.138, 78.593, 100, 100, 100, 100)
+#' curve <- "fn_piwise"
+#' sse_generic(params, t, y, curve)
+#' @noRd
+sse_generic <- function(params, t, y, curve) {
+  values <- paste(params, collapse = ", ")
+  string <- paste("sapply(t, FUN = ", curve, ", ", values, ")", sep = "")
+  y_hat <- eval(parse(text = string))
+  sse <- sum((y - y_hat)^2)
+  return(sse)
+}
 
-# params <- c(34.9, 61.8, 100)
-# t <- c(0, 29, 36, 42, 56, 76, 92, 100, 108)
-# y <- c(0, 0, 4.379, 26.138, 78.593, 100, 100, 100, 100)
-# fn <- fn_piwise
-#
-# sse_generic_can <- function(params, t, y, fn) {
-#   values <- paste(c(params[1:2], max(y, na.rm = TRUE)), collapse = ", ")
-#   values <- paste(c(params), collapse = ", ")
-#   string <- paste("sapply(t, FUN = fn, ", values, ")", sep = "")
-#   y_hat <- eval(parse(text = string))
-#   sse <- sum((y - y_hat)^2)
-#   return(sse)
-# }
-#
-# sse_generic_can(
-#   params <- c(34.9, 61.8, 100),
-#   t <- c(0, 29, 36, 42, 56, 76, 92, 100, 108),
-#   y <- c(0, 0, 4.379, 26.138, 78.593, 100, 100, 100, 100),
-#   fn <- fn_piwise
-# )
+#' @noRd
+create_call <- function(fn = "fn_piwise") {
+  arg <- formals(fn)
+  values <- paste(names(arg)[-1], collapse = ", ")
+  string <- paste(fn, "(time, ", values, ")", sep = "")
+  out <- rlang::parse_expr(string)
+  return(out)
+}

R/06_modeler.R

@@ -0,0 +1,179 @@
+#' Modeler HTP
+#'
+#' @param results An object of class \code{read_HTP}, containing the results of the \code{read_HTP()} function.
+#' @param index A string specifying the trait to be modeled. Default is \code{"GLI"}.
+#' @param plot_id An optional vector of plot IDs to filter the data. Default is \code{NULL}, meaning all plots are used.
+#' @param add_zero \code{TRUE} or \code{FALSE}. Add zero to the time series.\code{TRUE} by default.
+#' @param check_negative Logical. Convert negative values to zero. \code{TRUE} by default.
+#' @param method A vector of the methods to be used, each as a character string.
+#' See optimx package. c("subplex", "pracmanm", "anms") by default.
+#' @param return_method \code{TRUE} or \code{FALSE}. To return the method selected for the
+#' optimization in the table. \code{FALSE} by default.
+#' @param parameters Initial values for the parameters to be optimized over. c(t1 = 38.7, t2 = 62, t3 = 90, k = 0.32, beta = -0.01) by default.
+#' @param lower Bounds on the variables for methods such as "L-BFGS-B" that can handle box (or bounds) constraints.
+#' @param upper Bounds on the variables for methods such as "L-BFGS-B" that can handle box (or bounds) constraints.
+#' @param initial_vals data.frame with columns <plot, genotype, t1, t2, ... , and all the initial parameters>.
+#' @param fn String character with the name of the function "fn_lin_pl_lin".
+#' @param n_points Number of time points to approximate the AUC. 1000 by default.
+#' @param max_time Maximum time value for calculating the AUC. \code{NULL} by default takes the last time point.
+#' @return An object of class \code{modeler_HTP}, which is a list containing the following elements:
+#' \describe{
+#'   \item{\code{param}}{A data frame containing the optimized parameters and related information.}
+#'   \item{\code{dt}}{A data frame with data used.}
+#'   \item{\code{fn}}{The call used to calculate the AUC.}
+#'   \item{\code{max_time}}{Maximum time value used for calculating the AUC.}
+#' }
+#' @export
+#'
+#' @examples
+#' library(exploreHTP)
+#' data(dt_potato)
+#' dt_potato <- dt_potato
+#' results <- read_HTP(
+#'   data = dt_potato,
+#'   genotype = "Gen",
+#'   time = "DAP",
+#'   plot = "Plot",
+#'   traits = c("Canopy", "GLI_2"),
+#'   row = "Row",
+#'   range = "Range"
+#' )
+#' names(results)
+#' mat <- modeler_HTP(
+#'   results = results,
+#'   index = "GLI_2",
+#'   plot_id = c(195, 40),
+#'   parameters = c(t1 = 38.7, t2 = 62, t3 = 90, k = 0.32, beta = -0.01),
+#'   fn = "fn_lin_pl_lin",
+#' )
+#' plot(mat, plot_id = c(195, 40))
+#' mat$param
+#' can <- modeler_HTP(
+#'   results = results,
+#'   index = "Canopy",
+#'   plot_id = c(195, 40),
+#'   parameters = c(t1 = 45,t2 = 80, k = 0.9),
+#'   fn = "fn_piwise",
+#' )
+#' plot(can, plot_id = c(195, 40))
+#' can$param
+#' @import optimx
+#' @import tibble
+#' @import dplyr
+modeler_HTP <- function(results,
+                        index = "GLI",
+                        plot_id = NULL,
+                        check_negative = TRUE,
+                        add_zero = TRUE,
+                        method = c("subplex", "pracmanm", "anms"),
+                        return_method = FALSE,
+                        parameters = NULL,
+                        lower = -Inf,
+                        upper = Inf,
+                        initial_vals = NULL,
+                        fn = "fn_piwise",
+                        n_points = 1000,
+                        max_time = NULL) {
+  if (!inherits(results, "read_HTP")) {
+    stop("The object should be of read_HTP class")
+  }
+
+  dt <- results$dt_long |>
+    filter(trait %in% index) |>
+    filter(!is.na(value)) |>
+    droplevels()
+
+  if (check_negative) {
+    dt <- mutate(dt, value = ifelse(value < 0, 0, value))
+  }
+
+  if (add_zero) {
+    dt <- dt |>
+      mutate(time = 0, value = 0) |>
+      unique.data.frame() |>
+      rbind.data.frame(dt) |>
+      arrange(plot, time)
+  }
+
+  if (!is.null(initial_vals)) {
+    init <- initial_vals |>
+      pivot_longer(
+        cols = -c(plot, genotype),
+        names_to = "coef",
+        values_to = "val"
+      ) |>
+      nest_by(plot, genotype, .key = "param") |>
+      mutate(param = list(pull(param, val, coef)))
+  } else {
+    init <- dt |>
+      select(plot, genotype) |>
+      unique.data.frame() |>
+      cbind(data.frame(t(parameters))) |>
+      pivot_longer(
+        cols = -c(plot, genotype),
+        names_to = "coef",
+        values_to = "val"
+      ) |>
+      nest_by(plot, genotype, .key = "param") |>
+      mutate(param = list(pull(param, val, coef)))
+  }
+
+  if (!is.null(plot_id)) {
+    dt <- droplevels(filter(dt, plot %in% plot_id))
+    init <- droplevels(filter(init, plot %in% plot_id))
+  }
+
+  dt_nest <- dt |>
+    nest_by(plot, genotype, row, range) |>
+    full_join(init, by = c("plot", "genotype"))
+
+  param_mat <- dt_nest |>
+    summarise(
+      res = list(
+        opm(
+          par = param,
+          fn = sse_generic,
+          t = data$time,
+          y = data$value,
+          curve = fn,
+          method = method,
+          lower = lower,
+          upper = upper
+        ) |>
+          rownames_to_column(var = "method") |>
+          arrange(value) |>
+          rename(sse = value) |>
+          select(2:(length(parameters) + 1), method, sse) |>
+          slice(1)
+      ),
+      .groups = "drop"
+    ) |>
+    unnest(cols = res)
+
+  if (is.null(max_time)) {
+    max_time <- max(dt$time, na.rm = TRUE)
+  }
+
+  # AUC
+  density <- create_call(fn)
+  sq <- seq(0, max_time, length.out = n_points)
+  auc <- full_join(
+    x = expand.grid(time = sq, plot = unique(dt$plot)),
+    y = param_mat,
+    by = "plot"
+  ) |>
+    group_by(time, plot) |>
+    mutate(hat = !!density) |>
+    group_by(plot) |>
+    mutate(trapezoid_area = (lead(hat) + hat) / 2 * (lead(time) - time)) |>
+    filter(!is.na(trapezoid_area)) |>
+    summarise(total_area = sum(trapezoid_area))
+  param_mat <- full_join(param_mat, auc, by = "plot")
+
+  if (!return_method) {
+    param_mat <- select(param_mat, -method)
+  }
+  out <- list(param = param_mat, dt = dt, fn = density, max_time = max_time)
+  class(out) <- "modeler_HTP"
+  return(invisible(out))
+}

R/global.R

@@ -36,6 +36,7 @@ utils::globalVariables(
     "coef",
     "trapezoid_area",
     "hat",
-    "k"
+    "k",
+    "param"
   )
 )