## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set(echo = TRUE) library(hestia) library(dplyr) library(posterior) library(ggplot2) ## ----------------------------------------------------------------------------- head(sir) ## ----------------------------------------------------------------------------- # Basic SIR, fit recovery rate inf_process <- make_infection_model(transmit(from = "S", to = "I"), progress(from = "I", to = "R", gamma = NA)) ## ----------------------------------------------------------------------------- obs_process <- make_observation_model( pcr = c("S" = 0.05, "I" = 0.95, "R" = 0.05), igg = c("S" = 0.01, "I" = 0.01, "R" = 0.8) ) ## ----message=FALSE, warning=FALSE, eval=FALSE--------------------------------- # # # Note this is very computationally intensive so the block is not evaluated # # when knitting # # ## Inputs # # 1. Infection process model # # 2. Observation process model # # 3. Raw data # # 4. Covariates (optional) # # 5. Starting state probabilities # # sir_res <- run_model( # inf_model = inf_process, # obs_model = obs_process, # data = sir, # init_probs = c(1 - 2 * 1e-10, 1e-10, 1e-10), # stan_opts = stan_options(iter = 1000) # ) # ## ----------------------------------------------------------------------------- # The output from the model in the previous code chunk is available in the # sir_res package data object # Get summary statistics. Parameters are returned on the natural (model) scale, # so no back-transformation is needed. draws_sum <- summarise_draws(sir_res, "mean", "median", ~ quantile2(.x, probs = c(0.025, 0.975))) |> mutate(var_type = c(rep("Infection probability", 2), "Recovery rate")) # "True" values from underlying simulation draws_sum$true_value <- c(0.01, 0.05, 0.2) draws_sum ## ----fig.width=7, fig.height=4.5---------------------------------------------- draws_plot <- draws_sum |> mutate( label = factor( c("Extra-household infection prob (eh_prob)", "Intra-household infection prob (ih_prob)", "Recovery rate (gamma)"), levels = c("Extra-household infection prob (eh_prob)", "Intra-household infection prob (ih_prob)", "Recovery rate (gamma)") ) ) ggplot(draws_plot, aes(y = label)) + geom_errorbar(aes(xmin = q2.5, xmax = q97.5), width = 0, color = "#2b5c8f", linewidth = 1) + geom_point(aes(x = median, fill = "Posterior median"), size = 3, shape = 21, color = "#0f2540", stroke = 1.4) + geom_point(aes(x = true_value, fill = "True simulation value"), size = 3, shape = 23, color = "#a63603", stroke = 1.4) + facet_wrap(~var_type, scales = "free", ncol = 1) + scale_fill_manual( name = "", values = c("Posterior median" = "#ffffff", "True simulation value" = "#d95f02"), guide = guide_legend(override.aes = list( shape = c(21, 23), color = c("#0f2540", "#a63603"), fill = c("#ffffff", "#d95f02") )) ) + labs( title = "Posterior Parameter Estimates (SIR Model)", subtitle = "Circle = posterior median; bar = 95% credible interval; diamond = true simulation value", x = "Parameter Value (Natural Scale)", y = NULL ) + theme_minimal() + theme( legend.position = "top", legend.justification = "right", panel.grid.minor = element_blank(), panel.grid.major.y = element_line(color = "grey90"), plot.title = element_text(face = "bold"), strip.text = element_text(face = "bold", hjust = 0) ) ## ----message=FALSE, warning=FALSE, eval=FALSE--------------------------------- # # Simulated data from a basic SIR model for 500 households # head(sir_cov$observations) # # # Run model with x1 and x2 as covariates # sir_cov_res <- run_model( # inf_model = inf_process, # obs_model = obs_process, # data = sir_cov$observations, # ih_cov = sir_cov$covariates, # eh_cov = sir_cov$covariates, # init_probs = c(1 - 2 * 1e-10, 1e-10, 1e-10), # stan_opts = stan_options(iter = 1000) # ) # ## ----------------------------------------------------------------------------- # The output from the model in the previous code chunk is available in the # sir_cov_res package data object # Get summary statistics. Infection probabilities and rates are returned on the # natural (model) scale and coefficients on the natural (exponentiated) scale, # so no back-transformation is needed. draws_sum <- summarise_draws(sir_cov_res, "mean", "median", ~ quantile2(.x, probs = c(0.025, 0.975))) |> mutate(var_type = c(rep("Infection probability", 2), "Recovery rate", rep("Coefficient", 4))) # "True" values from underlying simulation draws_sum$true_vals <- c(0.01, 0.05, 0.2, exp(-0.4), exp(0.7), exp(0.8), exp(0.1)) draws_sum ## ----fig.width=7, fig.height=5.5---------------------------------------------- draws_cov_plot <- draws_sum |> mutate( category = factor( c("Baseline Probabilities", "Baseline Probabilities", "Recovery Rate", rep("Covariate Effects (Natural / Exponentiated Scale)", 4)), levels = c("Baseline Probabilities", "Recovery Rate", "Covariate Effects (Natural / Exponentiated Scale)") ), label = factor( c("Extra-household baseline (eh_prob)", "Intra-household baseline (ih_prob)", "Recovery rate (gamma)", "x1 on extra-household (x1_eh)", "x2 on extra-household (x2_eh)", "x1 on intra-household (x1_ih)", "x2 on intra-household (x2_ih)"), levels = rev(c("Extra-household baseline (eh_prob)", "Intra-household baseline (ih_prob)", "Recovery rate (gamma)", "x1 on extra-household (x1_eh)", "x2 on extra-household (x2_eh)", "x1 on intra-household (x1_ih)", "x2 on intra-household (x2_ih)")) ) ) ggplot(draws_cov_plot, aes(y = label)) + geom_vline(data = filter(draws_cov_plot, grepl("Covariate", category)), aes(xintercept = 1), linetype = "dashed", color = "grey60") + geom_errorbar(aes(xmin = q2.5, xmax = q97.5), width = 0, color = "#1b9e77", linewidth = 1) + geom_point(aes(x = median, fill = "Posterior median"), size = 3, shape = 21, color = "#0d664c", stroke = 1.4) + geom_point(aes(x = true_vals, fill = "True simulation value"), size = 3, shape = 23, color = "#a63603", stroke = 1.4) + facet_wrap(~category, scales = "free", ncol = 1) + scale_fill_manual( name = "", values = c("Posterior median" = "#ffffff", "True simulation value" = "#d95f02"), guide = guide_legend(override.aes = list( shape = c(21, 23), color = c("#0d664c", "#a63603"), fill = c("#ffffff", "#d95f02") )) ) + labs( title = "Posterior Estimates with Covariates (SIR Model)", subtitle = "Circle = posterior median; bar = 95% CrI; diamond = true simulation value", x = "Estimate", y = NULL ) + theme_minimal() + theme( legend.position = "top", legend.justification = "right", panel.grid.minor = element_blank(), panel.grid.major.y = element_line(color = "grey90"), plot.title = element_text(face = "bold"), strip.text = element_text(face = "bold", hjust = 0) )