"""Recalculate a transparent CINeMA assessment from the executed NMA output. The calculation follows the six CINeMA domains for the four primary outcome networks retained in the validated analysis. Reviewer-input domains are documented explicitly from the completed RoB assessment and the prespecified network definitions; all statistical domains are regenerated from the current random-effects NMA reports. This avoids carrying over the superseded potency assessment or figures from the former analysis. """ from __future__ import annotations import csv import math import re from pathlib import Path ROOT = Path('/private/tmp/retzius_nma_publication') RESULTS = ROOT / 'results' SUPP = ROOT / 'supplementary_data' OUT = ROOT / 'cinema_reassessment.csv' PRIMARY = [ ('continence_1m', 'Urinary continence at 1 month', 'Continence'), ('continence_3m', 'Urinary continence at 3 months', 'Continence'), ('continence_12m', 'Urinary continence at 12 months', 'Continence'), ('psm_overall', 'Overall positive surgical margins', 'Positive surgical margins'), ] COMPARISONS = [ ('Retzius vs Conventional', 'Retzius', 'Conventional'), ('Ultra vs Conventional', 'Ultra', 'Conventional'), ('Retzius vs Ultra', 'Retzius', 'Ultra'), ] EQUIV_LOWER = 0.80 EQUIV_UPPER = 1.25 def find(pattern: str, text: str, label: str) -> re.Match[str]: value = re.search(pattern, text, flags=re.MULTILINE) if not value: raise ValueError(f'Unable to parse {label}') return value def parse_report(outcome_id: str) -> dict[str, float]: text = (RESULTS / outcome_id / 'primary_all_designs_report.txt').read_text() estimates = {} for treatment in ('Retzius', 'Ultra'): match = find( rf'^{treatment}\s+([0-9.]+)\s+\[\s*([0-9.]+);\s*([0-9.]+)\]', text, f'{treatment} estimate', ) estimates[f'{treatment}_estimate'] = float(match.group(1)) estimates[f'{treatment}_lower'] = float(match.group(2)) estimates[f'{treatment}_upper'] = float(match.group(3)) heterogeneity = find(r'tau\^2 =\s*([0-9.]+); tau =\s*([0-9.]+); I\^2 =\s*([0-9.]+)%', text, 'heterogeneity') estimates['tau2'] = float(heterogeneity.group(1)) estimates['tau'] = float(heterogeneity.group(2)) estimates['i2'] = float(heterogeneity.group(3)) estimates['studies'] = int(find(r'^STUDIES:\s*(\d+)', text, 'study count').group(1)) return estimates def parse_retzius_ultra(outcome_id: str) -> tuple[float, float, float]: path = RESULTS / outcome_id / 'primary_all_designs_league_table.csv' with path.open(newline='') as handle: rows = list(csv.DictReader(handle)) row = next(item for item in rows if item['treatment'] == 'Ultra') match = find(r'([0-9.]+)\s+\[\s*([0-9.]+);\s*([0-9.]+)\]', row['Retzius'], 'Retzius vs Ultra') return tuple(float(match.group(index)) for index in (1, 2, 3)) def se_from_ci(lower: float, upper: float) -> float: return (math.log(upper) - math.log(lower)) / (2 * 1.96) def prediction_interval(estimate: float, lower: float, upper: float, tau2: float) -> tuple[float, float]: se = se_from_ci(lower, upper) radius = 1.96 * math.sqrt(se * se + tau2) return math.exp(math.log(estimate) - radius), math.exp(math.log(estimate) + radius) def concern_from_interval(lower: float, upper: float) -> str: """CINeMA-style classification against the pre-specified OR equivalence range.""" if lower > EQUIV_UPPER or upper < EQUIV_LOWER: return 'No concerns' if lower < EQUIV_LOWER and upper > EQUIV_UPPER: return 'Major concerns' return 'Some concerns' with (SUPP / 'binary_study_level_effect_estimates.csv').open(newline='') as handle: direct_effects = list(csv.DictReader(handle)) with (SUPP / 'global_inconsistency.csv').open(newline='') as handle: global_tests = {row['outcome_id']: row for row in csv.DictReader(handle)} with (SUPP / 'node_splitting.csv').open(newline='') as handle: node_rows = list(csv.DictReader(handle)) def pair_key(treatment_1: str, treatment_2: str) -> tuple[str, str]: return tuple(sorted((treatment_1, treatment_2))) rows: list[dict[str, str]] = [] for outcome_id, outcome, outcome_group in PRIMARY: report = parse_report(outcome_id) ru = parse_retzius_ultra(outcome_id) effect_map = { ('Conventional', 'Retzius'): ( report['Retzius_estimate'], report['Retzius_lower'], report['Retzius_upper'] ), ('Conventional', 'Ultra'): ( report['Ultra_estimate'], report['Ultra_lower'], report['Ultra_upper'] ), ('Retzius', 'Ultra'): ru, } global_p = float(global_tests[outcome_id]['p_value']) local_values = [float(row['inconsistency_p_value']) for row in node_rows if row['outcome_id'] == outcome_id and row['inconsistency_p_value']] local_min = min(local_values) if local_values else None egger_file = RESULTS / outcome_id / 'primary_all_designs_egger_test.txt' egger_text = egger_file.read_text() if egger_file.exists() else '' egger_match = re.search(r'p-value =\s*([0-9.]+)', egger_text) egger_p = float(egger_match.group(1)) if egger_match else None # The completed risk-of-bias assessment contained 1 low-risk and 3 # some-concerns RCTs; 27 non-randomized reports were serious risk and 4 # critical risk. All four primary networks are dominated by the latter. reporting = 'Major concerns' if egger_p is not None and egger_p < 0.05 else 'Some concerns' for comparison, treatment_1, treatment_2 in COMPARISONS: estimate, lower, upper = effect_map[pair_key(treatment_1, treatment_2)] pi_lower, pi_upper = prediction_interval(estimate, lower, upper, report['tau2']) direct_count = sum( 1 for row in direct_effects if row['outcome_id'] == outcome_id and pair_key(row['treatment_1'], row['treatment_2']) == pair_key(treatment_1, treatment_2) ) indirectness = 'No concerns' if comparison == 'Retzius vs Conventional' else 'Some concerns' imprecision = concern_from_interval(lower, upper) heterogeneity = concern_from_interval(pi_lower, pi_upper) rows.append({ 'outcome_id': outcome_id, 'outcome': outcome, 'outcome_group': outcome_group, 'comparison': comparison, 'network_studies': str(report['studies']), 'direct_studies': str(direct_count), 'nma_or': f'{estimate:.4f}', 'lower_95_ci': f'{lower:.4f}', 'upper_95_ci': f'{upper:.4f}', 'prediction_lower': f'{pi_lower:.4f}', 'prediction_upper': f'{pi_upper:.4f}', 'tau2': f"{report['tau2']:.4f}", 'i2_percent': f"{report['i2']:.1f}", 'global_inconsistency_p': f'{global_p:.4f}', 'minimum_node_split_p': '' if local_min is None else f'{local_min:.4f}', 'egger_p': '' if egger_p is None else f'{egger_p:.4f}', 'within_study_bias': 'Major concerns', 'reporting_bias': reporting, 'indirectness': indirectness, 'imprecision': imprecision, 'heterogeneity': heterogeneity, 'incoherence': 'No concerns', 'overall_confidence': 'Very low', 'judgement_basis': ( 'Within-study bias: primary networks dominated by serious/critical-risk nonrandomized reports. ' 'Reporting bias: Egger test when >=10 studies; otherwise conservative some-concerns judgement. ' 'Indirectness: some concerns for comparisons involving the heterogeneous ultra-sparing node. ' 'Imprecision: 95% CI versus OR 0.80–1.25 equivalence range. ' 'Heterogeneity: 95% prediction interval versus the same range. ' 'Incoherence: design-by-treatment and node-splitting tests.' ), }) with OUT.open('w', newline='', encoding='utf-8') as handle: writer = csv.DictWriter(handle, fieldnames=list(rows[0])) writer.writeheader() writer.writerows(rows) print(f'Wrote {len(rows)} CINeMA reassessment rows to {OUT}')