Adding mitchell model

NAMESPACE

@@ -48,6 +48,7 @@ export(kamykowski_1985)
 export(lactin2_1995)
 export(lobry_1991)
 export(lrf_1991)
+export(mitchell_2009)
 export(modifiedgaussian_2006)
 export(oneill_1972)
 export(pawar_2018)

R/get_model_names.R

@@ -13,53 +13,54 @@
 get_model_names <- function(model, returnall=FALSE){
 
   mod_names <- c("analytiskontodimas_2004",
-                   "ashrafi1_2018",
-                   "ashrafi2_2018",
-                   "ashrafi3_2018",
-                   "ashrafi4_2018",
-                   "ashrafi5_2018",
-                   "atkin_2005",
-                   "beta_2012",
-                   "betatypesimplified_2008",
-                   "boatman_2017",
-                   "briere1_1999",
-                   "briere1simplified_1999",
-                   "briere2_1999",
-                   "briere2simplified_1999",
-                   "briereextended_2021",
-                   "briereextendedsimplified_2021",
-                   "delong_2017",
-                   "deutsch_2008",
-                   "eubank_1973",
-                   "flextpc_2024",
-                   "flinn_1991",
-                   "gaussian_1987",
-                   "gaussianmodified_2006",
-                   "hinshelwood_1947",
-                   "janisch1_1925",
-                   "janisch2_1925",
-                   "joehnk_2008",
-                   "johnsonlewin_1946",
-                   "kamykowski_1985",
-                   "lactin2_1995",
-                   "lobry_1991",
-                   "oneill_1972",
-                   "pawar_2018",
-                   "quadratic_2008",
-                   "ratkowsky_1983",
-                   "rezende_2019",
-                   "rosso_1993",
-                   "sharpeschoolfull_1981",
-                   "sharpeschoolhigh_1981",
-                   "sharpeschoollow_1981",
-                   "spain_1982",
-                   "stinner_1974",
-                   "taylorsexton_1972",
-                   "thomas_2012",
-                   "thomas_2017",
-                   "tomlinsonphillips_2015",
-                   "warrendreyer_2006",
-                   "weibull_1995")
+                 "ashrafi1_2018",
+                 "ashrafi2_2018",
+                 "ashrafi3_2018",
+                 "ashrafi4_2018",
+                 "ashrafi5_2018",
+                 "atkin_2005",
+                 "beta_2012",
+                 "betatypesimplified_2008",
+                 "boatman_2017",
+                 "briere1_1999",
+                 "briere1simplified_1999",
+                 "briere2_1999",
+                 "briere2simplified_1999",
+                 "briereextended_2021",
+                 "briereextendedsimplified_2021",
+                 "delong_2017",
+                 "deutsch_2008",
+                 "eubank_1973",
+                 "flextpc_2024",
+                 "flinn_1991",
+                 "gaussian_1987",
+                 "gaussianmodified_2006",
+                 "hinshelwood_1947",
+                 "janisch1_1925",
+                 "janisch2_1925",
+                 "joehnk_2008",
+                 "johnsonlewin_1946",
+                 "kamykowski_1985",
+                 "lactin2_1995",
+                 "lobry_1991",
+                 "mitchell_2009",
+                 "oneill_1972",
+                 "pawar_2018",
+                 "quadratic_2008",
+                 "ratkowsky_1983",
+                 "rezende_2019",
+                 "rosso_1993",
+                 "sharpeschoolfull_1981",
+                 "sharpeschoolhigh_1981",
+                 "sharpeschoollow_1981",
+                 "spain_1982",
+                 "stinner_1974",
+                 "taylorsexton_1972",
+                 "thomas_2012",
+                 "thomas_2017",
+                 "tomlinsonphillips_2015",
+                 "warrendreyer_2006",
+                 "weibull_1995")
 
   mod_deprecated <- c('modifiedgaussian_2006', 'lrf_1991')
 

R/mitchell_2009.R

