diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..3a01feb
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,6 @@
+.Rproj.user
+.Rhistory
+.RData
+results/
+*.parquet
+.DS_Store
diff --git a/pvisim/DESCRIPTION b/pvisim/DESCRIPTION
new file mode 100644
index 0000000..231a0c3
--- /dev/null
+++ b/pvisim/DESCRIPTION
@@ -0,0 +1,24 @@
+Package: pvisim
+Title: Known-Truth Simulator for Target Trial Emulation of Early Peripheral Vascular Intervention
+Version: 0.1.0
+Authors@R: person("Chen", "Dun", email = "chendun1995@gmail.com", role = c("aut", "cre"))
+Description: Data-generating process (claudicants nested in physicians, latent
+    severity, physician skill and style), truths by simulating interventions,
+    target trial emulation designs (naive, landmark, sequential trials,
+    clone-censor-weight), simulated RCTs, physician-level variants, negative
+    controls and performance metrics. Not for CRAN.
+License: MIT
+Encoding: UTF-8
+Depends: R (>= 4.3)
+Imports:
+    stats,
+    utils,
+    yaml,
+    data.table
+Suggests:
+    testthat (>= 3.0.0),
+    fastglm,
+    sandwich,
+    glmmTMB,
+    pkgload
+Config/testthat/edition: 3
diff --git a/pvisim/NAMESPACE b/pvisim/NAMESPACE
new file mode 100644
index 0000000..5aeddab
--- /dev/null
+++ b/pvisim/NAMESPACE
@@ -0,0 +1,3 @@
+exportPattern("^[a-z]")
+import(data.table)
+importFrom(stats, plogis, qlogis, qnorm, qpois, rnorm, runif, rbinom, rlnorm, rgamma, setNames)
diff --git a/pvisim/R/params.R b/pvisim/R/params.R
new file mode 100644
index 0000000..aa8b76d
--- /dev/null
+++ b/pvisim/R/params.R
@@ -0,0 +1,51 @@
+# Parameters: read the YAML, then apply a world and a truth scenario as overrides.
+
+#' Recursively merge `override` into `base` (lists only; scalars replace).
+deep_merge <- function(base, override) {
+  for (nm in names(override)) {
+    if (is.list(override[[nm]]) && is.list(base[[nm]])) {
+      base[[nm]] <- deep_merge(base[[nm]], override[[nm]])
+    } else {
+      base[[nm]] <- override[[nm]]
+    }
+  }
+  base
+}
+
+#' Path of the default parameter file, inside the installed or loaded package.
+default_params_path <- function() {
+  p <- system.file("params", "dgp_default.yml", package = "pvisim")
+  if (p == "") stop("dgp_default.yml not found; load the package with pkgload::load_all('pvisim')")
+  p
+}
+
+#' Load parameters for one world and one truth scenario.
+#'
+#' @param world name in `worlds:` (W0, W1, W1_rho0, W2, W3, W4, W5, WR).
+#' @param scenario name in `scenarios:` (null, harm, benefit_harm).
+#' @param path YAML file; defaults to the package's dgp_default.yml.
+#' @param overrides optional list merged last (used by calibration and tests).
+#' @param calibrated apply the world's entry in `calibration_file` (intercepts solved by
+#'   simulations/01_calibrate.R under the null scenario), after the world and before the
+#'   scenario. Calibration itself runs with FALSE.
+load_params <- function(world = "W0", scenario = "null", path = default_params_path(),
+                        overrides = NULL, calibrated = TRUE,
+                        calibration_file = "simulations/config/calibrated.yml") {
+  raw <- yaml::read_yaml(path)
+  if (!world %in% names(raw$worlds)) stop("unknown world: ", world)
+  if (!scenario %in% names(raw$scenarios)) stop("unknown scenario: ", scenario)
+  par <- raw[setdiff(names(raw), c("worlds", "scenarios"))]
+  par <- deep_merge(par, raw$worlds[[world]])
+  if (calibrated && file.exists(calibration_file)) {
+    cal <- yaml::read_yaml(calibration_file)[[world]]
+    if (!is.null(cal)) par <- deep_merge(par, cal)
+    par$calibrated <- !is.null(cal)
+  } else {
+    par$calibrated <- FALSE
+  }
+  par <- deep_merge(par, raw$scenarios[[scenario]])
+  if (!is.null(overrides)) par <- deep_merge(par, overrides)
+  par$world <- world
+  par$scenario <- scenario
+  par
+}
diff --git a/pvisim/R/records.R b/pvisim/R/records.R
new file mode 100644
index 0000000..305db8d
--- /dev/null
+++ b/pvisim/R/records.R
@@ -0,0 +1,69 @@
+# What the analyst sees: simulated Medicare claims (plan D17).
+# Latent quantities (Z_t, true smoking, frailty, skill, style, progression rate) never
+# leave this function except in `hidden`, which analysis code must not touch.
+
+HIDDEN_COLS <- c("skill", "style", "frailty", "fi_noise", "smoking_true", "z0", "rate", "venous", "z")
+
+#' Turn a population and its simulated course into claims-like records.
+#'
+#' @return list(patients = one row per patient: measured baseline and coded event times;
+#'   pm = person-months (id, t, sev_obs) from month 0 to death or month 59;
+#'   hidden = latent traits, for oracle analyses and checks only).
+make_records <- function(pop, course, par) {
+  pat <- pop$pat; ev <- course$events; ms <- par$measurement
+  M <- par$size$months
+  smoking <- if (isTRUE(ms$smoking) || identical(ms$smoking, "true")) pat$smoking_true else pat$smoking_coded
+  rel <- ms$frailty_reliability
+  fi <- sqrt(rel) * pat$frailty + sqrt(1 - rel) * pat$fi_noise
+
+  patients <- data.table(id = pat$id, phys_id = pat$phys_id, obl = pat$obl, age = pat$age,
+                         female = pat$female, black = pat$black, year = pat$year,
+                         diabetes = pat$diabetes, ckd = pat$ckd, chf = pat$chf, copd = pat$copd,
+                         prior_cvd = pat$prior_cvd, smoking = smoking, frailty_index = fi)
+  patients <- cbind(patients, ev[, !"id"])
+  # t_clti_true is latent: only the coded date is in claims
+  patients[, t_clti_true := NULL]
+
+  # observed severity: exact Z_t (clean worlds) or a claims proxy (world W5, WR)
+  last <- pmin(ifelse(is.na(ev$t_death), M - 1L, ev$t_death), M - 1L)
+  nrep <- last + 1L
+  id <- rep(pat$id, nrep)
+  t <- sequence(nrep) - 1L
+  idx <- cbind(id, t + 1L)
+  if (identical(ms$severity, "exact")) {
+    sev <- course$z[idx]
+  } else {
+    w <- ms$proxy_window
+    cs <- t(apply(course$counts, 1, function(x) { x[is.na(x)] <- 0L; cumsum(x) }))
+    lagged <- cbind(matrix(0, nrow(cs), w), cs)[, seq_len(M)]
+    sev <- log1p(cs[idx] - lagged[idx])
+  }
+  pm <- data.table(id = id, t = t, sev_obs = sev)
+
+  hidden <- data.table(id = pat$id, skill = pat$skill, style = pat$style, frailty = pat$frailty,
+                       smoking_true = pat$smoking_true, rate = pat$rate,
+                       t_clti_true = ev$t_clti_true)
+  list(patients = patients, pm = pm, hidden = hidden, par = par)
+}
+
+#' Stop if an analysis table carries a latent column.
+assert_no_hidden <- function(dt) {
+  bad <- intersect(names(dt), HIDDEN_COLS)
+  if (length(bad)) stop("latent columns in analysis data: ", paste(bad, collapse = ", "))
+  invisible(TRUE)
+}
+
+#' One replicate's records: population, natural course, claims view.
+simulate_records <- function(par, rep = 1L, n = par$size$n_patients) {
+  seed <- replicate_seed(par$world, par$scenario, rep)
+  pop <- make_population(par, n = n, seed = seed)
+  course <- simulate_course(pop, par, policy("natural"), seed = seed + 104729L, record = TRUE)
+  make_records(pop, course, par)
+}
+
+#' Reproducible seed for a world x scenario x replicate.
+replicate_seed <- function(world, scenario, rep) {
+  key <- paste(world, scenario, rep, sep = "|")
+  h <- sum(utf8ToInt(key) * seq_along(utf8ToInt(key))^2)
+  as.integer((h * 2654435761) %% 2147483647)
+}
diff --git a/pvisim/R/simulate.R b/pvisim/R/simulate.R
new file mode 100644
index 0000000..c2d6f62
--- /dev/null
+++ b/pvisim/R/simulate.R
@@ -0,0 +1,217 @@
+# The simulated world: physicians, patients, and the monthly course of disease and care.
+# Analysis plan v1.0, section 5.1.
+#
+# Common random numbers. simulate_course() draws the same uniforms for every patient
+# in every month, in a fixed order, whatever the patient's state or the policy. So two
+# runs on the same population with the same seed differ only through the policy. This
+# keeps truths under different strategies paired, and makes the sharp null exact at the
+# patient level (tested).
+
+STREAMS <- c("eps", "proxy", "clti", "false_code", "amp", "bypass", "rr", "clti_revasc",
+             "nco1", "nco2", "nce", "pvi", "death")
+
+#' Physicians: latent skill k and style s (standard normal, correlated), setting, panel weight.
+make_physicians <- function(par, seed) {
+  set.seed(seed)
+  J <- par$size$n_physicians
+  rho <- par$physician$rho_skill_style
+  z1 <- rnorm(J); z2 <- rnorm(J)
+  s <- z1
+  k <- rho * z1 + sqrt(1 - rho^2) * z2
+  obl <- rbinom(J, 1, plogis(par$physician$obl_intercept + par$physician$obl_style * s))
+  w <- rlnorm(J, 0, par$size$panel_sdlog)
+  data.table(phys_id = seq_len(J), skill = k, style = s, obl = obl, weight = w / sum(w))
+}
+
+#' Patients at diagnosis: measured baseline, latent traits, and assigned physician.
+make_patients <- function(par, n, phys, seed) {
+  set.seed(seed)
+  pp <- par$patient
+  phys_id <- sample(phys$phys_id, n, replace = TRUE, prob = phys$weight)
+  age <- pp$age_min + rgamma(n, shape = pp$age_shape, scale = pp$age_scale)
+  frailty <- rnorm(n)
+  a10 <- (age - 76) / 10
+  comorb <- lapply(pp$comorbidity, function(cm)
+    rbinom(n, 1, plogis(qlogis(cm$prev) + cm$age_coef * a10 + cm$frailty_coef * frailty)))
+  smoking <- rbinom(n, 1, plogis(qlogis(pp$smoking$prev) + pp$smoking$age_coef * a10))
+  sv <- par$severity
+  z0 <- sv$z0_diabetes * comorb$diabetes + sv$z0_smoking * smoking + sv$z0_ckd * comorb$ckd +
+    rnorm(n, 0, sv$z0_sd)
+  rate <- sv$r0 + sv$r_diabetes * comorb$diabetes + sv$r_smoking * smoking + rnorm(n, 0, sv$r_sd)
+  venous <- rnorm(n)
+  year <- sample(pp$year_min:pp$year_max, n, replace = TRUE)
+  # what claims show for smoking and frailty (the analyst's versions; plan D17)
+  ms <- par$measurement
+  smoking_coded <- ifelse(smoking == 1, rbinom(n, 1, ms$smoking_sens), rbinom(n, 1, 1 - ms$smoking_spec))
+  fi_noise <- rnorm(n)
+  pat <- data.table(id = seq_len(n), phys_id = phys_id, age = age,
+                    female = rbinom(n, 1, pp$p_female), black = rbinom(n, 1, pp$p_black),
+                    year = year)
+  for (nm in names(comorb)) set(pat, j = nm, value = comorb[[nm]])
+  pat[, `:=`(smoking_true = smoking, smoking_coded = smoking_coded, frailty = frailty,
+             fi_noise = fi_noise, z0 = z0, rate = rate, venous = venous)]
+  pat <- merge(pat, phys[, .(phys_id, skill, style, obl)], by = "phys_id", sort = FALSE)
+  setorder(pat, id)
+  pat[]
+}
+
+#' A full population for one replicate (or for the truths).
+make_population <- function(par, n = par$size$n_patients, seed = 1L, phys = NULL) {
+  if (is.null(phys)) phys <- make_physicians(par, seed)
+  pat <- make_patients(par, n, phys, seed + 7919L)
+  list(phys = phys, pat = pat)
+}
+
+# ---- policies -----------------------------------------------------------------------
+# A policy overrides the natural elective-PVI decision among patients eligible in month t
+# (alive, no PVI yet, no MALE yet). Non-elective revascularization after CLTI is never
+# overridden: it happens after a MALE event, which every strategy allows.
+#   natural            observed practice
+#   A                  natural in months 0..(w-2); everyone still eligible gets PVI at w-1
+#   B1                 no PVI in months 0..(w-1); natural afterwards
+#   B2                 no elective PVI ever
+#   seq_force / seq_forbid (m)  natural except PVI forced / forbidden at month m
+
+policy <- function(type = "natural", m = NA_integer_) {
+  stopifnot(type %in% c("natural", "A", "B1", "B2", "seq_force", "seq_forbid"))
+  list(type = type, m = m)
+}
+
+apply_policy <- function(pol, t, natural, window) {
+  switch(pol$type,
+    natural = natural,
+    A = if (t < window - 1) natural else if (t == window - 1) rep(TRUE, length(natural)) else natural,
+    B1 = if (t < window) rep(FALSE, length(natural)) else natural,
+    B2 = rep(FALSE, length(natural)),
+    seq_force = if (t == pol$m) rep(TRUE, length(natural)) else natural,
+    seq_forbid = if (t == pol$m) rep(FALSE, length(natural)) else natural)
+}
+
+# ---- the monthly course -------------------------------------------------------------
+
+#' Simulate months 0..(months-1) for every patient under a policy.
+#'
+#' Within a month the order is: severity update; outcome events (CLTI, coding, amputation,
+#' bypass, repeat revascularization, revascularization after CLTI, negative controls);
+#' death; then the elective PVI decision among survivors still eligible. A MALE event
+#' and death in the same month count the MALE event first.
+#'
+#' @return list(pat = patient table with event times, z = N x months matrix of Z_t
+#'   (record = TRUE only), counts = N x months matrix of proxy counts, pvi_month = vector).
+simulate_course <- function(pop, par, pol = policy("natural"), seed = 1L, record = TRUE) {
+  pat <- pop$pat
+  N <- nrow(pat)
+  M <- par$size$months
+  W <- par$size$window
+  tr <- par$treatment; oc <- par$outcomes; nc <- par$negative_controls; ef <- par$effect
+  a10 <- (pat$age - 76) / 10
+
+  # time-fixed parts of each linear predictor
+  lp_win <- tr$window$intercept + tr$window$frailty * pat$frailty + tr$window$age10 * a10 +
+    tr$window$diabetes * pat$diabetes + tr$window$ckd * pat$ckd + tr$window$smoking * pat$smoking_true +
+    tr$window$style * pat$style + tr$window$obl * pat$obl + tr$window$year * (pat$year - par$patient$year_min)
+  lp_def <- tr$deferred$intercept + tr$deferred$style * pat$style + tr$deferred$obl * pat$obl
+  d <- oc$death
+  lp_death0 <- d$intercept + d$chf * pat$chf + d$ckd * pat$ckd + d$copd * pat$copd + d$frailty * pat$frailty
+  cl <- oc$clti
+  lp_clti0 <- cl$intercept + cl$diabetes * pat$diabetes + cl$ckd * pat$ckd + cl$smoking * pat$smoking_true +
+    cl$frailty * pat$frailty
+  lp_false <- oc$clti_false_code$intercept + oc$clti_false_code$style * pat$style
+  am <- oc$amputation
+  lp_amp0 <- am$intercept + am$diabetes * pat$diabetes + am$ckd * pat$ckd
+  rr <- oc$repeat_revasc
+  lp_rr0 <- rr$intercept + rr$skill * pat$skill + rr$style * pat$style
+  lp_nco1 <- nc$nco1$intercept + nc$nco1$frailty * pat$frailty + nc$nco1$smoking * pat$smoking_true +
+    nc$nco1$age10 * a10
+  lp_nco2 <- nc$nco2$intercept + nc$nco2$style * pat$style + nc$nco2$venous * pat$venous
+  lp_nce <- nc$nce$intercept + nc$nce$style * pat$style + nc$nce$frailty * pat$frailty
+  kappa_post <- ef$kappa * exp(-ef$skill_kappa * pat$skill)
+
+  # state
+  Z <- pat$z0
+  alive <- rep(TRUE, N); treated <- rep(FALSE, N); clti_true <- rep(FALSE, N)
+  clti_code <- rep(FALSE, N); male <- rep(FALSE, N)
+  t_pvi <- t_clti_true <- t_clti_code <- t_amp <- t_bypass <- t_rr <- t_male <- t_death <-
+    t_nco1 <- t_nco2 <- t_nce <- rep(NA_integer_, N)
+  pvi_elective <- rep(NA, N)
+  n_proc <- integer(N)
+  zmat <- if (record) matrix(NA_real_, N, M) else NULL
+  cmat <- if (record) matrix(NA_integer_, N, M) else NULL
+  ns <- length(STREAMS)
+
+  set.seed(seed)
+  for (t in 0:(M - 1)) {
+    U <- matrix(runif(N * ns), N, ns, dimnames = list(NULL, STREAMS))
+    a <- alive
+    if (t > 0) Z[a] <- Z[a] + pat$rate[a] + kappa_post[a] * treated[a] +
+      par$severity$eps_sd * qnorm(U[a, "eps"])
+    if (record) {
+      zmat[a, t + 1] <- Z[a]
+      cmat[a, t + 1] <- qpois(U[a, "proxy"], exp(par$measurement$proxy_intercept + par$measurement$proxy_z * Z[a]))
+    }
+
+    # true CLTI, and its coding (true onset always coded; false codes by physician style)
+    new_clti <- a & !clti_true & U[, "clti"] < plogis(lp_clti0 + cl$z * Z)
+    clti_true[new_clti] <- TRUE; t_clti_true[new_clti] <- t
+    new_code <- a & !clti_code & (new_clti | U[, "false_code"] < plogis(lp_false))
+    clti_code[new_code] <- TRUE; t_clti_code[new_code] <- t
+
+    new_amp <- a & is.na(t_amp) & U[, "amp"] < plogis(lp_amp0 + am$clti * clti_true + am$z * Z)
+    t_amp[new_amp] <- t
+    new_byp <- a & is.na(t_bypass) & U[, "bypass"] < plogis(oc$bypass$intercept + oc$bypass$z * Z)
+    t_bypass[new_byp] <- t
+    n_proc[new_byp] <- n_proc[new_byp] + 1L
+
+    # repeat revascularization: only after the index PVI (from the month after it)
+    new_rr <- a & treated & !is.na(t_pvi) & t > t_pvi & U[, "rr"] < plogis(lp_rr0 + rr$z * Z)
+    t_rr[new_rr & is.na(t_rr)] <- t
+    n_proc[new_rr] <- n_proc[new_rr] + 1L
+
+    # revascularization after a CLTI diagnosis (not elective; allowed by every strategy)
+    cr <- a & clti_code & U[, "clti_revasc"] < plogis(oc$clti_revasc$intercept)
+    first <- cr & !treated
+    rep_ <- cr & treated & !new_rr
+    t_rr[rep_ & is.na(t_rr)] <- t
+    n_proc[cr & !new_rr] <- n_proc[cr & !new_rr] + 1L
+    if (any(first)) {
+      treated[first] <- TRUE; t_pvi[first] <- t; pvi_elective[first] <- FALSE
+      Z[first] <- Z[first] - ef$relief
+    }
+
+    # negative controls (never depend on PVI)
+    x <- a & is.na(t_nco1) & U[, "nco1"] < plogis(lp_nco1); t_nco1[x] <- t
+    x <- a & is.na(t_nco2) & U[, "nco2"] < plogis(lp_nco2); t_nco2[x] <- t
+    if (t < W) { x <- a & is.na(t_nce) & U[, "nce"] < plogis(lp_nce); t_nce[x] <- t }
+
+    # MALE: coded CLTI, amputation, bypass, or repeat revascularization
+    new_male <- a & !male & (clti_code | !is.na(t_amp) | !is.na(t_bypass) | !is.na(t_rr))
+    male[new_male] <- TRUE; t_male[new_male] <- t
+
+    # death
+    dies <- a & U[, "death"] < plogis(lp_death0 + d$age10 * a10 + d$age10 * t / 120 + d$z * Z + d$clti * clti_true)
+    alive[dies] <- FALSE; t_death[dies] <- t
+
+    # elective PVI among survivors still eligible
+    elig <- alive & !treated & !male
+    if (any(elig)) {
+      lp <- if (t < W) lp_win[elig] + tr$window$z * Z[elig] else lp_def[elig] + tr$deferred$z * Z[elig]
+      natural <- U[elig, "pvi"] < plogis(lp)
+      go <- apply_policy(pol, t, natural, W)
+      idx <- which(elig)[go]
+      if (length(idx)) {
+        treated[idx] <- TRUE; t_pvi[idx] <- t; pvi_elective[idx] <- TRUE
+        n_proc[idx] <- n_proc[idx] + 1L
+        Z[idx] <- Z[idx] - ef$relief
+      }
+    }
+  }
+
+  # disease composite (D34): coded CLTI or major amputation
+  t_disease <- pmin(t_clti_code, t_amp, na.rm = TRUE)
+  out <- data.table(id = pat$id, t_pvi, pvi_elective, t_clti_true, t_clti_code, t_amp, t_bypass,
+                    t_rr, t_male, t_disease, t_death, t_nco1, t_nco2, t_nce, n_proc)
+  list(events = out, z = zmat, counts = cmat, policy = pol)
+}
+
+#' Cumulative incidence by month `h` (event in months 0..h-1, before death or in the same month).
+cum_inc <- function(t_event, h = 60L) mean(!is.na(t_event) & t_event < h)
diff --git a/pvisim/inst/params/dgp_default.yml b/pvisim/inst/params/dgp_default.yml
new file mode 100644
index 0000000..b6a74e3
--- /dev/null
+++ b/pvisim/inst/params/dgp_default.yml
@@ -0,0 +1,125 @@
+# Every number that defines the simulated world. Analysis plan v1.0, section 5.
+# Time unit: month. Hazards are monthly, on the logit scale.
+# Starting values are placeholders, to be tuned by calibration
+# (simulations/01_calibrate.R, targets in plan section 5.4).
+# Physician traits (skill k, style s) and patient frailty F are standard normal;
+# their influence is set by the coefficients that multiply them.
+
+size:
+  n_patients: 50000
+  n_physicians: 1500
+  panel_sdlog: 1.0          # lognormal panel weights; median panel about 20 at n = 50,000
+  months: 60                # follow-up months 0-59
+  window: 6                 # strategy A: PVI in months 0-5
+
+physician:
+  rho_skill_style: -0.3     # corr(k, s). Has no effect unless skill matters (world W1)
+  obl_intercept: 0.0        # logit P(office-based lab / ASC); about 50% of physicians
+  obl_style: 0.8            # high-style physicians more often work in OBL/ASC
+
+patient:
+  age_min: 66
+  age_shape: 2.5            # age = 66 + Gamma(shape, scale); mean about 76
+  age_scale: 4.0
+  p_female: 0.50
+  p_black: 0.10
+  # comorbidities: logit P = qlogis(prev) + age_coef * (age - 76) / 10 + frailty_coef * F
+  comorbidity:
+    diabetes:  {prev: 0.40, age_coef: -0.10, frailty_coef: 0.4}
+    ckd:       {prev: 0.20, age_coef:  0.30, frailty_coef: 0.5}
+    chf:       {prev: 0.20, age_coef:  0.40, frailty_coef: 0.6}
+    copd:      {prev: 0.25, age_coef:  0.10, frailty_coef: 0.4}
+    prior_cvd: {prev: 0.45, age_coef:  0.20, frailty_coef: 0.3}
+  smoking: {prev: 0.25, age_coef: -0.40}
+  year_min: 2015
+  year_max: 2019
+
+severity:                    # latent Z_t: walking limitation / ABI on a standard scale
+  z0_diabetes: 0.3
+  z0_smoking: 0.3
+  z0_ckd: 0.2
+  z0_sd: 1.0
+  r0: 0.02                   # mean drift per month
+  r_diabetes: 0.01
+  r_smoking: 0.01
+  r_sd: 0.01                 # patient-specific progression rate
+  eps_sd: 0.05               # month-to-month noise
+
+treatment:
+  window:                    # elective PVI hazard, months 0-5, among the eligible
+    intercept: -5.4
+    z: 0.9                   # confounding by indication
+    frailty: -0.3            # frail patients are less often offered PVI
+    age10: -0.3
+    diabetes: 0.1
+    ckd: -0.1
+    smoking: 0.1
+    style: 0.5               # physician style; world W2 sets it to the strong level
+    obl: 0.3
+    year: 0.05               # per calendar year after year_min
+  deferred:                  # elective PVI hazard from month 6 on (deferral / crossover)
+    intercept: -7.2
+    z: 0.0                   # world W4 makes deferral severity-driven
+    style: 0.5
+    obl: 0.3
+
+effect:                      # what a PVI does. Set by the truth scenario
+  relief: 0.0                # immediate drop in Z at PVI
+  kappa: 0.0                 # extra Z drift per month after PVI
+  skill_kappa: 0.0           # kappa * exp(-skill_kappa * k): skill slows post-PVI progression
+
+outcomes:
+  death:  {intercept: -5.0, age10: 0.7, chf: 0.4, ckd: 0.4, copd: 0.3, frailty: 0.4, z: 0.1, clti: 0.8}
+  clti:   {intercept: -6.6, z: 0.6, diabetes: 0.4, ckd: 0.3, smoking: 0.3, frailty: 0.2}
+  clti_false_code: {intercept: -9.0, style: 0.0}     # factor 3 sets style
+  amputation: {intercept: -8.5, clti: 3.0, diabetes: 0.4, ckd: 0.4, z: 0.2}
+  bypass: {intercept: -7.8, z: 0.5}
+  repeat_revasc: {intercept: -5.3, skill: 0.0, style: 0.0, z: 0.2}   # only after the index PVI
+  clti_revasc: {intercept: -4.5}                      # PVI after CLTI (not elective); most CLTI patients get bypass, tibial work, amputation or nothing
+
+negative_controls:           # never depend on PVI
+  nco1: {intercept: -6.4, frailty: 0.5, smoking: 0.4, age10: 0.3}   # hip fracture / pneumonia
+  nco2: {intercept: -7.1, style: 0.8, venous: 0.5}                  # varicose-vein ablation
+  nce:  {intercept: -5.7, style: 0.8, frailty: -0.2}                # low-value procedure, months 0-5
+
+measurement:                 # what the analyst sees. World W5 switches to claims proxies
+  severity: exact            # exact: analyst sees Z_t. proxy: only visit/testing counts
+  proxy_intercept: -1.0      # monthly count ~ Poisson(exp(proxy_intercept + proxy_z * Z_t))
+  proxy_z: 0.5
+  proxy_window: 3            # observed severity = log1p(sum of counts over the last 3 months)
+  frailty_reliability: 1.0   # frailty index = sqrt(rel) * F + sqrt(1 - rel) * noise
+  smoking: true              # true: analyst sees true smoking. coded: sensitivity / specificity below
+  smoking_sens: 0.6
+  smoking_spec: 0.98
+
+worlds:                      # one factor at a time (plan section 5.2, D24). Overrides on the default.
+                             # Intercepts are then re-solved per world by calibration, so every world
+                             # matches the same Medicare targets unless `calibration_targets` says otherwise.
+  W0: {}
+  W1:
+    effect: {skill_kappa: 0.5}
+    outcomes: {repeat_revasc: {skill: -0.6}}
+  W1_rho0:
+    physician: {rho_skill_style: 0.0}
+    effect: {skill_kappa: 0.5}
+    outcomes: {repeat_revasc: {skill: -0.6}}
+  W2:
+    treatment: {window: {style: 1.6}, deferred: {style: 1.6}}
+  W3:
+    treatment: {window: {style: 1.6}, deferred: {style: 1.6}}
+    outcomes: {clti_false_code: {style: 0.8}, repeat_revasc: {style: 0.6}}
+  W4:                        # deferral common (as in trial conservative arms) and severity-driven
+    treatment: {deferred: {z: 0.8}}
+    calibration_targets: {late_not: 0.15}
+  W5:
+    measurement: {severity: proxy, frailty_reliability: 0.4, smoking: coded}
+  WR:
+    effect: {skill_kappa: 0.5}
+    treatment: {window: {style: 1.6}, deferred: {style: 1.6, z: 0.8}}   # deferral at the Sorber rate, but severity-driven
+    outcomes: {repeat_revasc: {skill: -0.6, style: 0.6}, clti_false_code: {style: 0.8}}
+    measurement: {severity: proxy, frailty_reliability: 0.4, smoking: coded}
+
+scenarios:                   # truth scenarios (plan section 5.3, D34)
+  "null": {effect: {relief: 0.0, kappa: 0.0}}           # sharp disease null
+  harm: {effect: {relief: 0.3, kappa: 0.03}}          # tuned to RD about +5 pp (calibration)
+  benefit_harm: {effect: {relief: 1.2, kappa: 0.04}}  # curves cross; tuned in calibration
diff --git a/pvisim/tests/testthat/test-simulate.R b/pvisim/tests/testthat/test-simulate.R
new file mode 100644
index 0000000..2732ca6
--- /dev/null
+++ b/pvisim/tests/testthat/test-simulate.R
@@ -0,0 +1,55 @@
+test_that("every world and scenario loads", {
+  raw <- yaml::read_yaml(default_params_path())
+  for (w in names(raw$worlds)) for (s in names(raw$scenarios)) {
+    par <- load_params(w, s)
+    expect_equal(par$world, w)
+    expect_true(is.numeric(par$effect$relief))
+  }
+  expect_equal(load_params("W2")$treatment$window$style, 1.6)
+  expect_equal(load_params("W2")$treatment$window$z, 0.9)   # untouched by the override
+})
+
+small_world <- function(world = "WR", scenario = "null", n = 4000) {
+  par <- load_params(world, scenario, overrides = list(size = list(n_physicians = 150)))
+  list(par = par, pop = make_population(par, n = n, seed = 11L))
+}
+
+test_that("under the sharp null, disease, death and negative controls are identical under every policy", {
+  sw <- small_world("WR", "null")
+  runs <- lapply(list(policy("natural"), policy("A"), policy("B1"), policy("B2"),
+                      policy("seq_force", 3L)),
+                 function(p) simulate_course(sw$pop, sw$par, p, seed = 5L, record = FALSE)$events)
+  for (r in runs[-1]) {
+    for (col in c("t_disease", "t_clti_code", "t_amp", "t_death", "t_nco1", "t_nco2", "t_nce")) {
+      expect_identical(r[[col]], runs[[1]][[col]], info = col)
+    }
+  }
+  # MALE is not null: repeat revascularization needs an index PVI (D34)
+  expect_gt(cum_inc(runs[[2]]$t_male), cum_inc(runs[[3]]$t_male))
+})
+
+test_that("policies are followed", {
+  sw <- small_world("W0", "harm")
+  W <- sw$par$size$window
+  a <- simulate_course(sw$pop, sw$par, policy("A"), seed = 5L, record = FALSE)$events
+  # under A, everyone who survives month W-1 without a MALE event has had a PVI by then
+  # (death and MALE in a month come before that month's PVI decision)
+  free <- (is.na(a$t_death) | a$t_death > W - 1) & (is.na(a$t_male) | a$t_male > W - 1 | (!is.na(a$t_pvi) & a$t_pvi < a$t_male))
+  expect_true(all(!is.na(a$t_pvi[free]) & a$t_pvi[free] <= W - 1))
+  b1 <- simulate_course(sw$pop, sw$par, policy("B1"), seed = 5L, record = FALSE)$events
+  expect_false(any(b1$pvi_elective %in% TRUE & b1$t_pvi < W))
+  b2 <- simulate_course(sw$pop, sw$par, policy("B2"), seed = 5L, record = FALSE)$events
+  expect_false(any(b2$pvi_elective %in% TRUE))
+})
+
+test_that("records carry no latent columns and person-months end at death", {
+  par <- load_params("W5", "harm", overrides = list(size = list(n_physicians = 100)))
+  rec <- simulate_records(par, rep = 1L, n = 2000)
+  expect_true(assert_no_hidden(rec$patients))
+  expect_true(assert_no_hidden(rec$pm))
+  last <- rec$pm[, .(tmax = max(t)), by = id]
+  dead <- rec$patients[!is.na(t_death), .(id, t_death)]
+  chk <- merge(last, dead, by = "id")
+  expect_true(all(chk$tmax == chk$t_death))
+  expect_true(all(is.finite(rec$pm$sev_obs)))
+})
diff --git a/simulations/01_calibrate.R b/simulations/01_calibrate.R
new file mode 100644
index 0000000..61e29ad
--- /dev/null
+++ b/simulations/01_calibrate.R
@@ -0,0 +1,80 @@
+# Calibrate each world's natural course to the Medicare targets (plan section 5.4).
+#
+# For every world, under the null scenario, solve one intercept per target, cycling
+# until all targets hold. Every world gets the same targets, so worlds differ only in
+# the threat they switch on, not in base rates.
+#
+# Usage: Rscript simulations/01_calibrate.R [n] [world ...]
+# Writes simulations/config/calibrated.yml (merged over existing entries).
+
+pkgload::load_all("pvisim", quiet = TRUE)
+suppressPackageStartupMessages(library(data.table))
+
+args <- commandArgs(TRUE)
+n <- if (length(args) && grepl("^[0-9e.]+$", args[1])) as.numeric(args[1]) else 1e5
+worlds <- setdiff(args, as.character(args[1]))
+if (!length(worlds)) worlds <- names(yaml::read_yaml(default_params_path())$worlds)
+H <- 53L   # Sorber 2023 median follow-up, about 4.4 years
+
+# Natural-course statistics that the targets are defined on.
+stats_of <- function(par) {
+  pop <- make_population(par, n = n, seed = 20261027L)
+  ev <- simulate_course(pop, par, policy("natural"), seed = 20261028L, record = FALSE)$events
+  early <- ev$pvi_elective %in% TRUE & ev$t_pvi < par$size$window
+  late_early <- !is.na(ev$t_rr) & ev$t_rr < H
+  late_not <- !is.na(ev$t_pvi) & ev$t_pvi >= par$size$window & ev$t_pvi < H
+  list(early = mean(early),
+       late_not = mean(late_not[!early]),
+       late_early = mean(late_early[early]),
+       clti_not = mean((!is.na(ev$t_clti_code) & ev$t_clti_code < H)[!early]),
+       death = cum_inc(ev$t_death))
+}
+
+# target name -> parameter path and value (plan 5.4)
+targets <- list(
+  early      = list(path = c("treatment", "window", "intercept"), value = 0.032),
+  late_not   = list(path = c("treatment", "deferred", "intercept"), value = 0.036),
+  late_early = list(path = c("outcomes", "repeat_revasc", "intercept"), value = 0.225),
+  clti_not   = list(path = c("outcomes", "clti", "intercept"), value = 0.078),
+  death      = list(path = c("outcomes", "death", "intercept"), value = 0.35)
+)
+
+set_path <- function(lst, path, value) {
+  if (length(path) == 1) { lst[[path]] <- value; return(lst) }
+  lst[[path[1]]] <- set_path(if (is.null(lst[[path[1]]])) list() else lst[[path[1]]], path[-1], value)
+  lst
+}
+get_path <- function(lst, path) Reduce(function(x, p) x[[p]], path, lst)
+
+# Each target depends mainly on its own intercept, and for rare monthly events
+# logit(rate) moves about one-for-one with the intercept. So update all intercepts
+# together: b <- b - damp * (logit(stat) - logit(target)), until every target is within
+# `tol` on the logit scale. Common random numbers make each evaluation deterministic.
+calibrate_world <- function(world, max_iter = 25, tol = 0.02, damp = 0.9) {
+  base <- load_params(world, "null", calibrated = FALSE)
+  for (tg in names(base$calibration_targets)) targets[[tg]]$value <- base$calibration_targets[[tg]]
+  ov <- list()
+  for (tg in names(targets)) ov <- set_path(ov, targets[[tg]]$path, get_path(base, targets[[tg]]$path))
+  for (it in seq_len(max_iter)) {
+    s <- stats_of(load_params(world, "null", overrides = ov, calibrated = FALSE))
+    gap <- sapply(names(targets), function(tg) qlogis(s[[tg]]) - qlogis(targets[[tg]]$value))
+    cat(sprintf("%s iter %2d: %s\n", world, it,
+                paste(sprintf("%s %.3f", names(s), unlist(s)), collapse = "  ")))
+    if (all(abs(gap) < tol)) break
+    for (tg in names(targets)) {
+      b <- get_path(ov, targets[[tg]]$path)
+      ov <- set_path(ov, targets[[tg]]$path, round(b - damp * gap[[tg]], 4))
+    }
+  }
+  if (any(abs(gap) >= tol)) warning(world, ": not converged; largest logit gap ", round(max(abs(gap)), 3))
+  ov
+}
+
+out_path <- "simulations/config/calibrated.yml"
+cal <- if (file.exists(out_path)) yaml::read_yaml(out_path) else list()
+for (w in worlds) {
+  t0 <- Sys.time()
+  cal[[w]] <- calibrate_world(w)
+  cat(sprintf("%s done in %.1f min\n", w, as.numeric(Sys.time() - t0, units = "mins")))
+  yaml::write_yaml(cal, out_path)
+}
diff --git a/simulations/config/calibrated.yml b/simulations/config/calibrated.yml
new file mode 100644
index 0000000..17bad7e
--- /dev/null
+++ b/simulations/config/calibrated.yml
@@ -0,0 +1,104 @@
+W0:
+  treatment:
+    window:
+      intercept: -6.3552
+    deferred:
+      intercept: -7.7172
+  outcomes:
+    repeat_revasc:
+      intercept: -5.6375
+    clti:
+      intercept: -7.4903
+    death:
+      intercept: -5.6504
+W4:
+  treatment:
+    window:
+      intercept: -6.3535
+    deferred:
+      intercept: -6.7744
+  outcomes:
+    repeat_revasc:
+      intercept: -5.6429
+    clti:
+      intercept: -7.4903
+    death:
+      intercept: -5.6542
+W1:
+  treatment:
+    window:
+      intercept: -6.3552
+    deferred:
+      intercept: -7.7172
+  outcomes:
+    repeat_revasc:
+      intercept: -5.7951
+    clti:
+      intercept: -7.4903
+    death:
+      intercept: -5.6504
+W1_rho0:
+  treatment:
+    window:
+      intercept: -6.3552
+    deferred:
+      intercept: -7.7172
+  outcomes:
+    repeat_revasc:
+      intercept: -5.7187
+    clti:
+      intercept: -7.4903
+    death:
+      intercept: -5.6504
+W2:
+  treatment:
+    window:
+      intercept: -7.3097
+    deferred:
+      intercept: -8.5615
+  outcomes:
+    repeat_revasc:
+      intercept: -5.5855
+    clti:
+      intercept: -7.4913
+    death:
+      intercept: -5.6505
+W3:
+  treatment:
+    window:
+      intercept: -7.3134
+    deferred:
+      intercept: -8.5848
+  outcomes:
+    repeat_revasc:
+      intercept: -6.4032
+    clti:
+      intercept: -7.5232
+    death:
+      intercept: -5.649
+W5:
+  treatment:
+    window:
+      intercept: -6.3552
+    deferred:
+      intercept: -7.7172
+  outcomes:
+    repeat_revasc:
+      intercept: -5.6375
+    clti:
+      intercept: -7.4903
+    death:
+      intercept: -5.6504
+WR:
+  treatment:
+    window:
+      intercept: -7.3131
+    deferred:
+      intercept: -9.598
+  outcomes:
+    repeat_revasc:
+      intercept: -6.8121
+    clti:
+      intercept: -7.5233
+    death:
+      intercept: -5.6489
diff --git a/simulations/explore_targets.R b/simulations/explore_targets.R
new file mode 100644
index 0000000..e35ffd1
--- /dev/null
+++ b/simulations/explore_targets.R
@@ -0,0 +1,26 @@
+# Quick look: the natural course of each world against the calibration targets (plan 5.4).
+pkgload::load_all("pvisim", quiet = TRUE)
+library(data.table)
+args <- commandArgs(TRUE); worlds <- if (length(args)) args else c("W0", "WR")
+for (w in worlds) {
+  par <- load_params(w, "harm")
+  t0 <- Sys.time()
+  pop <- make_population(par, seed = 1L)
+  ev <- simulate_course(pop, par, policy("natural"), seed = 2L, record = FALSE)$events
+  ev <- merge(ev, pop$pat[, .(id, phys_id)], by = "id")
+  early <- !is.na(ev$t_pvi) & ev$t_pvi < 6 & ev$pvi_elective %in% TRUE
+  h <- 53  # Sorber median follow-up about 4.4 years
+  late <- !is.na(ev$t_pvi) & ev$t_pvi >= 6 | !is.na(ev$t_rr) & ev$t_rr >= 6
+  late <- late & pmin(ev$t_pvi * NA^(ev$t_pvi < 6), ifelse(!is.na(ev$t_rr) & ev$t_rr >= 6, ev$t_rr, NA), na.rm = TRUE) < h
+  clti <- !is.na(ev$t_clti_code) & ev$t_clti_code < h
+  phys <- ev[, .(n = .N, rate = mean(!is.na(t_pvi) & t_pvi < 6 & pvi_elective %in% TRUE)), by = phys_id]
+  sdr <- sd(phys$rate); out <- mean(phys$rate >= mean(phys$rate) + 2 * sdr)
+  cat(sprintf("\n%s (%.0fs)\n  early PVI %.1f%% (target 3.2)\n  physician rate: median %.1f%%, max %.0f%%, >=2SD outliers %.1f%% (0, 58, 5.6)\n",
+              w, as.numeric(Sys.time() - t0, units = "secs"), 100 * mean(early),
+              100 * median(phys$rate), 100 * max(phys$rate), 100 * out))
+  cat(sprintf("  late PVI by 4.4y: early %.1f%% vs not %.1f%% (22.5 vs 3.6)\n", 100 * mean(late[early]), 100 * mean(late[!early])))
+  cat(sprintf("  coded CLTI by 4.4y: early %.1f%% vs not %.1f%% (16.4 vs 7.8)\n", 100 * mean(clti[early]), 100 * mean(clti[!early])))
+  cat(sprintf("  death by 5y %.1f%% (about 35)   MALE by 5y %.1f%%   amputation %.2f%%   disease composite %.1f%%\n",
+              100 * cum_inc(ev$t_death), 100 * cum_inc(ev$t_male), 100 * cum_inc(ev$t_amp), 100 * cum_inc(ev$t_disease)))
+  cat(sprintf("  NCO1 %.1f%%  NCO2 %.1f%%  NCE %.1f%%\n", 100 * cum_inc(ev$t_nco1), 100 * cum_inc(ev$t_nco2), 100 * cum_inc(ev$t_nce)))
+}