diff --git a/NAMESPACE b/NAMESPACE index 937623b..0c0e395 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -10,6 +10,7 @@ export(bsl_run_diagnostics) export(bsl_summarise_posteriors) export(bsl_to_mcmc_list_multi) export(check_coverage) +export(compile_stan_model) export(compute_iqd) export(compute_median_bias) export(compute_pathogen_model_bayes_factors) diff --git a/R/utils.R b/R/utils.R index 9a78c39..d769485 100644 --- a/R/utils.R +++ b/R/utils.R @@ -944,6 +944,35 @@ gamma_type2_reliable <- function(datasets, } +# Compile Stan model ------------------------------------------------------- + +#' Compile a ddsynth Stan model +#' +#' Compiles and returns one of the Stan models shipped with the package. +#' Pass the returned object to [pre_inference_checks()], [fit_model()], or +#' [run_simulation_study_generalized()]. +#' +#' @param model Character string: `"factorised"` (default) uses the +#' order-statistic factorised likelihood +#' (`hierarchical_data_synthesis_summary_stats.stan`); `"joint"` uses the +#' joint parameterisation +#' (`hierarchical_data_synthesis_summary_stats_joint.stan`). +#' @return A compiled `stanmodel` object. +#' @export +compile_stan_model <- function(model = c("factorised", "joint")) { + model <- match.arg(model) + filename <- switch(model, + factorised = "hierarchical_data_synthesis_summary_stats.stan", + joint = "hierarchical_data_synthesis_summary_stats_joint.stan" + ) + stan_file <- system.file("stan", filename, package = "ddsynth") + if (!nzchar(stan_file)) { + stop("Stan model file '", filename, "' not found in ddsynth installation.") + } + rstan::stan_model(stan_file) +} + + # Pre-inference checks ----------------------------------------------------- #' Run pre-inference checks on a list of datasets diff --git a/analysis/compare_factorised_vs_joint.R b/analysis/compare_factorised_vs_joint.R new file mode 100644 index 0000000..5d48c39 --- /dev/null +++ b/analysis/compare_factorised_vs_joint.R @@ -0,0 +1,655 @@ +# ============================================================================= +# compare_factorised_vs_joint.R +# ----------------------------------------------------------------------------- +# Fits the joint-likelihood Stan model to the same real data used in +# main_analysis.R and produces a comparison table of posterior summaries +# under the two likelihoods. +# +# The likelihood differs only for summary_type 1 (median + range) and +# summary_type 2 (median + IQR), where the factorised model treats each +# reported order statistic independently while the joint model uses the +# true joint density of the three order statistics. Pathogens with no +# type-1/2 datasets are included in the table with a note that the two +# models are equivalent. +# +# Strategy +# -------- +# Rather than re-preparing data, the script reuses the stan_data objects +# stored in main_results.rds, ensuring an exact apples-to-apples comparison +# (same data, same priors, same MCMC settings). The "filtered" analysis is +# used throughout, matching the primary results in the paper. +# +# Output +# ------ +# results/comparison_factorised_vs_joint.rds — raw stanfit objects + summaries +# results/comparison_factorised_vs_joint.csv — formatted comparison table +# +# The RDS contains: +# $joint_results nested list: [[pathogen]][[dist]] with $fit, $stan_data +# $comparison_tbl data.frame with all summary columns for both models +# ============================================================================= + +devtools::load_all(here::here(), quiet = TRUE) +library(rstan) + +rstan_options(auto_write = TRUE) +options(mc.cores = parallel::detectCores()) + + +# ── 1. Settings ─────────────────────────────────────────────────────────────── + +CHAINS <- 4 +ITER <- 12000 +WARMUP <- 2000 +THIN <- 1 +SEED <- 123 +CONTROL <- list(adapt_delta = 0.999, max_treedepth = 12, stepsize = 0.01) + +# Analysis label to compare (must exist in main_results.rds). +ANALYSIS_LABEL <- "filtered" + +# Distributions to include. Set to NULL to use every distribution present +# in the existing results for each pathogen. +DIST_NAMES <- c("lognormal", "gamma", "weibull", "burr", "gengamma") + +# Credible interval width for formatted cells. +CI_PROBS <- c(0.025, 0.5, 0.975) + +OUTPUT_DIR <- here::here("results") +OUTPUT_FILE <- file.path(OUTPUT_DIR, "comparison_factorised_vs_joint.rds") +CSV_FILE <- file.path(OUTPUT_DIR, "comparison_factorised_vs_joint.csv") + +dir.create(OUTPUT_DIR, showWarnings = FALSE) + +# Wipe joint results and refit from scratch when TRUE. +FORCE_RERUN <- FALSE + + +# ── 2. Compile Stan models ──────────────────────────────────────────────────── + +stan_file_fact <- system.file("stan", + "hierarchical_data_synthesis_summary_stats.stan", + package = "ddsynth") +stan_file_joint <- system.file("stan", + "hierarchical_data_synthesis_summary_stats_joint.stan", + package = "ddsynth") + +message("Compiling factorised Stan model...") +stan_model_fact <- rstan::stan_model(stan_file_fact) + +message("Compiling joint Stan model...") +stan_model_joint <- rstan::stan_model(stan_file_joint) + + +# ── 3. Load main results ────────────────────────────────────────────────────── + +main_results_file <- file.path(OUTPUT_DIR, "main_results.rds") +if (!file.exists(main_results_file)) + stop("main_results.rds not found in ", OUTPUT_DIR, + " — run main_analysis.R first.") + +message("Loading main results from: ", main_results_file) +main_results <- readRDS(main_results_file) + + +# ── 4. Load or initialise joint results ────────────────────────────────────── + +if (!FORCE_RERUN && file.exists(OUTPUT_FILE)) { + message("Resuming from existing comparison file: ", OUTPUT_FILE) + stored <- readRDS(OUTPUT_FILE) + joint_results <- stored$joint_results +} else { + joint_results <- list() +} + + +# ── 5. Helper functions ─────────────────────────────────────────────────────── + +# Format a posterior sample vector as "median (lo, hi)". +.fmt <- function(x, digits = 2L) { + x <- x[is.finite(x)] + if (length(x) == 0L) return(NA_character_) + q <- quantile(x, CI_PROBS) + sprintf("%.*f (%.*f, %.*f)", digits, q[2L], digits, q[1L], digits, q[3L]) +} + +# Extract comparison row for one stanfit. +# Returns a named list of formatted strings and raw posterior medians. +.extract_summaries <- function(fit, dist_name, n_datasets) { + sims <- tryCatch(rstan::extract(fit), error = function(e) NULL) + if (is.null(sims)) return(NULL) + + tau_fmt <- if (n_datasets < 5L) + "— (n<5)" + else + .fmt(sims$tau, digits = 2L) + + kappa_fmt <- if (dist_name %in% c("burr", "gengamma") && !is.null(sims$kappa)) + .fmt(sims$kappa, digits = 2L) + else + NA_character_ + + list( + pred_median_fmt = .fmt(sims$pred_median, digits = 1L), + pred_q95_fmt = .fmt(sims$pred_q95, digits = 1L), + mu0_fmt = .fmt(sims$mu0, digits = 2L), + phi_fmt = .fmt(sims$phi, digits = 2L), + tau_fmt = tau_fmt, + kappa_fmt = kappa_fmt, + # Raw posterior medians for computing deltas + pred_median_med = median(sims$pred_median[is.finite(sims$pred_median)]), + pred_q95_med = median(sims$pred_q95[is.finite(sims$pred_q95)]), + mu0_med = median(sims$mu0), + phi_med = median(sims$phi), + tau_med = if (n_datasets >= 5L) median(sims$tau) else NA_real_ + ) +} + + +# ── 6. Fit joint model for each pathogen × distribution ────────────────────── + +pathogens <- names(main_results) + +for (pathogen in pathogens) { + + if (is.null(joint_results[[pathogen]])) + joint_results[[pathogen]] <- list() + + analysis_slot <- main_results[[pathogen]][[ANALYSIS_LABEL]] + if (is.null(analysis_slot)) { + message("\n [SKIP] ", pathogen, ": no '", ANALYSIS_LABEL, "' analysis in main results.") + next + } + + dists_to_run <- if (is.null(DIST_NAMES)) names(analysis_slot) else DIST_NAMES + + for (dist_name in dists_to_run) { + + existing_fact <- analysis_slot[[dist_name]] + + # Skip distributions that were skipped or failed in the main analysis. + if (is.null(existing_fact)) next + if (isTRUE(existing_fact$skipped)) { + message("\n [SKIP] ", pathogen, " / ", dist_name, + " — skipped in main analysis (", existing_fact$reason, ").") + next + } + if (is.null(existing_fact$fit) || is.null(existing_fact$stan_data)) { + message("\n [SKIP] ", pathogen, " / ", dist_name, " — no fit or stan_data.") + next + } + + # Skip if already fitted in a previous run, but retry prior failures + # (where the entry exists but fit = NULL due to a sampling error). + existing_joint <- joint_results[[pathogen]][[dist_name]] + if (!is.null(existing_joint)) { + if (isTRUE(existing_joint$identical_likelihood)) { + message("\n [SKIP] ", pathogen, " / ", dist_name, + " — identical-likelihood sentinel already in results.") + next + } + if (!is.null(existing_joint$fit)) { + message("\n [SKIP] ", pathogen, " / ", dist_name, + " — joint fit already in results.") + next + } + message("\n [RETRY] ", pathogen, " / ", dist_name, + " — previous entry has no fit; retrying.") + } + + stan_data <- existing_fact$stan_data + n_type12 <- sum(stan_data$summary_type %in% c(1L, 2L)) + + # When no type-1/2 datasets are present the two likelihoods are identical; + # store a sentinel rather than wasting compute. + if (n_type12 == 0L) { + message("\n [NOTE] ", pathogen, " / ", dist_name, + " — no type-1/2 datasets; likelihoods identical; storing sentinel.") + joint_results[[pathogen]][[dist_name]] <- list( + identical_likelihood = TRUE, + stan_data = stan_data + ) + next + } + + message("\n Fitting joint model: ", pathogen, " | ", dist_name, + " (n_datasets=", stan_data$n_datasets, + ", n_type12=", n_type12, ")") + + fit_joint <- tryCatch( + rstan::sampling( + stan_model_joint, + data = stan_data, + chains = CHAINS, + iter = ITER, + warmup = WARMUP, + thin = THIN, + control = CONTROL, + seed = SEED, + refresh = 200 + ), + error = function(e) { + message(" [FAIL] Joint sampling failed: ", conditionMessage(e)) + NULL + } + ) + + joint_results[[pathogen]][[dist_name]] <- list( + fit = fit_joint, + stan_data = stan_data, + identical_likelihood = FALSE + ) + + # Persist after each fit so a crash does not lose earlier work. + saveRDS(list(joint_results = joint_results), OUTPUT_FILE) + } +} + +message("\n All joint fits complete.") + + +# ── 7. Build comparison table ───────────────────────────────────────────────── + +dist_labels <- c( + lognormal = "Log-normal", + gamma = "Gamma", + weibull = "Weibull", + burr = "Burr XII", + gengamma = "Gen. gamma" +) + +rows <- list() + +for (pathogen in pathogens) { + + analysis_slot <- main_results[[pathogen]][[ANALYSIS_LABEL]] + if (is.null(analysis_slot)) next + + dists_to_run <- if (is.null(DIST_NAMES)) names(analysis_slot) else DIST_NAMES + + for (dist_name in dists_to_run) { + + fact_res <- analysis_slot[[dist_name]] + joint_res <- joint_results[[pathogen]][[dist_name]] + + # Skip entirely absent or failed entries. + if (is.null(fact_res)) next + if (isTRUE(fact_res$skipped)) next + if (is.null(fact_res$fit) || is.null(fact_res$stan_data)) next + + stan_data <- fact_res$stan_data + n_type12 <- sum(stan_data$summary_type %in% c(1L, 2L)) + + fact_summ <- .extract_summaries(fact_res$fit, dist_name, stan_data$n_datasets) + if (is.null(fact_summ)) next + + # Joint summaries: either extracted from fit, or copied from factorised + # when the likelihoods are identical (no type-1/2 data). + if (!is.null(joint_res) && isTRUE(joint_res$identical_likelihood)) { + joint_summ <- fact_summ # identical by construction + note <- "identical (no type-1/2 data)" + } else if (!is.null(joint_res) && !is.null(joint_res$fit)) { + joint_summ <- .extract_summaries(joint_res$fit, dist_name, stan_data$n_datasets) + note <- "" + } else { + joint_summ <- NULL + note <- "joint fit unavailable" + } + + row <- data.frame( + pathogen = pathogen, + dist = dist_labels[[dist_name]], + n_datasets = stan_data$n_datasets, + n_type12 = n_type12, + note = note, + + # Factorised model + fact_pred_median = if (!is.null(fact_summ)) fact_summ$pred_median_fmt else NA, + fact_pred_q95 = if (!is.null(fact_summ)) fact_summ$pred_q95_fmt else NA, + fact_mu0 = if (!is.null(fact_summ)) fact_summ$mu0_fmt else NA, + fact_phi = if (!is.null(fact_summ)) fact_summ$phi_fmt else NA, + fact_tau = if (!is.null(fact_summ)) fact_summ$tau_fmt else NA, + + # Joint model + joint_pred_median = if (!is.null(joint_summ)) joint_summ$pred_median_fmt else NA, + joint_pred_q95 = if (!is.null(joint_summ)) joint_summ$pred_q95_fmt else NA, + joint_mu0 = if (!is.null(joint_summ)) joint_summ$mu0_fmt else NA, + joint_phi = if (!is.null(joint_summ)) joint_summ$phi_fmt else NA, + joint_tau = if (!is.null(joint_summ)) joint_summ$tau_fmt else NA, + + # Absolute difference in posterior median point estimates (joint - factorised) + delta_pred_median = if (!is.null(fact_summ) && !is.null(joint_summ)) + round(joint_summ$pred_median_med - fact_summ$pred_median_med, 2L) + else NA_real_, + delta_pred_q95 = if (!is.null(fact_summ) && !is.null(joint_summ)) + round(joint_summ$pred_q95_med - fact_summ$pred_q95_med, 2L) + else NA_real_, + + stringsAsFactors = FALSE + ) + + rows[[length(rows) + 1L]] <- row + } +} + +comparison_tbl <- do.call(rbind, rows) +rownames(comparison_tbl) <- NULL + + +# ── 8. Save and print ─────────────────────────────────="────────────────────── + +saveRDS( + list(joint_results = joint_results, comparison_tbl = comparison_tbl), + OUTPUT_FILE +) +message("\nRDS saved to: ", OUTPUT_FILE) + +write.csv(comparison_tbl, CSV_FILE, row.names = FALSE) +message("CSV saved to: ", CSV_FILE) + + +# ── 9. Print formatted table ────────────────────────────────────────────────── + +message("\n", strrep("=", 80)) +message("COMPARISON: FACTORISED vs JOINT LIKELIHOOD (analysis = '", + ANALYSIS_LABEL, "')") +message(strrep("=", 80)) +message( + "\nColumns:\n", + " pathogen : pathogen name\n", + " dist : distribution\n", + " n_datasets : total number of datasets\n", + " n_type12 : datasets with type-1/2 (order-statistic) summaries\n", + " note : 'identical' when n_type12=0; blank otherwise\n", + " fact_*/joint_* : posterior median (2.5%, 97.5%) for each model\n", + " delta_pred_median : joint minus factorised posterior median of pred_median\n", + " delta_pred_q95 : joint minus factorised posterior median of pred_q95\n" +) + +# Print in two blocks: predicted quantities and parameters. +cols_pred <- c("pathogen", "dist", "n_datasets", "n_type12", + "fact_pred_median", "joint_pred_median", "delta_pred_median", + "fact_pred_q95", "joint_pred_q95", "delta_pred_q95", + "note") +cols_params <- c("pathogen", "dist", "n_datasets", "n_type12", + "fact_mu0", "joint_mu0", + "fact_phi", "joint_phi", + "fact_tau", "joint_tau", + "note") + +message("\n--- Predicted quantities (days) ---\n") +print(comparison_tbl[, intersect(cols_pred, names(comparison_tbl))], + row.names = FALSE, right = FALSE) + +message("\n--- Parameter estimates ---\n") +print(comparison_tbl[, intersect(cols_params, names(comparison_tbl))], + row.names = FALSE, right = FALSE) + +message("\n", strrep("=", 80)) +message("Done.") + + +# ── 10. LaTeX comparison table ──────────────────────────────────────────────── +# +# Layout: two rows per (pathogen × distribution) pair — factorised on top, +# joint below. This keeps the column count to 12, which fits comfortably in +# landscape A4 at \footnotesize with \tabcolsep = 4 pt. +# +# Columns: +# Pathogen | Distribution | N | n_12 | Model | Median | p_95 | +# mu_0 | phi | tau | Delta-Median | Delta-p_95 +# +# Distribution, N, and n_12 use \multirow{2} to span both rows of each pair. +# Pathogen is shown in bold on the first factorised row of each group; the +# cell is left blank for subsequent rows (avoids variable-count multirow spans). +# A thin \cmidrule separates distribution pairs within a pathogen group; +# \midrule separates pathogen groups. +# +# Requires in the LaTeX preamble: +# \usepackage{booktabs} +# \usepackage{longtable} +# \usepackage{multirow} +# \usepackage{pdflscape} +# \usepackage[table]{xcolor} % for \rowcolor % or lscape + +.generate_comparison_latex_table <- function( + tbl, + pathogen_order = NULL, + caption = paste0( + "Comparison of posterior summaries under the factorised and joint ", + "order-statistic likelihoods (\\textit{filtered} analysis). ", + "Each pair of rows corresponds to one pathogen--distribution combination; ", + "the upper row gives results under the factorised model and the lower row ", + "under the joint model. ", + "$N$: total number of datasets; ", + "$n_{12}$: datasets contributing a type-1 (median\\,+\\,range) or ", + "type-2 (median\\,+\\,IQR) summary, i.e.\\ those for which the two ", + "likelihoods differ. ", + "All posterior summaries are reported as median (2.5\\%ile, 97.5\\%ile). ", + "$\\tau$ is shown as \\texttt{---} when $N < 5$ (prior-dominated). ", + "$\\Delta$: joint minus factorised posterior median of the predictive ", + "quantity (days). ", + "$^{\\dagger}$\\,No type-1/2 datasets present; the two likelihoods are ", + "identical and the joint row repeats the factorised estimates." + ), + label = "tab:comparison_factorised_joint" +) { + + # ── Helpers ───────────────────────────────────────────────────────────────── + + .esc <- function(x) { + x <- gsub("_", "\\_", x, fixed = TRUE) + x <- gsub("&", "\\&", x, fixed = TRUE) + x <- gsub("%", "\\%", x, fixed = TRUE) + x + } + + .cell <- function(x) { + if (is.na(x) || x == "") return("---") + x + } + + .delta_fmt <- function(d) { + if (is.na(d)) return("---") + sprintf("%+.2f", d) + } + + # ── Pathogen display labels ────────────────────────────────────────────────── + + plabels <- c( + Nipah = "Nipah", MVD = "Marburg (MVD)", + EVD = "Ebola (EVD)", Lassa = "Lassa fever", + SARS = "SARS", MERS = "MERS", + Zika = "Zika", Measles = "Measles", + Mpox = "Mpox", Cholera = "Cholera", + RVF = "Rift Valley fever", CCHF = "CCHF", + COVID_19 = "COVID-19", Dengue = "Dengue", + YFV = "Yellow fever", Typhoid = "Typhoid", + Smallpox = "Smallpox", Flu = "Influenza" + ) + + if (is.null(pathogen_order)) + pathogen_order <- unique(tbl$pathogen) + + # ── Column spec ────────────────────────────────────────────────────────────── + # Total nominal content width ≈ 194 mm; with \tabcolsep=4pt fits A4 landscape. + + col_spec <- paste0( + "@{} ", + "p{2.5cm} ", # 1 Pathogen + "p{1.6cm} ", # 2 Distribution + "r ", # 3 N + "r ", # 4 n_12 + "l ", # 5 Model + "p{2.3cm} ", # 6 Median (days) + "p{2.3cm} ", # 7 p_95 (days) + "p{2.2cm} ", # 8 mu_0 + "p{3.0cm} ", # 9 phi + "p{2.0cm} ", # 10 tau + "r ", # 11 Delta Median + "r ", # 12 Delta p_95 + "@{}" + ) + + # ── Header rows ────────────────────────────────────────────────────────────── + + span_header <- paste0( + "\\rowcolor[HTML]{F0E68C}", + " & & & & ", + "& \\multicolumn{5}{c}{Posterior summary --- median (95\\% CI)} ", + "& \\multicolumn{2}{c}{$\\Delta$ (joint $-$ fact.)} \\\\" + ) + span_cmidrule <- "\\cmidrule(lr){6-10} \\cmidrule(lr){11-12}" + + col_header <- paste0( + "\\rowcolor[HTML]{F0E68C}", + "Pathogen & Distribution & $N$ & $n_{12}$ & Model ", + "& Median (days) & $p_{95}$ (days) & $\\mu_0$ & $\\phi$ & $\\tau$ ", + "& Median & $p_{95}$ \\\\" + ) + + make_head <- function(is_first) { + c( + if (!is_first) + "\\multicolumn{12}{l}{\\small\\textit{continued from previous page}} \\\\", + "\\toprule", + span_header, + span_cmidrule, + col_header, + "\\midrule" + ) + } + + # ── Table body ─────────────────────────────────────────────────────────────── + + body <- character(0) + first_pathogen <- TRUE + + for (pathogen in pathogen_order) { + + sub <- tbl[tbl$pathogen == pathogen, , drop = FALSE] + if (nrow(sub) == 0L) next + + plabel <- if (!is.na(plabels[pathogen])) plabels[[pathogen]] else pathogen + + if (!first_pathogen) + body <- c(body, "\\midrule") + first_pathogen <- FALSE + + for (i in seq_len(nrow(sub))) { + + row <- sub[i, ] + + is_identical <- grepl("identical", row$note, fixed = FALSE) + joint_unavail <- grepl("unavailable", row$note, fixed = FALSE) + + # ── Factorised row ────────────────────────────────────────────────────── + + # Pathogen label: bold on first dist row only; blank thereafter. + path_cell <- if (i == 1L) + paste0("\\textbf{", .esc(plabel), "}") + else + "" + + fact_cells <- paste( + path_cell, + paste0("\\multirow{2}{*}{", .esc(row$dist), "}"), + paste0("\\multirow{2}{*}{", row$n_datasets, "}"), + paste0("\\multirow{2}{*}{", row$n_type12, "}"), + "Factorised", + .cell(row$fact_pred_median), + .cell(row$fact_pred_q95), + .cell(row$fact_mu0), + .cell(row$fact_phi), + .cell(row$fact_tau), + "", # Delta Median: blank on factorised row + "", # Delta p_95: blank on factorised row + sep = " & " + ) + body <- c(body, paste0(fact_cells, " \\\\")) + + # ── Joint row ────────────────────────────────────────────────────────── + + model_lbl <- if (is_identical) "Joint$^{\\dagger}$" + else if (joint_unavail) "Joint (failed)" + else "Joint" + + j_median <- if (joint_unavail) "---" else .cell(row$joint_pred_median) + j_q95 <- if (joint_unavail) "---" else .cell(row$joint_pred_q95) + j_mu0 <- if (joint_unavail) "---" else .cell(row$joint_mu0) + j_phi <- if (joint_unavail) "---" else .cell(row$joint_phi) + j_tau <- if (joint_unavail) "---" else .cell(row$joint_tau) + + d_median <- if (joint_unavail || is_identical) "---" + else .delta_fmt(row$delta_pred_median) + d_q95 <- if (joint_unavail || is_identical) "---" + else .delta_fmt(row$delta_pred_q95) + + joint_cells <- paste( + "", # Pathogen: blank (multirow from fact row covers it) + "", # Dist: blank + "", # N: blank + "", # n_12: blank + model_lbl, + j_median, j_q95, j_mu0, j_phi, j_tau, + d_median, d_q95, + sep = " & " + ) + body <- c(body, paste0(joint_cells, " \\\\")) + + # Thin rule between distribution pairs within the same pathogen group, + # but not after the last distribution for this pathogen. + if (i < nrow(sub)) + body <- c(body, "\\cmidrule(l{4pt}r{0pt}){2-12}") + + } + } + + # ── Assemble full table ────────────────────────────────────────────────────── + + lines <- c( + "% ============================================================", + "% Comparison table: factorised vs joint order-statistic likelihood", + "% Required LaTeX packages: booktabs, longtable, multirow, pdflscape, xcolor (table)", + "% ============================================================", + "{\\setlength{\\tabcolsep}{4pt}", + "\\begin{landscape}", + "\\footnotesize", + paste0("\\begin{longtable}{", col_spec, "}"), + paste0("\\caption{", caption, "}\\label{", label, "} \\\\"), + make_head(is_first = TRUE), + "\\endfirsthead", + make_head(is_first = FALSE), + "\\endhead", + "\\multicolumn{12}{r}{\\small\\textit{continued on next page}} \\\\", + "\\endfoot", + "\\bottomrule", + paste0( + "\\multicolumn{12}{l}{\\footnotesize", + " $^{\\dagger}$\\,No type-1/2 datasets; joint and factorised", + " likelihoods are identical.} \\\\" + ), + "\\endlastfoot", + body, + "\\end{longtable}", + "\\end{landscape}", + "} % end \\setlength{\\tabcolsep}" + ) + + paste(lines, collapse = "\n") +} + + +# ── Generate, write, and echo the table ────────────────────────────────────── + +latex_table <- .generate_comparison_latex_table(comparison_tbl) + +latex_file <- file.path(OUTPUT_DIR, "comparison_factorised_vs_joint.tex") +writeLines(latex_table, latex_file) +message("\nLaTeX table written to: ", latex_file) + +# Echo to console so it can be copied directly. +message("\n", strrep("-", 80)) +message("LaTeX table (copy into document):\n") +cat(latex_table, "\n") +message(strrep("-", 80)) diff --git a/analysis/main_analysis.R b/analysis/main_analysis.R index a59ce83..9c5d2ab 100644 --- a/analysis/main_analysis.R +++ b/analysis/main_analysis.R @@ -195,9 +195,18 @@ subgroup_config <- list( # ── 5. Compile Stan model (once) ────────────────────────────────────────────── - -stan_file <- system.file("stan", "hierarchical_data_synthesis_summary_stats.stan", - package = "ddsynth") +# Set to "factorised" (default) or "joint" to switch between Stan models. +STAN_MODEL <- "factorised" +if (!STAN_MODEL %in% c("factorised", "joint")) + stop('STAN_MODEL must be "factorised" or "joint", got: "', STAN_MODEL, '"') + +stan_model_file <- switch(STAN_MODEL, + factorised = "hierarchical_data_synthesis_summary_stats.stan", + joint = "hierarchical_data_synthesis_summary_stats_joint.stan" +) +stan_file <- system.file("stan", stan_model_file, package = "ddsynth") +if (!nzchar(stan_file)) + stop("Stan file '", stan_model_file, "' not found in ddsynth installation.") stan_model <- rstan::stan_model(stan_file) diff --git a/analysis/new_simulation_study.R b/analysis/new_simulation_study.R index f2daa59..47933ca 100644 --- a/analysis/new_simulation_study.R +++ b/analysis/new_simulation_study.R @@ -38,6 +38,11 @@ RESULTS_DIR <- file.path(MAIN_FOLDER, "results", "new_scenarios") N_SIM <- 100 # replications per scenario SEED <- 123 # global RNG seed (per-sim seeds derived from this) +# STAN_MODEL: "factorised" (default) or "joint". +STAN_MODEL <- "factorised" +if (!STAN_MODEL %in% c("factorised", "joint")) + stop('STAN_MODEL must be "factorised" or "joint", got: "', STAN_MODEL, '"') + # FORCE_RERUN: FALSE = resume; TRUE = rerun all; character vector = rerun named. FORCE_RERUN <- FALSE @@ -64,9 +69,12 @@ dir.create(RESULTS_DIR, recursive = TRUE, showWarnings = F dir.create(file.path(MAIN_FOLDER, "results"), recursive = TRUE, showWarnings = FALSE) cat("Compiling Stan model...\n") +stan_model_file <- switch(STAN_MODEL, + factorised = "hierarchical_data_synthesis_summary_stats.stan", + joint = "hierarchical_data_synthesis_summary_stats_joint.stan" +) stan_model <- rstan::stan_model( - file.path(MAIN_FOLDER, "inst", "stan", - "hierarchical_data_synthesis_summary_stats.stan") + file.path(MAIN_FOLDER, "inst", "stan", stan_model_file) ) cat("Done.\n\n") diff --git a/inst/stan/hierarchical_data_synthesis_summary_stats_joint.stan b/inst/stan/hierarchical_data_synthesis_summary_stats_joint.stan new file mode 100644 index 0000000..3b3070e --- /dev/null +++ b/inst/stan/hierarchical_data_synthesis_summary_stats_joint.stan @@ -0,0 +1,712 @@ +// hierarchical_data_synthesis_summary_stats_joint.stan +// +// Identical to hierarchical_data_synthesis_summary_stats.stan except that +// summary types 1 (median+range) and 2 (median+IQR) use the joint density of +// three order statistics rather than the product of their marginals. +// +// The joint density of order statistics X_(i1) < X_(i2) < X_(i3) from n IID +// observations with CDF F and density f is: +// +// f(x1,x2,x3) = n! / [(i1-1)! (i2-i1-1)! (i3-i2-1)! (n-i3)!] +// * F(x1)^{i1-1} * f(x1) +// * [F(x2)-F(x1)]^{i2-i1-1} * f(x2) +// * [F(x3)-F(x2)]^{i3-i2-1} * f(x3) +// * [1-F(x3)]^{n-i3} +// +// Gap terms with zero exponent are omitted to avoid 0*(-Inf) = NaN. +// Gaps are computed via log_diff_exp for numerical stability. + +functions { + // CDF for each distribution + real dist_cdf_fun(real x, int dist_type, real loc, real phi, real kappa) { + if (dist_type == 1) { // lognormal + return lognormal_cdf(x | loc, phi); + } else if (dist_type == 2) { // gamma + real mean_d = exp(loc); + real shape = phi; + real rate = shape / mean_d; + return gamma_cdf(x | shape, rate); + } else if (dist_type == 3) { // weibull + real scale = exp(loc); + real shape = phi; + return weibull_cdf(x | shape, scale); + } else if (dist_type == 4) { // burr XII: lambda=exp(loc), c=phi, k=kappa + real u = log(x) - loc; + return -expm1(-kappa * log1p_exp(phi * u)); + } else if (dist_type == 5) { // generalised gamma (Prentice): mu=loc, sigma=phi, Q=kappa + real gamma_shape = 1.0 / (kappa * kappa); + real w = (log(x) - loc) / phi; + return gamma_cdf(gamma_shape * exp(kappa * w) | gamma_shape, 1); + } + return 0; + } + + // "log pdf" for each distribution + real dist_logpdf_fun(real x, int dist_type, real loc, real phi, real kappa) { + if (dist_type == 1) { // lognormal + return lognormal_lpdf(x | loc, phi); + } else if (dist_type == 2) { // gamma + real mean_d = exp(loc); + real shape = phi; + real rate = shape / mean_d; + return gamma_lpdf(x | shape, rate); + } else if (dist_type == 3) { // weibull + real scale = exp(loc); + real shape = phi; + return weibull_lpdf(x | shape, scale); + } else if (dist_type == 4) { // burr XII: lambda=exp(loc), c=phi, k=kappa + real u = log(x) - loc; + real log_term = log1p_exp(phi * u); // log(1 + (x/lambda)^c) + return log(phi) + log(kappa) + (phi - 1) * u - loc - (kappa + 1) * log_term; + } else if (dist_type == 5) { // generalised gamma (Prentice): mu=loc, sigma=phi, Q=kappa + real gamma_shape = 1.0 / (kappa * kappa); + real w = (log(x) - loc) / phi; + return log(kappa) - log(phi) - log(x) + + gamma_shape * log(gamma_shape) + + gamma_shape * kappa * w + - gamma_shape * exp(kappa * w) + - lgamma(gamma_shape); + } + return 0; + } + + // Log-scale CDF for each distribution. + real dist_log_cdf_fun(real x, int dist_type, real loc, real phi, real kappa) { + if (dist_type == 1) { // lognormal + return lognormal_lcdf(x | loc, phi); + } else if (dist_type == 2) { // gamma + real mean_d = exp(loc); + real rate = phi / mean_d; + return gamma_lcdf(x | phi, rate); + } else if (dist_type == 3) { // weibull + real scale = exp(loc); + return weibull_lcdf(x | phi, scale); + } else if (dist_type == 4) { // burr XII: lambda=exp(loc), c=phi, k=kappa + real u = log(x) - loc; + return log1m_exp(-kappa * log1p_exp(phi * u)); + } else if (dist_type == 5) { // generalised gamma (Prentice): mu=loc, sigma=phi, Q=kappa + real gamma_shape = 1.0 / (kappa * kappa); + real w = (log(x) - loc) / phi; + return gamma_lcdf(gamma_shape * exp(kappa * w) | gamma_shape, 1); + } + return negative_infinity(); + } + + // Log-scale CCDF (log survivor function) for each distribution. + real dist_log_ccdf_fun(real x, int dist_type, real loc, real phi, real kappa) { + if (dist_type == 1) { // lognormal + return lognormal_lccdf(x | loc, phi); + } else if (dist_type == 2) { // gamma + real mean_d = exp(loc); + real rate = phi / mean_d; + return gamma_lccdf(x | phi, rate); + } else if (dist_type == 3) { // weibull + real scale = exp(loc); + return weibull_lccdf(x | phi, scale); + } else if (dist_type == 4) { // burr XII: log CCDF = -k * log(1 + (x/lambda)^c) + real u = log(x) - loc; + return -kappa * log1p_exp(phi * u); + } else if (dist_type == 5) { // generalised gamma (Prentice) + real gamma_shape = 1.0 / (kappa * kappa); + real w = (log(x) - loc) / phi; + return gamma_lccdf(gamma_shape * exp(kappa * w) | gamma_shape, 1); + } + return negative_infinity(); + } + + // Log likelihood for order statistic: + // k-th order statistic out of n observations. + real order_stat_logpdf_fun(real x, int n, int k, int dist_type, real loc, real phi, real kappa) { + real log_f; + real log_F = 0.0; + real log_1mF = 0.0; + + if (dist_type == 1) { // lognormal + log_f = lognormal_lpdf(x | loc, phi); + log_F = lognormal_lcdf(x | loc, phi); + log_1mF = lognormal_lccdf(x | loc, phi); + + } else if (dist_type == 2) { // gamma + real mean_d = exp(loc); + real rate = phi / mean_d; + log_f = gamma_lpdf(x | phi, rate); + log_F = gamma_lcdf(x | phi, rate); + log_1mF = gamma_lccdf(x | phi, rate); + + } else if (dist_type == 3) { // weibull + real scale = exp(loc); + log_f = weibull_lpdf(x | phi, scale); + log_F = weibull_lcdf(x | phi, scale); + log_1mF = weibull_lccdf(x | phi, scale); + + } else if (dist_type == 4) { // burr XII: lambda=exp(loc), c=phi, k=kappa + real u = log(x) - loc; + real log_term = log1p_exp(phi * u); + log_f = log(phi) + log(kappa) + (phi - 1) * u - loc - (kappa + 1) * log_term; + log_1mF = -kappa * log_term; + log_F = log1m_exp(log_1mF); + + } else if (dist_type == 5) { // generalised gamma (Prentice) + real gamma_shape = 1.0 / (kappa * kappa); + real w = (log(x) - loc) / phi; + real arg = gamma_shape * exp(kappa * w); + log_f = log(kappa) - log(phi) - log(x) + + gamma_shape * log(gamma_shape) + + gamma_shape * kappa * w + - arg + - lgamma(gamma_shape); + if (k > 1) log_F = gamma_lcdf(arg | gamma_shape, 1); + if (k < n) log_1mF = gamma_lccdf(arg | gamma_shape, 1); + } + + real log_dens = lchoose(n, k) + log_f; + if (k > 1) log_dens += (k - 1) * log_F; + if (k < n) log_dens += (n - k) * log_1mF; + return log_dens; + } + + // Joint log density of three order statistics x1 < x2 < x3 with ranks + // i1 < i2 < i3 from a sample of size n. + // + // log f(x1,x2,x3) = log_nc + // + (i1-1)*log_F(x1) [omitted when i1==1] + // + log_f(x1) + // + (i2-i1-1)*log[F(x2)-F(x1)] [omitted when i2-i1==1] + // + log_f(x2) + // + (i3-i2-1)*log[F(x3)-F(x2)] [omitted when i3-i2==1] + // + log_f(x3) + // + (n-i3)*log[1-F(x3)] [omitted when i3==n] + // + // where log_nc = lgamma(n+1) - lgamma(i1) - lgamma(i2-i1) - lgamma(i3-i2) + // - lgamma(n-i3+1) + real joint_3_order_stat_logpdf( + real x1, real x2, real x3, + int i1, int i2, int i3, + int n, + int dist_type, real loc, real phi, real kappa + ) { + if (x1 >= x2 || x2 >= x3) return negative_infinity(); + + real log_nc = lgamma(n + 1) + - lgamma(i1) + - lgamma(i2 - i1) + - lgamma(i3 - i2) + - lgamma(n - i3 + 1); + + real log_f1 = dist_logpdf_fun(x1, dist_type, loc, phi, kappa); + real log_f2 = dist_logpdf_fun(x2, dist_type, loc, phi, kappa); + real log_f3 = dist_logpdf_fun(x3, dist_type, loc, phi, kappa); + + real log_F1 = dist_log_cdf_fun(x1, dist_type, loc, phi, kappa); + real log_F2 = dist_log_cdf_fun(x2, dist_type, loc, phi, kappa); + real log_F3 = dist_log_cdf_fun(x3, dist_type, loc, phi, kappa); + real log_1mF3 = dist_log_ccdf_fun(x3, dist_type, loc, phi, kappa); + + real lp = log_nc + log_f1 + log_f2 + log_f3; + + if (i1 > 1) lp += (i1 - 1) * log_F1; + if (i2 - i1 > 1) lp += (i2 - i1 - 1) * log_diff_exp(log_F2, log_F1); + if (i3 - i2 > 1) lp += (i3 - i2 - 1) * log_diff_exp(log_F3, log_F2); + if (i3 < n) lp += (n - i3) * log_1mF3; + + return lp; + } + + // Helper function to compute gamma quantile approximation + real gamma_quantile_approx(real p, real shape, real scale) { + real z = inv_Phi(p); + if (shape > 1) { + return shape * scale * pow(1 - 1.0/(9*shape) + z/(3*sqrt(shape)), 3); + } else { + return shape * scale * (1 + z / sqrt(shape)); + } + } +} + +data { + int n_datasets; + array[n_datasets] int n_obs; + array[n_datasets] int summary_type; + int dist_type; + + array[n_datasets] real obs_stat1; + array[n_datasets] real obs_stat2; + array[n_datasets] real obs_stat3; + + int n_freq_total; + array[n_freq_total] real freq_value; + array[n_freq_total] int freq_count; + array[n_datasets] int freq_start; + array[n_datasets] int freq_len; + + array[n_freq_total] real freq_lower; + array[n_freq_total] real freq_upper; + + real mu0_mean; + real mu0_sd; + + real log_tau_mean; + real log_tau_sd; + + real log_phi_mean; + real log_phi_sd; + + real log_kappa_mean; + real log_kappa_sd; +} + +transformed data { + real z_q25 = -0.6745; + real z_q75 = 0.6745; + real z_q90 = 1.28155; + real z_q95 = 1.64485; + + for (d in 1:n_datasets) { + if (summary_type[d] == 1) { + if (obs_stat2[d] >= obs_stat1[d] || obs_stat1[d] >= obs_stat3[d]) { + reject("For summary_type=1, must have min < median < max"); + } + } else if (summary_type[d] == 2) { + if (obs_stat2[d] >= obs_stat1[d] || obs_stat1[d] >= obs_stat3[d]) { + reject("For summary_type=2, must have q25 < median < q75"); + } + } else if (summary_type[d] == 4) { + if (freq_len[d] == 0) { + reject("For summary_type=4, freq_len must be > 0"); + } + } else if (summary_type[d] == 5) { + if (freq_len[d] == 0) { + reject("For summary_type=5, freq_len must be > 0"); + } + } + } +} + +parameters { + real mu0; + real log_tau; + real log_phi; + real log_kappa; + vector[n_datasets] loc_d_raw; +} + +transformed parameters { + real tau = exp(log_tau); + real phi = exp(log_phi); + real kappa = exp(log_kappa); + vector[n_datasets] loc_d = mu0 + tau * loc_d_raw; +} + +model { + mu0 ~ normal(mu0_mean, mu0_sd); + log_tau ~ normal(log_tau_mean, log_tau_sd); + log_phi ~ normal(log_phi_mean, log_phi_sd); + log_kappa ~ normal(log_kappa_mean, log_kappa_sd); + + loc_d_raw ~ std_normal(); + + for (d in 1:n_datasets) { + real loc = loc_d[d]; + int n = n_obs[d]; + + if (summary_type[d] == 1) { // median + range: joint density of (min, median, max) + int k_median = (n + 1) %/% 2; + target += joint_3_order_stat_logpdf( + obs_stat2[d], obs_stat1[d], obs_stat3[d], // x1=min, x2=median, x3=max + 1, k_median, n, // i1=1, i2=k_median, i3=n + n, dist_type, loc, phi, kappa); + } + + else if (summary_type[d] == 2) { // median + IQR: joint density of (q25, median, q75) + int k_median = (n + 1) %/% 2; + int k_q25 = (n + 1) %/% 4; + if (k_q25 < 1) k_q25 = 1; + int k_q75 = (3 * (n + 1)) %/% 4; + if (k_q75 <= k_q25) k_q75 = k_q25 + 1; + if (k_q75 > n) k_q75 = n; + target += joint_3_order_stat_logpdf( + obs_stat2[d], obs_stat1[d], obs_stat3[d], // x1=q25, x2=median, x3=q75 + k_q25, k_median, k_q75, + n, dist_type, loc, phi, kappa); + } + + else if (summary_type[d] == 3) { // mean + sd (unchanged from factorised model) + real expected_mean; + real expected_sd; + real se_mean; + real se_sd; + + if (dist_type == 1) { // lognormal + expected_mean = exp(loc + phi^2 / 2); + real var_ = (exp(phi^2) - 1) * exp(2 * loc + phi^2); + expected_sd = sqrt(var_); + se_mean = expected_sd / sqrt(n); + se_sd = expected_sd / sqrt(2 * (n - 1)); + + } else if (dist_type == 2) { // gamma + real mean_d = exp(loc); + real shape = phi; + real scale_param = mean_d / shape; + expected_mean = mean_d; + expected_sd = sqrt(shape * scale_param^2); + se_mean = expected_sd / sqrt(n); + se_sd = expected_sd / sqrt(2 * (n - 1)); + + } else if (dist_type == 3) { // weibull + real scale = exp(loc); + real shape = phi; + expected_mean = scale * tgamma(1 + 1.0 / shape); + real var_ = scale^2 * (tgamma(1 + 2.0 / shape) - pow(tgamma(1 + 1.0 / shape), 2)); + expected_sd = sqrt(var_); + se_mean = expected_sd / sqrt(n); + se_sd = expected_sd / sqrt(2 * (n - 1)); + + } else if (dist_type == 4) { // burr XII + if (kappa * phi <= 2.0) { + target += negative_infinity(); + } else { + real lam = exp(loc); + expected_mean = lam * kappa * exp(lbeta(kappa - 1.0/phi, 1.0 + 1.0/phi)); + real e2 = lam^2 * kappa * exp(lbeta(kappa - 2.0/phi, 1.0 + 2.0/phi)); + expected_sd = sqrt(fabs(e2 - expected_mean^2)); + se_mean = expected_sd / sqrt(n); + se_sd = expected_sd / sqrt(2 * (n - 1)); + obs_stat1[d] ~ normal(expected_mean, se_mean); + obs_stat2[d] ~ normal(expected_sd, se_sd); + } + + } else if (dist_type == 5) { // generalised gamma + real gamma_shape = 1.0 / (kappa * kappa); + real log_ET = loc + 2.0*phi/kappa * log(kappa) + + lgamma(gamma_shape + phi/kappa) - lgamma(gamma_shape); + real log_ET2 = 2.0*loc + 4.0*phi/kappa * log(kappa) + + lgamma(gamma_shape + 2.0*phi/kappa) - lgamma(gamma_shape); + expected_mean = exp(log_ET); + expected_sd = sqrt(fabs(exp(log_ET2) - expected_mean^2)); + se_mean = expected_sd / sqrt(n); + se_sd = expected_sd / sqrt(2 * (n - 1)); + obs_stat1[d] ~ normal(expected_mean, se_mean); + obs_stat2[d] ~ normal(expected_sd, se_sd); + } + + if (dist_type <= 3) { + obs_stat1[d] ~ normal(expected_mean, se_mean); + obs_stat2[d] ~ normal(expected_sd, se_sd); + } + } + + else if (summary_type[d] == 4) { // raw frequency table (unchanged) + int s = freq_start[d]; + int len = freq_len[d]; + for (i in s:(s + len - 1)) { + target += freq_count[i] * dist_logpdf_fun(freq_value[i], dist_type, loc, phi, kappa); + } + } + + else if (summary_type[d] == 5) { // interval-censored frequency table (unchanged) + int s = freq_start[d]; + int len = freq_len[d]; + for (i in s:(s + len - 1)) { + if (freq_lower[i] == freq_upper[i]) { + target += freq_count[i] * dist_logpdf_fun(freq_lower[i], dist_type, loc, phi, kappa); + } else { + real log_cdf_u = dist_log_cdf_fun(freq_upper[i], dist_type, loc, phi, kappa); + real log_cdf_l = dist_log_cdf_fun(freq_lower[i], dist_type, loc, phi, kappa); + target += freq_count[i] * log_diff_exp(log_cdf_u, log_cdf_l); + } + } + } + } +} + +generated quantities { + real pred_mean; + real pred_median; + real pred_q25; + real pred_q75; + real pred_q90; + real pred_q95; + real pred_sd; + // Joint log-likelihood per dataset (3 slots each for interface compatibility; + // for types 1 and 2 the joint value is stored in slot idx and slots idx+1, + // idx+2 are set to 0, since the joint density cannot be decomposed per statistic). + vector[n_datasets * 3] log_lik; + + { + if (n_datasets < 5) { + real loc_pred = mean(loc_d); + + if (dist_type == 1) { // lognormal + pred_mean = exp(loc_pred + phi^2 / 2); + pred_median = exp(loc_pred); + pred_q25 = exp(loc_pred - 0.6745 * phi); + pred_q75 = exp(loc_pred + 0.6745 * phi); + pred_q90 = exp(loc_pred + 1.28155 * phi); + pred_q95 = exp(loc_pred + 1.64485 * phi); + pred_sd = sqrt((exp(phi^2) - 1) * exp(2 * loc_pred + phi^2)); + + } else if (dist_type == 2) { // gamma + real mean_d = exp(loc_pred); + real scale_param = mean_d / phi; + pred_mean = mean_d; + pred_sd = sqrt(mean_d * scale_param); + pred_median = gamma_quantile_approx(0.5, phi, scale_param); + pred_q25 = gamma_quantile_approx(0.25, phi, scale_param); + pred_q75 = gamma_quantile_approx(0.75, phi, scale_param); + pred_q90 = gamma_quantile_approx(0.90, phi, scale_param); + pred_q95 = gamma_quantile_approx(0.95, phi, scale_param); + + } else if (dist_type == 3) { // weibull + real scale = exp(loc_pred); + pred_mean = scale * tgamma(1 + 1.0 / phi); + pred_median = scale * pow(log(2), 1.0 / phi); + pred_q25 = scale * pow(log(4.0 / 3.0), 1.0 / phi); + pred_q75 = scale * pow(log(4.0), 1.0 / phi); + pred_q90 = scale * pow(log(10.0), 1.0 / phi); + pred_q95 = scale * pow(log(20.0), 1.0 / phi); + pred_sd = sqrt(scale^2 * (tgamma(1 + 2.0/phi) - pow(tgamma(1 + 1.0/phi), 2))); + + } else if (dist_type == 4) { // burr XII + real lam = exp(loc_pred); + pred_median = lam * pow(pow(0.5, -1.0/kappa) - 1.0, 1.0/phi); + pred_q25 = lam * pow(pow(0.75, -1.0/kappa) - 1.0, 1.0/phi); + pred_q75 = lam * pow(pow(0.25, -1.0/kappa) - 1.0, 1.0/phi); + pred_q90 = lam * pow(pow(0.10, -1.0/kappa) - 1.0, 1.0/phi); + pred_q95 = lam * pow(pow(0.05, -1.0/kappa) - 1.0, 1.0/phi); + if (kappa * phi > 1.0) { + pred_mean = lam * kappa * exp(lbeta(kappa - 1.0/phi, 1.0 + 1.0/phi)); + } else { + pred_mean = positive_infinity(); + } + if (kappa * phi > 2.0) { + real e2 = lam^2 * kappa * exp(lbeta(kappa - 2.0/phi, 1.0 + 2.0/phi)); + pred_sd = sqrt(fabs(e2 - pred_mean^2)); + } else { + pred_sd = positive_infinity(); + } + + } else if (dist_type == 5) { // generalised gamma: Monte Carlo via gamma_rng + int L_gg = 200; + vector[L_gg] gg_samples; + real gs = 1.0 / (kappa * kappa); + for (l in 1:L_gg) { + real y = gamma_rng(gs, 1); + gg_samples[l] = exp(loc_pred + phi / kappa * log(kappa * kappa * y)); + } + gg_samples = sort_asc(gg_samples); + pred_mean = mean(gg_samples); + pred_sd = sd(gg_samples); + pred_median = gg_samples[100]; + pred_q25 = gg_samples[50]; + pred_q75 = gg_samples[150]; + pred_q90 = gg_samples[180]; + pred_q95 = gg_samples[190]; + } + + } else { + int L = 2000; + + if (dist_type == 5) { + int L_gg = 200; + vector[L_gg] gg_t_samples; + real gs = 1.0 / (kappa * kappa); + for (l in 1:L_gg) { + real loc_sample = normal_rng(mu0, tau); + real y = gamma_rng(gs, 1); + gg_t_samples[l] = exp(loc_sample + phi / kappa * log(kappa * kappa * y)); + } + gg_t_samples = sort_asc(gg_t_samples); + pred_mean = mean(gg_t_samples); + pred_sd = sd(gg_t_samples); + pred_median = gg_t_samples[100]; + pred_q25 = gg_t_samples[50]; + pred_q75 = gg_t_samples[150]; + pred_q90 = gg_t_samples[180]; + pred_q95 = gg_t_samples[190]; + + } else { + vector[L] means; + vector[L] medians; + vector[L] q25s; + vector[L] q75s; + vector[L] q90s; + vector[L] q95s; + vector[L] sds; + + for (l in 1:L) { + real loc_sample = normal_rng(mu0, tau); + + if (dist_type == 1) { // lognormal + means[l] = exp(loc_sample + phi^2 / 2); + medians[l] = exp(loc_sample); + q25s[l] = exp(loc_sample - 0.6745 * phi); + q75s[l] = exp(loc_sample + 0.6745 * phi); + q90s[l] = exp(loc_sample + 1.28155 * phi); + q95s[l] = exp(loc_sample + 1.64485 * phi); + sds[l] = sqrt((exp(phi^2) - 1) * exp(2 * loc_sample + phi^2)); + + } else if (dist_type == 2) { // gamma + real mean_d = exp(loc_sample); + real scale_param = mean_d / phi; + means[l] = mean_d; + sds[l] = sqrt(mean_d * scale_param); + medians[l] = gamma_quantile_approx(0.5, phi, scale_param); + q25s[l] = gamma_quantile_approx(0.25, phi, scale_param); + q75s[l] = gamma_quantile_approx(0.75, phi, scale_param); + q90s[l] = gamma_quantile_approx(0.90, phi, scale_param); + q95s[l] = gamma_quantile_approx(0.95, phi, scale_param); + + } else if (dist_type == 3) { // weibull + real scale = exp(loc_sample); + means[l] = scale * tgamma(1 + 1.0 / phi); + medians[l] = scale * pow(log(2), 1.0 / phi); + q25s[l] = scale * pow(log(4.0 / 3.0), 1.0 / phi); + q75s[l] = scale * pow(log(4.0), 1.0 / phi); + q90s[l] = scale * pow(log(10.0), 1.0 / phi); + q95s[l] = scale * pow(log(20.0), 1.0 / phi); + real var_weib = scale^2 * (tgamma(1 + 2.0/phi) - pow(tgamma(1 + 1.0/phi), 2)); + sds[l] = sqrt(var_weib); + + } else if (dist_type == 4) { // burr XII + real lam_l = exp(loc_sample); + medians[l] = lam_l * pow(pow(0.5, -1.0/kappa) - 1.0, 1.0/phi); + q25s[l] = lam_l * pow(pow(0.75, -1.0/kappa) - 1.0, 1.0/phi); + q75s[l] = lam_l * pow(pow(0.25, -1.0/kappa) - 1.0, 1.0/phi); + q90s[l] = lam_l * pow(pow(0.10, -1.0/kappa) - 1.0, 1.0/phi); + q95s[l] = lam_l * pow(pow(0.05, -1.0/kappa) - 1.0, 1.0/phi); + if (kappa * phi > 1.0) { + means[l] = lam_l * kappa * exp(lbeta(kappa - 1.0/phi, 1.0 + 1.0/phi)); + } else { + means[l] = positive_infinity(); + } + if (kappa * phi > 2.0) { + real e2_l = lam_l^2 * kappa * exp(lbeta(kappa - 2.0/phi, 1.0 + 2.0/phi)); + sds[l] = (kappa * phi > 1.0) ? sqrt(fabs(e2_l - means[l]^2)) : positive_infinity(); + } else { + sds[l] = positive_infinity(); + } + } + } + + pred_mean = mean(means); + pred_median = mean(medians); + pred_q25 = mean(q25s); + pred_q75 = mean(q75s); + pred_q90 = mean(q90s); + pred_q95 = mean(q95s); + pred_sd = mean(sds); + } + } + } + + // Log likelihood + { + int idx = 1; + + for (d in 1:n_datasets) { + real loc = loc_d[d]; + int n = n_obs[d]; + + if (summary_type[d] == 1) { // joint density of (min, median, max) + int k_median = (n + 1) %/% 2; + log_lik[idx] = joint_3_order_stat_logpdf( + obs_stat2[d], obs_stat1[d], obs_stat3[d], + 1, k_median, n, + n, dist_type, loc, phi, kappa); + log_lik[idx + 1] = 0; + log_lik[idx + 2] = 0; + + } else if (summary_type[d] == 2) { // joint density of (q25, median, q75) + int k_median = (n + 1) %/% 2; + int k_q25 = (n + 1) %/% 4; + if (k_q25 < 1) k_q25 = 1; + int k_q75 = (3 * (n + 1)) %/% 4; + if (k_q75 <= k_q25) k_q75 = k_q25 + 1; + if (k_q75 > n) k_q75 = n; + log_lik[idx] = joint_3_order_stat_logpdf( + obs_stat2[d], obs_stat1[d], obs_stat3[d], + k_q25, k_median, k_q75, + n, dist_type, loc, phi, kappa); + log_lik[idx + 1] = 0; + log_lik[idx + 2] = 0; + + } else if (summary_type[d] == 3) { // mean + sd (unchanged) + real expected_mean; + real expected_sd; + + if (dist_type == 1) { + expected_mean = exp(loc + phi^2 / 2); + real var_ = (exp(phi^2) - 1) * exp(2 * loc + phi^2); + expected_sd = sqrt(var_); + } else if (dist_type == 2) { + real mean_d = exp(loc); + real shape = phi; + real scale_param = mean_d / shape; + expected_mean = mean_d; + expected_sd = sqrt(mean_d * scale_param); + } else if (dist_type == 3) { + real scale = exp(loc); + real shape = phi; + expected_mean = scale * tgamma(1 + 1.0 / shape); + real var_ = scale^2 * (tgamma(1 + 2.0 / shape) - pow(tgamma(1 + 1.0 / shape), 2)); + expected_sd = sqrt(var_); + } else if (dist_type == 4) { + if (kappa * phi > 2.0) { + real lam = exp(loc); + expected_mean = lam * kappa * exp(lbeta(kappa - 1.0/phi, 1.0 + 1.0/phi)); + real e2 = lam^2 * kappa * exp(lbeta(kappa - 2.0/phi, 1.0 + 2.0/phi)); + expected_sd = sqrt(fabs(e2 - expected_mean^2)); + } else { + log_lik[idx] = negative_infinity(); + log_lik[idx + 1] = negative_infinity(); + log_lik[idx + 2] = 0; + } + } else if (dist_type == 5) { + real gamma_shape = 1.0 / (kappa * kappa); + real log_ET = loc + 2.0*phi/kappa * log(kappa) + + lgamma(gamma_shape + phi/kappa) - lgamma(gamma_shape); + real log_ET2 = 2.0*loc + 4.0*phi/kappa * log(kappa) + + lgamma(gamma_shape + 2.0*phi/kappa) - lgamma(gamma_shape); + expected_mean = exp(log_ET); + expected_sd = sqrt(fabs(exp(log_ET2) - expected_mean^2)); + } + + if (dist_type != 4 || kappa * phi > 2.0) { + log_lik[idx] = normal_lpdf(obs_stat1[d] | expected_mean, expected_sd / sqrt(n)); + log_lik[idx + 1] = normal_lpdf(obs_stat2[d] | expected_sd, expected_sd / sqrt(2 * (n - 1))); + log_lik[idx + 2] = 0; + } + + } else if (summary_type[d] == 4) { // raw frequency table + real ll_type4 = 0; + int s = freq_start[d]; + int len = freq_len[d]; + for (i in s:(s + len - 1)) { + ll_type4 += freq_count[i] * dist_logpdf_fun(freq_value[i], dist_type, loc, phi, kappa); + } + log_lik[idx] = ll_type4; + log_lik[idx + 1] = 0; + log_lik[idx + 2] = 0; + + } else if (summary_type[d] == 5) { // interval-censored frequency table + real ll_type5 = 0; + int s = freq_start[d]; + int len = freq_len[d]; + for (i in s:(s + len - 1)) { + if (freq_lower[i] == freq_upper[i]) { + ll_type5 += freq_count[i] * dist_logpdf_fun(freq_lower[i], dist_type, loc, phi, kappa); + } else { + real log_cdf_u = dist_log_cdf_fun(freq_upper[i], dist_type, loc, phi, kappa); + real log_cdf_l = dist_log_cdf_fun(freq_lower[i], dist_type, loc, phi, kappa); + ll_type5 += freq_count[i] * log_diff_exp(log_cdf_u, log_cdf_l); + } + } + log_lik[idx] = ll_type5; + log_lik[idx + 1] = 0; + log_lik[idx + 2] = 0; + } + + idx += 3; + } + } +} diff --git a/man/compile_stan_model.Rd b/man/compile_stan_model.Rd new file mode 100644 index 0000000..051dd81 --- /dev/null +++ b/man/compile_stan_model.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils.R +\name{compile_stan_model} +\alias{compile_stan_model} +\title{Compile a ddsynth Stan model} +\usage{ +compile_stan_model(model = c("factorised", "joint")) +} +\arguments{ +\item{model}{Character string: \code{"factorised"} (default) uses the +order-statistic factorised likelihood +(\code{hierarchical_data_synthesis_summary_stats.stan}); \code{"joint"} uses the +joint parameterisation +(\code{hierarchical_data_synthesis_summary_stats_joint.stan}).} +} +\value{ +A compiled \code{stanmodel} object. +} +\description{ +Compiles and returns one of the Stan models shipped with the package. +Pass the returned object to \code{\link[=pre_inference_checks]{pre_inference_checks()}}, \code{\link[=fit_model]{fit_model()}}, or +\code{\link[=run_simulation_study_generalized]{run_simulation_study_generalized()}}. +} diff --git a/tests/testthat/test-utils.R b/tests/testthat/test-utils.R index 6e968b5..28ababd 100644 --- a/tests/testthat/test-utils.R +++ b/tests/testthat/test-utils.R @@ -66,3 +66,17 @@ test_that("make_stan_init_fn sets loc_d_raw to zero vector of length n_datasets" init <- make_stan_init_fn(sd)() expect_equal(init$loc_d_raw, rep(0.0, 3L)) }) + +# compile_stan_model ---------------------------------------------------------- + +test_that("compile_stan_model rejects invalid model names", { + expect_error(compile_stan_model("invalid"), "should be one of") + expect_error(compile_stan_model("FACTORISED"), "should be one of") +}) + +test_that("compile_stan_model accepts valid model names without error on arg check", { + # match.arg() resolves partial and exact matches before any file I/O; + # test that argument parsing succeeds for both valid choices. + expect_true(match.arg("factorised", c("factorised", "joint")) == "factorised") + expect_true(match.arg("joint", c("factorised", "joint")) == "joint") +}) diff --git a/vignettes/new_simulation_results_all.rds b/vignettes/new_simulation_results_all.rds new file mode 100644 index 0000000..3917333 Binary files /dev/null and b/vignettes/new_simulation_results_all.rds differ