Bootstrap fits to estimate uncertainty

Code

library(tidyverse)
library(targets)
library(GGally)
library(kableExtra)

# load "R/*" scripts and saved R objects from the targets pipeline
tar_source()
tar_load(c(exp1_data, exp2_data, exp1_data_agg, exp2_data_agg, exp3_data_agg))
set.seed(213)

# make sure that if the model fails to converge, it doesn't stop the whole process
safe_boot_est <- purrr::safely(boot_est)

Bootstrap the original model

Code

res_asy <- run_or_load(
  t(replicate(1000, safe_boot_est(exp1_data)$result))[1, , ],
  "output/res_boot1000_asy.rds")
res_asy_exp2 <- run_or_load(
  t(replicate(1000, safe_boot_est(exp2_data)$result))[1, , ],
  "output/res_boot1000_asy_exp2.rds"
  )
res_lin <- run_or_load(
  t(replicate(1000, safe_boot_est(exp1_data, model = "linear")$result))[1, , ],
  "output/res_boot1000_lin.rds"
  )
res_lin_exp2 <- run_or_load(
  t(replicate(1000, safe_boot_est(exp2_data, model = "linear")$result))[1, , ],
  "output/res_boot1000_lin_exp2.rds"
  )

plot the results

Code

plot_bootstrap_results(res_asy)

Code

plot_bootstrap_results(res_asy_exp2)

Code

plot_bootstrap_results(res_lin)

Code

plot_bootstrap_results(res_lin_exp2)

Code

library(HDInterval)
hdi(unlist(res_asy[,'rate']))

     lower      upper 
0.05588434 0.21303652 
attr(,"credMass")
[1] 0.95

Code

hdi(unlist(res_asy_exp2[,'rate']))

     lower      upper 
0.05242345 0.21291100 
attr(,"credMass")
[1] 0.95

Code

hdi(unlist(res_lin[,'rate']))

      lower       upper 
0.005521282 0.031903235 
attr(,"credMass")
[1] 0.95

Code

hdi(unlist(res_lin_exp2[,'rate']))

      lower       upper 
0.003978507 0.061735364 
attr(,"credMass")
[1] 0.95