Fitting many curves using rTPC

Daniel Padfield

2023-08-16

A brief example of how to fit many models to multiple TPCs using rTPC, nls.multstart and the tidyverse.

# load packages
library(rTPC)
library(nls.multstart)
library(broom)
library(tidyverse)

In the final part of the general pipeline, we demonstrate how multiple models can be fitted to multiple TPCs. Instead of picking all 24 model formulations to demonstrate this approach, we will use only 2 models in this example: gaussian_1987() and sharpeschoolhigh_1981().

We can demonstrate the fitting of multiple curves by modelling all 60 TPCs from the example dataset of rTPC curve from the example dataset rTPC. These TPCs are of respiration and photosynthesis of the aquatic algae, Chlorella vulgaris. These algae differed in their growth temperature, growth_temp, and how long they had been grown at that temperature, process, either for ~100 or ~10 generations.

Using a similar approach to vignette('fit_many_models'), models can be fitted to each curve using list columns and purrr::map() to fit and store multiple models in a data frame.

When fitting lots of models at once, it is useful to know the progress the code as it may take a long time to run.The R package purrr does now has progress bars built in; you can add a progress bar to the map() function using the argument .progress. However, when fitting multiple models with multiple map() calls a new progress bar each time, so currently I have not found a way to have a single progress bar for multiple map() calls, so we will keep with what we have at the moment.

# load in data
data("chlorella_tpc")
d <- chlorella_tpc

# when scaling up our code to fit hundreds of models, its nice to have a progress bar
# edit nls_multstart to allow for a progress bar
nls_multstart_progress <- function(formula, data = parent.frame(), iter, start_lower, 
                                   start_upper, supp_errors = c("Y", "N"), convergence_count = 100, 
                                   control, modelweights, ...){
  if(!is.null(pb)){
    pb$tick()
  }
  nls_multstart(formula = formula, data = data, iter = iter, start_lower = start_lower, 
                start_upper = start_upper, supp_errors = supp_errors, convergence_count = convergence_count, 
                control = control, modelweights = modelweights, ...)
}

# start progress bar and estimate time it will take
number_of_models <- 2
number_of_curves <- length(unique(d$curve_id))

# setup progress bar
pb <- progress::progress_bar$new(total = number_of_curves*number_of_models,
                                 clear = FALSE,
                                 format ="[:bar] :percent :elapsedfull")

# fit two chosen model formulation in rTPC
d_fits <- nest(d, data = c(temp, rate)) %>%
  mutate(gaussian = map(data, ~nls_multstart_progress(rate~gaussian_1987(temp = temp, rmax,topt,a),
                        data = .x,
                        iter = c(3,3,3),
                        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)),
         sharpeschoolhigh = map(data, ~nls_multstart_progress(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 15),
                        data = .x,
                        iter = c(3,3,3,3),
                        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)))
[=========================================================] 100% 00:00:40

Like previous vignettes, the predictions of each model can be estimated using broom::augment(). To do this, we first create a new list column containing high resolution temperature values by taking the min and max of each curve. Next we stack the models and finally we get the new predictions using the map2(), which allows us to use both fit and new_data list columns. After unnesting the preds column, we are then left with high resolution predictions for each curve. As this code covers a lot of steps, each line of the code is commented.

# create new list column of for high resolution data
d_preds <- mutate(d_fits, new_data = map(data, ~tibble(temp = seq(min(.x$temp), max(.x$temp), length.out = 100)))) %>%
  # get rid of original data column
  select(., -data) %>%
  # stack models into a single column, with an id column for model_name
  pivot_longer(., names_to = 'model_name', values_to = 'fit', c(gaussian,sharpeschoolhigh)) %>%
  # create new list column containing the predictions
  # this uses both fit and new_data list columns
  mutate(preds = map2(fit, new_data, ~augment(.x, newdata = .y))) %>%
  # select only the columns we want to keep
  select(curve_id, growth_temp, process, flux, model_name, preds) %>%
  # unlist the preds list column
  unnest(preds)

glimpse(d_preds)
#> Rows: 12,000
#> Columns: 7
#> $ curve_id    <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
#> $ growth_temp <dbl> 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20…
#> $ process     <chr> "acclimation", "acclimation", "acclimation", "acclimation"…
#> $ flux        <chr> "respiration", "respiration", "respiration", "respiration"…
#> $ model_name  <chr> "gaussian", "gaussian", "gaussian", "gaussian", "gaussian"…
#> $ temp        <dbl> 16.00000, 16.33333, 16.66667, 17.00000, 17.33333, 17.66667…
#> $ .fitted     <dbl> 0.02790076, 0.03175769, 0.03607053, 0.04088150, 0.04623512…

We can then plot the predictions of each curve using ggplot2.

# plot
ggplot(d_preds) +
  geom_line(aes(temp, .fitted, col = model_name)) +
  geom_point(aes(temp, rate), d) +
  facet_wrap(~curve_id, scales = 'free_y', ncol = 6) +
  theme_bw() +
  theme(legend.position = 'none') +
  scale_color_brewer(type = 'qual', palette = 2) +
  labs(x = 'Temperature (ºC)',
       y = 'Metabolic rate',
       title = 'All fitted thermal performance curves',
       subtitle = 'gaussian in green; sharpeschoolhigh in orange')

The traits of each thermal performance curve can also easily be calculated.

# stack models and calculate extra params
d_params <- pivot_longer(d_fits, names_to = 'model_name', values_to = 'fit', c(gaussian,sharpeschoolhigh)) %>%
  mutate(params = map(fit, calc_params)) %>%
  select(curve_id, growth_temp, process, flux, model_name, params) %>%
  unnest(params)

glimpse(d_params)
#> Rows: 120
#> Columns: 16
#> $ curve_id              <dbl> 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, …
#> $ growth_temp           <dbl> 20, 20, 20, 20, 23, 23, 27, 27, 27, 27, 30, 30, …
#> $ process               <chr> "acclimation", "acclimation", "acclimation", "ac…
#> $ flux                  <chr> "respiration", "respiration", "respiration", "re…
#> $ model_name            <chr> "gaussian", "sharpeschoolhigh", "gaussian", "sha…
#> $ rmax                  <dbl> 1.4972473, 1.8127063, 1.7935416, 1.9798352, 0.94…
#> $ topt                  <dbl> 36.338, 41.646, 35.554, 39.007, 34.615, 35.170, …
#> $ ctmin                 <dbl> 18.695, 2.539, 18.835, 6.086, 17.368, 11.970, 15…
#> $ ctmax                 <dbl> 53.981, 45.558, 52.272, 47.680, 51.861, 53.064, …
#> $ e                     <dbl> 0.6758954, 0.5802496, 0.8427748, 0.6621701, 1.02…
#> $ eh                    <dbl> 1.0767825, 11.4840311, 1.0954099, 2.4774826, 1.1…
#> $ q10                   <dbl> 2.373766, 2.063147, 2.927011, 2.305129, 3.722119…
#> $ thermal_safety_margin <dbl> 17.643, 3.912, 16.718, 8.673, 17.246, 17.894, 18…
#> $ thermal_tolerance     <dbl> 35.286, 43.019, 33.437, 41.594, 34.493, 41.094, …
#> $ breadth               <dbl> 9.651, 5.371, 9.165, 6.495, 9.442, 7.714, 10.354…
#> $ skewness              <dbl> -4.008871e-01, -1.090378e+01, -2.526351e-01, -1.…