Quick bootstrap of a nls model using dplyr and broom

nls_bootstrap.R

# quick bootstrap of a TPC model

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

# write function for schoolfield high
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)
}

# load in data
data(Chlorella_TRC)

# filter for one curve
d_1 <- filter(Chlorella_TRC, curve_id == 1)
d_8 <- filter(Chlorella_TRC, curve_id == 8)
d_32 <- filter(Chlorella_TRC, curve_id == 32)

d_boots <- d_32 %>% bootstrap(m = 500) %>%
do(model1 = nls_multstart(ln.rate ~ schoolfield_high(lnc, E, Eh, Th, temp = K, Tc = 20),
data = .,
iter = 10,
param_bds = c(-10, 10, 0.1, 2, 0.5, 5, 285, 330),
lower = c(lnc=-10, E=0, Eh=0, Th=0),
upper = c(lnc = 5, E = 10, Eh = 30, Th = 350),
supp_errors = 'Y')
)


# get parameters ####
d_param_boot <- tidy(d_boots, model1) %>%
ungroup()

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

# get predictions ####
d_pred_boot <- augment(d_boots, model1) %>%
ungroup()

# get mean and upper and lower quantiles of each prediction ####
d_mean_boot <- group_by(d_pred_boot, K) %>%
summarise(., mu = mean(.fitted),
lwr_CI = quantile(.fitted, 0.025),
upr_CI = quantile(.fitted, 0.975)) %>%
ungroup()

# plot
p2 <- ggplot() +
geom_point(aes(K, ln.rate), d_pred_boot, alpha = .01) +
geom_line(aes(K, .fitted, group = factor(replicate)), d_pred_boot, alpha = .01) +
geom_line(aes(K, mu), d_mean_boot) +
geom_line(aes(K, lwr_CI), d_mean_boot, linetype = 2) +
geom_line(aes(K, upr_CI), d_mean_boot, linetype = 2)

# plot both
p3 <- p1 + p2
#ggsave('~/Desktop/boot_plot.pdf', p3, height = 6, width = 12)