Bootstrap non-linear regression with purrr and modelr

This post was updated on the 02/03/2018 to reflect changes to nls.multstart.

Introduction

For my first academic publication, a reviewer asked for the \(r^{2}\) values of the thermal performance curves I fitted using non-linear regression. I bowed to the request as is often the case with reviewer comments, but would now resist as the \(r^{2}\) is not necessarily an effective goodness of fit measure for non-linear regression (see this SO answer). It does raise the question of how to determine how well a biologically meaningful model fits the data it is fitted to. I generally just plot every curve to its data, but it tells me nothing of the uncertainty around the curve.

Step forward the bootstrap! Bootstrapping involes simulating “new” datasets produced from the existing data by sampling with replacement. The same model is then fitted separately on each individual bootstrapped dataset. Doing this over and over allows us to visualise uncertainty of predictions and produce confidence intervals of estimated parameters. This blog post was inspired by posts by Andrew MacDonald and Hadley Wickham, as well as a broom vignette which use this approach. I have taken their approaches and again applied them to thermal performance curves. The broom approach in these blog posts has since been replaced by modelr::bootstrap(), another package of the tidyverse.

Bootstrapping predictions for a single curve

I will demonstrate this approach by using the thermal performance curves for phytoplankton metabolism that I used in a previous post. The Sharpe-Schoolfield equation and meaning of the parameters can be found in more detail in the previous post.

Firstly lets load the packages used in the script and write the Sharpe-Schoolfield model as a function.

# load packages
library(nls.multstart) # devtools::install_github('padpadpadpad/nls.multstart')
library(patchwork) # devtools::install_github('thomasp85/patchwork')
library(ggplot2)
library(broom)
library(purrr)
library(dplyr)
library(tidyr)
library(nlstools)
library(modelr)

# write function for sharpe schoolfield model
schoolfield_high <- function(lnc, E, Eh, Th, temp, Tc) {
  Tc <- 273.15 + Tc
  k <- 8.62e-5
  boltzmann.term <- lnc + log(exp(E/k*(1/Tc - 1/temp)))
  inactivation.term <- log(1/(1 + exp(Eh/k*(1/Th - 1/temp))))
  return(boltzmann.term + inactivation.term)
}

Then load in the data and have a look at the its structure using glimpse().

# load in data
data(Chlorella_TRC)

# look at data
glimpse(Chlorella_TRC)
## Observations: 649
## Variables: 7
## $ curve_id    <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2,...
## $ growth.temp <dbl> 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20...
## $ process     <chr> "acclimation", "acclimation", "acclimation", "accl...
## $ flux        <chr> "respiration", "respiration", "respiration", "resp...
## $ temp        <dbl> 16, 19, 22, 25, 28, 31, 34, 37, 40, 43, 46, 49, 16...
## $ K           <dbl> 289.15, 292.15, 295.15, 298.15, 301.15, 304.15, 30...
## $ ln.rate     <dbl> -2.06257833, -1.32437939, -0.95416807, -0.79443675...

There are 60 curves here, 30 each for photosynthesis and respiration. The treatments are growth temperature (20, 23, 27, 30, 33 ºC) and adaptive process (acclimation or adaptation) that reflects the number of generations cultures were grown at each temperature. Bootstrapping the uncertainty on each individual curve is difficult for thermal performance curves because rates generally rapidly decrease after the optimum temperature, \(T_{opt}\), making data collection difficult.

This means it is likely some of the bootstraps will not include points after the optimum, making the unimodal model formulation unsuitable. Because of this I will pool the replicates curves within treatments together to demonstrate the bootstrapping approach, giving 20 curves in total. This post therefore ignores the non-independence of data points within replicates (a little naughty!).

To bootstrap a single curve, we can filter the dataset for a single flux (photosynthesis) at a singe growth temperature (20 ºC) and generations of growth (~ 100) and plot the data.

# filter for one curve
d_examp <- filter(Chlorella_TRC, growth.temp == 20, flux == 'photosynthesis', process == 'adaptation')

# plot 
ggplot(d_examp, aes(K - 273.15, ln.rate)) +
  geom_point(col = 'green4') +
  ylab('log Metabolic rate') +
  xlab('Assay temperature (ºC)') +
  theme_bw(base_size = 12, base_family = 'Helvetica')

We can use nls_multstart(), that allows for multiple start parameters, to fit a single model to the data. We can then use tidy() and augment() from broom to get the parameters and predictions of the model.

# run nls_multstart
fit <- nls_multstart(ln.rate ~ schoolfield_high(lnc, E, Eh, Th, temp = K, Tc = 20),
                     data = d_examp,
                     iter = 500,
                     start_lower = c(lnc = -10, E = 0.1, Eh = 0.2, Th = 285),
                     start_upper = c(lnc = 10, E = 2, Eh = 5, Th = 330),
                     supp_errors = 'Y',
                     na.action = na.omit,
                     lower = c(lnc = -10, E = 0, Eh = 0, Th = 0))

