Under construction
In this vignette, we’ll demonstrate how to use primarycensoreddist in conjunction with Stan for Bayesian inference of epidemiological delay distributions. We’ll cover the following key points:
primarycensoreddist
This vignette builds on the concepts introduced in the Getting Started with primarycensoreddist vignette and assumes familiarity with using Stan tools as covered in the How to use primarycensoreddist with Stan vignette.
Alongside the primarycensoreddist package we will use the cmdstanr package for interfacing with cmdstan. We will also use the ggplot2 package for plotting and dplyr for data manipulation.
cmdstanr
ggplot2
dplyr
+library(primarycensoreddist) +library(cmdstanr) +library(ggplot2) +library(dplyr)
library(primarycensoreddist) +library(cmdstanr) +library(ggplot2) +library(dplyr)
We’ll start by simulating some censored and truncated delay distribution data. We’ll use the rprimarycensoreddist function (actually we will use the rpcens alias for brevity).
rprimarycensoreddist
rpcens
+set.seed(123) # For reproducibility + +# Define the true distribution parameters +n <- 1000 +meanlog <- 1.5 +sdlog <- 0.75 + +# Generate varying pwindow, swindow, and obs_time lengths +pwindows <- sample(1, n, replace = TRUE) +swindows <- sample(1, n, replace = TRUE) +obs_times <- 10 + +# Function to generate a single sample +generate_sample <- function(pwindow, swindow, obs_time) { + rpcens( + 1, rlnorm, + meanlog = meanlog, sdlog = sdlog, + pwindow = pwindow, swindow = swindow, D = obs_time + ) +} + +# Generate samples +samples <- mapply(generate_sample, pwindows, swindows, obs_times) + +# Create initial data frame +delay_data <- data.frame( + pwindow = pwindows, + swindow = swindows, + obs_time = obs_times, + observed_delay = samples +) + +head(delay_data)
set.seed(123) # For reproducibility + +# Define the true distribution parameters +n <- 1000 +meanlog <- 1.5 +sdlog <- 0.75 + +# Generate varying pwindow, swindow, and obs_time lengths +pwindows <- sample(1, n, replace = TRUE) +swindows <- sample(1, n, replace = TRUE) +obs_times <- 10 + +# Function to generate a single sample +generate_sample <- function(pwindow, swindow, obs_time) { + rpcens( + 1, rlnorm, + meanlog = meanlog, sdlog = sdlog, + pwindow = pwindow, swindow = swindow, D = obs_time + ) +} + +# Generate samples +samples <- mapply(generate_sample, pwindows, swindows, obs_times) + +# Create initial data frame +delay_data <- data.frame( + pwindow = pwindows, + swindow = swindows, + obs_time = obs_times, + observed_delay = samples +) + +head(delay_data)
## pwindow swindow obs_time observed_delay +## 1 1 1 10 2 +## 2 1 1 10 1 +## 3 1 1 10 8 +## 4 1 1 10 6 +## 5 1 1 10 5 +## 6 1 1 10 4
+# Aggregate to unique combinations and count occurrences +# Aggregate to unique combinations and count occurrences +delay_counts <- delay_data |> + summarise(n = n(), .by = c(pwindow, swindow, obs_time, observed_delay)) + +head(delay_counts)
# Aggregate to unique combinations and count occurrences +# Aggregate to unique combinations and count occurrences +delay_counts <- delay_data |> + summarise(n = n(), .by = c(pwindow, swindow, obs_time, observed_delay)) + +head(delay_counts)
## pwindow swindow obs_time observed_delay n +## 1 1 1 10 2 160 +## 2 1 1 10 1 96 +## 3 1 1 10 8 59 +## 4 1 1 10 6 113 +## 5 1 1 10 5 119 +## 6 1 1 10 4 150
+# Compare the samples with and without secondary censoring to the true +# distribution +# Calculate empirical CDF +empirical_cdf <- ecdf(samples) + +# Create a sequence of x values for the theoretical CDF +x_seq <- seq(min(samples), max(samples), length.out = 100) + +# Calculate theoretical CDF +theoretical_cdf <- plnorm(x_seq, meanlog = meanlog, sdlog = sdlog) + +# Create a long format data frame for plotting +cdf_data <- data.frame( + x = rep(x_seq, 2), + probability = c(empirical_cdf(x_seq), theoretical_cdf), + type = rep(c("Observed", "Theoretical"), each = length(x_seq)), + stringsAsFactors = FALSE +) + +# Plot +ggplot(cdf_data, aes(x = x, y = probability, color = type)) + + geom_step(linewidth = 1) + + scale_color_manual( + values = c(Observed = "#4292C6", Theoretical = "#252525") + ) + + geom_vline( + aes(xintercept = mean(samples), color = "Observed"), + linetype = "dashed", linewidth = 1 + ) + + geom_vline( + aes(xintercept = exp(meanlog + sdlog^2 / 2), color = "Theoretical"), + linetype = "dashed", linewidth = 1 + ) + + labs( + title = "Comparison of Observed vs Theoretical CDF", + x = "Delay", + y = "Cumulative Probability", + color = "CDF Type" + ) + + theme_minimal() + + theme( + panel.grid.minor = element_blank(), + plot.title = element_text(hjust = 0.5), + legend.position = "bottom" + )
# Compare the samples with and without secondary censoring to the true +# distribution +# Calculate empirical CDF +empirical_cdf <- ecdf(samples) + +# Create a sequence of x values for the theoretical CDF +x_seq <- seq(min(samples), max(samples), length.out = 100) + +# Calculate theoretical CDF +theoretical_cdf <- plnorm(x_seq, meanlog = meanlog, sdlog = sdlog) + +# Create a long format data frame for plotting +cdf_data <- data.frame( + x = rep(x_seq, 2), + probability = c(empirical_cdf(x_seq), theoretical_cdf), + type = rep(c("Observed", "Theoretical"), each = length(x_seq)), + stringsAsFactors = FALSE +) + +# Plot +ggplot(cdf_data, aes(x = x, y = probability, color = type)) + + geom_step(linewidth = 1) + + scale_color_manual( + values = c(Observed = "#4292C6", Theoretical = "#252525") + ) + + geom_vline( + aes(xintercept = mean(samples), color = "Observed"), + linetype = "dashed", linewidth = 1 + ) + + geom_vline( + aes(xintercept = exp(meanlog + sdlog^2 / 2), color = "Theoretical"), + linetype = "dashed", linewidth = 1 + ) + + labs( + title = "Comparison of Observed vs Theoretical CDF", + x = "Delay", + y = "Cumulative Probability", + color = "CDF Type" + ) + + theme_minimal() + + theme( + panel.grid.minor = element_blank(), + plot.title = element_text(hjust = 0.5), + legend.position = "bottom" + )
We’ve aggregated the data to unique combinations of pwindow, swindow, and obs_time and counted the number of occurrences of each observed_delay for each combination. This is the data we will use to fit our model.
pwindow
swindow
obs_time
observed_delay
We’ll start by fitting a naive model using cmdstan. We’ll use the cmdstanr package to interface with cmdstan. We define the model in a string and then write it to a file as in the How to use primarycensoreddist with Stan vignette.
+writeLines( + "data { + int<lower=0> N; // number of observations + vector[N] y; // observed delays + vector[N] n; // number of occurrences for each delay + } + parameters { + real mu; + real<lower=0> sigma; + } + model { + // Priors + mu ~ normal(1, 1); + sigma ~ normal(0.5, 1); + + // Likelihood + target += n .* lognormal_lpdf(y | mu, sigma); + }", + con = file.path(tempdir(), "naive_lognormal.stan") +)
writeLines( + "data { + int<lower=0> N; // number of observations + vector[N] y; // observed delays + vector[N] n; // number of occurrences for each delay + } + parameters { + real mu; + real<lower=0> sigma; + } + model { + // Priors + mu ~ normal(1, 1); + sigma ~ normal(0.5, 1); + + // Likelihood + target += n .* lognormal_lpdf(y | mu, sigma); + }", + con = file.path(tempdir(), "naive_lognormal.stan") +)
Now let’s compile the model
+naive_model <- cmdstan_model(file.path(tempdir(), "naive_lognormal.stan"))
naive_model <- cmdstan_model(file.path(tempdir(), "naive_lognormal.stan"))
and now let’s fit the compiled model.
+naive_fit <- naive_model$sample( + data = list( + # Add a small constant to avoid log(0) + y = delay_counts$observed_delay + 1e-6, + n = delay_counts$n, + N = nrow(delay_counts) + ), + chains = 4, + parallel_chains = 4, + show_messages = FALSE, + refresh = 0 +) +naive_fit
naive_fit <- naive_model$sample( + data = list( + # Add a small constant to avoid log(0) + y = delay_counts$observed_delay + 1e-6, + n = delay_counts$n, + N = nrow(delay_counts) + ), + chains = 4, + parallel_chains = 4, + show_messages = FALSE, + refresh = 0 +) +naive_fit
## Warning: NAs introduced by coercion
## variable mean median sd mad q5 q95 rhat ess_bulk +## lp__ -28480.45 -28480.10 1.00 0.74 -28482.40 -28479.50 1.00 1999 +## mu -0.10 -0.10 0.05 0.04 -0.18 -0.02 1.00 2691 +## sigma 4.61 4.61 0.03 0.03 4.56 4.67 1.00 3820 +## ess_tail +## NA +## 2310 +## 2869
We see that the model has converged and the diagnostics look good. However, just from the model posterior summary we see that we might not be very happy with the fit. mu is smaller than the target 1.5 and sigma is larger than the target 0.75.
mu
sigma
We’ll now fit an improved model using the primarycensoreddist package. The main improvement is that we will use the primary_censored_dist_lpdf function to fit the model. This is the Stan version of the pcens() function and accounts for the primary and secondary censoring windows as well as the truncation. We encode that our primary distribution is a lognormal distribution by passing 1 as the dist_id parameter and that our primary event distribution is uniform by passing 1 as the primary_dist_id parameter. See the Stan documentation for more details on the primary_censored_dist_lpdf function.
primary_censored_dist_lpdf
pcens()
dist_id
primary_dist_id
+writeLines( + " + functions { + #include primary_censored_dist.stan + #include expgrowth.stan + } + data { + int<lower=0> N; // number of observations + array[N] int<lower=0> y; // observed delays + array[N] int<lower=0> n; // number of occurrences for each delay + array[N] int<lower=0> pwindow; // primary censoring window + array[N] int<lower=0> swindow; // secondary censoring window + array[N] int<lower=0> D; // maximum delay + } + transformed data { + array[0] real primary_params; + } + parameters { + real mu; + real<lower=0> sigma; + } + model { + // Priors + mu ~ normal(1, 1); + sigma ~ normal(0.5, 0.5); + + // Likelihood + for (i in 1:N) { + target += n[i] * primary_censored_dist_lpmf( + y[i] | 1, {mu, sigma}, + pwindow[i], swindow[i], D[i], + 1, primary_params + ); + } + }", + con = file.path(tempdir(), "primarycensoreddist_lognormal.stan") +)
writeLines( + " + functions { + #include primary_censored_dist.stan + #include expgrowth.stan + } + data { + int<lower=0> N; // number of observations + array[N] int<lower=0> y; // observed delays + array[N] int<lower=0> n; // number of occurrences for each delay + array[N] int<lower=0> pwindow; // primary censoring window + array[N] int<lower=0> swindow; // secondary censoring window + array[N] int<lower=0> D; // maximum delay + } + transformed data { + array[0] real primary_params; + } + parameters { + real mu; + real<lower=0> sigma; + } + model { + // Priors + mu ~ normal(1, 1); + sigma ~ normal(0.5, 0.5); + + // Likelihood + for (i in 1:N) { + target += n[i] * primary_censored_dist_lpmf( + y[i] | 1, {mu, sigma}, + pwindow[i], swindow[i], D[i], + 1, primary_params + ); + } + }", + con = file.path(tempdir(), "primarycensoreddist_lognormal.stan") +)
+primarycensoreddist_model <- cmdstan_model( + file.path(tempdir(), "primarycensoreddist_lognormal.stan"), + include_paths = pcd_stan_path() +)
primarycensoreddist_model <- cmdstan_model( + file.path(tempdir(), "primarycensoreddist_lognormal.stan"), + include_paths = pcd_stan_path() +)
Now let’s fit the compiled model.
+primarycensoreddist_fit <- primarycensoreddist_model$sample( + data = list( + y = delay_counts$observed_delay, + n = delay_counts$n, + pwindow = delay_counts$pwindow, + swindow = delay_counts$swindow, + D = delay_counts$obs_time, + N = nrow(delay_counts) + ), + chains = 4, + parallel_chains = 4, + init = list( # we use this to resolve initialisation issues + list(mu = 1.5, sigma = 0.6), + list(mu = 1.5, sigma = 0.4), + list(mu = 1.5, sigma = 0.3), + list(mu = 1.5, sigma = 0.55) + ), + refresh = 0, + show_messages = FALSE +) +primarycensoreddist_fit
primarycensoreddist_fit <- primarycensoreddist_model$sample( + data = list( + y = delay_counts$observed_delay, + n = delay_counts$n, + pwindow = delay_counts$pwindow, + swindow = delay_counts$swindow, + D = delay_counts$obs_time, + N = nrow(delay_counts) + ), + chains = 4, + parallel_chains = 4, + init = list( # we use this to resolve initialisation issues + list(mu = 1.5, sigma = 0.6), + list(mu = 1.5, sigma = 0.4), + list(mu = 1.5, sigma = 0.3), + list(mu = 1.5, sigma = 0.55) + ), + refresh = 0, + show_messages = FALSE +) +primarycensoreddist_fit
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail +## lp__ -2134.42 -2134.13 0.97 0.70 -2136.43 -2133.49 1.00 732 952 +## mu 1.53 1.53 0.05 0.05 1.46 1.61 1.01 620 686 +## sigma 0.77 0.77 0.04 0.04 0.71 0.83 1.00 633 641
We see that the model has converged and the diagnostics look good. We also see that the posterior means are very near the true parameters and the 90% credible intervals include the true parameters.
Stan is a probabilistic programming language for statistical inference. It provides a flexible and efficient platform for Bayesian modeling and is widely used in various fields of data science and statistics. In this vignette, we’ll use Stan in conjunction with primarycensoreddist to perform Bayesian inference on censored data.
For more information on Stan: -- Stan’s official website -- Stan documentation -- Stan forums for community support and discussions
For more information on Stan:
Stan functions are accessed using the pcd_load_stan_functions() function. This function takes the name of the function as an argument and returns the function as a string. It can additionally write the functions to a file and wrap them in a functions{} block.
pcd_load_stan_functions()
functions{}
pcd_load_stan_functions("primary_censored_dist_lpmf")
## [1] "// Stan functions from primarycensoreddist version 0.1.0.1000\nreal primary_censored_dist_lpmf(int d, int dist_id, array[] real params,\n data real pwindow, real swindow, data real D,\n int primary_dist_id, array[] real primary_params) {\n\n real d_upper = d + swindow;\n if (d_upper > D) {\n reject(\"Upper truncation point is greater than D. It is \", d_upper,\n \" and D is \", D, \". Resolve this by increasing D to be greater or equal to d + swindow.\");\n }\n real log_cdf_upper = primary_censored_dist_lcdf(\n d_upper | dist_id, params, pwindow, D, primary_dist_id, primary_params\n );\n real log_cdf_lower = primary_censored_dist_lcdf(\n d | dist_id, params, pwindow, D, primary_dist_id, primary_params\n );\n return log_diff_exp(log_cdf_upper, log_cdf_lower);\n}"
## [1] "// Stan functions from primarycensoreddist version 0.2.0.1000\nreal primary_censored_dist_lpmf(int d, int dist_id, array[] real params,\n data real pwindow, real swindow, data real D,\n int primary_dist_id, array[] real primary_params) {\n\n real d_upper = d + swindow;\n if (d_upper > D) {\n reject(\"Upper truncation point is greater than D. It is \", d_upper,\n \" and D is \", D, \". Resolve this by increasing D to be greater or equal to d + swindow.\");\n }\n real log_cdf_upper = primary_censored_dist_lcdf(\n d_upper | dist_id, params, pwindow, D, primary_dist_id, primary_params\n );\n real log_cdf_lower = primary_censored_dist_lcdf(\n d | dist_id, params, pwindow, D, primary_dist_id, primary_params\n );\n return log_diff_exp(log_cdf_upper, log_cdf_lower);\n}"
## Stan functions written to: /tmp/RtmpaPRJO3/expgrowth_rng.stan
## Stan functions written to: /tmp/RtmpaPzbyH/expgrowth_rng.stan
This can now be compiled and used in the same way as any other cmdstanr model.
model <- cmdstan_model(expgrowth_rng_file) model
## functions { -## // Stan functions from primarycensoreddist version 0.1.0.1000 +## // Stan functions from primarycensoreddist version 0.2.0.1000 ## real expgrowth_rng(real min, real max, real r) { ## real u = uniform_rng(0, 1); ## if (abs(r) < 1e-10) { @@ -203,20 +206,44 @@ text = c( "functions {", "#include expgrowth.stan", + "}", + "generated quantities {", + " real y = expgrowth_rng(0, 1, 0.4);", "}" ), con = expgrowth_stan_file )
}
We can now use this file to compile a model. **Note** that we need to include the path to the `primarycensoreddist` Stan functions using the `include_paths` argument to `cmdstan_model()`. - - -``` r -model <- cmdstan_model(expgrowth_stan_file, include_paths = pcd_stan_path()) -model
We can now use this file to compile a model. Note that we need to include the path to the primarycensoreddist Stan functions using the include_paths argument to cmdstan_model().
include_paths
cmdstan_model()
+model <- cmdstan_model(expgrowth_stan_file, include_paths = pcd_stan_path()) +model
model <- cmdstan_model(expgrowth_stan_file, include_paths = pcd_stan_path()) +model
## functions { ## #include expgrowth.stan +## } +## generated quantities { +## real y = expgrowth_rng(0, 1, 0.4); ## }
We can then sample from the model (we set fixed_param = TRUE here as our toy example doesn’t require MCMC sampling).
fixed_param = TRUE
+samples <- model$sample(chains = 1, fixed_param = TRUE)
samples <- model$sample(chains = 1, fixed_param = TRUE)
## Running MCMC with 1 chain... +## +## Chain 1 Iteration: 1 / 1000 [ 0%] (Sampling) +## Chain 1 Iteration: 100 / 1000 [ 10%] (Sampling) +## Chain 1 Iteration: 200 / 1000 [ 20%] (Sampling) +## Chain 1 Iteration: 300 / 1000 [ 30%] (Sampling) +## Chain 1 Iteration: 400 / 1000 [ 40%] (Sampling) +## Chain 1 Iteration: 500 / 1000 [ 50%] (Sampling) +## Chain 1 Iteration: 600 / 1000 [ 60%] (Sampling) +## Chain 1 Iteration: 700 / 1000 [ 70%] (Sampling) +## Chain 1 Iteration: 800 / 1000 [ 80%] (Sampling) +## Chain 1 Iteration: 900 / 1000 [ 90%] (Sampling) +## Chain 1 Iteration: 1000 / 1000 [100%] (Sampling) +## Chain 1 finished in 0.0 seconds.
+samples
samples
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail +## y 0.51 0.52 0.29 0.38 0.06 0.95 1.00 1065 983
Whilst it is possible to use Stan functions directly in R it is not recommended for most use cases (use the R functions in primarycensoreddist instead). However, it can be useful to understand what is going on under the hood or for exploration (indeed we use this internally in primarycensoreddist to check our functions against the R implementations). To do this we use the expose_functions() method on our already compiled model.
expose_functions()
+ model$expose_functions(global = TRUE) -## ar: creating stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_nvecserial.a -## ar: creating stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_cvodes.a -## ar: creating stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_idas.a -## ar: creating stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_kinsol.a -## /home/runner/.cmdstan/cmdstan-2.35.0/stan/lib/stan_math/lib/tbb_2020.3/build/Makefile.tbb:28: CONFIG: cfg=release arch=intel64 compiler=gcc target=linux runtime=cc11.4.0_libc2.35_kernel6.5.0 We can now use the function in R. Note that this may get slightly more complicated if our stan function depends on other stan functions (i.e. you need to have those included in your compiled model as well). - + expgrowth_rng(0, 1, 0.4) -## [1] 0.5643954 +## [1] 0.6524848
model$expose_functions(global = TRUE)
## ar: creating stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_nvecserial.a
## ar: creating stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_cvodes.a
## ar: creating stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_idas.a
## ar: creating stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_kinsol.a
## /home/runner/.cmdstan/cmdstan-2.35.0/stan/lib/stan_math/lib/tbb_2020.3/build/Makefile.tbb:28: CONFIG: cfg=release arch=intel64 compiler=gcc target=linux runtime=cc11.4.0_libc2.35_kernel6.5.0
We can now use the function in R. Note that this may get slightly more complicated if our stan function depends on other stan functions (i.e. you need to have those included in your compiled model as well).
+ expgrowth_rng(0, 1, 0.4) -## [1] 0.5643954 +## [1] 0.6524848
expgrowth_rng(0, 1, 0.4)
## [1] 0.5643954
## [1] 0.6524848
Development version.
This release puts in place initial documentation and vignettes. It also includes a new primary censored distribution interface to allow for non-secondary event censored distributions. Development of this release as identified some numerical issues in the gradient evaluations for the primary censored distributions which may lead to breaking interface changes in 0.3.0 for the Stan code.
0.3.0
swindow = 0