Case resampling: Resampling the original data with replacement
Bootstrapping involves 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.
First, we will demonstrate this case resampling approach using data
from a recent paper by Padfield et al. (2020), that measures
the thermal performance of the bacteria, Pseudomonas
fluorescens, in the presence and absence of its phage, \(\phi 2\). In this study, each single growth
rate estimate is a technical replicate, coming from an isogenic strain
of bacteria either inoculated with, or without, the phage. As such, all
the data points within each phage treatment can be used to estimate the
same curve and there is no rep column as in the
chlorella_tpc dataset. We can visualise one of the curves
(bacteria in the absence of phage), using ggplot2.
# load in data
data("bacteria_tpc")
# keep just a single curve
d <- filter(bacteria_tpc, phage == 'nophage')
# show the data
ggplot(d, aes(temp, rate)) +
geom_point(size = 2, alpha = 0.5) +
theme_bw(base_size = 12) +
labs(x = 'Temperature (ºC)',
y = 'Growth rate',
title = 'Growth rate across temperatures')
As in the study, we can fit the Sharpe-Schoolfield model to the data
and plot the predictions using the approaches in
vignette(rTPC) and
vignette(fit_many_models).
# fit Sharpe-Schoolfield model
d_fit <- nest(d, data = c(temp, rate)) %>%
mutate(sharpeschoolhigh = map(data, ~nls_multstart(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)),
# create new temperature data
new_data = map(data, ~tibble(temp = seq(min(.x$temp), max(.x$temp), length.out = 100))),
# predict over that data,
preds = map2(sharpeschoolhigh, new_data, ~augment(.x, newdata = .y)))
# unnest predictions
d_preds <- select(d_fit, preds) %>%
unnest(preds)
# plot data and predictions
ggplot() +
geom_line(aes(temp, .fitted), d_preds, col = 'blue') +
geom_point(aes(temp, rate), d, size = 2, alpha = 0.5) +
theme_bw(base_size = 12) +
labs(x = 'Temperature (ºC)',
y = 'Growth rate',
title = 'Growth rate across temperatures')
Here we have the best fit to the data. If we want confidence bands
around this prediction, we can get those by resampling the data a number
of times. The R package car contains the function
Boot() that provides a wrapper for the widely used
function boot::boot() that is tailored to bootstrapping
regression models.
nls_multstart() is designed to fit models across a wide possible parameter space, but as it samples multiple start parameters for each model, using it with bootstrapping becomes computationally expensive. Instead, we refit the model using minpack.lm::nlsLM(), using the coefficients of nls_multstart() as the start values. The Boot() function then refits the model 999 times and stores the model coefficients.
# refit model using nlsLM
fit_nlsLM <- minpack.lm::nlsLM(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 15),
data = d,
start = coef(d_fit$sharpeschoolhigh[[1]]),
lower = get_lower_lims(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981'),
upper = get_upper_lims(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981'),
weights = rep(1, times = nrow(d)))
# bootstrap using case resampling
boot1 <- Boot(fit_nlsLM, method = 'case')
# look at the data
head(boot1$t)
#> r_tref e eh th
#> [1,] 0.2873690 0.7577894 2.414822 30.52012
#> [2,] 0.2966582 0.7187752 2.415575 30.95682
#> [3,] 0.2950011 0.7072504 2.447862 31.17755
#> [4,] 0.2767081 0.7695552 2.400215 30.63889
#> [5,] 0.2515705 0.8235706 2.412467 30.39302
#> [6,] 0.2267187 0.9857155 2.299607 28.98983The parameters of each bootstrapped refit are returned. All methods that are available for boot() and Boot() are supported for these objects. This includes the hist.boot() to look at the distribution of each parameter.
We can easily create predictions for each of these models and through this confidence intervals around the original fitted predictions. We can then plot (1) the bootstrapped fits and (2) the confidence regions around the model predictions.
# create predictions of each bootstrapped model
boot1_preds <- boot1$t %>%
as.data.frame() %>%
drop_na() %>%
mutate(iter = 1:n()) %>%
group_by_all() %>%
do(data.frame(temp = seq(min(d$temp), max(d$temp), length.out = 100))) %>%
ungroup() %>%
mutate(pred = sharpeschoolhigh_1981(temp, r_tref, e, eh, th, tref = 15))
# calculate bootstrapped confidence intervals
boot1_conf_preds <- group_by(boot1_preds, temp) %>%
summarise(conf_lower = quantile(pred, 0.025),
conf_upper = quantile(pred, 0.975)) %>%
ungroup()
# plot bootstrapped CIs
p1 <- ggplot() +
geom_line(aes(temp, .fitted), d_preds, col = 'blue') +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), boot1_conf_preds, fill = 'blue', alpha = 0.3) +
geom_point(aes(temp, rate), d, size = 2, alpha = 0.5) +
theme_bw(base_size = 12) +
labs(x = 'Temperature (ºC)',
y = 'Growth rate',
title = 'Growth rate across temperatures')
# plot bootstrapped predictions
p2 <- ggplot() +
geom_line(aes(temp, .fitted), d_preds, col = 'blue') +
geom_line(aes(temp, pred, group = iter), boot1_preds, col = 'blue', alpha = 0.007) +
geom_point(aes(temp, rate), d, size = 2, alpha = 0.5) +
theme_bw(base_size = 12) +
labs(x = 'Temperature (ºC)',
y = 'Growth rate',
title = 'Growth rate across temperatures')
p1 + p2
This method works well here, because there are many points beyond the optimum temperature and multiple independent points at each temperature.
This method becomes more problematic when there is a small sample
size and the coverage of temperature values beyond the optimum
temperature is small. This means that many of the bootstrapped models
will not have any points beyond the optimum, which is problematic for
mathematical models that expect a unimodal shape. The effect of this can
be seen by case resampling a curve from the chlorella_tpc
dataset used throughout these vignettes. Here we again fit the model
using nls_multstart(), refit the model using
nlsLM(), then bootstrap the model using
Boot().
# load in chlorella data
data('chlorella_tpc')
d2 <- filter(chlorella_tpc, curve_id == 1)
# fit Sharpe-Schoolfield model to raw data
d_fit <- nest(d2, data = c(temp, rate)) %>%
mutate(sharpeschoolhigh = map(data, ~nls_multstart(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)),
# create new temperature data
new_data = map(data, ~tibble(temp = seq(min(.x$temp), max(.x$temp), length.out = 100))),
# predict over that data,
preds = map2(sharpeschoolhigh, new_data, ~augment(.x, newdata = .y)))
# refit model using nlsLM
fit_nlsLM2 <- nlsLM(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 15),
data = d2,
start = coef(d_fit$sharpeschoolhigh[[1]]),
lower = get_lower_lims(d2$temp, d2$rate, model_name = 'sharpeschoolhigh_1981'),
upper = get_upper_lims(d2$temp, d2$rate, model_name = 'sharpeschoolhigh_1981'),
control = nls.lm.control(maxiter=500),
weights = rep(1, times = nrow(d2)))
# bootstrap using case resampling
boot2 <- Boot(fit_nlsLM2, method = 'case')
#>
#> Number of bootstraps was 990 out of 999 attemptedNumber of bootstraps was 994 out of 999 attempted We can then create predictions for each bootstrapped model and
calculate 95% confidence intervals around the predictions. Models that
don’t fit and return NA for the parameter estimates are
dropped.
# unnest predictions of original model fit
d_preds <- select(d_fit, preds) %>%
unnest(preds)
# predict over new data
boot2_preds <- boot2$t %>%
as.data.frame() %>%
drop_na() %>%
mutate(iter = 1:n()) %>%
group_by_all() %>%
do(data.frame(temp = seq(min(d2$temp), max(d2$temp), length.out = 100))) %>%
ungroup() %>%
mutate(pred = sharpeschoolhigh_1981(temp, r_tref, e, eh, th, tref = 15))
# calculate bootstrapped confidence intervals
boot2_conf_preds <- group_by(boot2_preds, temp) %>%
summarise(conf_lower = quantile(pred, 0.025),
conf_upper = quantile(pred, 0.975)) %>%
ungroup()
# plot bootstrapped CIs
p1 <- ggplot() +
geom_line(aes(temp, .fitted), d_preds, col = 'blue') +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), boot2_conf_preds, fill = 'blue', alpha = 0.3) +
geom_point(aes(temp, rate), d2, size = 2) +
theme_bw(base_size = 12) +
labs(x = 'Temperature (ºC)',
y = 'Growth rate',
title = 'Growth rate across temperatures')
# plot bootstrapped predictions
p2 <- ggplot() +
geom_line(aes(temp, .fitted), d_preds, col = 'blue') +
geom_line(aes(temp, pred, group = iter), boot2_preds, col = 'blue', alpha = 0.007) +
geom_point(aes(temp, rate), d2, size = 2) +
theme_bw(base_size = 12) +
labs(x = 'Temperature (ºC)',
y = 'Growth rate',
title = 'Growth rate across temperatures')
p1 + p2
As can be seen, bootstrapping-with-replacement with only a single point at each temperature can lead to a large variety of fits. In the second panel, we can see the variation of the curve fits, clustering around 4 possible paths for the decrease in rate beyond the optimum temperature. This occurs because in many instances there are no points sampled at the very high temperatures, leading to this clustering in curve fits.
For this curve, the profiling method failed and the asymptotic method
has very wide confidence intervals for 
You can see that the confidence intervals around certain parameters,
such as