# broom functions to tidy up model
params <- tidy(fit)
preds <- augment(fit)

# plot with predictions
ggplot(d_examp, aes(K - 273.15, ln.rate)) +
  geom_point(col = 'green4') +
  geom_line(aes(K - 273.15, .fitted), preds) +
  ylab('log Metabolic rate') +
  xlab('Assay temperature (ºC)') +
  theme_bw(base_size = 12, base_family = 'Helvetica') +
  ggtitle('Single TPC with fitted model')

This fit provides p-values and confidence intervals can be calculated using nlstools::confint2(). However, bootstrapping can provide confidence intervals around predictions and for estimated parameters.

The bootstrap() function in modelr samples bootstrap replicates (here we do 200), each of which is randomly sampled with replacement. This creates a list column in our tibble called strap which contains the bootsrapped dataset, and a new column called boot_num that is the number of that bootstrap (from 1 to 200).

We can then create a new list column of the fit for each strap using purrr::map().

fit_boots <- d_examp %>% 
  modelr::bootstrap(n = 200, id = 'boot_num') %>%
  group_by(boot_num) %>%
  mutate(fit = map(strap, ~nls_multstart(ln.rate ~ schoolfield_high(lnc, E, Eh, Th, temp = K, Tc = 20),
                        data = data.frame(.),
                        iter = 100,
                        start_lower = c(lnc = -10, E = 0.1, Eh = 0.2, Th = 285),
                        start_upper = c(lnc = 10, E = 2, Eh = 5, Th = 330),
                        lower = c(lnc=-10, E=0, Eh=0, Th=0),
                        supp_errors = 'Y')
  ))

fit_boots
## # A tibble: 200 x 3
## # Groups:   boot_num [200]
##    strap          boot_num fit      
##    <list>         <chr>    <list>   
##  1 <S3: resample> 001      <S3: nls>
##  2 <S3: resample> 002      <S3: nls>
##  3 <S3: resample> 003      <S3: nls>
##  4 <S3: resample> 004      <S3: nls>
##  5 <S3: resample> 005      <S3: nls>
##  6 <S3: resample> 006      <S3: nls>
##  7 <S3: resample> 007      <S3: nls>
##  8 <S3: resample> 008      <S3: nls>
##  9 <S3: resample> 009      <S3: nls>
## 10 <S3: resample> 010      <S3: nls>
## # ... with 190 more rows

Each bootstrap replicate is stored in a list column within the tibble . This then allows us to apply the tidy() and augment() functions used earlier, using unnest() to combine the list column into a dataframe.

# get parameters ####
params_boot <- fit_boots %>%
  unnest(fit %>% map(tidy)) %>%
  ungroup()

# get predictions
preds_boot <- fit_boots %>%
  unnest(fit %>% map(augment)) %>%
  ungroup()

Using these two dataframes, we can plot each set of bootstrapped predictions alongside the fit of the original data, preds, and plot the distribution of each estimated parameter.

The relatively new package patchwork by Thomas Lin Pedersen can help add multiple graphs together simply by saying plot_1 + plot_2.

# plot distribution of estimated parameters
p1 <- ggplot(params_boot, aes(estimate)) +
  geom_histogram(col = 'black', fill = 'white') +
  facet_wrap(~ term, scales = 'free_x')

# plot points with predictions
p2 <- ggplot() +
  geom_line(aes(K - 273.15, .fitted, group = boot_num), preds_boot, alpha = .03) +
  geom_line(aes(K - 273.15, .fitted), preds) +
  geom_point(aes(K - 273.15, ln.rate), d_examp, col = 'green4') +
  ylab('log Metabolic rate') +
  xlab('Assay temperature (ºC)') +
  theme_bw(base_size = 12, base_family = 'Helvetica')
  
# plot both
p1 + p2

We can smooth our predictions over smaller increments of our predictor variable by passing a new dataset to augment(). Alongside this, for every value of the predictor we can calculate the 2.5% and 97.5% quantiles which gives confidence bands around the predictions.

Personally I prefer this approach rather than plotting each bootstrapped replicate.

# new data frame of predictions
new_preds <- d_examp %>%
  do(., data.frame(K = seq(min(.$K), max(.$K), length.out = 250), stringsAsFactors = FALSE))

# create smooth predictions for best fit
preds <- augment(fit, newdata = new_preds)

