inst/doc/fit_many_models.R

In rTPC: Fitting and Analysing Thermal Performance Curves

## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  #tidy.opts=list(width.cutoff=60),
  #tidy=TRUE,
  fig.align = 'center',
  warning=FALSE
)

## ----setup, message=FALSE-----------------------------------------------------
# load packages
library(rTPC)
library(nls.multstart)
library(broom)
library(tidyverse)

# write function to label ggplot2 panels
label_facets_num <- function(string){
  len <- length(string)
  string = paste('(', 1:len, ') ', string, sep = '')
  return(string)
}

## ----first_plot, fig.width=6, fig.height = 4----------------------------------
# load in data
data("chlorella_tpc")

# keep just a single curve
d <- filter(chlorella_tpc, curve_id == 1)

# show the data
ggplot(d, aes(temp, rate)) +
  geom_point() +
  theme_bw(base_size = 12) +
  labs(x = 'Temperature (ºC)',
       y = 'Metabolic rate',
       title = 'Respiration across temperatures')


## ----many_mods, message=FALSE-------------------------------------------------
# fit every model formulation in rTPC
d_fits <- nest(d, data = c(temp, rate)) %>%
  mutate(beta = map(data, ~nls_multstart(rate~beta_2012(temp = temp, a, b, c, d, e),
                        data = .x,
                        iter = c(6,6,6,6,6),
                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'beta_2012') - 10,
                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'beta_2012') + 10,
                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'beta_2012'),
                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'beta_2012'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)),
         boatman = map(data, ~nls_multstart(rate~boatman_2017(temp = temp, rmax, tmin, tmax, a,b),
                        data = .x,
                        iter = c(4,4,4,4,4),
                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'boatman_2017') - 10,
                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'boatman_2017') + 10,
                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'boatman_2017'),
                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'boatman_2017'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)),
         briere2 = map(data, ~nls_multstart(rate~briere2_1999(temp = temp, tmin, tmax, a,b),
                        data = .x,
                        iter = c(4,4,4,4),
                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'briere2_1999') - 10,
                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'briere2_1999') + 10,
                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'briere2_1999'),
                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'briere2_1999'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)),
         delong = map(data, ~nls_multstart(rate~delong_2017(temp = temp, c, eb, ef, tm, ehc),
                        data = .x,
                        iter = c(4,4,4,4,4),
                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'delong_2017') - 10,
                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'delong_2017') + 10,
                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'delong_2017'),
                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'delong_2017'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)),
         deutsch = map(data, ~nls_multstart(rate~deutsch_2008(temp = temp, rmax, topt, ctmax, a),
                        data = .x,
                        iter = c(4,4,4,4),
                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'deutsch_2008') - 10,
                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'deutsch_2008') + 10,
                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'deutsch_2008'),
                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'deutsch_2008'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)),
         flinn = map(data, ~nls_multstart(rate~flinn_1991(temp = temp, a, b, c),
                        data = .x,
                        iter = c(5,5,5),
                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'flinn_1991') - 10,
                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'flinn_1991') + 10,
                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'flinn_1991'),
                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'flinn_1991'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)),
         gaussian = map(data, ~nls_multstart(rate~gaussian_1987(temp = temp, rmax, topt, a),
                        data = .x,
                        iter = c(4,4,4),
                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'gaussian_1987') - 10,
                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'gaussian_1987') + 10,
                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'gaussian_1987'),
                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'gaussian_1987'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)),
         hinshelwood = map(data, ~nls_multstart(rate~hinshelwood_1947(temp = temp, a, e, b, eh),
                        data = .x,
                        iter = c(5,5,5,5),
                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'hinshelwood_1947') - 1,
                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'hinshelwood_1947') + 1,
                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'hinshelwood_1947'),
                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'hinshelwood_1947'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)),
         joehnk = map(data, ~nls_multstart(rate~joehnk_2008(temp = temp, rmax, topt, a, b, c),
                        data = .x,
                        iter = c(4,4,4,4, 4),
                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'joehnk_2008') - 10,
                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'joehnk_2008') + 10,
                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'joehnk_2008'),
                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'joehnk_2008'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)),
         johnson_lewin = map(data, ~suppressWarnings(nls_multstart(rate~ johnsonlewin_1946(temp = temp, r0, e, eh, topt),
                        data = .x,
                        iter = c(4,4,4,4),
                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'johnsonlewin_1946') - 1,
                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'johnsonlewin_1946') + 1,
                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'johnsonlewin_1946'),
                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'johnsonlewin_1946'),
                        supp_errors = 'Y',
                        convergence_count = FALSE))),
         kamykowski = map(data, ~nls_multstart(rate~kamykowski_1985(temp = temp, tmin, tmax, a, b, c),
                        data = .x,
                        iter = c(4,4,4,4,4),
                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'kamykowski_1985') - 10,
                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'kamykowski_1985') + 10,
                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'kamykowski_1985'),
                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'kamykowski_1985'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)),
         lactin2 = map(data, ~nls_multstart(rate~lactin2_1995(temp = temp, a, b, tmax, delta_t),
                        data = .x,
                        iter = c(4,4,4,4),
                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'lactin2_1995') - 10,
                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'lactin2_1995') + 10,
                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'lactin2_1995'),
                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'lactin2_1995'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)),
         lrf = map(data, ~nls_multstart(rate~lrf_1991(temp = temp, rmax, topt, tmin, tmax),
                  data = d,
                  iter = c(3,3,3,3),
                  start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'lrf_1991') - 10,
                  start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'lrf_1991') + 10,
                  lower = get_lower_lims(.x$temp, .x$rate, model_name = 'lrf_1991'),
                  upper = get_upper_lims(.x$temp, .x$rate, model_name = 'lrf_1991'),
                  supp_errors = 'Y',
                  convergence_count = FALSE)),
         modifiedgaussian = map(data, ~nls_multstart(rate~modifiedgaussian_2006(temp = temp, rmax, topt, a, b),
                        data = .x,
                        iter = c(4,4,4,4),
                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'modifiedgaussian_2006') - 10,
                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'modifiedgaussian_2006') + 10,
                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'modifiedgaussian_2006'),
                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'modifiedgaussian_2006'),
                        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),
                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'oneill_1972') - 10,
                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'oneill_1972') + 10,
                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'oneill_1972'),
                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'oneill_1972'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)),
         pawar = map(data, ~nls_multstart(rate~pawar_2018(temp = temp, r_tref, e, eh, topt, tref = 15),
                        data = .x,
                        iter = c(4,4,4,4),
                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'pawar_2018') - 10,
                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'pawar_2018') + 10,
                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'pawar_2018'),
                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'pawar_2018'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)),
         quadratic = map(data, ~nls_multstart(rate~quadratic_2008(temp = temp, a, b, c),
                        data = .x,
                        iter = c(4,4,4),
                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'quadratic_2008') - 0.5,
                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'quadratic_2008') + 0.5,
                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'quadratic_2008'),
                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'quadratic_2008'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)),
         ratkowsky = map(data, ~nls_multstart(rate~ratkowsky_1983(temp = temp, tmin, tmax, a, b),
                        data = .x,
                        iter = c(4,4,4,4),
                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'ratkowsky_1983') - 10,
                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'ratkowsky_1983') + 10,
                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'ratkowsky_1983'),
                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'ratkowsky_1983'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)),
         rezende = map(data, ~nls_multstart(rate~rezende_2019(temp = temp, q10, a,b,c),
                        data = .x,
                        iter = c(4,4,4,4),
                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'rezende_2019') - 10,
                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'rezende_2019') + 10,
                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'rezende_2019'),
                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'rezende_2019'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)),
         sharpeschoolfull = map(data, ~nls_multstart(rate~sharpeschoolfull_1981(temp = temp, r_tref,e,el,tl,eh,th, tref = 15),
                        data = .x,
                        iter = c(4,4,4,4,4,4),
                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'sharpeschoolfull_1981') - 10,
                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'sharpeschoolfull_1981') + 10,
                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'sharpeschoolfull_1981'),
                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'sharpeschoolfull_1981'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)),
         sharpeschoolhigh = map(data, ~nls_multstart(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 15),
                        data = .x,
                        iter = c(4,4,4,4),
                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'sharpeschoolhigh_1981') - 10,
                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'sharpeschoolhigh_1981') + 10,
                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'sharpeschoolhigh_1981'),
                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'sharpeschoolhigh_1981'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)),
         sharpeschoollow = map(data, ~nls_multstart(rate~sharpeschoollow_1981(temp = temp, r_tref,e,el,tl, tref = 15),
                        data = .x,
                        iter = c(4,4,4,4),
                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'sharpeschoollow_1981') - 10,
                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'sharpeschoollow_1981') + 10,
                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'sharpeschoollow_1981'),
                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'sharpeschoollow_1981'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)),
         spain = map(data, ~nls_multstart(rate~spain_1982(temp = temp, a,b,c,r0),
                        data = .x,
                        iter = c(4,4,4,4),
                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'spain_1982') - 1,
                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'spain_1982') + 1,
                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'spain_1982'),
                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'spain_1982'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)),
         thomas1 = map(data, ~nls_multstart(rate~thomas_2012(temp = temp, a,b,c,topt),
                        data = .x,
                        iter = c(4,4,4,4),
                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'thomas_2012') - 1,
                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'thomas_2012') + 2,
                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'thomas_2012'),
                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'thomas_2012'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)),
         thomas2 = map(data, ~nls_multstart(rate~thomas_2017(temp = temp, a,b,c,d,e),
                        data = .x,
                        iter = c(3,3,3,3,3),
                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'thomas_2017') - 10,
                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'thomas_2017') + 10,
                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'thomas_2017'),
                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'thomas_2017'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)),
         weibull = map(data, ~nls_multstart(rate~weibull_1995(temp = temp, a,topt,b,c),
                        data = .x,
                        iter = c(4,4,4,4),
                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'weibull_1995') - 10,
                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'weibull_1995') + 10,
                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'weibull_1995'),
                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'weibull_1995'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)))