@@ -0,0 +1,91 @@
+#' Mitchell Angilletta model for fitting thermal performance curves
+#'
+#' @param temp temperature in degrees centigrade
+#' @param topt optimum temperature (ºC) where rate is maximal
+#' @param a scale parameter to convert the value of the cosine density to the appropriate magnitude
+#' @param b parameter dictating the performance breadth
+#' @author Daniel Padfield
+#' @return a numeric vector of rate values based on the temperatures and parameter values provided to the function
+#' @references Mitchell, W. A., & Angilletta Jr, M. J. (2009). Thermal games: frequency-dependent models of thermal adaptation. Functional Ecology, 510-520.
+#' @details Equation:
+#' \deqn{rate=\frac{1}{2 \cdot b} \cdot (1 + cos(\frac{temp - t_{opt}}{b} \cdot \pi)) \cdot a }{%
+#' rate = (1/(2. b)).(1 + cos(((temp - topt)/b).pi)).a}
+#' 
+#' When temperatures fall below topt - b or above topt + b, rates are set to 0 to prevent multimodality.
+#'
+#' Start values in \code{get_start_vals} are derived from the data or sensible values from the literature.
+#'
+#' Limits in \code{get_lower_lims} and \code{get_upper_lims} are derived from the data or based extreme values that are unlikely to occur in ecological settings.
+#' 
+#' @note Generally we found this model easy to fit.
+#' @concept model
+#' @examples
+#' # load in ggplot
+#' library(ggplot2)
+#'
+#' # subset for the first TPC curve
+#' data('chlorella_tpc')
+#' d <- subset(chlorella_tpc, curve_id == 1)
+#'
+#' # get start values and fit model
+#' start_vals <- get_start_vals(d$temp, d$rate, model_name = 'mitchell_2009')
+#' # fit model
+#' mod <- nls.multstart::nls_multstart(rate~mitchell_2009(temp = temp, topt, a, b),
+#' data = d,
+#' iter = c(3,3,3),
+#' start_lower = start_vals - 10,
+#' start_upper = start_vals + 10,
+#' lower = get_lower_lims(d$temp, d$rate, model_name = 'mitchell_2009'),
+#' upper = get_upper_lims(d$temp, d$rate, model_name = 'mitchell_2009'),
+#' supp_errors = 'Y',
+#' convergence_count = FALSE)
+#'
+#' # look at model fit
+#' summary(mod)
+#'
+#' # get predictions
+#' preds <- data.frame(temp = seq(min(d$temp), max(d$temp), length.out = 100))
+#' preds <- broom::augment(mod, newdata = preds)
+#'
+#' # plot
+#' ggplot(preds) +
+#' geom_point(aes(temp, rate), d) +
+#' geom_line(aes(temp, .fitted), col = 'blue') +
+#' theme_bw()
+#' @export mitchell_2009
+
+mitchell_2009 <- function(temp, topt, a, b){
+
+  result <- ( 1 / (2 * b) ) * ( 1 + cos( ( (temp - topt) / b ) * pi ) ) * a
+
+  # Make sure that all values are positive and that the resulting curve is not multimodal
+
+  # if results are negative, make them zero
+  result[which(result <= 0)] <- 0
+  # in this function b is 1/2 the breadth, so topt - b and topt + b are the inflection points of the curve where the line hits 0
+  result[which(temp < topt - b)] <- 0
+  result[which(temp > topt + b)] <- 0
+
+  return(result)
+}  
+
+mitchell_2009.starting_vals <- function(d){
+  topt = mean(d$x[d$y == max(d$y, na.rm = TRUE)])
+  a = 2 * 10^-4
+  b = max(d$x, na.rm = TRUE) - min(d$x, na.rm = TRUE)
+  return(c(topt = topt, a = a, b=b))
+}
+
+mitchell_2009.lower_lims <- function(d){
+  topt = min(d$x, na.rm = TRUE)
+  a = 0
+  b = 0
+  return(c(topt = topt, a = a, b=b))
+}
+
+mitchell_2009.upper_lims <- function(d){
+  topt = max(d$x, na.rm = TRUE) *10
+  a = Inf
+  b = Inf
+  return(c(topt = topt, a = a, b=b))
+}

README.md

@@ -14,7 +14,7 @@ mirror](https://cranlogs.r-pkg.org/badges/grand-total/rTPC)](https://www.r-pkg.o
 <!-- badges: end -->
 
 **rTPC** is an R package that helps fit thermal performance curves
-(TPCs) in R. **rTPC** contains 48 model formulations previously used to
+(TPCs) in R. **rTPC** contains 49 model formulations previously used to
 fit TPCs and has helper functions to help set sensible start parameters,
 upper and lower parameter limits and estimate parameters useful in
 downstream analyses, such as cardinal temperatures, maximum rate and

tests/testthat/test-mitchell_2009.R

@@ -0,0 +1,43 @@
+# do not run the test on CRAN as they take too long
+testthat::skip_on_cran()
+
+# method: fit model and get predictions. Check these against others.
+
+# load in ggplot
+library(ggplot2)
+
+# subset for the first TPC curve
+data('chlorella_tpc')
+d <- subset(chlorella_tpc, curve_id == 1)
+modelname <- "mitchell_2009"
+
+# get start values and fit model
+start_vals <- get_start_vals(d$temp, d$rate, model_name = modelname)
+
+# fit model
+mod <- nls.multstart::nls_multstart(rate~mitchell_2009(temp = temp, topt, a, b),
+                                    data = d,
+                                    iter = c(3,3,3),
+                                    start_lower = start_vals - 10,
+                                    start_upper = start_vals + 10,
+                                    lower = get_lower_lims(d$temp, d$rate, model_name = modelname),
+                                    upper = get_upper_lims(d$temp, d$rate, model_name = modelname),
+                                    supp_errors = 'Y',
+                                    convergence_count = FALSE)
+
+# get predictions
+preds <- broom::augment(mod)
+# dput(round(preds$.fitted, 1))
+
+# plot
+ggplot(preds) +
+  geom_point(aes(temp, rate)) +
+  geom_line(aes(temp, .fitted)) +
+  theme_bw()
+
+# run test
+testthat::test_that(paste(modelname, "function works"), {
+  testthat::expect_equal(
+    round(preds$.fitted, 1),
+    c(0, 0, 0.1, 0.4, 0.8, 1.2, 1.4, 1.5, 1.3, 1, 0.6, 0.2))
+})

vignettes/fit_many_models.Rmd

@@ -355,6 +355,15 @@ d_fits <- nest(d, data = c(temp, rate)) %>%
                         upper = get_upper_lims(.x$temp, .x$rate, model_name = 'lobry_1991'),
                         supp_errors = 'Y',
                         convergence_count = FALSE)),
+         mitchell = map(data, ~nls_multstart(rate~mitchell_2009(temp = temp, topt, a, b),
+                        data = .x,
+                        iter = c(4,4,4),
+                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'mitchell_2009') - 10,
+                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'mitchell_2009') + 10,
+                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'mitchell_2009'),
+                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'mitchell_2009'),
+                        supp_errors = 'Y',
+                        convergence_count = FALSE)),
          oneill = map(data, ~nls_multstart(rate~oneill_1972(temp = temp, rmax, ctmax, topt, q10),
                         data = .x,
                         iter = c(4,4,4,4),