# create smoother predictions for bootstrapped replicate
preds <- fit_boots %>%
  unnest(fit %>% map(augment, newdata = new_preds)) %>%
  # group by each value of K and get quantiles
  group_by(., K) %>%
  summarise(lwr_CI = quantile(.fitted, 0.025),
            upr_CI = quantile(.fitted, 0.975)) %>%
  ungroup() %>%
  merge(., preds, by = 'K')

# plot
ggplot() +
  geom_point(aes(K - 273.15, ln.rate), d_examp) +
  geom_line(aes(K - 273.15, .fitted), preds) +
  geom_ribbon(aes(K - 273.15, ymin = lwr_CI, ymax = upr_CI), fill = 'green4', preds, alpha = .2) +
  ylab('log Metabolic rate') +
  xlab('Assay temperature (ºC)') +
  theme_bw(base_size = 12, base_family = 'Helvetica') +
  ggtitle('Single TPC with confidence intervals')

This is starting to look pretty nice, and is a great approach to visualising uncertainty of non-linear regressions for many types of data.

From params_boot we can calculate confidence intervals of each estimated parameter by taking the desired quantiles of the data. This can be compared the output from confint2() from nlstools.

# calculate confidence intervals of estimated parameters using confint2()
confint_1 <- confint2(fit) %>%
  data.frame() %>%
  rename(., conf_low = X2.5.., conf_high = X97.5..) %>%
  mutate(method = 'nlstools') %>%
  cbind(., select(params, term, estimate)) %>%
  select(., term, estimate, conf_low, conf_high, method)

# calculate confidence intervals using bootstraps
confint_2 <- group_by(params_boot, term) %>%
  summarise(.,
            conf_low = quantile(estimate, 0.025),
            conf_high = quantile(estimate, 0.975),
            estimate = quantile(estimate, 0.5)) %>%
  ungroup() %>%
  mutate(., method = 'boot')

# bind both methods
confint <- bind_rows(confint_1, confint_2)

# plot each method side by side
ggplot(confint, aes(method, estimate, col = method)) +
  geom_point(size = 3) +
  geom_linerange(aes(ymin = conf_low, ymax = conf_high)) +
  facet_wrap(~ term, scales = 'free_y') +
  theme_bw() +
  theme(legend.position = 'none') +
  ggtitle('Comparison of confidence interval calculation for estimated parameters')

Bootstrapping gives similar mean estimates, but gives wider, asymmetric, confidence intervals compared to those calculated using nlstools::confint2(). In this instance this could be because I am only running 200 bootstrap replicates, whereas the number of bootstraps done in published analyses are commonly around 10,000.

Crucially, bootstrapping allows the calculation of confidence intervals for parameters derived from the model that were not present in the initial fitting process. For example, the optimum temperature of a thermal performance curve, \(T_{opt}\) is calculated as:

\[T_{opt} = \frac{E_{h}T_{h}}{E_{h} + k T_{h} ln(\frac{E_{h}}{E} - 1)}\] We can calculate \(T_{opt}\) for each iteration of the bootstrap and then plot the distribution of derived parameters. This can be done by using tidyr::spread() to have a column for each estimated parameter, and then calculate Topt for each boot_num.

# function for calculating Topt
Topt <- function(E, Th, Eh){
  return((Eh*Th)/(Eh + (8.62e-05 *Th*log((Eh/E) - 1))))
}

Topt_boot <- select(params_boot, boot_num, term, estimate) %>%
  spread(., term, estimate) %>%
  mutate(., Topt = Topt(E, Th, Eh))

# plot distribution of Topt
ggplot(Topt_boot, aes(Topt - 273.15)) +
  geom_histogram(col = 'black', fill = 'white') +
  xlab('Optimum Temperature (ºC)') +
  ggtitle('Distribution of optimum temperatures via bootstrapping')

Bootstrapping allows us to get uncertainty estimates for parameters outside of the original curve fitting process!

Bootstrapping multiple curves

Bootstrapping over each curve can be done by combining functions from the tidyverse to the bootstrap() call. To fit a single model to each curve, I use nest(), mutate() and map() as shown previously. I searched for a way of using the same workflow for bootstrapping, and finally came across the answer. Each element of strap is not strictly a dataframe (more of a promise to be a dataframe), so the only difference to fitting multiple non-linear regressions is the need to specify the data using dataframe(.) within map().

After grouping the dataframe, new datasets are bootstrapped using modelr::bootstrap(). Using unnest() gives a new column called boot_num which represents the bootstrap replicate within each group. The tibble is then re-grouped to include boot_num and the model can finally be fitted to each bootstrapped dataset. Obviously the total number of models your code fits is increases as you up the number of bootstraps, so be aware this code may take a fair while to run!

