Add functions for quickly fitting tpcs

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: rTPC
 Title: Fitting and Analysing Thermal Performance Curves
-Version: 1.0.7
+Version: 1.1.0
 Authors@R: 
     c(person(given = "Daniel",
              family = "Padfield",
@@ -16,34 +16,40 @@ Authors@R:
              role = "aut",
              email = "francis.windram17@imperial.ac.uk"))
 Maintainer: Daniel Padfield <d.padfield@exeter.ac.uk>
-Description: Helps to fit thermal performance curves (TPCs). 'rTPC' contains 26 model formulations previously used to fit TPCs and has helper functions to set sensible start parameters, upper and lower parameter limits and estimate parameters useful in downstream analyses, such as cardinal temperatures, maximum rate and optimum temperature. See Padfield et al. (2021) <doi:10.1111/2041-210X.13585>.
+Description: Helps to fit thermal performance curves (TPCs). 'rTPC' contains 49 model formulations previously used to fit TPCs and has helper functions to set sensible start parameters, upper and lower parameter limits and estimate parameters useful in downstream analyses, such as cardinal temperatures, maximum rate and optimum temperature. See Padfield et al. (2021) <doi:10.1111/2041-210X.13585>.
 License: GPL (>= 3)
 URL: https://github.com/padpadpadpad/rTPC, https://padpadpadpad.github.io/rTPC/
 BugReports: https://github.com/padpadpadpad/rTPC/issues
 Depends: 
     R (>= 2.10)
 Imports: 
     cli,
+    glue,
     rlang,
     stats
 Suggests:
     boot,
     broom,
     car,
+    dplyr,
     forcats,
-    ggplot2,
     ggrepel,
     knitr,
+    lubridate,
     minpack.lm,
+    mirai,
     MuMIn,
     nls.multstart,
     nlstools,
     patchwork,
     progress,
+    purrr,
     RColorBrewer,
     rmarkdown,
     stringr,
     testthat,
+    tibble,
+    tidyr,
     tidyverse
 VignetteBuilder: 
     knitr

NAMESPACE

@@ -38,6 +38,7 @@ export(get_start_vals)
 export(get_thermalsafetymargin)
 export(get_thermaltolerance)
 export(get_topt)
+export(get_tpc_as_formula)
 export(get_upper_lims)
 export(hinshelwood_1947)
 export(janisch1_1925)
@@ -53,6 +54,8 @@ export(modifiedgaussian_2006)
 export(oneill_1972)
 export(pawar_2018)
 export(quadratic_2008)
+export(quickfit_tpc)
+export(quickfit_tpc_multi)
 export(ratkowsky_1983)
 export(rezende_2019)
 export(rosso_1993)

R/get_tpc_as_formula.R

@@ -0,0 +1,47 @@
+#' Get a formula object for calling a TPC
+#'
+#' @param model_name the name of the model being fitted
+#' @param temp the name of the temperature column
+#' @param trait the name of the trait column
+#' @param explicit whether to return the formula constructed using the explicit form of the tpc function (e.g. \code{rTPC::briere1_1999()})
+#' 
+#' @author Francis Windram
+#' @return A formula calling the expected TPC
+#' @concept helper
+#' 
+#' @examples
+#' get_tpc_as_formula("briere1_1999", "temperature", "rate")
+#' # > rate ~ briere1_1999(temp = temperature, tmin, tmax, a)
+#' 
+#' @export get_tpc_as_formula
+
+get_tpc_as_formula <- function(model_name, temp, trait, explicit = FALSE){
+
+  if (length(temp) != 1) {
+    cli::cli_abort(c("x"="Supplied {.arg temp} is not a column name (contains {.val {length(temp)}} element{?s})"))
+  }
+
+  if (length(trait) != 1) {
+    cli::cli_abort(c("x"="Supplied {.arg trait} is not a column name (contains {.val {length(trait)}} element{?s})"))
+  }
+
+  mod_names <- rTPC::get_model_names(returnall = TRUE)
+  model_name <- tryCatch(rlang::arg_match(model_name, mod_names), error = function(e){
+    cli::cli_abort(c("x"="Supplied {.arg model_name} ({.val {model_name}}) is not an available model in rTPC.",
+                     "!"="Please check the spelling of {.arg model_name}.",
+                     " "="(run {.fn rTPC::get_model_names} to see all valid names.)",
+                     ""), call=rlang::caller_env(n=4))
+  })
+
+  tpcfunc <- get(model_name, envir = asNamespace("rTPC"))
+  tpcargs <- names(formals(tpcfunc))
+
+  if (explicit) {
+    model <- paste0("rTPC::", model_name) 
+  } else {
+    model <- model_name
+  }
+
+  formulastr <- glue::glue("{trait}~{model}(temp = {temp}, {paste(tpcargs[tpcargs != 'temp'], collapse = ', ')})")
+  return(stats::as.formula(formulastr))
+}

