Model Diagnostics • bayescoveragemodel

This article describes how to assess model fit and convergence for Bayesian transition models. Results can also be analysed using the bayescoveragemodelchecks package.

Setup

library(bayescoveragemodel)
library(dplyr)
library(ggplot2)
library(localhierarchy)

Load Fitted Models

Load previously fitted models for diagnostics:

# Choose indicator and model type
indicator <- "anc4"
model_name <- "spline"  # or "rw2"

# Load step 1ab fit (main global model)
fit <- get_fit(indicator, runstep = "step1ab", folder_suffix = "run_workflow")

# Load validation fit
fit_val <- get_fit(indicator, runstep = "step1ab", folder_suffix = "val_workflow")

Visualize Estimates

Create plots for all countries to visually assess model fit:

# Plot estimates from main fit
plots_main <- plot_estimates_local_all(
  results = fit,
  save_plots = TRUE,
  indicator_name = indicator
)

# Compare main fit vs validation
plots_comparison <- plot_estimates_local_all(
  results = fit,
  results2 = fit_val,
  modelnames = c("Main", "Validation"),
  save_plots = TRUE,
  indicator_name = indicator
)

Convergence Diagnostics

Check MCMC convergence using standard diagnostics:

# Get convergence diagnostics
diag <- get_convergence_diagnostics(fit)
print(diag)

Key diagnostics to check:

Rhat: Should be close to 1.0 (< 1.01 is good)
ESS (Effective Sample Size): Should be sufficiently large

Hierarchical Model Checks

Some model parameters are modeled hierarchically. We can check these parameters as follows:

# Plot hierarchical checks
plot_hierchecks(fit)

These plots help assess:

Prior-posterior comparison for hierarchical hyperparameters
Shrinkage behavior across groups

Compare Model Fits

Compare estimates from different model runs (e.g., spline vs RW2):

# Load both model types
fit_spline <- get_fit(indicator, runstep = "step1ab", folder_suffix = "run_spline")
fit_rw2 <- get_fit(indicator, runstep = "step1ab", folder_suffix = "run_rw2")

# Compare spline and RW2 models
plots_model_comparison <- plot_estimates_local_all(
  results = fit_spline,
  results2 = fit_rw2,
  modelnames = c("Spline", "RW2"),
  save_plots = TRUE,
  indicator_name = indicator
)