## ----glimpse_fits-------------------------------------------------------------
glimpse(select(d_fits, 1:7))

## ----look_at_beta-------------------------------------------------------------
d_fits$beta[[1]]

## ----preds_and_plot, fig.width = 10, fig.height = 10--------------------------
# stack models
d_stack <- select(d_fits, -data) %>%
  pivot_longer(., names_to = 'model_name', values_to = 'fit', beta:weibull)

# get parameters using tidy
params <- d_stack %>%
  mutate(., est = map(fit, tidy)) %>%
  select(-fit) %>%
  unnest(est)

# get predictions using augment
newdata <- tibble(temp = seq(min(d$temp), max(d$temp), length.out = 100))
d_preds <- d_stack %>%
  mutate(., preds = map(fit, augment, newdata = newdata)) %>%
  select(-fit) %>%
  unnest(preds)

# plot
ggplot(d_preds, aes(temp, rate)) +
  geom_point(aes(temp, rate), d) +
  geom_line(aes(temp, .fitted), col = 'blue') +
  facet_wrap(~model_name, labeller = labeller(model_name = label_facets_num), scales = 'free', ncol = 5) +
  theme_bw(base_size = 12) +
  theme(legend.position = 'none',
        strip.text = element_text(hjust = 0),
        strip.background = element_blank()) +
  labs(x = 'Temperature (ºC)',
       y = 'Metabolic rate',
       title = 'Fits of every model available in rTPC') +
  geom_hline(aes(yintercept = 0), linetype = 2)