From d14cc4c325d1c11ffc47b7c38c6d4aec9450441c Mon Sep 17 00:00:00 2001 From: cm401 Date: Mon, 15 Jun 2026 10:04:30 +0100 Subject: [PATCH 01/10] small plot changes --- analysis/plot_data_format_ablation.R | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/analysis/plot_data_format_ablation.R b/analysis/plot_data_format_ablation.R index 1a2ca03..d5ad60a 100644 --- a/analysis/plot_data_format_ablation.R +++ b/analysis/plot_data_format_ablation.R @@ -116,7 +116,7 @@ names(comparison_tbl) <- .rename_arm_cols(names(comparison_tbl)) } # Consistent theme for all figures. -theme_ablation <- function(base_size = 10) { +theme_ablation <- function(base_size = 18) { theme_minimal(base_size = base_size) + theme( panel.grid.major.y = element_blank(), @@ -447,7 +447,7 @@ if (length(top_pathogens) == 0L) { geom_text(data = js_labels, aes(x = x_pos, y = y_pos, label = label), inherit.aes = FALSE, - size = 2.6, hjust = 0, parse = TRUE, colour = "grey30") + + size = 4.6, hjust = 0, parse = TRUE, colour = "grey30") + scale_fill_aaas(labels = ARM_LABELS) + scale_colour_aaas(labels = ARM_LABELS) + scale_linewidth_manual(values = line_lwd, labels = ARM_LABELS) + @@ -887,7 +887,7 @@ pA <- ggplot(forest_both, scale_size_manual( values = ARM_SIZES_C_FULL) + scale_linewidth_manual(values = ARM_EBW_C_FULL) + labs(x = "Days (log scale)", - subtitle = "Best-fitting distribution per pathogen (main analysis model weights)") + + subtitle = "Best-fitting distribution per pathogen\n(main analysis)") + theme_ablation() + guides(colour = guide_legend(override.aes = list(size = 2.5)), size = "none", @@ -906,8 +906,8 @@ pC <- ggplot(tau_long, geom_point(size = 3.5, position = position_dodge(width = COMB_DODGE)) + scale_colour_aaas() + scale_shape_manual( values = ARM_SHAPES_FULL) + - labs(x = expression(tau~"(heterogeneity SD)"), - subtitle = "Between-study heterogeneity") + + labs(x = expression(atop(tau, "(heterogeneity SD)")), + subtitle = "Between-study\nheterogeneity") + y_shared # Panel D: predictive CrI ratio — the combined (confounded) signal. @@ -950,7 +950,7 @@ pD <- ggplot(pred_ratio, geom_point(size = 3.5) + scale_colour_aaas() + scale_shape_manual( values = COMP_SHAPES) + - labs(x = "Predictive P50 CrI ratio (arm / F)", + labs(x = "Predictive P50 CrI\nratio (arm / F)", subtitle = "Combined effect (P50)") + y_shared @@ -962,7 +962,7 @@ pE2 <- ggplot(pred_q95_ratio, geom_point(size = 3.5) + scale_colour_aaas() + scale_shape_manual(values = COMP_SHAPES) + - labs(x = "Predictive P95 CrI ratio (arm / F)", + labs(x = "Predictive P95 CrI\nratio (arm / F)", subtitle = "Combined effect (P95)") + y_shared @@ -989,7 +989,7 @@ pTop <- ggplot(mu0_ratio_all, theme( axis.title.y = element_text(), axis.title.x = element_blank(), - axis.text.x = element_text(angle = 45, hjust = 1, size = 8) + axis.text.x = element_text(angle = 45, hjust = 1, size = 14) ) pTop2 <- ggplot(pred_q95_ratio_all, @@ -1007,7 +1007,7 @@ pTop2 <- ggplot(pred_q95_ratio_all, theme( axis.title.y = element_text(), axis.title.x = element_blank(), - axis.text.x = element_text(angle = 45, hjust = 1, size = 8) + axis.text.x = element_text(angle = 45, hjust = 1, size = 14) ) # ── 8i. Assemble and save ───────────────────────────────────────────────────── @@ -1019,7 +1019,7 @@ fig4 <- (pTop / pTop2 / bottom_row) + plot_layout(heights = c(2, 2, 5), guides = "collect") + plot_annotation(tag_levels = "A") & theme(legend.position = "bottom", - plot.tag = element_text(face = "bold", size = 10)) + plot.tag = element_text(face = "bold", size = 18)) fig4_height <- max(12, n_tau_pathogens * FIG_HEIGHT_ROW + 9) From cb63ee6b77c4dccd89226dfd1ff7d43fcd30f7b6 Mon Sep 17 00:00:00 2001 From: cm401 Date: Mon, 15 Jun 2026 10:04:54 +0100 Subject: [PATCH 02/10] make plot more readible --- R/ploting_utils.R | 41 ++++------------------------------------- 1 file changed, 4 insertions(+), 37 deletions(-) diff --git a/R/ploting_utils.R b/R/ploting_utils.R index 1d5602d..0e6eba8 100644 --- a/R/ploting_utils.R +++ b/R/ploting_utils.R @@ -1112,54 +1112,21 @@ generate_matched_moments_plot <- function(burr_kappas = 2, colour = NULL ) + - # Distribution label + parameters — bottom-right, colour-coded, italic. - # Built as a plotmath expression (parse = TRUE) so Greek letters (phi, kappa) - # render correctly in every PDF device, not just cairo_pdf. + # Distribution name — bottom-right, colour-coded, italic. + # Parameters (phi, kappa, mean) are omitted here; see the parameter table. { - sims <- tryCatch(rstan::extract(result_filtered$fit), error = function(e) NULL) - # Sanitise the distribution label: replace non-breaking spaces (\u00a0) # with ordinary spaces so the plotmath parser accepts them inside quotes. dl_safe <- gsub("\u00a0", " ", dist_label) - - if (!is.null(sims)) { - # Line 1: distribution name (quoted plain text in plotmath) - l1 <- sprintf("'%s'", dl_safe) - - # Line 2: mean (optional; quoted plain text) - l2 <- if (!is.null(sims$pred_mean)) - sprintf("'Mean %.1f d'", stats::median(sims$pred_mean)) - else NULL - - # Line 3: shape parameters using plotmath Greek symbols. - # phi → φ (all distributions) - # kappa → κ (Burr XII and GG only) - has_phi <- !is.null(sims$phi) - has_kappa <- best_dist %in% c("burr", "gengamma") && !is.null(sims$kappa) - l3 <- if (has_phi && has_kappa) - sprintf("phi==%.2f~','~kappa==%.2f", - stats::median(sims$phi), stats::median(sims$kappa)) - else if (has_phi) - sprintf("phi==%.2f", stats::median(sims$phi)) - else NULL - - # Stack non-NULL lines using nested atop(), then wrap in italic(). - # Reduce builds: atop(atop(l1, l2), l3) for 3 lines, atop(l1, l2) for 2, etc. - active <- Filter(Negate(is.null), list(l1, l2, l3)) - stacked <- Reduce(function(acc, x) sprintf("atop(%s, %s)", acc, x), active) - annot_label <- sprintf("italic(%s)", stacked) - } else { - annot_label <- sprintf("italic('%s')", dl_safe) - } + annot_label <- sprintf("italic('%s')", dl_safe) ggplot2::annotate( "text", x = x_max * 0.97, y = 0.02, label = annot_label, colour = dist_col, - size = 2.1, + size = 3.5, hjust = 1, vjust = 0, - lineheight = 1.1, parse = TRUE ) } From 4adb02d87339f31e6571aba2fdb1681741e0009f Mon Sep 17 00:00:00 2001 From: cm401 Date: Mon, 15 Jun 2026 10:05:33 +0100 Subject: [PATCH 03/10] Add type=6: double censored intervals --- R/utils.R | 130 +++++++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 109 insertions(+), 21 deletions(-) diff --git a/R/utils.R b/R/utils.R index d769485..ae86c1e 100644 --- a/R/utils.R +++ b/R/utils.R @@ -256,6 +256,13 @@ extract_quantiles <- function(cdf_summary, probs = c(0.5, 0.95), cdf_mat = NULL) #' in that interval. When `freq_lower[i] == freq_upper[i]` the #' observation is treated as exact. `n` is optional and defaults to #' `sum(freq_count)`.} +#' \item{`expo_lower`, `expo_upper`, `event_lower`, `event_upper`, +#' `freq_count`}{Double interval-censored frequency table (summary +#' type 6). Both the exposure window (`expo_lower`, `expo_upper`) and +#' the event window (`event_lower`, `event_upper`) are interval-censored. +#' All five vectors must have the same length. Requires +#' `expo_upper[i] <= event_lower[i]` for all i (positive delays). `n` +#' is optional and defaults to `sum(freq_count)`.} #' } #' Each element may also contain an optional `source` field — a free-text #' character string recording the bibliographic reference for that dataset @@ -292,14 +299,16 @@ extract_quantiles <- function(cdf_summary, probs = c(0.5, 0.95), cdf_mat = NULL) #' fields (populated with zeros for non-type-5 datasets), as these are #' required by the Stan model regardless of which summary types are present. #' -#' @note **Backward compatibility:** The Stan model requires `freq_lower` and -#' `freq_upper` to be present in the data list for all runs, including those -#' that contain only type 1--4 datasets. This is handled automatically when -#' using this function. If you construct the Stan data list manually (rather -#' than via this function), you must include these fields explicitly, e.g.: +#' @note **Backward compatibility:** The Stan model requires `freq_lower`, +#' `freq_upper`, `event_lower`, and `event_upper` to be present in the data +#' list for all runs. These are populated automatically by this function. +#' If constructing the Stan data list manually, fill all four arrays with +#' zeros for non-applicable datasets, e.g.: #' ```r -#' stan_data$freq_lower <- rep(0, stan_data$n_freq_total) -#' stan_data$freq_upper <- rep(0, stan_data$n_freq_total) +#' stan_data$freq_lower <- rep(0, stan_data$n_freq_total) +#' stan_data$freq_upper <- rep(0, stan_data$n_freq_total) +#' stan_data$event_lower <- rep(0, stan_data$n_freq_total) +#' stan_data$event_upper <- rep(0, stan_data$n_freq_total) #' ``` #' @importFrom utils modifyList #' @export @@ -369,11 +378,13 @@ prepare_stan_data_from_datasets <- function(datasets, dist_type = 1, obs_stat1 <- numeric(n_datasets) obs_stat2 <- numeric(n_datasets) obs_stat3 <- numeric(n_datasets) - # Frequency table flat arrays (for summary_type == 4 and 5) - freq_value_all <- numeric(0) - freq_lower_all <- numeric(0) - freq_upper_all <- numeric(0) - freq_count_all <- integer(0) + # Frequency table flat arrays (for summary_type == 4, 5, and 6) + freq_value_all <- numeric(0) + freq_lower_all <- numeric(0) + freq_upper_all <- numeric(0) + event_lower_all <- numeric(0) # event window lower bounds for type 6; 0 elsewhere + event_upper_all <- numeric(0) # event window upper bounds for type 6; 0 elsewhere + freq_count_all <- integer(0) freq_start_vec <- integer(n_datasets) freq_len_vec <- integer(n_datasets) running_start <- 1L @@ -416,10 +427,12 @@ prepare_stan_data_from_datasets <- function(datasets, dist_type = 1, obs_stat3[i] <- 0 # placeholder freq_start_vec[i] <- running_start freq_len_vec[i] <- length(d$freq_value) - freq_value_all <- c(freq_value_all, as.numeric(d$freq_value)) - freq_lower_all <- c(freq_lower_all, rep(0, length(d$freq_value))) # unused for type 4 - freq_upper_all <- c(freq_upper_all, rep(0, length(d$freq_value))) # unused for type 4 - freq_count_all <- c(freq_count_all, as.integer(d$freq_count)) + freq_value_all <- c(freq_value_all, as.numeric(d$freq_value)) + freq_lower_all <- c(freq_lower_all, rep(0, length(d$freq_value))) # unused for type 4 + freq_upper_all <- c(freq_upper_all, rep(0, length(d$freq_value))) # unused for type 4 + event_lower_all <- c(event_lower_all, rep(0, length(d$freq_value))) # unused for type 4 + event_upper_all <- c(event_upper_all, rep(0, length(d$freq_value))) # unused for type 4 + freq_count_all <- c(freq_count_all, as.integer(d$freq_count)) running_start <- running_start + freq_len_vec[i] } else if (!is.null(d$freq_lower) && !is.null(d$freq_upper) && !is.null(d$freq_count)) { @@ -438,12 +451,57 @@ prepare_stan_data_from_datasets <- function(datasets, dist_type = 1, obs_stat3[i] <- 0 # placeholder freq_start_vec[i] <- running_start freq_len_vec[i] <- length(d$freq_lower) - freq_value_all <- c(freq_value_all, rep(0, length(d$freq_lower))) # unused for type 5 - freq_lower_all <- c(freq_lower_all, as.numeric(d$freq_lower)) - freq_upper_all <- c(freq_upper_all, as.numeric(d$freq_upper)) - freq_count_all <- c(freq_count_all, as.integer(d$freq_count)) + freq_value_all <- c(freq_value_all, rep(0, length(d$freq_lower))) # unused for type 5 + freq_lower_all <- c(freq_lower_all, as.numeric(d$freq_lower)) + freq_upper_all <- c(freq_upper_all, as.numeric(d$freq_upper)) + event_lower_all <- c(event_lower_all, rep(0, length(d$freq_lower))) # unused for type 5 + event_upper_all <- c(event_upper_all, rep(0, length(d$freq_lower))) # unused for type 5 + freq_count_all <- c(freq_count_all, as.integer(d$freq_count)) running_start <- running_start + freq_len_vec[i] + } else if (!is.null(d$expo_lower) && !is.null(d$expo_upper) && + !is.null(d$event_lower) && !is.null(d$event_upper) && + !is.null(d$freq_count)) { + # Type 6: double interval-censored frequency table. + # expo_lower / expo_upper: exposure window bounds. + # event_lower / event_upper: event window bounds. + # In Stan, freq_lower / freq_upper carry the exposure bounds; the new + # event_lower / event_upper arrays carry the event bounds. + n_len <- length(d$expo_lower) + if (length(d$expo_upper) != n_len || length(d$event_lower) != n_len || + length(d$event_upper) != n_len || length(d$freq_count) != n_len) { + stop(paste("Dataset", i, + ": expo_lower, expo_upper, event_lower, event_upper and freq_count", + "must all have the same length")) + } + if (any(d$expo_lower > d$expo_upper)) { + stop(paste("Dataset", i, + ": all expo_lower values must be <= their corresponding expo_upper values")) + } + if (any(d$event_lower > d$event_upper)) { + stop(paste("Dataset", i, + ": all event_lower values must be <= their corresponding event_upper values")) + } + if (any(d$expo_upper > d$event_lower)) { + stop(paste("Dataset", i, + ": all expo_upper values must be <= their corresponding event_lower values", + "(delays must be non-negative)")) + } + n_obs_vec[i] <- if (!is.null(d$n)) d$n else sum(d$freq_count) + summary_type[i] <- 6L + obs_stat1[i] <- 0 # placeholder + obs_stat2[i] <- 0 # placeholder + obs_stat3[i] <- 0 # placeholder + freq_start_vec[i] <- running_start + freq_len_vec[i] <- n_len + freq_value_all <- c(freq_value_all, rep(0, n_len)) # unused for type 6 + freq_lower_all <- c(freq_lower_all, as.numeric(d$expo_lower)) # exposure lower bound + freq_upper_all <- c(freq_upper_all, as.numeric(d$expo_upper)) # exposure upper bound + event_lower_all <- c(event_lower_all, as.numeric(d$event_lower)) # event lower bound + event_upper_all <- c(event_upper_all, as.numeric(d$event_upper)) # event upper bound + freq_count_all <- c(freq_count_all, as.integer(d$freq_count)) + running_start <- running_start + n_len + } else { stop(paste("Dataset", i, "does not have recognized summary statistics")) } @@ -468,6 +526,16 @@ prepare_stan_data_from_datasets <- function(datasets, dist_type = 1, fc <- freq_count_all[s:(s + ln - 1)] mid <- (fl + fu) / 2 central_estimates[i] <- sum(mid * fc) / sum(fc) + } else if (summary_type[i] == 6L && freq_len_vec[i] > 0) { + s <- freq_start_vec[i] + ln <- freq_len_vec[i] + el <- freq_lower_all[s:(s + ln - 1)] # expo_lower + er <- freq_upper_all[s:(s + ln - 1)] # expo_upper + evl <- event_lower_all[s:(s + ln - 1)] # event_lower + evu <- event_upper_all[s:(s + ln - 1)] # event_upper + fc <- freq_count_all[s:(s + ln - 1)] + delay_mid <- ((evl + evu) / 2) - ((el + er) / 2) + central_estimates[i] <- sum(delay_mid * fc) / sum(fc) } } valid_centrals <- central_estimates[central_estimates > 0] @@ -485,6 +553,8 @@ prepare_stan_data_from_datasets <- function(datasets, dist_type = 1, freq_value = freq_value_all, freq_lower = freq_lower_all, freq_upper = freq_upper_all, + event_lower = event_lower_all, + event_upper = event_upper_all, freq_count = freq_count_all, freq_start = as.array(freq_start_vec), freq_len = as.array(freq_len_vec) @@ -593,6 +663,14 @@ update_phi_prior <- function(stan_data, datasets) { w <- d$freq_count / sum(d$freq_count) mean_est <- sum(mid * w) sd_est <- sqrt(sum(w * (mid - mean_est)^2)) + } else if (!is.null(d$expo_lower) && !is.null(d$expo_upper) && + !is.null(d$event_lower) && !is.null(d$event_upper) && + !is.null(d$freq_count)) { + delay_mid <- ((d$event_lower + d$event_upper) / 2) - + ((d$expo_lower + d$expo_upper) / 2) + w <- d$freq_count / sum(d$freq_count) + mean_est <- sum(delay_mid * w) + sd_est <- sqrt(sum(w * (delay_mid - mean_est)^2)) } if (is.na(mean_est) || is.na(sd_est) || sd_est <= 0) return(NA_real_) @@ -777,6 +855,14 @@ should_attempt_gg <- function(datasets, w <- d$freq_count / sum(d$freq_count) mean_est <- sum(mid * w) sd_est <- sqrt(sum(w * (mid - mean_est)^2)) + } else if (!is.null(d$expo_lower) && !is.null(d$expo_upper) && + !is.null(d$event_lower) && !is.null(d$event_upper) && + !is.null(d$freq_count)) { + delay_mid <- ((d$event_lower + d$event_upper) / 2) - + ((d$expo_lower + d$expo_upper) / 2) + w <- d$freq_count / sum(d$freq_count) + mean_est <- sum(delay_mid * w) + sd_est <- sqrt(sum(w * (delay_mid - mean_est)^2)) } if (is.na(mean_est) || is.na(sd_est) || mean_est <= 0 || sd_est <= 0) @@ -800,7 +886,9 @@ should_attempt_gg <- function(datasets, # ── Check 2: information richness ────────────────────────────────────────── is_rich <- vapply(datasets, function(d) { - !is.null(d$freq_value) || (!is.null(d$freq_lower) && !is.null(d$freq_upper)) + !is.null(d$freq_value) || + (!is.null(d$freq_lower) && !is.null(d$freq_upper)) || + (!is.null(d$expo_lower) && !is.null(d$event_lower)) }, logical(1)) rich_frac <- mean(is_rich) From a6d034f2bd7c6edc50f0bef6aed3e43f56bd1e2d Mon Sep 17 00:00:00 2001 From: cm401 Date: Mon, 15 Jun 2026 10:06:10 +0100 Subject: [PATCH 04/10] Add double-censored intervals to stan model --- ...archical_data_synthesis_summary_stats.stan | 101 +++++++++++++++++- 1 file changed, 100 insertions(+), 1 deletion(-) diff --git a/inst/stan/hierarchical_data_synthesis_summary_stats.stan b/inst/stan/hierarchical_data_synthesis_summary_stats.stan index 874b193..0486b4d 100644 --- a/inst/stan/hierarchical_data_synthesis_summary_stats.stan +++ b/inst/stan/hierarchical_data_synthesis_summary_stats.stan @@ -149,12 +149,70 @@ functions { return shape * scale * (1 + z / sqrt(shape)); // Fallback for small shape } } + + // Double-censored log-likelihood for one observation (summary type 6). + // + // Computes log P(E in [ex_l, ex_r], O in [ev_l, ev_r]) proportional to theta, + // where D = O - E ~ f_D(d; theta), assuming E uniform on [ex_l, ex_r] and + // O uniform on [ev_l, ev_r]. + // + // The normalising constant (ex_r - ex_l)(ev_r - ev_l) is constant in theta + // and is dropped, consistent with how type 5 drops 1/(ex_r - ex_l). + // + // Four cases: + // Both endpoints point-observed -> log f_D(ev_l - ex_l) + // Point exposure, interval event -> log[F(ev_r - ex_l) - F(ev_l - ex_l)] + // Interval exposure, point event (= type 5) -> log[F(ev_l - ex_l) - F(ev_l - ex_r)] + // Both interval-censored (general) -> 7-point Gauss-Legendre quadrature + real dc_log_lik(real ex_l, real ex_r, real ev_l, real ev_r, + int dist_type, real loc, real phi, real kappa) { + // Case 1: both endpoints point-observed + if (ex_l == ex_r && ev_l == ev_r) { + return dist_logpdf_fun(ev_l - ex_l, dist_type, loc, phi, kappa); + } + // Case 2: point exposure, interval event + if (ex_l == ex_r) { + return log_diff_exp( + dist_log_cdf_fun(ev_r - ex_l, dist_type, loc, phi, kappa), + dist_log_cdf_fun(ev_l - ex_l, dist_type, loc, phi, kappa) + ); + } + // Case 3: interval exposure, point event — reduces to type 5 + if (ev_l == ev_r) { + return log_diff_exp( + dist_log_cdf_fun(ev_l - ex_l, dist_type, loc, phi, kappa), + dist_log_cdf_fun(ev_l - ex_r, dist_type, loc, phi, kappa) + ); + } + // Case 4: both interval-censored — 7-point Gauss-Legendre quadrature on [ex_l, ex_r]. + // Nodes and weights on [-1, 1] (Abramowitz & Stegun table 25.4). + // The Jacobian factor (ex_r - ex_l)/2 is constant in theta and is dropped. + array[7] real t = {-0.9491079123427585, -0.7415311855993945, + -0.4058451513773832, 0.0, + 0.4058451513773832, 0.7415311855993945, + 0.9491079123427585}; + array[7] real w = { 0.1294849661688697, 0.2797053914892767, + 0.3818300505051189, 0.4179591836734694, + 0.3818300505051189, 0.2797053914892767, + 0.1294849661688697}; + real mid = 0.5 * (ex_r + ex_l); + real half = 0.5 * (ex_r - ex_l); + array[7] real log_terms; + for (k in 1:7) { + real e_k = mid + half * t[k]; + log_terms[k] = log(w[k]) + log_diff_exp( + dist_log_cdf_fun(ev_r - e_k, dist_type, loc, phi, kappa), + dist_log_cdf_fun(ev_l - e_k, dist_type, loc, phi, kappa) + ); + } + return log_sum_exp(log_terms); + } } data { int n_datasets; // Number of datasets array[n_datasets] int n_obs; // Sample sizes for each dataset - array[n_datasets] int summary_type; // 1=median+range, 2=median+IQR, 3=mean+sd, 4=raw freq table, 5=interval-censored freq table + array[n_datasets] int summary_type; // 1=median+range, 2=median+IQR, 3=mean+sd, 4=raw freq table, 5=interval-censored freq table, 6=double interval-censored freq table int dist_type; // 1=lognormal, 2=gamma, 3=weibull, 4=burr XII, 5=gen. gamma // Observed summaries - organized by dataset (used for summary_type 1, 2, 3) @@ -177,6 +235,12 @@ data { array[n_freq_total] real freq_lower; // lower bound of censoring interval array[n_freq_total] real freq_upper; // upper bound of censoring interval + // Event window bounds for summary_type == 6 (double interval-censored). + // For all other summary types these arrays are populated with zeros. + // For type 6, freq_lower / freq_upper carry the exposure window (expo_lower, expo_upper). + array[n_freq_total] real event_lower; // lower bound of event window + array[n_freq_total] real event_upper; // upper bound of event window + // Prior hyperparameters for mu0 ~ normal(mu0_mean, mu0_sd) real mu0_mean; real mu0_sd; @@ -221,6 +285,10 @@ transformed data { if (freq_len[d] == 0) { reject("For summary_type=5, freq_len must be > 0"); } + } else if (summary_type[d] == 6) { + if (freq_len[d] == 0) { + reject("For summary_type=6, freq_len must be > 0"); + } } } } @@ -375,6 +443,22 @@ model { } } } + + else if (summary_type[d] == 6) { // double interval-censored frequency table + // Likelihood: count * dc_log_lik(expo_lower, expo_upper, event_lower, event_upper). + // freq_lower / freq_upper carry the exposure window bounds for type 6. + // Uses 7-point GL quadrature; degenerates correctly to type 5 when + // event_lower[i] == event_upper[i] (point event observation). + int s = freq_start[d]; + int len = freq_len[d]; + for (i in s:(s + len - 1)) { + target += freq_count[i] * dc_log_lik( + freq_lower[i], freq_upper[i], + event_lower[i], event_upper[i], + dist_type, loc, phi, kappa + ); + } + } } } @@ -677,6 +761,21 @@ generated quantities { log_lik[idx] = ll_type5; log_lik[idx + 1] = 0; // unused log_lik[idx + 2] = 0; // unused + + } else if (summary_type[d] == 6) { // double interval-censored frequency table + real ll_type6 = 0; + int s = freq_start[d]; + int len = freq_len[d]; + for (i in s:(s + len - 1)) { + ll_type6 += freq_count[i] * dc_log_lik( + freq_lower[i], freq_upper[i], + event_lower[i], event_upper[i], + dist_type, loc, phi, kappa + ); + } + log_lik[idx] = ll_type6; + log_lik[idx + 1] = 0; // unused + log_lik[idx + 2] = 0; // unused } idx += 3; From 2331e1ef6a5e477c2c2d274f35bf326a9b800186 Mon Sep 17 00:00:00 2001 From: cm401 Date: Mon, 15 Jun 2026 10:06:17 +0100 Subject: [PATCH 05/10] Add double-censored intervals to stan model --- ...al_data_synthesis_summary_stats_joint.stan | 79 ++++++++++++++++++- 1 file changed, 78 insertions(+), 1 deletion(-) diff --git a/inst/stan/hierarchical_data_synthesis_summary_stats_joint.stan b/inst/stan/hierarchical_data_synthesis_summary_stats_joint.stan index 3b3070e..3643c4e 100644 --- a/inst/stan/hierarchical_data_synthesis_summary_stats_joint.stan +++ b/inst/stan/hierarchical_data_synthesis_summary_stats_joint.stan @@ -221,12 +221,52 @@ functions { return shape * scale * (1 + z / sqrt(shape)); } } + + // Double-censored log-likelihood for one observation (summary type 6). + // Identical to the version in hierarchical_data_synthesis_summary_stats.stan. + real dc_log_lik(real ex_l, real ex_r, real ev_l, real ev_r, + int dist_type, real loc, real phi, real kappa) { + if (ex_l == ex_r && ev_l == ev_r) { + return dist_logpdf_fun(ev_l - ex_l, dist_type, loc, phi, kappa); + } + if (ex_l == ex_r) { + return log_diff_exp( + dist_log_cdf_fun(ev_r - ex_l, dist_type, loc, phi, kappa), + dist_log_cdf_fun(ev_l - ex_l, dist_type, loc, phi, kappa) + ); + } + if (ev_l == ev_r) { + return log_diff_exp( + dist_log_cdf_fun(ev_l - ex_l, dist_type, loc, phi, kappa), + dist_log_cdf_fun(ev_l - ex_r, dist_type, loc, phi, kappa) + ); + } + array[7] real t = {-0.9491079123427585, -0.7415311855993945, + -0.4058451513773832, 0.0, + 0.4058451513773832, 0.7415311855993945, + 0.9491079123427585}; + array[7] real w = { 0.1294849661688697, 0.2797053914892767, + 0.3818300505051189, 0.4179591836734694, + 0.3818300505051189, 0.2797053914892767, + 0.1294849661688697}; + real mid = 0.5 * (ex_r + ex_l); + real half = 0.5 * (ex_r - ex_l); + array[7] real log_terms; + for (k in 1:7) { + real e_k = mid + half * t[k]; + log_terms[k] = log(w[k]) + log_diff_exp( + dist_log_cdf_fun(ev_r - e_k, dist_type, loc, phi, kappa), + dist_log_cdf_fun(ev_l - e_k, dist_type, loc, phi, kappa) + ); + } + return log_sum_exp(log_terms); + } } data { int n_datasets; array[n_datasets] int n_obs; - array[n_datasets] int summary_type; + array[n_datasets] int summary_type; int dist_type; array[n_datasets] real obs_stat1; @@ -242,6 +282,12 @@ data { array[n_freq_total] real freq_lower; array[n_freq_total] real freq_upper; + // Event window bounds for summary_type == 6 (double interval-censored). + // For all other summary types these arrays are populated with zeros. + // For type 6, freq_lower / freq_upper carry the exposure window (expo_lower, expo_upper). + array[n_freq_total] real event_lower; // lower bound of event window + array[n_freq_total] real event_upper; // upper bound of event window + real mu0_mean; real mu0_sd; @@ -278,6 +324,10 @@ transformed data { if (freq_len[d] == 0) { reject("For summary_type=5, freq_len must be > 0"); } + } else if (summary_type[d] == 6) { + if (freq_len[d] == 0) { + reject("For summary_type=6, freq_len must be > 0"); + } } } } @@ -416,6 +466,18 @@ model { } } } + + else if (summary_type[d] == 6) { // double interval-censored frequency table + int s = freq_start[d]; + int len = freq_len[d]; + for (i in s:(s + len - 1)) { + target += freq_count[i] * dc_log_lik( + freq_lower[i], freq_upper[i], + event_lower[i], event_upper[i], + dist_type, loc, phi, kappa + ); + } + } } } @@ -704,6 +766,21 @@ generated quantities { log_lik[idx] = ll_type5; log_lik[idx + 1] = 0; log_lik[idx + 2] = 0; + + } else if (summary_type[d] == 6) { // double interval-censored frequency table + real ll_type6 = 0; + int s = freq_start[d]; + int len = freq_len[d]; + for (i in s:(s + len - 1)) { + ll_type6 += freq_count[i] * dc_log_lik( + freq_lower[i], freq_upper[i], + event_lower[i], event_upper[i], + dist_type, loc, phi, kappa + ); + } + log_lik[idx] = ll_type6; + log_lik[idx + 1] = 0; + log_lik[idx + 2] = 0; } idx += 3; From f4c8411cb10ee7cb94b75eb66f435422235b19ee Mon Sep 17 00:00:00 2001 From: cm401 Date: Tue, 16 Jun 2026 13:34:46 +0100 Subject: [PATCH 06/10] Add data_type 7 and associated functionality for right truncation adjustment. --- ...archical_data_synthesis_summary_stats.stan | 130 +++++++++++++++++- ...al_data_synthesis_summary_stats_joint.stan | 123 ++++++++++++++++- 2 files changed, 243 insertions(+), 10 deletions(-) diff --git a/inst/stan/hierarchical_data_synthesis_summary_stats.stan b/inst/stan/hierarchical_data_synthesis_summary_stats.stan index 0486b4d..ca9fab2 100644 --- a/inst/stan/hierarchical_data_synthesis_summary_stats.stan +++ b/inst/stan/hierarchical_data_synthesis_summary_stats.stan @@ -150,6 +150,59 @@ functions { } } + // Log of integral_{ex_l}^{ex_r} F_D(T - e) de — the truncation denominator + // for type 7 (right-truncated, doubly censored). Uses 7-point GL quadrature. + // Degenerates to log F_D(T - ex_l) when ex_l == ex_r (point exposure). + real trunc_denom_log(real ex_l, real ex_r, real T, + int dist_type, real loc, real phi, real kappa) { + if (ex_l == ex_r) { + return dist_log_cdf_fun(T - ex_l, dist_type, loc, phi, kappa); + } + array[7] real t = {-0.9491079123427585, -0.7415311855993945, + -0.4058451513773832, 0.0, + 0.4058451513773832, 0.7415311855993945, + 0.9491079123427585}; + array[7] real w = { 0.1294849661688697, 0.2797053914892767, + 0.3818300505051189, 0.4179591836734694, + 0.3818300505051189, 0.2797053914892767, + 0.1294849661688697}; + real mid = 0.5 * (ex_r + ex_l); + real half = 0.5 * (ex_r - ex_l); + array[7] real log_terms; + for (k in 1:7) { + real e_k = mid + half * t[k]; + real max_del = T - e_k; + log_terms[k] = log(w[k]) + dist_log_cdf_fun(max_del, dist_type, loc, phi, kappa); + } + return log_sum_exp(log_terms); + } + + // Log of integral_{ex_l}^{ex_r} [1 - F_D(T - e)] de — right-censored contribution + // for type 7 (onset not yet observed by T). Uses 7-point GL quadrature. + real right_censor_log_lik(real ex_l, real ex_r, real T, + int dist_type, real loc, real phi, real kappa) { + if (ex_l == ex_r) { + return log1m_exp(dist_log_cdf_fun(T - ex_l, dist_type, loc, phi, kappa)); + } + array[7] real t = {-0.9491079123427585, -0.7415311855993945, + -0.4058451513773832, 0.0, + 0.4058451513773832, 0.7415311855993945, + 0.9491079123427585}; + array[7] real w = { 0.1294849661688697, 0.2797053914892767, + 0.3818300505051189, 0.4179591836734694, + 0.3818300505051189, 0.2797053914892767, + 0.1294849661688697}; + real mid = 0.5 * (ex_r + ex_l); + real half = 0.5 * (ex_r - ex_l); + array[7] real log_terms; + for (k in 1:7) { + real e_k = mid + half * t[k]; + real max_del = T - e_k; + log_terms[k] = log(w[k]) + log1m_exp(dist_log_cdf_fun(max_del, dist_type, loc, phi, kappa)); + } + return log_sum_exp(log_terms); + } + // Double-censored log-likelihood for one observation (summary type 6). // // Computes log P(E in [ex_l, ex_r], O in [ev_l, ev_r]) proportional to theta, @@ -212,7 +265,7 @@ functions { data { int n_datasets; // Number of datasets array[n_datasets] int n_obs; // Sample sizes for each dataset - array[n_datasets] int summary_type; // 1=median+range, 2=median+IQR, 3=mean+sd, 4=raw freq table, 5=interval-censored freq table, 6=double interval-censored freq table + array[n_datasets] int summary_type; // 1=median+range, 2=median+IQR, 3=mean+sd, 4=raw freq table, 5=interval-censored freq table, 6=double interval-censored freq table, 7=double interval-censored with right truncation/censoring int dist_type; // 1=lognormal, 2=gamma, 3=weibull, 4=burr XII, 5=gen. gamma // Observed summaries - organized by dataset (used for summary_type 1, 2, 3) @@ -235,11 +288,17 @@ data { array[n_freq_total] real freq_lower; // lower bound of censoring interval array[n_freq_total] real freq_upper; // upper bound of censoring interval - // Event window bounds for summary_type == 6 (double interval-censored). + // Event window bounds for summary_type == 6 and 7 (double interval-censored). // For all other summary types these arrays are populated with zeros. - // For type 6, freq_lower / freq_upper carry the exposure window (expo_lower, expo_upper). - array[n_freq_total] real event_lower; // lower bound of event window - array[n_freq_total] real event_upper; // upper bound of event window + // For types 6 and 7, freq_lower / freq_upper carry the exposure window (expo_lower, expo_upper). + array[n_freq_total] real event_lower; // lower bound of event window; 0 for right-censored (type 7) + array[n_freq_total] real event_upper; // upper bound of event window; 0 for right-censored (type 7) + + // Type 7: right truncation / right censoring. + // event_observed[i] = 1 if onset was recorded, 0 if right-censored (onset not yet seen). + // truncation_time[d] = analysis date T in days from the reference date; ignored for other types. + array[n_freq_total] int event_observed; // 1=onset seen, 0=right-censored + array[n_datasets] real truncation_time; // analysis date T per dataset // Prior hyperparameters for mu0 ~ normal(mu0_mean, mu0_sd) real mu0_mean; @@ -289,6 +348,13 @@ transformed data { if (freq_len[d] == 0) { reject("For summary_type=6, freq_len must be > 0"); } + } else if (summary_type[d] == 7) { + if (freq_len[d] == 0) { + reject("For summary_type=7, freq_len must be > 0"); + } + if (truncation_time[d] <= 0) { + reject("For summary_type=7, truncation_time must be > 0"); + } } } } @@ -459,6 +525,35 @@ model { ); } } + + else if (summary_type[d] == 7) { // double interval-censored + right truncation/censoring + // For observed individuals (event_observed[i] == 1): + // log L = dc_log_lik(...) - trunc_denom_log(...) + // i.e. numerator = integral F(O_R-e) - F(O_L-e) de over exposure window + // denominator = integral F(T-e) de over exposure window + // For right-censored individuals (event_observed[i] == 0): + // log L = right_censor_log_lik(...) + // i.e. integral [1 - F(T-e)] de over exposure window + int s = freq_start[d]; + int len = freq_len[d]; + real T = truncation_time[d]; + for (i in s:(s + len - 1)) { + if (event_observed[i] == 1) { + target += freq_count[i] * ( + dc_log_lik(freq_lower[i], freq_upper[i], + event_lower[i], event_upper[i], + dist_type, loc, phi, kappa) + - trunc_denom_log(freq_lower[i], freq_upper[i], T, + dist_type, loc, phi, kappa) + ); + } else { + target += freq_count[i] * right_censor_log_lik( + freq_lower[i], freq_upper[i], T, + dist_type, loc, phi, kappa + ); + } + } + } } } @@ -776,6 +871,31 @@ generated quantities { log_lik[idx] = ll_type6; log_lik[idx + 1] = 0; // unused log_lik[idx + 2] = 0; // unused + + } else if (summary_type[d] == 7) { // double interval-censored + right truncation/censoring + real ll_type7 = 0; + int s = freq_start[d]; + int len = freq_len[d]; + real T = truncation_time[d]; + for (i in s:(s + len - 1)) { + if (event_observed[i] == 1) { + ll_type7 += freq_count[i] * ( + dc_log_lik(freq_lower[i], freq_upper[i], + event_lower[i], event_upper[i], + dist_type, loc, phi, kappa) + - trunc_denom_log(freq_lower[i], freq_upper[i], T, + dist_type, loc, phi, kappa) + ); + } else { + ll_type7 += freq_count[i] * right_censor_log_lik( + freq_lower[i], freq_upper[i], T, + dist_type, loc, phi, kappa + ); + } + } + log_lik[idx] = ll_type7; + log_lik[idx + 1] = 0; // unused + log_lik[idx + 2] = 0; // unused } idx += 3; diff --git a/inst/stan/hierarchical_data_synthesis_summary_stats_joint.stan b/inst/stan/hierarchical_data_synthesis_summary_stats_joint.stan index 3643c4e..1b64ee5 100644 --- a/inst/stan/hierarchical_data_synthesis_summary_stats_joint.stan +++ b/inst/stan/hierarchical_data_synthesis_summary_stats_joint.stan @@ -222,6 +222,59 @@ functions { } } + // Log of integral_{ex_l}^{ex_r} F_D(T - e) de — truncation denominator + // for type 7 (right-truncated, doubly censored). Uses 7-point GL quadrature. + // Degenerates to log F_D(T - ex_l) when ex_l == ex_r (point exposure). + real trunc_denom_log(real ex_l, real ex_r, real T, + int dist_type, real loc, real phi, real kappa) { + if (ex_l == ex_r) { + return dist_log_cdf_fun(T - ex_l, dist_type, loc, phi, kappa); + } + array[7] real t = {-0.9491079123427585, -0.7415311855993945, + -0.4058451513773832, 0.0, + 0.4058451513773832, 0.7415311855993945, + 0.9491079123427585}; + array[7] real w = { 0.1294849661688697, 0.2797053914892767, + 0.3818300505051189, 0.4179591836734694, + 0.3818300505051189, 0.2797053914892767, + 0.1294849661688697}; + real mid = 0.5 * (ex_r + ex_l); + real half = 0.5 * (ex_r - ex_l); + array[7] real log_terms; + for (k in 1:7) { + real e_k = mid + half * t[k]; + real max_del = T - e_k; + log_terms[k] = log(w[k]) + dist_log_cdf_fun(max_del, dist_type, loc, phi, kappa); + } + return log_sum_exp(log_terms); + } + + // Log of integral_{ex_l}^{ex_r} [1 - F_D(T - e)] de — right-censored contribution + // for type 7 (onset not yet observed by T). Uses 7-point GL quadrature. + real right_censor_log_lik(real ex_l, real ex_r, real T, + int dist_type, real loc, real phi, real kappa) { + if (ex_l == ex_r) { + return log1m_exp(dist_log_cdf_fun(T - ex_l, dist_type, loc, phi, kappa)); + } + array[7] real t = {-0.9491079123427585, -0.7415311855993945, + -0.4058451513773832, 0.0, + 0.4058451513773832, 0.7415311855993945, + 0.9491079123427585}; + array[7] real w = { 0.1294849661688697, 0.2797053914892767, + 0.3818300505051189, 0.4179591836734694, + 0.3818300505051189, 0.2797053914892767, + 0.1294849661688697}; + real mid = 0.5 * (ex_r + ex_l); + real half = 0.5 * (ex_r - ex_l); + array[7] real log_terms; + for (k in 1:7) { + real e_k = mid + half * t[k]; + real max_del = T - e_k; + log_terms[k] = log(w[k]) + log1m_exp(dist_log_cdf_fun(max_del, dist_type, loc, phi, kappa)); + } + return log_sum_exp(log_terms); + } + // Double-censored log-likelihood for one observation (summary type 6). // Identical to the version in hierarchical_data_synthesis_summary_stats.stan. real dc_log_lik(real ex_l, real ex_r, real ev_l, real ev_r, @@ -266,7 +319,7 @@ functions { data { int n_datasets; array[n_datasets] int n_obs; - array[n_datasets] int summary_type; + array[n_datasets] int summary_type; // 1=median+range, 2=median+IQR, 3=mean+sd, 4=raw freq table, 5=interval-censored freq table, 6=double interval-censored freq table, 7=double interval-censored with right truncation/censoring int dist_type; array[n_datasets] real obs_stat1; @@ -282,11 +335,17 @@ data { array[n_freq_total] real freq_lower; array[n_freq_total] real freq_upper; - // Event window bounds for summary_type == 6 (double interval-censored). + // Event window bounds for summary_type == 6 and 7 (double interval-censored). // For all other summary types these arrays are populated with zeros. - // For type 6, freq_lower / freq_upper carry the exposure window (expo_lower, expo_upper). - array[n_freq_total] real event_lower; // lower bound of event window - array[n_freq_total] real event_upper; // upper bound of event window + // For types 6 and 7, freq_lower / freq_upper carry the exposure window (expo_lower, expo_upper). + array[n_freq_total] real event_lower; // lower bound of event window; 0 for right-censored (type 7) + array[n_freq_total] real event_upper; // upper bound of event window; 0 for right-censored (type 7) + + // Type 7: right truncation / right censoring. + // event_observed[i] = 1 if onset was recorded, 0 if right-censored (onset not yet seen). + // truncation_time[d] = analysis date T in days from the reference date; ignored for other types. + array[n_freq_total] int event_observed; // 1=onset seen, 0=right-censored + array[n_datasets] real truncation_time; // analysis date T per dataset real mu0_mean; real mu0_sd; @@ -328,6 +387,13 @@ transformed data { if (freq_len[d] == 0) { reject("For summary_type=6, freq_len must be > 0"); } + } else if (summary_type[d] == 7) { + if (freq_len[d] == 0) { + reject("For summary_type=7, freq_len must be > 0"); + } + if (truncation_time[d] <= 0) { + reject("For summary_type=7, truncation_time must be > 0"); + } } } } @@ -478,6 +544,28 @@ model { ); } } + + else if (summary_type[d] == 7) { // double interval-censored + right truncation/censoring + int s = freq_start[d]; + int len = freq_len[d]; + real T = truncation_time[d]; + for (i in s:(s + len - 1)) { + if (event_observed[i] == 1) { + target += freq_count[i] * ( + dc_log_lik(freq_lower[i], freq_upper[i], + event_lower[i], event_upper[i], + dist_type, loc, phi, kappa) + - trunc_denom_log(freq_lower[i], freq_upper[i], T, + dist_type, loc, phi, kappa) + ); + } else { + target += freq_count[i] * right_censor_log_lik( + freq_lower[i], freq_upper[i], T, + dist_type, loc, phi, kappa + ); + } + } + } } } @@ -781,6 +869,31 @@ generated quantities { log_lik[idx] = ll_type6; log_lik[idx + 1] = 0; log_lik[idx + 2] = 0; + + } else if (summary_type[d] == 7) { // double interval-censored + right truncation/censoring + real ll_type7 = 0; + int s = freq_start[d]; + int len = freq_len[d]; + real T = truncation_time[d]; + for (i in s:(s + len - 1)) { + if (event_observed[i] == 1) { + ll_type7 += freq_count[i] * ( + dc_log_lik(freq_lower[i], freq_upper[i], + event_lower[i], event_upper[i], + dist_type, loc, phi, kappa) + - trunc_denom_log(freq_lower[i], freq_upper[i], T, + dist_type, loc, phi, kappa) + ); + } else { + ll_type7 += freq_count[i] * right_censor_log_lik( + freq_lower[i], freq_upper[i], T, + dist_type, loc, phi, kappa + ); + } + } + log_lik[idx] = ll_type7; + log_lik[idx + 1] = 0; // unused + log_lik[idx + 2] = 0; // unused } idx += 3; From f812555e1e6006360924f1cb2f72c50b9ffccb35 Mon Sep 17 00:00:00 2001 From: cm401 Date: Tue, 16 Jun 2026 13:35:18 +0100 Subject: [PATCH 07/10] Initial revision --- vignettes/linelist_analysis.Rmd | 333 ++++++++++++++++++++++++++++++++ 1 file changed, 333 insertions(+) create mode 100644 vignettes/linelist_analysis.Rmd diff --git a/vignettes/linelist_analysis.Rmd b/vignettes/linelist_analysis.Rmd new file mode 100644 index 0000000..ec47d22 --- /dev/null +++ b/vignettes/linelist_analysis.Rmd @@ -0,0 +1,333 @@ +--- +title: "Real-Time Line List Analysis with Right Truncation Correction" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Real-Time Line List Analysis with Right Truncation Correction} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r chunk-opts, include=FALSE, purl=FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + eval = TRUE, + purl = FALSE +) +``` + +## Overview + +During an active outbreak, a line list collects cases as they are reported. +Cases with long delays between exposure and symptom onset are less likely to +have been reported by the analysis date -- they appear to be absent, not +because they don't exist but because they haven't happened yet or haven't been +recorded. This is **right truncation**: the observed sample is a biased subset +of all cases, enriched for short delays. + +A second complication is **double interval censoring**: both the exposure date +(e.g. start of a contact window) and the onset date are often known only to the +day, not the exact hour. + +The **type 7** data format in `ddsynth` handles both simultaneously: + +- Doubly interval-censored delays (both exposure and onset known to the day). +- Right truncation correction: only cases with onset before the analysis date + $T$ are observed; the denominator of the likelihood accounts for this. +- Right-censored cases: individuals with a known exposure window but no onset + recorded yet contribute the survival function $1 - F_D(T - x_i)$. + +## The likelihood + +Let $e_i \sim \text{Uniform}(E^L_i, E^R_i)$ be the latent exposure time and +$s_i \sim \text{Uniform}(S^L_i, S^R_i)$ the latent onset time. Marginalising +over both latent times gives: + +**Observed cases** (onset recorded, $s_i < T$): + +$$ +L_i \;\propto\; \frac{ + \int_{E^L_i}^{E^R_i} \bigl[F_D(S^R_i - e) - F_D(S^L_i - e)\bigr]\, de +}{ + \int_{E^L_i}^{E^R_i} F_D(T - e)\, de +} +$$ + +**Right-censored cases** (exposure known, onset not yet seen): + +$$ +L_i \;\propto\; \int_{E^L_i}^{E^R_i} \bigl[1 - F_D(T - e)\bigr]\, de +$$ + +Both integrals are evaluated via 7-point Gauss-Legendre quadrature, requiring +only repeated CDF evaluations and adding negligible overhead over type 6. + +## Setup + +```{r setup, purl=FALSE} +library(ddsynth) +library(rstan) +library(dplyr) +library(ggplot2) +library(patchwork) + +options(mc.cores = parallel::detectCores()) +rstan_options(auto_write = TRUE) +``` + +## Simulating a line list + +We simulate an outbreak of 300 cases. The true incubation period is +log-normal with $\mu_0 = \log(6)$ and $\phi = 0.45$ (mean $\approx$ 7 days). +The analysis date is day 40 from outbreak start. Cases with exposures close to +day 40 may not yet have developed symptoms -- these become right-censored +observations. + +```{r simulate, purl=FALSE} +set.seed(42) + +reference_date <- as.Date("2024-01-01") +T_cal <- as.Date("2024-02-10") # analysis date +T_num <- as.numeric(T_cal - reference_date) # calendar day 40 + +n_cases <- 300 +true_mu0 <- log(6) +true_phi <- 0.45 + +# Exposure start uniform over outbreak period; 1-day window +expo_start_num <- runif(n_cases, 0, 35) +expo_end_num <- expo_start_num + 1 + +# True latent exposure and incubation delay +expo_true <- runif(n_cases, expo_start_num, expo_end_num) +delay_true <- rlnorm(n_cases, meanlog = true_mu0, sdlog = true_phi) +onset_true <- expo_true + delay_true + +# Only cases with onset before T are observed +observed_mask <- onset_true < T_num + +# For observed cases: onset known to the day (floor to 1-day window) +onset_lower <- floor(onset_true) +onset_upper <- onset_lower + 1 + +# Build line list +linelist <- data.frame( + exposure_start = reference_date + floor(expo_start_num), + exposure_end = reference_date + floor(expo_start_num) + 1, + onset_lower = ifelse(observed_mask, reference_date + onset_lower, NA), + onset_upper = ifelse(observed_mask, reference_date + onset_lower + 1, NA) +) +linelist$onset_lower <- as.Date(linelist$onset_lower, origin = "1970-01-01") +linelist$onset_upper <- as.Date(linelist$onset_upper, origin = "1970-01-01") + +cat("Total cases in line list:", nrow(linelist), "\n") +cat("Observed (onset recorded):", sum(observed_mask), "\n") +cat("Right-censored (no onset yet):", sum(!observed_mask), "\n") +``` + +## Preparing the dataset + +Convert dates to numeric (days from `reference_date`). All bounds in a type 7 +dataset must share the same calendar origin. + +```{r prepare-data, purl=FALSE} +to_num <- function(d) as.numeric(d - reference_date) + +dataset_ll <- list( + expo_lower = to_num(linelist$exposure_start), + expo_upper = to_num(linelist$exposure_end), + event_lower = ifelse(!is.na(linelist$onset_lower), + to_num(linelist$onset_lower), NA_real_), + event_upper = ifelse(!is.na(linelist$onset_upper), + to_num(linelist$onset_upper), NA_real_), + freq_count = rep(1L, nrow(linelist)), + truncation_time = T_num +) + +# Sanity checks +stopifnot(all(dataset_ll$event_upper[!is.na(dataset_ll$event_upper)] <= T_num)) +stopifnot(all(dataset_ll$expo_upper < T_num)) +``` + +The table below summarises the fields: + +| Field | Description | +|---|---| +| `expo_lower`, `expo_upper` | Exposure window (calendar days from reference) | +| `event_lower`, `event_upper` | Onset window; `NA` for right-censored rows | +| `freq_count` | Number of cases with this combination | +| `truncation_time` | Analysis date $T$ (calendar days from reference) | + +## Compiling the model and preparing Stan data + +```{r compile, purl=FALSE} +stan_model_code <- compile_stan_model(model = "factorised") +``` + +```{r stan-data, purl=FALSE} +stan_data_ll <- prepare_stan_data_from_datasets( + datasets = list(linelist = dataset_ll), + dist_type = 1, # log-normal + custom_priors = list( + mu0_mean = log(5), + mu0_sd = 1.0, + log_tau_mean = 0.0, + log_tau_sd = 0.5, + log_phi_mean = -0.7, + log_phi_sd = 0.5 + ) +) +``` + +## Fitting the model + +```{r fit, purl=FALSE} +fit_ll <- rstan::sampling( + stan_model_code, + data = stan_data_ll, + chains = 4, + iter = 4000, + warmup = 1000, + control = list(adapt_delta = 0.9, max_treedepth = 10), + seed = 42 +) + +print(fit_ll, pars = c("mu0", "phi")) +``` + +## Summarising the posterior + +```{r posterior-summary, purl=FALSE} +draws <- rstan::extract(fit_ll) + +mean_delay <- exp(draws$mu0 + draws$phi^2 / 2) + +cat(sprintf( + "Mean incubation period: %.1f days (95%% CI: %.1f -- %.1f)\n", + median(mean_delay), + quantile(mean_delay, 0.025), + quantile(mean_delay, 0.975) +)) +cat("True mean:", round(exp(true_mu0 + true_phi^2 / 2), 1), "days\n") +``` + +## Temporal trend diagnostics + +Even with truncation correction, it is important to verify that the apparent +shortening of delays near the analysis date is fully accounted for by the +model. Two diagnostics address this. + +**Panel A -- Truncation envelope**: Individual observed delays are plotted +against onset date. The coloured curves show model-predicted quantiles of the +delay distribution *truncated* at the maximum observable delay for that onset +date ($T - \text{exposure midpoint}$). If the observed delays track the +envelope, truncation alone explains the pattern. If delays fall systematically +below the envelope, there is genuine shortening beyond what truncation predicts. + +**Panel B -- PIT over time**: The probability integral transform value +$F_D(d_i) / F_D(T - x_i)$ should be approximately Uniform(0, 1) with no +temporal trend if truncation is the only source of apparent delay shortening. +A downward loess trend reveals a genuine delay reduction over time. + +```{r temporal-diagnostics, purl=FALSE} +posterior <- extract_posterior_draws(fit_ll, stan_data_ll) + +plot_temporal_trend( + dataset = dataset_ll, + posterior = posterior, + reference_date = reference_date, + probs = c(0.5, 0.75, 0.95), + n_draws = 200L +) +``` + +### Interpreting the plots + +| Pattern | Interpretation | +|---|---| +| Points track envelope; PIT flat and centred | Truncation fully explains apparent shortening | +| Points below envelope; PIT trends downward | Genuine delay reduction -- consider a time-varying model | +| Points above envelope; PIT trends upward | Delays genuinely increasing or reporting lag | +| PIT multimodal or U-shaped | Parametric family mis-specified | + +## Combining with historic data + +A major advantage of the `ddsynth` framework is the ability to pool line list +estimates with estimates from earlier studies that used different data formats. +Six hypothetical historic studies are defined below, mixing three summary types: + +- **Type 1** (median + range): Studies C and F +- **Type 2** (median + IQR): Studies A, B, and E +- **Type 3** (mean + SD): Study D + +All values are consistent with a log-normal distribution with phi (log-SD) +around 0.45 but slightly varying population means, reflecting genuine +between-study heterogeneity. This is important: `phi` is a **shared** +parameter across all datasets, so summary statistics must be compatible with +the same distributional shape. + +With 7 datasets total (6 historic + 1 line list), the between-study +heterogeneity parameter `tau` becomes identifiable from data, and the +predicted quantities are averaged over the full predictive distribution. + +```{r combine-historic, purl=FALSE} +# --- Type 2: median + IQR --- +# Study A: log-normal(log(5.5), 0.45) +dataset_A <- list(n = 60L, median = 5.50, Q1 = 4.06, Q3 = 7.45) + +# Study B: log-normal(log(6.5), 0.45) +dataset_B <- list(n = 80L, median = 6.50, Q1 = 4.80, Q3 = 8.81) + +# Study E: log-normal(log(6.2), 0.45) +dataset_E <- list(n = 120L, median = 6.20, Q1 = 4.58, Q3 = 8.40) + +# --- Type 1: median + range (min/max as 1st and nth order statistics) --- +# Study C: log-normal(log(7.0), 0.45), n=45 +dataset_C <- list(n = 45L, median = 7.00, min = 2.82, max = 17.37) + +# Study F: log-normal(log(5.8), 0.45), n=35 +dataset_F <- list(n = 35L, median = 5.80, min = 2.45, max = 13.73) + +# --- Type 3: mean + SD --- +# Study D: log-normal(log(6.0), 0.45) -> mean = 6.64, SD = 3.15 +dataset_D <- list(n = 100L, mean = 6.64, sd = 3.15) + +stan_data_joint <- prepare_stan_data_from_datasets( + datasets = list( + study_A = dataset_A, + study_B = dataset_B, + study_C = dataset_C, + study_D = dataset_D, + study_E = dataset_E, + study_F = dataset_F, + linelist = dataset_ll + ), + dist_type = 1, + custom_priors = list( + mu0_mean = log(6), + mu0_sd = 1.0, + log_tau_mean = -1.0, # prior centred on tau ~ 0.37; realistic between-study SD + log_tau_sd = 0.5, + log_phi_mean = -0.7, + log_phi_sd = 0.5 + ) +) + +fit_joint <- rstan::sampling( + stan_model_code, + data = stan_data_joint, + chains = 4, + iter = 4000, + warmup = 1000, + control = list(adapt_delta = 0.95, max_treedepth = 12), + seed = 42 +) + +print(fit_joint, pars = c("mu0", "phi", "tau")) +``` + +With more than 5 datasets, `tau` is identifiable from data and the reported +`pred_mean`, `pred_median` etc. are averaged over the predictive distribution +for a new study, incorporating both within-study uncertainty and +between-study heterogeneity. The type 7 likelihood for the line list +automatically corrects for truncation and right-censoring. From 3312b186ec1f0b76e5ba006bd19645584eb35315 Mon Sep 17 00:00:00 2001 From: cm401 Date: Tue, 16 Jun 2026 13:35:36 +0100 Subject: [PATCH 08/10] Allow for data_type 7 now --- R/utils.R | 155 ++++++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 116 insertions(+), 39 deletions(-) diff --git a/R/utils.R b/R/utils.R index ae86c1e..5673dbc 100644 --- a/R/utils.R +++ b/R/utils.R @@ -378,15 +378,17 @@ prepare_stan_data_from_datasets <- function(datasets, dist_type = 1, obs_stat1 <- numeric(n_datasets) obs_stat2 <- numeric(n_datasets) obs_stat3 <- numeric(n_datasets) - # Frequency table flat arrays (for summary_type == 4, 5, and 6) - freq_value_all <- numeric(0) - freq_lower_all <- numeric(0) - freq_upper_all <- numeric(0) - event_lower_all <- numeric(0) # event window lower bounds for type 6; 0 elsewhere - event_upper_all <- numeric(0) # event window upper bounds for type 6; 0 elsewhere - freq_count_all <- integer(0) - freq_start_vec <- integer(n_datasets) - freq_len_vec <- integer(n_datasets) + # Frequency table flat arrays (for summary_type == 4, 5, 6, and 7) + freq_value_all <- numeric(0) + freq_lower_all <- numeric(0) + freq_upper_all <- numeric(0) + event_lower_all <- numeric(0) # event window lower bounds for types 6/7; 0 elsewhere + event_upper_all <- numeric(0) # event window upper bounds for types 6/7; 0 elsewhere + event_observed_all <- integer(0) # 1=onset seen, 0=right-censored; type 7 only (1 elsewhere) + freq_count_all <- integer(0) + freq_start_vec <- integer(n_datasets) + freq_len_vec <- integer(n_datasets) + truncation_time_vec <- numeric(n_datasets) # analysis date T for type 7; 0 elsewhere running_start <- 1L # Process each dataset @@ -427,13 +429,14 @@ prepare_stan_data_from_datasets <- function(datasets, dist_type = 1, obs_stat3[i] <- 0 # placeholder freq_start_vec[i] <- running_start freq_len_vec[i] <- length(d$freq_value) - freq_value_all <- c(freq_value_all, as.numeric(d$freq_value)) - freq_lower_all <- c(freq_lower_all, rep(0, length(d$freq_value))) # unused for type 4 - freq_upper_all <- c(freq_upper_all, rep(0, length(d$freq_value))) # unused for type 4 - event_lower_all <- c(event_lower_all, rep(0, length(d$freq_value))) # unused for type 4 - event_upper_all <- c(event_upper_all, rep(0, length(d$freq_value))) # unused for type 4 - freq_count_all <- c(freq_count_all, as.integer(d$freq_count)) - running_start <- running_start + freq_len_vec[i] + freq_value_all <- c(freq_value_all, as.numeric(d$freq_value)) + freq_lower_all <- c(freq_lower_all, rep(0, length(d$freq_value))) # unused for type 4 + freq_upper_all <- c(freq_upper_all, rep(0, length(d$freq_value))) # unused for type 4 + event_lower_all <- c(event_lower_all, rep(0, length(d$freq_value))) # unused for type 4 + event_upper_all <- c(event_upper_all, rep(0, length(d$freq_value))) # unused for type 4 + event_observed_all <- c(event_observed_all, rep(1L, length(d$freq_value))) # unused for type 4 + freq_count_all <- c(freq_count_all, as.integer(d$freq_count)) + running_start <- running_start + freq_len_vec[i] } else if (!is.null(d$freq_lower) && !is.null(d$freq_upper) && !is.null(d$freq_count)) { # Type 5: interval-censored frequency table @@ -451,17 +454,18 @@ prepare_stan_data_from_datasets <- function(datasets, dist_type = 1, obs_stat3[i] <- 0 # placeholder freq_start_vec[i] <- running_start freq_len_vec[i] <- length(d$freq_lower) - freq_value_all <- c(freq_value_all, rep(0, length(d$freq_lower))) # unused for type 5 - freq_lower_all <- c(freq_lower_all, as.numeric(d$freq_lower)) - freq_upper_all <- c(freq_upper_all, as.numeric(d$freq_upper)) - event_lower_all <- c(event_lower_all, rep(0, length(d$freq_lower))) # unused for type 5 - event_upper_all <- c(event_upper_all, rep(0, length(d$freq_lower))) # unused for type 5 - freq_count_all <- c(freq_count_all, as.integer(d$freq_count)) + freq_value_all <- c(freq_value_all, rep(0, length(d$freq_lower))) # unused for type 5 + freq_lower_all <- c(freq_lower_all, as.numeric(d$freq_lower)) + freq_upper_all <- c(freq_upper_all, as.numeric(d$freq_upper)) + event_lower_all <- c(event_lower_all, rep(0, length(d$freq_lower))) # unused for type 5 + event_upper_all <- c(event_upper_all, rep(0, length(d$freq_lower))) # unused for type 5 + event_observed_all <- c(event_observed_all, rep(1L, length(d$freq_lower))) # unused for type 5 + freq_count_all <- c(freq_count_all, as.integer(d$freq_count)) running_start <- running_start + freq_len_vec[i] } else if (!is.null(d$expo_lower) && !is.null(d$expo_upper) && !is.null(d$event_lower) && !is.null(d$event_upper) && - !is.null(d$freq_count)) { + !is.null(d$freq_count) && is.null(d$truncation_time)) { # Type 6: double interval-censored frequency table. # expo_lower / expo_upper: exposure window bounds. # event_lower / event_upper: event window bounds. @@ -494,13 +498,69 @@ prepare_stan_data_from_datasets <- function(datasets, dist_type = 1, obs_stat3[i] <- 0 # placeholder freq_start_vec[i] <- running_start freq_len_vec[i] <- n_len - freq_value_all <- c(freq_value_all, rep(0, n_len)) # unused for type 6 - freq_lower_all <- c(freq_lower_all, as.numeric(d$expo_lower)) # exposure lower bound - freq_upper_all <- c(freq_upper_all, as.numeric(d$expo_upper)) # exposure upper bound - event_lower_all <- c(event_lower_all, as.numeric(d$event_lower)) # event lower bound - event_upper_all <- c(event_upper_all, as.numeric(d$event_upper)) # event upper bound - freq_count_all <- c(freq_count_all, as.integer(d$freq_count)) - running_start <- running_start + n_len + freq_value_all <- c(freq_value_all, rep(0, n_len)) # unused for type 6 + freq_lower_all <- c(freq_lower_all, as.numeric(d$expo_lower)) # exposure lower bound + freq_upper_all <- c(freq_upper_all, as.numeric(d$expo_upper)) # exposure upper bound + event_lower_all <- c(event_lower_all, as.numeric(d$event_lower)) # event lower bound + event_upper_all <- c(event_upper_all, as.numeric(d$event_upper)) # event upper bound + event_observed_all <- c(event_observed_all, rep(1L, n_len)) # unused for type 6 + freq_count_all <- c(freq_count_all, as.integer(d$freq_count)) + running_start <- running_start + n_len + + } else if (!is.null(d$expo_lower) && !is.null(d$expo_upper) && + !is.null(d$freq_count) && !is.null(d$truncation_time)) { + # Type 7: doubly interval-censored with right truncation/censoring. + # event_lower / event_upper may contain NAs for right-censored individuals + # (onset not yet observed by the analysis date T). + n_len <- length(d$expo_lower) + evl_raw <- if (!is.null(d$event_lower)) d$event_lower else rep(NA_real_, n_len) + evu_raw <- if (!is.null(d$event_upper)) d$event_upper else rep(NA_real_, n_len) + + if (length(d$expo_upper) != n_len || length(d$freq_count) != n_len || + length(evl_raw) != n_len || length(evu_raw) != n_len) { + stop(paste("Dataset", i, + ": expo_lower, expo_upper, event_lower, event_upper, and freq_count", + "must all have the same length")) + } + if (any(d$expo_lower > d$expo_upper)) { + stop(paste("Dataset", i, + ": all expo_lower values must be <= their corresponding expo_upper values")) + } + obs_mask <- !is.na(evl_raw) + if (any(obs_mask & (evl_raw > evu_raw), na.rm = TRUE)) { + stop(paste("Dataset", i, + ": all event_lower values must be <= their corresponding event_upper values")) + } + if (any(obs_mask & (d$expo_upper > evl_raw), na.rm = TRUE)) { + stop(paste("Dataset", i, + ": all expo_upper values must be <= their corresponding event_lower values", + "(delays must be non-negative)")) + } + if (any(d$expo_upper >= d$truncation_time)) { + stop(paste("Dataset", i, + ": all expo_upper values must be < truncation_time")) + } + if (any(obs_mask & (evu_raw > d$truncation_time), na.rm = TRUE)) { + stop(paste("Dataset", i, + ": all event_upper values must be <= truncation_time")) + } + n_obs_vec[i] <- if (!is.null(d$n)) d$n else sum(d$freq_count) + summary_type[i] <- 7L + truncation_time_vec[i] <- d$truncation_time + obs_stat1[i] <- 0 + obs_stat2[i] <- 0 + obs_stat3[i] <- 0 + freq_start_vec[i] <- running_start + freq_len_vec[i] <- n_len + freq_value_all <- c(freq_value_all, rep(0, n_len)) + freq_lower_all <- c(freq_lower_all, as.numeric(d$expo_lower)) + freq_upper_all <- c(freq_upper_all, as.numeric(d$expo_upper)) + # Replace NA event bounds with 0 for right-censored rows (ignored in Stan) + event_lower_all <- c(event_lower_all, as.numeric(ifelse(obs_mask, evl_raw, 0))) + event_upper_all <- c(event_upper_all, as.numeric(ifelse(obs_mask, evu_raw, 0))) + event_observed_all <- c(event_observed_all, as.integer(obs_mask)) + freq_count_all <- c(freq_count_all, as.integer(d$freq_count)) + running_start <- running_start + n_len } else { stop(paste("Dataset", i, "does not have recognized summary statistics")) @@ -536,6 +596,21 @@ prepare_stan_data_from_datasets <- function(datasets, dist_type = 1, fc <- freq_count_all[s:(s + ln - 1)] delay_mid <- ((evl + evu) / 2) - ((el + er) / 2) central_estimates[i] <- sum(delay_mid * fc) / sum(fc) + } else if (summary_type[i] == 7L && freq_len_vec[i] > 0) { + s <- freq_start_vec[i] + ln <- freq_len_vec[i] + el <- freq_lower_all[s:(s + ln - 1)] + er <- freq_upper_all[s:(s + ln - 1)] + evl <- event_lower_all[s:(s + ln - 1)] + evu <- event_upper_all[s:(s + ln - 1)] + fc <- freq_count_all[s:(s + ln - 1)] + eobs <- event_observed_all[s:(s + ln - 1)] + # Use only fully observed rows for the central estimate + obs_idx <- eobs == 1L + if (any(obs_idx)) { + delay_mid <- ((evl[obs_idx] + evu[obs_idx]) / 2) - ((el[obs_idx] + er[obs_idx]) / 2) + central_estimates[i] <- sum(delay_mid * fc[obs_idx]) / sum(fc[obs_idx]) + } } } valid_centrals <- central_estimates[central_estimates > 0] @@ -549,15 +624,17 @@ prepare_stan_data_from_datasets <- function(datasets, dist_type = 1, obs_stat1 = as.array(obs_stat1), obs_stat2 = as.array(obs_stat2), obs_stat3 = as.array(obs_stat3), - n_freq_total = length(freq_value_all), - freq_value = freq_value_all, - freq_lower = freq_lower_all, - freq_upper = freq_upper_all, - event_lower = event_lower_all, - event_upper = event_upper_all, - freq_count = freq_count_all, - freq_start = as.array(freq_start_vec), - freq_len = as.array(freq_len_vec) + n_freq_total = length(freq_value_all), + freq_value = freq_value_all, + freq_lower = freq_lower_all, + freq_upper = freq_upper_all, + event_lower = event_lower_all, + event_upper = event_upper_all, + event_observed = event_observed_all, + freq_count = freq_count_all, + freq_start = as.array(freq_start_vec), + freq_len = as.array(freq_len_vec), + truncation_time = as.array(truncation_time_vec) ) stan_data$mu0_mean <- if (length(valid_centrals) > 0) log(mean(valid_centrals)) else 0 From 44ab5c4f9d3066c03015a12d686837bd595b1656 Mon Sep 17 00:00:00 2001 From: cm401 Date: Tue, 16 Jun 2026 13:36:11 +0100 Subject: [PATCH 09/10] Temporal trend diagnostics for data type 7 line list fits --- R/ploting_utils.R | 212 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 212 insertions(+) diff --git a/R/ploting_utils.R b/R/ploting_utils.R index 0e6eba8..4c2e25f 100644 --- a/R/ploting_utils.R +++ b/R/ploting_utils.R @@ -1896,3 +1896,215 @@ plot_simulation_study_sens_figure <- function(summary_res, base_size = 9) { fig & guides(color = guide_legend(nrow = 2, byrow = FALSE), shape = guide_legend(nrow = 2, byrow = FALSE)) } + + +# ── Temporal trend diagnostics for type-7 line list fits ────────────────────── + +#' Extract posterior draws of distribution parameters from a stanfit object +#' +#' Returns a list with elements \code{loc}, \code{phi}, and \code{kappa} +#' (each a vector of posterior draws), plus \code{dist_type} (integer). +#' Designed for use with the type-7 temporal diagnostic functions below. +#' +#' @param fit A \code{stanfit} object returned by \code{rstan::sampling()}. +#' @param stan_data The Stan data list passed to fit (needed for \code{dist_type}). +#' @return A list with draws of loc (mu0), phi, kappa, and the dist_type integer. +#' @export +extract_posterior_draws <- function(fit, stan_data) { + draws <- rstan::extract(fit) + list( + loc = as.numeric(draws[["mu0"]]), + phi = as.numeric(draws[["phi"]]), + kappa = as.numeric(draws[["kappa"]]), + dist_type = stan_data$dist_type + ) +} + +# Internal: CDF of the fitted distribution evaluated at scalar x +.dist_cdf <- function(x, dist_type, loc, phi, kappa) { + if (x <= 0) return(0) + if (dist_type == 1L) return(plnorm(x, meanlog = loc, sdlog = phi)) + if (dist_type == 2L) { rate <- phi / exp(loc); return(pgamma(x, shape = phi, rate = rate)) } + if (dist_type == 3L) return(pweibull(x, shape = phi, scale = exp(loc))) + stop("dist_type ", dist_type, " not supported in temporal diagnostics") +} + +# Internal: quantile of the delay distribution *truncated* at max_delay. +# Returns the p-th quantile of F_D(d) / F_D(max_delay). +.trunc_quantile <- function(p, max_delay, dist_type, loc, phi, kappa) { + F_max <- .dist_cdf(max_delay, dist_type, loc, phi, kappa) + if (F_max <= 0) return(NA_real_) + target_p <- p * F_max + if (dist_type == 1L) return(qlnorm(target_p, meanlog = loc, sdlog = phi)) + if (dist_type == 2L) { rate <- phi / exp(loc); return(qgamma(target_p, shape = phi, rate = rate)) } + if (dist_type == 3L) return(qweibull(target_p, shape = phi, scale = exp(loc))) + stop("dist_type ", dist_type, " not supported in temporal diagnostics") +} + +# Internal: PIT value = F_D(delay) / F_D(max_delay) for one observation. +.pit_value <- function(delay, max_delay, dist_type, loc, phi, kappa) { + F_max <- .dist_cdf(max_delay, dist_type, loc, phi, kappa) + if (!is.finite(F_max) || F_max <= 0) return(NA_real_) + .dist_cdf(delay, dist_type, loc, phi, kappa) / F_max +} + +#' Temporal trend diagnostics for a type-7 line list fit +#' +#' Produces a two-panel diagnostic plot: +#' +#' \describe{ +#' \item{Panel A -- Truncation envelope}{Individual observed delays vs onset +#' date, overlaid with posterior median quantile curves (50th, 75th, 95th +#' percentile) of the delay distribution *truncated* at each individual's +#' maximum observable delay (\code{T - expo_mid}). If points track the +#' envelope, truncation accounts for the apparent shortening of delays +#' over time. If points fall systematically below the envelope, delays +#' are genuinely shrinking faster than truncation predicts.} +#' \item{Panel B -- PIT over time}{Probability integral transform values +#' \eqn{F_D(d_i) / F_D(T - x_i)} for each observed case plotted against +#' onset date. Under the null (truncation only, no trend) PIT values are +#' approximately Uniform(0, 1) with no temporal pattern. A downward loess +#' trend indicates a genuine reduction in delays over time.} +#' } +#' +#' @param dataset A single type-7 dataset list (the element passed to +#' \code{prepare_stan_data_from_datasets()}). +#' @param posterior A list returned by \code{extract_posterior_draws()}. +#' @param reference_date A \code{Date} giving day 0 of the calendar axis. +#' If \code{NULL}, the x-axis is labelled in days from reference. +#' @param probs Numeric vector of quantile levels for Panel A. +#' Default: \code{c(0.5, 0.75, 0.95)}. +#' @param n_draws Number of posterior draws to use for envelope uncertainty. +#' Default: 200. +#' @return A patchwork of two ggplot objects. +#' @export +plot_temporal_trend <- function(dataset, posterior, + reference_date = NULL, + probs = c(0.5, 0.75, 0.95), + n_draws = 200L) { + + if (!requireNamespace("patchwork", quietly = TRUE)) + stop("Package 'patchwork' is required for plot_temporal_trend().") + + # ── extract observed rows only ───────────────────────────────────────────── + obs_mask <- !is.na(dataset$event_lower) + expo_lower <- dataset$expo_lower[obs_mask] + expo_upper <- dataset$expo_upper[obs_mask] + evl <- dataset$event_lower[obs_mask] + evu <- dataset$event_upper[obs_mask] + counts <- dataset$freq_count[obs_mask] + T_val <- dataset$truncation_time + + if (sum(obs_mask) == 0L) + stop("No observed (non-censored) rows found in dataset.") + + expo_mid <- (expo_lower + expo_upper) / 2 + event_mid <- (evl + evu) / 2 + delay_obs <- event_mid - expo_mid + max_delay <- T_val - expo_mid + + # Expand counts to individual rows + idx <- rep(seq_along(delay_obs), times = counts) + df <- data.frame( + expo_mid = expo_mid[idx], + event_mid = event_mid[idx], + delay_obs = delay_obs[idx], + max_delay = max_delay[idx] + ) + + # Convert to dates if reference_date is provided + if (!is.null(reference_date)) { + df$onset_date <- reference_date + df$event_mid + x_lab <- "Onset date" + } else { + df$onset_date <- df$event_mid + x_lab <- "Days from reference" + } + + # ── sample posterior draws ───────────────────────────────────────────────── + n_post <- length(posterior$loc) + draw_idx <- sample.int(n_post, min(n_draws, n_post)) + dist_type <- posterior$dist_type + + # ── Panel A: truncation envelope ────────────────────────────────────────── + # For a grid of onset dates, compute posterior-median truncated quantiles. + onset_grid <- seq(min(df$onset_date), max(df$onset_date), length.out = 60) + + if (!is.null(reference_date)) { + max_del_grid <- T_val - as.numeric(onset_grid - reference_date) + } else { + max_del_grid <- T_val - onset_grid + } + + env_list <- lapply(seq_along(probs), function(pi) { + p <- probs[pi] + draws_mat <- vapply(draw_idx, function(j) { + vapply(seq_along(onset_grid), function(g) { + .trunc_quantile(p, max_del_grid[g], + dist_type, posterior$loc[j], + posterior$phi[j], posterior$kappa[j]) + }, numeric(1)) + }, numeric(length(onset_grid))) + data.frame( + onset_date = onset_grid, + q_med = apply(draws_mat, 1, median, na.rm = TRUE), + q_lo = apply(draws_mat, 1, quantile, 0.1, na.rm = TRUE), + q_hi = apply(draws_mat, 1, quantile, 0.9, na.rm = TRUE), + prob_label = paste0(round(p * 100), "th percentile") + ) + }) + env_df <- do.call(rbind, env_list) + + pA <- ggplot2::ggplot(df, ggplot2::aes(x = onset_date, y = delay_obs)) + + ggplot2::geom_point(alpha = 0.4, size = 1.5, colour = "steelblue") + + ggplot2::geom_ribbon( + data = env_df, + ggplot2::aes(x = onset_date, ymin = q_lo, ymax = q_hi, + fill = prob_label, y = NULL), + alpha = 0.15, inherit.aes = FALSE + ) + + ggplot2::geom_line( + data = env_df, + ggplot2::aes(x = onset_date, y = q_med, colour = prob_label), + linewidth = 0.8, inherit.aes = FALSE + ) + + ggplot2::scale_colour_brewer(palette = "Dark2", name = "Truncation envelope") + + ggplot2::scale_fill_brewer(palette = "Dark2", name = "Truncation envelope") + + ggplot2::labs(x = x_lab, y = "Observed delay (days)", + title = "A: Truncation envelope") + + ggplot2::theme_bw(base_size = 11) + + # ── Panel B: PIT over time ───────────────────────────────────────────────── + # Compute posterior-median PIT per individual. + pit_vals <- vapply(seq_len(nrow(df)), function(i) { + pit_draws <- vapply(draw_idx, function(j) { + .pit_value(df$delay_obs[i], df$max_delay[i], + dist_type, posterior$loc[j], + posterior$phi[j], posterior$kappa[j]) + }, numeric(1)) + median(pit_draws, na.rm = TRUE) + }, numeric(1)) + + df$pit <- pit_vals + + pB <- ggplot2::ggplot(df, ggplot2::aes(x = onset_date, y = pit)) + + ggplot2::geom_jitter(alpha = 0.4, size = 1.5, colour = "steelblue", + height = 0.01, width = 0) + + ggplot2::geom_smooth(method = "loess", formula = y ~ x, + se = TRUE, colour = "firebrick", + fill = "firebrick", alpha = 0.15, linewidth = 0.8) + + ggplot2::geom_hline(yintercept = 0.5, linetype = "dashed", + colour = "grey40", linewidth = 0.5) + + ggplot2::scale_y_continuous(limits = c(0, 1), + labels = scales::label_number(accuracy = 0.1)) + + ggplot2::labs(x = x_lab, + y = expression(F[D](d[i]) / F[D](T - x[i])), + title = "B: PIT over time", + caption = paste0( + "Dashed line = 0.5 (expected under no trend).\n", + "Downward loess trend indicates genuine delay reduction." + )) + + ggplot2::theme_bw(base_size = 11) + + patchwork::wrap_plots(pA, pB, ncol = 1) +} From db7dbbaaef9cc57c9c6da843bcaa4b794b50cba9 Mon Sep 17 00:00:00 2001 From: cm401 Date: Tue, 16 Jun 2026 13:51:04 +0100 Subject: [PATCH 10/10] initial revision --- vignettes/double_censoring.Rmd | 447 +++++++++++++++++++++++++++++++++ 1 file changed, 447 insertions(+) create mode 100644 vignettes/double_censoring.Rmd diff --git a/vignettes/double_censoring.Rmd b/vignettes/double_censoring.Rmd new file mode 100644 index 0000000..7694142 --- /dev/null +++ b/vignettes/double_censoring.Rmd @@ -0,0 +1,447 @@ +--- +title: "Double Interval Censoring for General Delays" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Double Interval Censoring for General Delays} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r chunk-opts, include=FALSE, purl=FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + eval = FALSE, # Stan models are too heavy to run in CI; code is shown for reference + purl = FALSE +) +``` + +## Overview + +The package's existing **type 5** data format handles delays where the +**exposure** time is interval-censored (known only within a window) but the +**event** time is observed exactly. This is the typical setting for incubation +period estimation: the date of symptom onset is known but the precise moment of +infection is not. + +For delays other than the incubation period -- such as the **serial interval** +(symptom onset to symptom onset), **generation time** (infection to infection), +or **time from onset to hospitalisation** -- *both* endpoints may be +interval-censored. This vignette demonstrates the **type 6** data format, which +handles this general double-censoring case correctly. + +## The problem with the naive approach + +A tempting shortcut is to form the outer delay bounds +$[\,O_L - E_R,\; O_R - E_L\,]$ and treat the result as a standard +interval-censored delay (type 5). This is equivalent to assigning uniform +probability to all delays in that outer interval. However, the correct +weight function is **trapezoidal**: delays near the centre of the outer +interval are more consistent with the observed windows than delays near the +edges. Using outer bounds therefore over-disperses the likelihood and biases +parameter estimates, with bias growing as censoring windows widen relative to +the mean delay. + +The type 6 likelihood evaluates the correct 1-D integral + +$$ +L \;\propto\; \int_{E_L}^{E_R} + \bigl[F_D(O_R - e) - F_D(O_L - e)\bigr]\, de +$$ + +via Gauss-Legendre quadrature, which requires only CDF evaluations and adds +negligible computational cost. + +## Setup + +```{r setup, purl=FALSE} +library(ddsynth) +library(rstan) +library(dplyr) +library(tidyr) +library(ggplot2) +library(ggsci) +library(patchwork) + +options(mc.cores = parallel::detectCores()) +rstan_options(auto_write = TRUE) +``` + +## Simulating doubly-censored delay data + +We simulate a hypothetical study of the **serial interval** -- the time between +symptom onset in a primary case and symptom onset in a secondary case. The true +distribution is log-normal with population mean $\mu_0 = \log(6)$ and +log-standard-deviation $\phi = 0.45$ (mean serial interval approximately 7 +days). + +Censoring structure: + +- **Exposure window** (`expo_lower`, `expo_upper`): symptom onset of the + primary case is known only to within a 1-3 day window (e.g. recalled + approximately). +- **Event window** (`event_lower`, `event_upper`): symptom onset of the + secondary case is recorded to the nearest day, yielding a 1-day interval. + +```{r simulate, purl=FALSE} +simulate_double_censored <- function(n, + mu = log(6), + sigma = 0.45, + seed = 101) { + set.seed(seed) + + # True serial intervals + D_true <- rlnorm(n, meanlog = mu, sdlog = sigma) + + # Exposure window: onset of primary case known within 1-3 days + expo_width <- sample(1:3, n, replace = TRUE) + expo_lower <- floor(runif(n, 0, 20)) # arbitrary calendar anchor + expo_upper <- expo_lower + expo_width + + # True exposure time: uniform within window + E_true <- runif(n, expo_lower, expo_upper) + + # True event time and daily-rounded event window + O_true <- E_true + D_true + event_lower <- floor(O_true) + event_upper <- event_lower + 1L + + # Ensure positive delay bounds (drop any negative outer-bound rows) + naive_lower <- event_lower - expo_upper + naive_upper <- event_upper - expo_lower + keep <- naive_lower > 0 + + data.frame( + expo_lower = expo_lower[keep], + expo_upper = expo_upper[keep], + event_lower = event_lower[keep], + event_upper = event_upper[keep], + naive_lower = naive_lower[keep], + naive_upper = naive_upper[keep], + D_true = D_true[keep] + ) +} + +dat <- simulate_double_censored(n = 200) + +cat("Simulated", nrow(dat), "observations\n") +cat("True delay summary (days):\n") +print(summary(dat$D_true)) +``` + +## Formatting data for the model + +### Type 6: correct double-censored likelihood + +Each observation contributes one row with its four bounds. Identical window +combinations are collapsed into a frequency table to reduce the size of the +Stan data object. + +```{r format-type6, purl=FALSE} +dat_freq6 <- dat |> + group_by(expo_lower, expo_upper, event_lower, event_upper) |> + summarise(freq_count = n(), .groups = "drop") + +dataset_dc <- list( + expo_lower = dat_freq6$expo_lower, + expo_upper = dat_freq6$expo_upper, + event_lower = dat_freq6$event_lower, + event_upper = dat_freq6$event_upper, + freq_count = as.integer(dat_freq6$freq_count) +) + +cat("Type 6 frequency table: ", nrow(dat_freq6), "distinct interval combinations\n") +``` + +### Type 5 (naive): outer-bounds approximation + +For comparison, we also prepare the same data using the naive approach: each +observation's outer delay bounds $[O_L - E_R,\; O_R - E_L]$ are passed as a +type 5 interval-censored dataset. + +```{r format-type5, purl=FALSE} +dat_freq5 <- dat |> + group_by(naive_lower, naive_upper) |> + summarise(freq_count = n(), .groups = "drop") + +dataset_naive <- list( + freq_lower = dat_freq5$naive_lower, + freq_upper = dat_freq5$naive_upper, + freq_count = as.integer(dat_freq5$freq_count) +) + +cat("Type 5 (naive) frequency table:", nrow(dat_freq5), "distinct intervals\n") +``` + +## Preparing Stan data and compiling the model + +```{r compile, purl=FALSE} +stan_model_code <- compile_stan_model(model = "factorised") +``` + +```{r stan-data, purl=FALSE} +# True parameter values for reference +true_mu <- log(6) +true_sigma <- 0.45 + +# --- Type 6 (correct) --- +stan_data_dc <- prepare_stan_data_from_datasets( + datasets = list(serial_interval = dataset_dc), + dist_type = 1, # log-normal + custom_priors = list( + mu0_sd = 1.0, + log_tau_mean = 0.0, + log_tau_sd = 0.5, + log_phi_mean = -0.7, + log_phi_sd = 0.5 + ) +) + +# --- Type 5 naive --- +stan_data_naive <- prepare_stan_data_from_datasets( + datasets = list(serial_interval = dataset_naive), + dist_type = 1, + custom_priors = list( + mu0_sd = 1.0, + log_tau_mean = 0.0, + log_tau_sd = 0.5, + log_phi_mean = -0.7, + log_phi_sd = 0.5 + ) +) +``` + +## Fitting both models + +```{r fit, purl=FALSE} +sampling_args <- list( + chains = 4, + iter = 4000, + warmup = 1000, + thin = 1, + control = list(adapt_delta = 0.9, max_treedepth = 10), + seed = 42 +) + +cat("Fitting type 6 (double-censored) model...\n") +fit_dc <- do.call( + rstan::sampling, + c(list(object = stan_model_code, data = stan_data_dc), sampling_args) +) + +cat("Fitting type 5 (naive outer-bounds) model...\n") +fit_naive <- do.call( + rstan::sampling, + c(list(object = stan_model_code, data = stan_data_naive), sampling_args) +) +``` + +## Convergence diagnostics + +```{r diagnostics, purl=FALSE} +check_convergence <- function(fit, label) { + s <- summary(fit)$summary + rhat <- s[, "Rhat"] + n_eff <- s[, "n_eff"] + cat(sprintf( + "[%s] Max Rhat: %.3f | Min n_eff: %.0f | Rhat > 1.05: %d\n", + label, + max(rhat, na.rm = TRUE), + min(n_eff, na.rm = TRUE), + sum(rhat > 1.05, na.rm = TRUE) + )) +} + +check_convergence(fit_dc, "Type 6 (correct)") +check_convergence(fit_naive, "Type 5 (naive) ") +``` + +## Parameter recovery + +We compare posterior estimates of $\mu_0$ (population log-mean) and +$\phi$ (log-standard-deviation) against the true simulation values. + +```{r recovery, purl=FALSE} +extract_summary <- function(fit, label) { + sims <- rstan::extract(fit) + data.frame( + model = label, + parameter = c("mu0 (log-mean)", "phi (log-SD)"), + true_val = c(true_mu, true_sigma), + post_mean = c(mean(sims$mu0), mean(sims$phi)), + post_med = c(median(sims$mu0), median(sims$phi)), + ci_low = c(quantile(sims$mu0, 0.025), quantile(sims$phi, 0.025)), + ci_high = c(quantile(sims$mu0, 0.975), quantile(sims$phi, 0.975)) + ) +} + +recovery_df <- bind_rows( + extract_summary(fit_dc, "Type 6 (correct)"), + extract_summary(fit_naive, "Type 5 (naive)") +) + +print(recovery_df, digits = 3) +``` + +The type 6 model should recover both parameters well. The naive type 5 model +will tend to over-estimate $\phi$ (attributing the extra dispersion introduced +by the wide outer bounds to the delay distribution itself) and may show a small +bias in $\mu_0$. + +## Visualisation + +### Posterior distributions for key parameters + +```{r plot-posteriors, fig.width=9, fig.height=4, purl=FALSE} +draws_dc <- rstan::extract(fit_dc) +draws_naive <- rstan::extract(fit_naive) + +posterior_df <- bind_rows( + data.frame(model = "Type 6 (correct)", mu0 = draws_dc$mu0, phi = draws_dc$phi), + data.frame(model = "Type 5 (naive)", mu0 = draws_naive$mu0, phi = draws_naive$phi) +) + +p_mu0 <- ggplot(posterior_df, aes(x = mu0, fill = model, colour = model)) + + geom_density(alpha = 0.35, linewidth = 0.8) + + geom_vline(xintercept = true_mu, linetype = "dashed", linewidth = 0.8) + + scale_fill_npg() + + scale_colour_npg() + + labs(title = expression(paste("Posterior: ", mu[0], " (log-mean)")), + subtitle = "Dashed line = true value", + x = expression(mu[0]), + y = "Density", + fill = NULL, colour = NULL) + + theme_bw(base_size = 12) + + theme(legend.position = "bottom") + +p_phi <- ggplot(posterior_df, aes(x = phi, fill = model, colour = model)) + + geom_density(alpha = 0.35, linewidth = 0.8) + + geom_vline(xintercept = true_sigma, linetype = "dashed", linewidth = 0.8) + + scale_fill_npg() + + scale_colour_npg() + + labs(title = expression(paste("Posterior: ", phi, " (log-SD)")), + subtitle = "Dashed line = true value", + x = expression(phi), + y = "Density", + fill = NULL, colour = NULL) + + theme_bw(base_size = 12) + + theme(legend.position = "bottom") + +p_mu0 + p_phi + plot_layout(guides = "collect") & + theme(legend.position = "bottom") +``` + +### Posterior predictive CDFs + +```{r plot-cdfs, fig.width=8, fig.height=5, purl=FALSE} +x_seq <- seq(0.1, 25, length.out = 400) + +cdf_dc <- compute_predictive_cdf(fit_dc, "lognormal", x_seq = x_seq, n_draws = 500) +cdf_naive <- compute_predictive_cdf(fit_naive, "lognormal", x_seq = x_seq, n_draws = 500) + +# True CDF for comparison +true_cdf_df <- data.frame( + x = x_seq, + cdf = plnorm(x_seq, meanlog = true_mu, sdlog = true_sigma) +) + +cdf_plot_df <- bind_rows( + cdf_dc$summary |> mutate(model = "Type 6 (correct)"), + cdf_naive$summary |> mutate(model = "Type 5 (naive)") +) + +model_colours <- c("Type 6 (correct)" = "#0072B5", "Type 5 (naive)" = "#BC3C29") + +ggplot() + + geom_ribbon( + data = cdf_plot_df, + aes(x = x, ymin = low, ymax = high, fill = model), + alpha = 0.20 + ) + + geom_line( + data = cdf_plot_df, + aes(x = x, y = median, colour = model), + linewidth = 1.0 + ) + + geom_line( + data = true_cdf_df, + aes(x = x, y = cdf), + colour = "black", linetype = "dashed", linewidth = 0.9 + ) + + scale_fill_manual(values = model_colours) + + scale_colour_manual(values = model_colours) + + scale_y_continuous(breaks = seq(0, 1, 0.2), limits = c(0, 1)) + + labs( + title = "Posterior predictive CDFs: type 6 vs naive type 5", + subtitle = "Dashed black line = true log-normal CDF | Shaded bands = 95% credible intervals", + x = "Delay (days)", + y = "Cumulative probability", + fill = NULL, colour = NULL + ) + + theme_bw(base_size = 13) + + theme(legend.position = "bottom") +``` + +### Illustrating the naive over-dispersion + +The plot below shows the distribution of observed delay intervals under the two +approaches for a random sample of 50 observations. The outer-bounds approach +produces wider, overlapping intervals that inflate apparent delay uncertainty. + +```{r plot-intervals, fig.width=8, fig.height=5, purl=FALSE} +set.seed(7) +idx <- sample(nrow(dat), 50) +samp <- dat[idx, ] |> + arrange(D_true) |> + mutate(obs = row_number()) + +interval_df <- bind_rows( + samp |> transmute(obs, approach = "Type 6 (correct)", + lower = event_lower - expo_upper, + upper = event_upper - expo_lower, + mid = (lower + upper) / 2), + samp |> transmute(obs, approach = "Type 5 (naive)", + lower = naive_lower, upper = naive_upper, + mid = (lower + upper) / 2) +) + +ggplot(interval_df, aes(y = obs, colour = approach)) + + geom_segment(aes(x = lower, xend = upper, yend = obs), + linewidth = 0.6, alpha = 0.7) + + geom_point(data = samp |> mutate(obs = row_number()), + aes(x = D_true, y = obs), + colour = "black", shape = 4, size = 1.5, inherit.aes = FALSE) + + scale_colour_manual(values = model_colours) + + labs( + title = "Observed delay intervals for 50 sampled cases", + subtitle = "Crosses mark the true simulated delay", + x = "Delay (days)", + y = "Observation (sorted by true delay)", + colour = NULL + ) + + theme_bw(base_size = 12) + + theme(legend.position = "bottom", + axis.text.y = element_blank(), + axis.ticks.y = element_blank()) +``` + +## When to use type 6 vs type 5 + +| Situation | Recommended format | +|---|---| +| Event time is recorded to the exact day, exposure window is 0-1 days wide | Type 5 (or type 4) | +| Event time is exact; exposure window spans several days | Type 5 | +| Both endpoints interval-censored; windows narrow relative to the mean delay | Type 5 with outer bounds is acceptable | +| Both endpoints interval-censored; windows wide (e.g. > 30% of mean delay) | **Type 6** | +| Serial interval, generation time, or time-to-event with uncertain start | **Type 6** | + +As a rule of thumb, if the sum of the exposure and event window widths +exceeds roughly a third of the expected delay, the naive outer-bounds +approximation will produce a meaningfully biased estimate of $\phi$ and +should be replaced by type 6. + +## Session Info + +```{r session-info, eval=TRUE, purl=FALSE} +sessionInfo() +```