R/quickfit_tpc.R

@@ -0,0 +1,83 @@
+#' Perform a quick tpc fit
+#'
+#' @description Performs a simple TPC fit using \code{\link[nls.multstart]{nls_multstart}}. This function tries to use a sensible default configuration, 
+#' however if you need to use the more esoteric elements of \code{\link[nls.multstart]{nls_multstart}} then you will need to construct your own.
+#' @param data the data to fit a model to
+#' @param model_name the model name as a string
+#' @param temp the column name (as a string) containing the temperature data
+#' @param trait the column name (as a string) containing the temperature data
+#' @param start_adjusts any adjustments to make to the lower and upper starting bounds. If \code{0 < start_adjusts < 1}, this will be interpreted as a proportion of the base starting values.
+#' @param iter number of combinations of starting parameters which will be tried (as in \code{\link[nls.multstart]{nls_multstart}})
+#' @param lhstype method to use for Latin Hypercube Sampling using \code{\link[lhs]{lhs}} (as in \code{\link[nls.multstart]{nls_multstart}})
+#' @param gridstart whether to run a gridstart approach (interpreting iter as the number of samples to take across each parameter, 
+#' so \code{3} will become \code{c(3,3,3)} for a 3-parameter model)
+#' @param force whether to force a gridstart even with very large numbers of iterations
+#' @author Francis Windram
+#' @return The nls model object of the fit model
+#' @concept helper
+#' 
+#' @examples
+#' \dontrun{
+#' data("chlorella_tpc")
+#' subs <- subset(chlorella_tpc, curve_id == 1)
+#' quickfit_tpc(subs, "briere1_1999", "temp", "rate")
+#' 
+#' quickfit_tpc(
+#'   subs,
+#'   "briere1_1999",
+#'   "temp",
+#'   "rate",
+#'   start_adjusts = 10,
+#'   iter = 150,
+#'   lhstype = "maximin"
+#' )
+#' 
+#' quickfit_tpc(
+#'   subs,
+#'   "briere1_1999",
+#'   "temp",
+#'   "rate",
+#'   start_adjusts = 10,
+#'   iter = 5,
+#'   gridstart = TRUE
+#' )
+#' }
+#' 
+#' @export quickfit_tpc
+
+
+quickfit_tpc <- function(data, model_name, temp, trait, start_adjusts = 0, iter = 150, lhstype = "none", gridstart=FALSE, force=FALSE) {
+
+  rlang::check_installed("nls.multstart")
+
+  # The form of column extraction used below only works properly on dfs.
+  data <- as.data.frame(data)
+
+  start_vals <- rTPC::get_start_vals(data[,temp], data[,trait], model_name = model_name)
+
+  if (0 < start_adjusts && start_adjusts < 1) {
+    start_adjusts <- start_vals * start_adjusts
+  }
+
+  if (gridstart){
+    iter <- rep(iter, length(start_vals))
+    total_iterations <- prod(iter)
+    if (total_iterations > 1500 && !force) {
+      cli::cli_abort(c(
+        "x"="Massive gridstart detected! ({.val {total_iterations}} iterations)",
+        "!"="If this is intended, set force={.val {TRUE}}, else set a lower {.arg iter} value such as {.val {3}}."))
+    }
+  }
+
+  mod <- nls.multstart::nls_multstart(rTPC::get_tpc_as_formula(model_name, temp, trait),
+                       data = data,
+                       iter = iter,
+                       start_lower = start_vals - start_adjusts,
+                       start_upper = start_vals + start_adjusts,
+                       lower = rTPC::get_lower_lims(data[, temp], data[, trait], model_name = model_name),
+                       upper = rTPC::get_upper_lims(data[, temp], data[, trait], model_name = model_name),
+                       supp_errors = 'Y',
+                       lhstype = lhstype,
+                       convergence_count = FALSE)
+  return(mod)
+}

R/quickfit_tpc_multi.R

