""" bluebelle_case_sim.py ===================================================================== Reproducible illustrative computations for the manuscript "Value of Information with Imprecise Probabilities for Decision Modeling: Research Prioritization When Evidence Does Not Identify a Single Precise Probability Measure" (BMC Medical Research Methodology) This single script regenerates EVERY numeric value reported in the applied wound-dressing example, together with Monte Carlo uncertainty for the EVPI / EVPPI / EVSI endpoints, and writes the result tables used by the manuscript. Outputs (written to the working directory): bluebelle_main_results.csv per-measure EVPI / EVPPI / ENB / pbest bluebelle_main_evsi.csv per-measure EVSI (mean over replicates) bluebelle_voi_uncertainty.csv reference + envelope endpoints with MC SE bluebelle_validation.csv validation / simulation-study checks Reproducibility --------------- All randomness is driven by NumPy Generator objects with fixed seeds, so the script is deterministic: rerunning it on the same NumPy version reproduces the CSV files byte-for-byte. The EVSI estimator is a common-random-number (CRN) nested Monte Carlo estimator on the normal-normal log-odds-ratio update; see compute_evsi() for details. Run: python3 bluebelle_case_sim.py Tested with numpy>=1.24, pandas>=1.5, scipy>=1.10, scikit-learn>=1.1. """ import numpy as np import pandas as pd from scipy.special import expit, logit from numpy.linalg import inv, cholesky from sklearn.preprocessing import PolynomialFeatures from sklearn.linear_model import RidgeCV from sklearn.pipeline import make_pipeline # --------------------------------------------------------------------------- # Global settings and fixed model inputs # --------------------------------------------------------------------------- MASTER_SEED = 20260619 WTP = 20000.0 strategies = np.array(['Exposed', 'Simple', 'Glue', 'Complex']) dressing_cost = np.array([0.0, 5.25, 13.86, 21.39]) ssi_cost_meanlog = 8.972237608 ssi_cost_sdlog = 0.163148238 Nops_per_year = 1200000 wounds_per_operation = 1.84 time_horizon_years = 4.673079209 population_wounds = Nops_per_year * wounds_per_operation * time_horizon_years risk_sources = { 'PHE_base_case': (0.137984898, 0.001772214), 'Jenks_sensitivity': (0.089390519, 0.006062126), } qaly_loss_grid = [0.06, 0.12] rho_grid = [-0.4, 0, 0.4] # Monte Carlo budgets ------------------------------------------------------ N_PSA = 20000 # inner PSA draws for EVPI / EVPPI N_REP_VOI = 8 # independent replicates for MC uncertainty (EVPI/EVPPI) EVSI_N_TOTAL = 2000 # planned future-trial total sample size EVSI_N_OUTER = 6000 # preposterior (outer) draws EVSI_N_INNER = 8000 # shared inner pool size (common random numbers) N_REP_EVSI = 8 # independent replicates for EVSI MC uncertainty # --------------------------------------------------------------------------- # Evidence inputs: odds ratios -> log-odds-ratio means / standard errors # --------------------------------------------------------------------------- def invert_or(or_, l, u): """Invert an odds ratio and its CI (effect expressed vs. Simple dressing).""" return 1 / or_, 1 / u, 1 / l def logor_params(or_, l, u, ci=0.95): """Return (log-OR mean, log-OR SE) from a point estimate and 95% CI.""" z = 1.959963984540054 return np.log(or_), (np.log(u) - np.log(l)) / (2 * z) raw = { 'Omitting_Bluebelle': dict(ex=(0.979, 0.561, 1.546), glue=(1.049, 0.371, 2.413), complex=(0.858, 0.535, 1.263)), 'Including_Bluebelle_ITT': dict(ex=(1.019, 0.641, 1.494), glue=(0.983, 0.480, 1.752), complex=(0.868, 0.563, 1.228)), 'Including_Bluebelle_PP': dict(ex=(1.003, 0.627, 1.478), glue=(0.999, 0.482, 1.828), complex=(0.869, 0.566, 1.240)), } effect_models = {} for k, r in raw.items(): ex = invert_or(*r['ex']) gl = invert_or(*r['glue']) cx = r['complex'] mus = [] ses = [] for vals in [ex, gl, cx]: mu, se = logor_params(*vals) mus.append(mu) ses.append(se) effect_models[k] = (np.array(mus), np.array(ses)) # order: Exposed, Glue, Complex # --------------------------------------------------------------------------- # Distribution helpers # --------------------------------------------------------------------------- def beta_ab_from_mean_se(mean, se): """Method-of-moments Beta(a, b) from a mean and standard error.""" var = se ** 2 tmp = mean * (1 - mean) / var - 1 return mean * tmp, (1 - mean) * tmp def cov_from_se_rho(se, rho): """Build a 3x3 covariance from marginal SEs and a common correlation rho.""" R = np.full((3, 3), rho) np.fill_diagonal(R, 1.0) return np.diag(se) @ R @ np.diag(se) def profile_to_pars(row): """Translate one admissible-measure profile into model parameters.""" mean_p0, se_p0 = risk_sources[row['risk_source']] mu, se = effect_models[row['effect_source']] shifts = np.array([row.get('logor_shift_exposed', 0.0), row.get('logor_shift_glue', 0.0), row.get('logor_shift_complex', 0.0)]) mult = row.get('effect_se_multiplier', 1.0) return dict(mean_p0=mean_p0, se_p0=se_p0, qaly_loss=row['qaly_loss'], logor_mu=mu + shifts, logor_cov=cov_from_se_rho(se * mult, row['rho_logor']), ssi_cost_multiplier=row.get('ssi_cost_multiplier', 1.0)) # --------------------------------------------------------------------------- # Core decision model: net benefit for each strategy # --------------------------------------------------------------------------- def simulate(row, n=20000, seed=1, posterior=None): """Simulate n PSA draws and return (nmb, draws). nmb is an n x 4 array of net benefits (Exposed, Simple, Glue, Complex); larger is preferred. If `posterior` is supplied, the log-OR mean/cov are replaced (used by the EVSI inner conditional model). """ rng = np.random.default_rng(seed) pars = profile_to_pars(row) if posterior is not None: pars['logor_mu'] = posterior['mu'] pars['logor_cov'] = posterior['cov'] a, b = beta_ab_from_mean_se(pars['mean_p0'], pars['se_p0']) p0 = rng.beta(a, b, size=n) c_ssi = rng.lognormal(ssi_cost_meanlog, ssi_cost_sdlog, size=n) * pars['ssi_cost_multiplier'] logor3 = rng.multivariate_normal(pars['logor_mu'], pars['logor_cov'], size=n) logor4 = np.column_stack([logor3[:, 0], np.zeros(n), logor3[:, 1], logor3[:, 2]]) p_ssi = expit(logit(p0)[:, None] + logor4) costs = p_ssi * c_ssi[:, None] + dressing_cost[None, :] qalys = -p_ssi * pars['qaly_loss'] nmb = WTP * qalys - costs draws = pd.DataFrame({'p0': p0, 'c_ssi': c_ssi, 'logor_exposed': logor3[:, 0], 'logor_glue': logor3[:, 1], 'logor_complex': logor3[:, 2]}) return nmb, draws def summarize(nmb): """Per-measure EVPI and decision summaries from a PSA sample.""" enb = nmb.mean(axis=0) evpi = np.mean(np.max(nmb, axis=1)) - np.max(enb) opt = strategies[np.argmax(enb)] gap = np.sort(enb)[-1] - np.sort(enb)[-2] pbest = np.bincount(np.argmax(nmb, axis=1), minlength=4) / nmb.shape[0] return enb, evpi, opt, gap, pbest def evppi_grouping(nmb, x, nbins=30, seed=1): """Approximate EVPPI. One-dimensional inputs use quantile grouping; multivariate inputs use a cross-fitted polynomial-ridge regression estimator for the conditional mean net benefit. This is an illustrative estimator; manuscript-grade analyses should use a validated EVPPI method appropriate for the model. """ n = nmb.shape[0] x = np.asarray(x) cur = np.max(nmb.mean(axis=0)) if x.ndim == 1: qs = np.quantile(x, np.linspace(0, 1, nbins + 1)) qs = np.unique(qs) if len(qs) <= 2: return 0.0 g = np.digitize(x, qs[1:-1], right=True) val = 0.0 for gr in np.unique(g): idx = (g == gr) val += idx.mean() * np.max(nmb[idx].mean(axis=0)) return max(0.0, float(val - cur)) xs = (x - x.mean(axis=0)) / x.std(axis=0) rng = np.random.default_rng(seed) fold = rng.integers(0, 2, size=n) pred = np.zeros_like(nmb) for k in range(nmb.shape[1]): for f in [0, 1]: train = fold != f test = fold == f model = make_pipeline(PolynomialFeatures(degree=2, include_bias=False), RidgeCV(alphas=np.logspace(-3, 3, 13))) model.fit(xs[train], nmb[train, k]) pred[test, k] = model.predict(xs[test]) val = np.mean(np.max(pred, axis=1)) return max(0.0, float(val - cur)) # --------------------------------------------------------------------------- # EVSI: common-random-number nested Monte Carlo on the normal-normal update # --------------------------------------------------------------------------- def approx_lor_likelihood_variance(row, n_total): """Approximate likelihood covariance for log odds ratios vs. a common control. Off-diagonal terms reflect the shared Simple-dressing control arm in a four-arm trial that allocates n_total/4 wounds per arm. """ pars = profile_to_pars(row) n_per = n_total / 4.0 p0 = pars['mean_p0'] p_treat = expit(logit(p0) + pars['logor_mu']) control_var = 1 / (n_per * p0) + 1 / (n_per * (1 - p0)) treat_var = 1 / (n_per * p_treat) + 1 / (n_per * (1 - p_treat)) V = np.full((3, 3), control_var) np.fill_diagonal(V, treat_var + control_var) return V def compute_evsi(row, n_total, n_outer=EVSI_N_OUTER, n_inner=EVSI_N_INNER, seed=1): """Common-random-number (CRN) nested-MC EVSI for a four-arm trial. The future trial updates the three relative-effect log-odds-ratios under a normal-normal model. The posterior covariance Spost is data-independent, so the only quantity that varies across simulated future datasets is the posterior MEAN. We therefore: 1. Draw ONE fixed inner pool of baseline risk p0, SSI cost c_ssi, and standardized normal effect draws Zstd, and reuse it for the current value and for every preposterior evaluation (common random numbers). 2. Draw outer preposterior data summaries z ~ N(mu0, S0 + V) and map each to its analytic posterior mean. 3. Evaluate the conditional expected net benefit of each strategy by shifting/scaling the SAME standardized effect pool to the posterior, take the row-wise best, and average. Sharing random numbers removes the large between-draw sampling noise that a naive nested estimator (fresh inner draws per posterior sample) introduces, so the estimator is low-variance and respects EVSI <= EVPPI <= EVPI. Returns a scalar EVSI estimate (per wound). """ rng = np.random.default_rng(seed) pars = profile_to_pars(row) mu0 = pars['logor_mu'] S0 = pars['logor_cov'] V = approx_lor_likelihood_variance(row, n_total) P0 = inv(S0) Pl = inv(V) Spost = inv(P0 + Pl) Lpost = cholesky(Spost) L0 = cholesky(S0) # Shared inner pool (common random numbers) ----------------------------- a, b = beta_ab_from_mean_se(pars['mean_p0'], pars['se_p0']) p0_base = rng.beta(a, b, size=n_inner) c_ssi_base = rng.lognormal(ssi_cost_meanlog, ssi_cost_sdlog, size=n_inner) * pars['ssi_cost_multiplier'] Zstd = rng.standard_normal((n_inner, 3)) lp0 = logit(p0_base) def inner_enb(mu, L): logor3 = mu + Zstd @ L.T logor4 = np.column_stack([logor3[:, 0], np.zeros(n_inner), logor3[:, 1], logor3[:, 2]]) p_ssi = expit(lp0[:, None] + logor4) nmb = WTP * (-p_ssi * pars['qaly_loss']) - (p_ssi * c_ssi_base[:, None] + dressing_cost[None, :]) return nmb.mean(axis=0) # Current value: effects from the prior, same inner pool current_value = float(np.max(inner_enb(mu0, L0))) # Preposterior outer loop ---------------------------------------------- zpred = rng.multivariate_normal(mu0, S0 + V, size=n_outer) post_means = ((P0 @ mu0)[None, :] + zpred @ Pl.T) @ Spost.T # n_outer x 3 vals = np.array([np.max(inner_enb(m, Lpost)) for m in post_means]) return float(vals.mean() - current_value) def compute_evsi_replicates(row, n_total, n_rep=N_REP_EVSI, base_seed=1): """Return (mean, mc_se) of EVSI over n_rep independent CRN replicates.""" vals = np.array([compute_evsi(row, n_total, seed=base_seed + 1000 * k) for k in range(n_rep)]) return float(vals.mean()), float(vals.std(ddof=1) / np.sqrt(n_rep)) # --------------------------------------------------------------------------- # Build the 36 admissible-measure profiles # --------------------------------------------------------------------------- def build_profiles(): profiles = [] pid = 1 for risk in risk_sources: for qloss in qaly_loss_grid: for eff in effect_models: for rho in rho_grid: profiles.append(dict(profile_id=pid, risk_source=risk, qaly_loss=qloss, effect_source=eff, rho_logor=rho, effect_se_multiplier=1.0, ssi_cost_multiplier=1.0, logor_shift_exposed=0.0, logor_shift_glue=0.0, logor_shift_complex=0.0, scenario_label='Main profile')) pid += 1 return pd.DataFrame(profiles) # --------------------------------------------------------------------------- # Main EVPI / EVPPI computation with Monte Carlo uncertainty # --------------------------------------------------------------------------- def run_evpi_evppi(profile_df): """Per-measure EVPI / EVPPI (point estimates) plus replicate MC SEs.""" rows = [] mc_rows = [] for _, r in profile_df.iterrows(): row = r.to_dict() pid = int(row['profile_id']) # Replicates for MC uncertainty on EVPI and EVPPI(effects) evpi_reps = [] evppi_eff_reps = [] for k in range(N_REP_VOI): nmb_k, draws_k = simulate(row, N_PSA, seed=3000 + pid + 7919 * k) _, evpi_k, _, _, _ = summarize(nmb_k) evpi_reps.append(evpi_k) evppi_eff_reps.append( evppi_grouping(nmb_k, draws_k[['logor_exposed', 'logor_glue', 'logor_complex']].values)) evpi_reps = np.array(evpi_reps) evppi_eff_reps = np.array(evppi_eff_reps) # Primary point estimate uses the first replicate's PSA sample nmb, draws = simulate(row, N_PSA, seed=3000 + pid) enb, evpi, opt, gap, pbest = summarize(nmb) evppi_baseline = evppi_grouping(nmb, draws.p0.values) evppi_cost = evppi_grouping(nmb, draws.c_ssi.values) evppi_effects = evppi_grouping(nmb, draws[['logor_exposed', 'logor_glue', 'logor_complex']].values) rows.append(dict(profile_id=pid, risk_source=row['risk_source'], qaly_loss=row['qaly_loss'], effect_source=row['effect_source'], rho_logor=row['rho_logor'], optimal_strategy=opt, EVPI=evpi, population_EVPI=evpi * population_wounds, decision_gap=gap, ENB_Exposed=enb[0], ENB_Simple=enb[1], ENB_Glue=enb[2], ENB_Complex=enb[3], pbest_Exposed=pbest[0], pbest_Simple=pbest[1], pbest_Glue=pbest[2], pbest_Complex=pbest[3], EVPPI_baseline=evppi_baseline, EVPPI_cost=evppi_cost, EVPPI_effects=evppi_effects)) mc_rows.append(dict(profile_id=pid, EVPI_mean=evpi_reps.mean(), EVPI_se=evpi_reps.std(ddof=1) / np.sqrt(N_REP_VOI), EVPPI_effects_mean=evppi_eff_reps.mean(), EVPPI_effects_se=evppi_eff_reps.std(ddof=1) / np.sqrt(N_REP_VOI))) return pd.DataFrame(rows), pd.DataFrame(mc_rows) # --------------------------------------------------------------------------- # Validation / simulation study (Priority 3) # --------------------------------------------------------------------------- def validation_study(): """Independent checks that the estimators behave as theory requires. Returns a tidy DataFrame of (check, quantity, value, reference, note). """ from scipy.stats import norm out = [] rng = np.random.default_rng(MASTER_SEED + 99) # --- Check 1: closed-form EVPI for a two-strategy normal model ---------- # Strategy A: NB = 0. Strategy B: NB = X, X ~ N(mu, s^2). # EVPI = E[max(0, X)] - max(0, mu) = s*phi(mu/s) + mu*Phi(mu/s) - max(0, mu). for mu, s in [(0.0, 1.0), (0.5, 1.0), (-0.5, 2.0)]: n = 4_000_000 X = rng.normal(mu, s, size=n) nmb = np.column_stack([np.zeros(n), X]) evpi_mc = np.mean(np.max(nmb, axis=1)) - np.max(nmb.mean(axis=0)) z = mu / s evpi_cf = s * norm.pdf(z) + mu * norm.cdf(z) - max(0.0, mu) out.append(dict(check='Closed-form EVPI (2-strategy normal)', quantity=f'mu={mu}, s={s}', value=round(evpi_mc, 4), reference=round(evpi_cf, 4), note=f'abs diff {abs(evpi_mc - evpi_cf):.4f}')) # --- Check 2: EVSI limits on the reference profile ---------------------- ref = dict(profile_id=14, risk_source='PHE_base_case', qaly_loss=0.12, effect_source='Including_Bluebelle_ITT', rho_logor=0.0, effect_se_multiplier=1.0, ssi_cost_multiplier=1.0, logor_shift_exposed=0.0, logor_shift_glue=0.0, logor_shift_complex=0.0) nmb_ref, draws_ref = simulate(ref, N_PSA, seed=4242) evppi_eff = evppi_grouping(nmb_ref, draws_ref[['logor_exposed', 'logor_glue', 'logor_complex']].values) _, evpi_ref, _, _, _ = summarize(nmb_ref) evsi_small = compute_evsi(ref, 20, seed=4243) # tiny trial -> ~0 evsi_mid = compute_evsi(ref, 2000, seed=4244) # the reported design evsi_large = compute_evsi(ref, 5_000_000, seed=4245) # huge trial -> EVPPI(effects) out.append(dict(check='EVSI lower limit (n_total=20)', quantity='EVSI', value=round(evsi_small, 2), reference='~0 (small)', note='small trials carry little information')) out.append(dict(check='EVSI at reported design (n_total=2000)', quantity='EVSI', value=round(evsi_mid, 2), reference=f'<= EVPPI_eff {evppi_eff:.1f}', note='must not exceed EVPPI(effects)')) out.append(dict(check='EVSI upper limit (n_total=5e6)', quantity='EVSI', value=round(evsi_large, 2), reference=f'-> EVPPI_eff {evppi_eff:.1f}', note='large trials approach EVPPI(effects)')) ordering_ok = (evsi_mid <= evppi_eff + 2) and (evppi_eff <= evpi_ref + 2) out.append(dict(check='Ordering at reference measure', quantity='EVSI<=EVPPI<=EVPI', value=f'{evsi_mid:.1f}<={evppi_eff:.1f}<={evpi_ref:.1f}', reference='True' if ordering_ok else 'False', note='within MC tolerance')) return pd.DataFrame(out) # --------------------------------------------------------------------------- # Driver # --------------------------------------------------------------------------- def main(): profile_df = build_profiles() base = profile_df[(profile_df.risk_source == 'PHE_base_case') & (profile_df.qaly_loss == 0.12) & (profile_df.effect_source == 'Including_Bluebelle_ITT') & (profile_df.rho_logor == 0)].iloc[0].to_dict() base_id = int(base['profile_id']) print('population_wounds', population_wounds) print('base_profile_id', base_id) # --- EVPI / EVPPI (with MC uncertainty) -------------------------------- res, mc_evpi = run_evpi_evppi(profile_df) base_row = res[res.profile_id == base_id].iloc[0] # --- EVSI across all 36 admissible measures ---------------------------- # One CRN estimate per measure builds the envelope (single-run MC SE is # small, ~GBP 2). Replicated MC uncertainty is then reported for the three # decision-relevant measures (reference and the two envelope endpoints). evsi_records = [] for _, r in profile_df.iterrows(): row = r.to_dict() pid = int(row['profile_id']) evsi_val = compute_evsi(row, EVSI_N_TOTAL, seed=25000 + pid) evsi_records.append(dict(profile_id=pid, risk_source=row['risk_source'], qaly_loss=row['qaly_loss'], effect_source=row['effect_source'], rho_logor=row['rho_logor'], EVSI=evsi_val, population_EVSI=evsi_val * population_wounds)) main_evsi = pd.DataFrame(evsi_records) # Replicated MC uncertainty for reference + envelope-endpoint measures evsi_lo_pid = int(main_evsi.loc[main_evsi.EVSI.idxmin()].profile_id) evsi_hi_pid = int(main_evsi.loc[main_evsi.EVSI.idxmax()].profile_id) evsi_se_map = {} for pid in {base_id, evsi_lo_pid, evsi_hi_pid}: row = profile_df[profile_df.profile_id == pid].iloc[0].to_dict() mean_pid, se_pid = compute_evsi_replicates(row, EVSI_N_TOTAL, base_seed=80000 + pid) evsi_se_map[pid] = (mean_pid, se_pid) # --- Assemble endpoint / reference uncertainty table ------------------- def endpoint(df, col, which): idx = df[col].idxmax() if which == 'max' else df[col].idxmin() return df.loc[idx] evpi_ref = base_row evpi_lo = endpoint(res, 'EVPI', 'min') evpi_hi = endpoint(res, 'EVPI', 'max') evsi_ref = main_evsi[main_evsi.profile_id == base_id].iloc[0] evsi_lo = endpoint(main_evsi, 'EVSI', 'min') evsi_hi = endpoint(main_evsi, 'EVSI', 'max') def mc_se_for(pid, col): r = mc_evpi[mc_evpi.profile_id == int(pid)].iloc[0] return r[col] def evsi_se_for(pid): return evsi_se_map.get(int(pid), (np.nan, np.nan))[1] unc_rows = [ dict(target='EVPI', role='reference (P0)', profile_id=int(evpi_ref.profile_id), per_wound=evpi_ref.EVPI, mc_se=mc_se_for(evpi_ref.profile_id, 'EVPI_se'), population=evpi_ref.EVPI * population_wounds), dict(target='EVPI', role='envelope lower', profile_id=int(evpi_lo.profile_id), per_wound=evpi_lo.EVPI, mc_se=mc_se_for(evpi_lo.profile_id, 'EVPI_se'), population=evpi_lo.EVPI * population_wounds), dict(target='EVPI', role='envelope upper', profile_id=int(evpi_hi.profile_id), per_wound=evpi_hi.EVPI, mc_se=mc_se_for(evpi_hi.profile_id, 'EVPI_se'), population=evpi_hi.EVPI * population_wounds), dict(target='EVPPI_effects', role='reference (P0)', profile_id=int(evpi_ref.profile_id), per_wound=evpi_ref.EVPPI_effects, mc_se=mc_se_for(evpi_ref.profile_id, 'EVPPI_effects_se'), population=evpi_ref.EVPPI_effects * population_wounds), dict(target='EVPPI_effects', role='envelope lower', profile_id=int(res.loc[res.EVPPI_effects.idxmin()].profile_id), per_wound=res.EVPPI_effects.min(), mc_se=np.nan, population=res.EVPPI_effects.min() * population_wounds), dict(target='EVPPI_effects', role='envelope upper', profile_id=int(res.loc[res.EVPPI_effects.idxmax()].profile_id), per_wound=res.EVPPI_effects.max(), mc_se=np.nan, population=res.EVPPI_effects.max() * population_wounds), dict(target='EVSI', role='reference (P0)', profile_id=int(evsi_ref.profile_id), per_wound=evsi_ref.EVSI, mc_se=evsi_se_for(evsi_ref.profile_id), population=evsi_ref.EVSI * population_wounds), dict(target='EVSI', role='envelope lower', profile_id=int(evsi_lo.profile_id), per_wound=evsi_lo.EVSI, mc_se=evsi_se_for(evsi_lo.profile_id), population=evsi_lo.EVSI * population_wounds), dict(target='EVSI', role='envelope upper', profile_id=int(evsi_hi.profile_id), per_wound=evsi_hi.EVSI, mc_se=evsi_se_for(evsi_hi.profile_id), population=evsi_hi.EVSI * population_wounds), ] unc = pd.DataFrame(unc_rows) # --- Validation study -------------------------------------------------- validation = validation_study() # --- Console summary --------------------------------------------------- pd.set_option('display.width', 160) print('\nBase EVPI', base_row[['optimal_strategy', 'EVPI', 'population_EVPI', 'decision_gap']].to_dict()) print('\nEVPI envelope per wound: [%.1f, %.1f]' % (res.EVPI.min(), res.EVPI.max())) print('EVPPI(effects) envelope per wound: [%.1f, %.1f]' % (res.EVPPI_effects.min(), res.EVPPI_effects.max())) print('decision gap range: [%.1f, %.1f]' % (res.decision_gap.min(), res.decision_gap.max())) print('\nEVSI envelope per wound: [%.1f, %.1f]; reference %.1f' % (main_evsi.EVSI.min(), main_evsi.EVSI.max(), evsi_ref.EVSI)) print('\n=== Reference / envelope endpoints with Monte Carlo SE ===') print(unc.to_string(index=False)) print('\n=== Validation / simulation study ===') print(validation.to_string(index=False)) # --- Write outputs ----------------------------------------------------- res.to_csv('bluebelle_main_results.csv', index=False) main_evsi.to_csv('bluebelle_main_evsi.csv', index=False) unc.to_csv('bluebelle_voi_uncertainty.csv', index=False) validation.to_csv('bluebelle_validation.csv', index=False) print('\nWrote: bluebelle_main_results.csv, bluebelle_main_evsi.csv, ' 'bluebelle_voi_uncertainty.csv, bluebelle_validation.csv') if __name__ == '__main__': main()