# run nls.multstart on each curve of the original data ####
fit_many <- group_by(Chlorella_TRC, growth.temp, process, flux) %>%
  nest() %>%
  mutate(., fit = purrr::map(data, ~nls_multstart(ln.rate ~ schoolfield_high(lnc, E, Eh, Th, temp = K, Tc = 20),
                                   data = .x,
                                   iter = 500,
                                   start_lower = c(lnc = -10, E = 0.1, Eh = 0.2, Th = 285),
                                   start_upper = c(lnc = 10, E = 2, Eh = 5, Th = 330),
                                   supp_errors = 'Y',
                                   na.action = na.omit,
                                   lower = c(lnc = -10, E = 0, Eh = 0, Th = 0))))

# run bootstrap over many curves ####
boot_many <- group_by(Chlorella_TRC, growth.temp, process, flux) %>%
  # create 200 bootstrap replicates per curve
  do(., boot = modelr::bootstrap(., n = 200, id = 'boot_num')) %>%
  # unnest to show bootstrap number, .id
  unnest() %>%
  # regroup to include the boot_num
  group_by(., growth.temp, process, flux, boot_num) %>%
  # run the model using map()
  mutate(fit = map(strap, ~nls_multstart(ln.rate ~ schoolfield_high(lnc, E, Eh, Th, temp = K, Tc = 20),
                        data = data.frame(.),
                        iter = 50,
                        start_lower = c(lnc = -10, E = 0.1, Eh = 0.2, Th = 285),
                        start_upper = c(lnc = 10, E = 2, Eh = 5, Th = 330),
                        lower = c(lnc=-10, E=0, Eh=0, Th=0),
                        supp_errors = 'Y')
  ))

Smooth predictions can then be calculated from fit_many and boot_many and plotted to demonstrate the uncertainty of multiple curves. I do some wrangling to get the max and min temperature of each curve so that I only plot predictions over the range of each specific curve (some have measurements up to 46 ºC, some to 49 ºC).

# new data frame for smooth predictions
new_preds <- Chlorella_TRC %>%
  do(., data.frame(K = seq(min(.$K), max(.$K), length.out = 150), stringsAsFactors = FALSE))

# get max and min for each curve
max_min <- group_by(Chlorella_TRC, growth.temp, flux, process) %>%
  summarise(., min_K = min(K), max_K = max(K)) %>%
  ungroup()

# create smoother predictions for unbootstrapped models
preds_many <- fit_many %>%
  unnest(fit %>% map(augment, newdata = new_preds))

# create smoother predictions for bootstrapped replicates
preds_many <- boot_many %>%
  unnest(fit %>% map(augment, newdata = new_preds)) %>%
  ungroup() %>%
  # group by each value of K and get quantiles
  group_by(., growth.temp, process, flux, K) %>%
  summarise(lwr_CI = quantile(.fitted, 0.025),
            upr_CI = quantile(.fitted, 0.975)) %>%
  ungroup() %>%
  merge(., preds_many, by = c('K', 'growth.temp', 'flux', 'process')) %>%
  # merge with max_min to delete predictions outside of the max and min temperatures of each curve
  merge(., max_min, by = c('growth.temp', 'flux', 'process')) %>%
  group_by(., growth.temp, flux, process) %>%
  filter(., K >= unique(min_K) & K <= unique(max_K)) %>%
  rename(., ln.rate = .fitted) %>%
  ungroup()

# plot predictions 
ggplot(Chlorella_TRC, aes(K - 273.15, ln.rate, group = flux)) +
  geom_point(alpha = 0.5, size = 0.5) +
  geom_line(data = preds_many) +
  geom_ribbon(aes(ymin = lwr_CI, ymax = upr_CI, fill = flux), data = preds_many, alpha = .2) + 
  scale_fill_manual(values = c('green4', 'black')) +
  theme_bw(base_size = 12, base_family = 'Helvetica') +
  ylab('log Metabolic rate') +
  xlab('Assay temperature (ºC)') +
  facet_wrap(~ process + growth.temp, labeller = labeller(.multi_line = FALSE)) +
  theme_bw(base_size = 12, base_family = 'Helvetica') +
  theme(legend.position = c(0.9, 0.15)) +
  ggtitle('Multiple TPCs with confidence intervals')

When in doubt, it seems that purrr has the answer. Bootstrapping is a super useful method for visualising the uncertainty of predictions for non-linear regressions, and better allow us to understand how well a particular model fits the data. The tidyverse, as usual, provides a set of tools that makes this method easy to understand and implement, sort of.

Next steps

Another way of visualising uncertainty in predictions is by using a Bayesian approach. When using Stan, the generated quantities{} block is a great way to create new predictions from the data. This allows for models that can account for the structure of the data (i.e. temperatures within replicates) and visualise the uncertainty of the model fits. The package brms looks to be a great help in fitting non-linear mixed models in a Bayesian framework. Just need to get round to experimenting with it!

comments powered by Disqus