@@ -0,0 +1,117 @@
+#' Perform a parallelised quick tpc fit across models and curves
+#'
+#' @description Performs a parallelised TPC fit using \code{\link[nls.multstart]{nls_multstart}} and \code{\link[purrr]{map}}. This function tries to use a sensible default configuration, 
+#' however if you need to use the more esoteric elements of \code{\link[nls.multstart]{nls_multstart}} then you may need to construct your own running script.
+#' @param d the data to fit a model to
+#' @param model_names a vector of model names to fit as strings
+#' @param temp the column name (as a string) containing the temperature data
+#' @param trait the column name (as a string) containing the temperature data
+#' @param start_adjusts any adjustments to make to the lower and upper starting bounds. If \code{0 < start_adjusts < 1}, this will be interpreted as a proportion of the base starting values.
+#' @param iter number of combinations of starting parameters which will be tried (as in \code{\link[nls.multstart]{nls_multstart}})
+#' @param lhstype method to use for Latin Hypercube Sampling using \code{\link[lhs]{lhs}} (as in \code{\link[nls.multstart]{nls_multstart}})
+#' @param gridstart whether to run a gridstart approach (interpreting iter as the number of samples to take across each parameter, 
+#' so \code{3} will become \code{c(3,3,3)} for a 3-parameter model)
+#' @param force whether to force a gridstart even with very large numbers of iterations
+#' 
+#' @note The parameters \code{temp}, \code{trait}, \code{start_adjusts}, \code{iter}, \code{lhstype}, \code{gridstart}, or \code{force} can be specified per-model
+#'  by providing a vector of values of a length equal to the number of models to be fit.
+#' 
+#' @author Francis Windram
+#' @return A tibble of model fits
+#' @concept helper
+#' 
+#' @examples
+#' \dontrun{
+#' data("chlorella_tpc")
+#' d2 <- subset(chlorella_tpc, curve_id <= 60)
+#' 
+#' # Set up daemons for parallelism
+#' mirai::daemons(2)
+#' 
+#' quickfit_tpc_multi(d2, c("briere1_1999", "briere2_1999"), "temp", "rate")
+#' 
+#' quickfit_tpc_multi(d2, c("briere1_1999", "briere2_1999"), "temp", "rate", start_adjusts = 10)
+#' 
+#' quickfit_tpc_multi(d2, c("briere1_1999", "briere2_1999"), "temp", "rate", iter = c(100, 150))
+#' 
+#' mirai::daemons(0)
+#' }
+#' @export quickfit_tpc_multi
+
+quickfit_tpc_multi <- function(d, model_names, temp, trait, start_adjusts = 0, iter = 150, lhstype = "none", gridstart=FALSE, force=FALSE) {
+
+  rlang::check_installed(c("purrr", "mirai", "dplyr", "tidyr"), version = c("1.1.0", NA, NA, NA))
+
+  if (!mirai::daemons_set()) {
+    cli::cli_abort(c(
+      "x" = "To perform parallel fitting, you must initialise mirai daemons.",
+      ">" = "Please run {.fn mirai::daemons} with the number of cores you wish to use.",
+      ""))
+  }
+
+  rlang::check_installed(c("nls.multstart", "mirai", "carrier"))
+
+  # Check model_names to make sure models are in the list (rather than doing this on each process)
+  mod_names <- rTPC::get_model_names(returnall = TRUE)
+  to_fit <- intersect(model_names, mod_names)
+  unknown <- setdiff(model_names, mod_names)
+
+  if (length(to_fit) < 1) {
+    cli::cli_abort(c("x" = "No valid models found! Got: {.val {model_names}}",
+                     "!"="Please check the spellings within {.arg model_names}.",
+                     " "="(run {.fn rTPC::get_model_names} to see all valid names.)",
+                     ""))
+  } else if (length(unknown) > 0) {
+    cli::cli_warn(c("!" = "Unknown models specified. Ignoring {.val {unknown}}"))
+  }
+
+  nesting_cols <- c(temp, trait)
+
+  fit_models <- d |> 
+    tidyr::nest(data = dplyr::all_of(nesting_cols)) |> 
+    dplyr::mutate(
+      mods = purrr::map(
+        data,
+        purrr::in_parallel(
+          \(x) .quickfit_multi(x, 
+                               model_names = to_fit,
+                               temp = temp,
+                               trait = trait,
+                               start_adjusts = start_adjusts,
+                               iter = iter,
+                               lhstype = lhstype,
+                               gridstart = gridstart,
+                               force = force
+          ),
+          .quickfit_multi = .quickfit_multi,
+          to_fit = to_fit,
+          temp = temp,
+          trait = trait,
+          start_adjusts = start_adjusts,
+          iter = iter,
+          lhstype = lhstype,
+          gridstart = gridstart,
+          force = force
+        ),
+        .progress = TRUE
+      )
+    ) |>
+    tidyr::unnest(mods)
+
+  return(fit_models)
+
+}
+
+.quickfit_multi <- function(data, model_names, temp, trait, start_adjusts=0, iter=150, lhstype="none", gridstart=FALSE, force=FALSE) {
+  # Coerce args (basically make df of params)
+  df <- tibble::tibble(model=model_names, temp, trait, start_adjusts, iter, lhstype , gridstart, force)
+
+  # Map tpcs using quickfit
+  out <- purrr::pmap(df, \(...) list(rTPC::quickfit_tpc(data, ...)))
+
+  # Get some lovely names
+  names(out) <- model_names
+
+  # Cast to tibble and return
+  return(tibble::as.tibble(out))
+}