#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Supplementary File S6 - Reproducible analysis code Manuscript: "Serum iron and total Kt/V in peritoneal dialysis-associated peritonitis: a cross-sectional study" Input : Data File S5 (de-identified dataset, 170 observations; deposited at Zenodo, DOI: 10.5281/zenodo.22816756, restricted access) Output: all summary statistics and regression estimates reported in Table 1, Table 2, Table 3 and Supplementary Tables S1, S3, and S4. Also outputs: VIF for fully adjusted model, iron-CRP correlations, robust Poisson regression (RR), iron and Kt/V tipping-point analyses, cluster bootstrap by patient, and OR-based E-values. Also outputs: Iron x log-CRP interaction (Wald + LRT). Also outputs: Firth penalised logistic regression (unadjusted and fully adjusted specifications) and the comparison of case characteristics by total Kt/V measurement availability (Supplementary Table S3b). Also outputs: the finite-sample correction comparison for cluster-robust SEs (Section 3f) and the pre-specified, data-gated pre-admission Kt/V restriction (Section 3g). The pathogen spectrum (Table S2, Fig. 4) is summarised from clinical microbiology records and is not contained in the S5 dataset. Python 3.10+, pandas, numpy, scipy, scikit-learn, statsmodels. No random seed ambiguity: multiple imputation uses random_state = 2024 + imputation index. """ import numpy as np, pandas as pd from scipy import stats from scipy.optimize import minimize from sklearn.experimental import enable_iterative_imputer # noqa: F401 from sklearn.impute import IterativeImputer from sklearn.linear_model import BayesianRidge import statsmodels.api as sm_api # also used under alias sm_glm below sm_glm = sm_api # ---------------------------------------------------------------------- # 1. Load data # ---------------------------------------------------------------------- # NOTE: Data File S5 contains the data-quality corrections documented in the # Supplementary Codebook (File S7): three EF values re-expressed as # percentages (obs_id 14/35/36) and five implausible values set to missing # (TIBC, phosphorus, sdLDL, KtV, BMI). df = pd.read_csv("09_Supplementary_Data_S5.csv") df["y"] = (df["group"] == "case").astype(float) # pooled standard deviations used for per-1-SD exposure scaling SD_IRON = df["iron"].std(ddof=1) # 6.367 umol/L SD_KTV = df["KtV"].std(ddof=1) # 0.672 df["iron_sd"] = df["iron"] / SD_IRON df["ktv_sd"] = df["KtV"] / SD_KTV df["male"] = df["sex"] # ---------------------------------------------------------------------- # 2. Descriptive statistics (Table 1 / Supplementary Table S1) # mean +/- SD (sample SD, ddof=1) when Shapiro-Wilk p>0.05, else # median (IQR) with numpy/pandas linear-interpolation quantiles. # ---------------------------------------------------------------------- def describe(col): out = {} for g, label in [("case", "cases"), ("ctrl", "controls")]: x = df.loc[df.group == g, col].dropna() W, pn = stats.shapiro(x) if len(x) >= 3 else (np.nan, 1.0) if pn > 0.05: out[label] = f"{x.mean():.2f} \u00b1 {x.std(ddof=1):.2f}" else: out[label] = (f"{x.median():.2f} " f"({x.quantile(0.25):.2f}\u2013{x.quantile(0.75):.2f})") out[label + "_n"] = len(x) a = df.loc[df.y == 1, col].dropna() b = df.loc[df.y == 0, col].dropna() _, p = stats.mannwhitneyu(a, b, alternative="two-sided") out["p"] = p return out # FDR screening: Mann-Whitney U + Benjamini-Hochberg over 57 candidates EXCLUDE = {"group", "is_case", "y", "patient_code", "obs_id", "bacteria", "antibiotic", "LOS", "effluent_WBC", "effluent_RBC", "patient_status", "group_binary", "log_crp", "log_crp_sd", "iron_sd", "ktv_sd", "male"} cands = [c for c in df.columns if c not in EXCLUDE and df[c].dtype.kind in "if" and df[c].notna().sum() > 20] rows = [] for v in cands: a = df.loc[df.y == 1, v].dropna() b = df.loc[df.y == 0, v].dropna() _, p = stats.mannwhitneyu(a, b, alternative="two-sided") rows.append((v, len(a), len(b), p)) scr = pd.DataFrame(rows, columns=["var", "n_case", "n_ctrl", "p"]) scr = scr.sort_values("p").reset_index(drop=True) m = len(scr) q = (scr["p"].values * m / np.arange(1, m + 1)) q = np.minimum.accumulate(q[::-1])[::-1] scr["q"] = np.minimum(q, 1) scr["sig"] = scr["q"] < 0.05 # 27 variables significant print("\n=== Univariate screening (Mann-Whitney U, BH FDR over all candidates) ===") with pd.option_context("display.max_rows", None): print(scr[["var", "n_case", "n_ctrl", "p", "q", "sig"]].to_string(index=False)) # Output descriptive statistics for Table 1 / S1 print("\n=== Descriptive Statistics (Table 1 / S1) ===") TABLE1_VARS = ["age", "sex", "hemoglobin", "WBC", "platelet", "CRP", "PCT", "iron", "TIBC", "TSAT", "ferritin", "albumin", "prealbumin", "calcium", "phosphorus", "PTH", "KtV", "peritoneal_KtV", "residual_KtV", "residual_Ccr", "urine_vol", "total_Ccr", "potassium", "sodium", "chloride", "BMI", "ESR", "IL6"] # Only clearance variables carry a unit; Kt/V (total/residual/peritoneal) is a # dimensionless ratio and is intentionally left unlabelled. UNITS = {"residual_Ccr": "L/week/1.73 m²", "total_Ccr": "L/week/1.73 m²"} for v in TABLE1_VARS: if v in df.columns: desc = describe(v) unit_tag = f" [{UNITS[v]}]" if v in UNITS else "" print(f"\n{v}{unit_tag}:") print(f" Cases: {desc.get('cases', 'N/A')} (n={desc.get('cases_n', 'N/A')})") print(f" Controls: {desc.get('controls', 'N/A')} (n={desc.get('controls_n', 'N/A')})") print(f" P value (MWU): {desc['p']:.3f}") # ---------------------------------------------------------------------- # 3. Logistic regression with cluster-robust (sandwich) standard errors # ---------------------------------------------------------------------- def fit_logit(d, xcols, cluster_col="patient_code"): X1 = np.column_stack([np.ones(len(d))] + [d[c].values for c in xcols]) y = d["y"].values def nll(b): z = np.clip(X1 @ b, -30, 30) p = np.clip(1 / (1 + np.exp(-z)), 1e-10, 1 - 1e-10) return -np.sum(y * np.log(p) + (1 - y) * np.log(1 - p)) def grad(b): z = np.clip(X1 @ b, -30, 30) p = 1 / (1 + np.exp(-z)) return X1.T @ (p - y) r = minimize(nll, np.zeros(X1.shape[1]), jac=grad, method="BFGS", options={"maxiter": 3000, "ftol": 1e-12}) b = r.x z = np.clip(X1 @ b, -30, 30) p = 1 / (1 + np.exp(-z)) W = p * (1 - p) H = X1.T @ (W[:, None] * X1) Hinv = np.linalg.inv(H) cov_model = Hinv # model-based (iid) SE if cluster_col is not None and cluster_col in d.columns: # No finite-sample correction factor (e.g. G/(G-1)) is applied to the # clustered sandwich SE; with 102-111 patient clusters per model its # effect is modest, and finite-sample-corrected variants of the fully # adjusted estimates are reported in Section 3f. meat = np.zeros_like(H) # clustered sandwich SE for g_ in np.unique(d[cluster_col].values): idx = np.where(d[cluster_col].values == g_)[0] sc = X1[idx].T @ ((y - p)[idx]) meat += np.outer(sc, sc) cov_cluster = Hinv @ meat @ Hinv else: cov_cluster = cov_model # independent obs (imputed/dedup sets) return b, cov_model, cov_cluster def or_row(b, cov, i): se = np.sqrt(cov[i, i]) z = b[i] / se pv = 2 * stats.norm.sf(abs(z)) return dict(OR=np.exp(b[i]), lo=np.exp(b[i] - 1.96 * se), hi=np.exp(b[i] + 1.96 * se), P=pv) SPECS = { "unadjusted": (["iron_sd", "ktv_sd"], ["iron", "KtV"]), "adjusted": (["iron_sd", "ktv_sd", "male", "age", "hemoglobin", "albumin"], ["iron", "KtV", "sex", "age", "hemoglobin", "albumin"]), "fully": (["iron_sd", "ktv_sd", "male", "age", "hemoglobin", "albumin", "log_crp"], ["iron", "KtV", "sex", "age", "hemoglobin", "albumin", "log_crp"]), } for name, (xcols, req) in SPECS.items(): d = df.dropna(subset=req) b, cov_model, cov_cl = fit_logit(d, xcols) print(f"\n{name}: n={len(d)} clusters={d.patient_code.nunique()}") for i, c in enumerate(xcols, start=1): rc = or_row(b, cov_cl, i) rm = or_row(b, cov_model, i) print(f" {c:10s} cluster OR={rc['OR']:.3f} " f"({rc['lo']:.3f}-{rc['hi']:.3f}) P={rc['P']:.4f} | " f"model OR={rm['OR']:.3f} ({rm['lo']:.3f}-{rm['hi']:.3f})") # ---------------------------------------------------------------------- # 3b. VIF calculation for fully adjusted model # ---------------------------------------------------------------------- print("\n=== VIF (Fully Adjusted Model) ===") d_vif = df.dropna(subset=["iron", "KtV", "sex", "age", "hemoglobin", "albumin", "CRP"]).copy() # log_crp is the CSV column ln(CRP+1), identical transform to the primary regression X_vif = d_vif[["iron_sd", "ktv_sd", "male", "age", "hemoglobin", "albumin", "log_crp"]] # Calculate VIF for each predictor from numpy.linalg import inv corr_matrix = X_vif.corr().values inv_corr = inv(corr_matrix) vif_values = np.diag(inv_corr) var_names = ["iron_sd", "ktv_sd", "male", "age", "hemoglobin", "albumin", "log_crp"] for var, vif in zip(var_names, vif_values): print(f" {var}: VIF = {vif:.2f}") # ---------------------------------------------------------------------- # 3c. Iron × log-CRP interaction (Wald + likelihood ratio test) # Same fully adjusted 7-variable specification and same cluster-robust # SEs as the primary analysis (Section 3, "fully" model): # iron_sd, ktv_sd, sex, age, hemoglobin, albumin, log_crp [= ln(CRP+1), # the CSV column], plus the iron_sd × log_crp product. Complete-case # set matches the primary fully adjusted model (n = 109). # ---------------------------------------------------------------------- df_cc = df.dropna(subset=["iron_sd", "ktv_sd", "sex", "age", "hemoglobin", "albumin", "log_crp", "patient_code"]).copy() df_cc["ix"] = df_cc["iron_sd"] * df_cc["log_crp"] base = ["iron_sd", "ktv_sd", "sex", "age", "hemoglobin", "albumin", "log_crp"] # Reduced (no interaction) and full (with interaction) nested models X0 = sm_api.add_constant(df_cc[base]) X1 = sm_api.add_constant(df_cc[base + ["ix"]]) m0 = sm_api.Logit(df_cc["y"], X0).fit(disp=0) m1 = sm_api.Logit(df_cc["y"], X1).fit( disp=0, cov_type="cluster", cov_kwds={"groups": df_cc["patient_code"], "use_correction": False}) b, se = m1.params["ix"], m1.bse["ix"] print(f"\n=== Iron x log-CRP Interaction (fully adjusted, n={len(df_cc)}) ===") print(f" interaction OR = {np.exp(b):.2f} " f"({np.exp(b - 1.96*se):.2f}-{np.exp(b + 1.96*se):.2f}), " f"Wald P = {m1.pvalues['ix']:.3f}") lrt = 2 * (m1.llf - m0.llf) print(f" LRT chi2 = {lrt:.3f}, df = 1, " f"P = {1 - stats.chi2.cdf(lrt, 1):.3f}") # -> interaction OR = 1.04 (0.50-2.17), Wald P = 0.913; # LRT chi2 = 0.013, P = 0.910 # ---------------------------------------------------------------------- # 3d. Firth penalised logistic regression (Jeffreys-prior penalised # likelihood; Heinze-Schemper modified-score iteration). Small-event # sensitivity analysis for the unadjusted and fully adjusted # specifications; Wald-type 95% CIs from the penalised information. # ---------------------------------------------------------------------- print("\n=== Firth Penalised Logistic Regression ===") def firth_logit(X1, y, tol=1e-9, maxiter=500): X1 = np.asarray(X1, float); y = np.asarray(y, float) b = np.zeros(X1.shape[1]) for it in range(maxiter): eta = np.clip(X1 @ b, -30, 30) p = 1.0 / (1.0 + np.exp(-eta)) W = p * (1 - p) I = X1.T @ (W[:, None] * X1) Iinv = np.linalg.pinv(I) h = np.einsum("ij,jk,ik->i", X1, Iinv, X1) * W # hat-matrix diagonal U = X1.T @ (y - p + h * (0.5 - p)) # modified score step = Iinv @ U b = b + step if np.max(np.abs(step)) < tol: break eta = np.clip(X1 @ b, -30, 30) p = 1.0 / (1.0 + np.exp(-eta)) W = p * (1 - p) I = X1.T @ (W[:, None] * X1) cov = np.linalg.pinv(I) return b, cov, it for _name, (_xcols, _req) in SPECS.items(): if _name not in ("unadjusted", "fully"): continue _d = df.dropna(subset=_req) _X1 = np.column_stack([np.ones(len(_d))] + [_d[c].values for c in _xcols]) _b, _cov, _it = firth_logit(_X1, _d["y"].values) print(f"{_name}: n={len(_d)}, events={int(_d.y.sum())}") for _i, _c in enumerate(_xcols[:2], start=1): _se = np.sqrt(_cov[_i, _i]); _z = _b[_i] / _se print(f" {_c:8s} Firth OR={np.exp(_b[_i]):.3f} " f"({np.exp(_b[_i]-1.96*_se):.3f}-{np.exp(_b[_i]+1.96*_se):.3f}) " f"P={2*stats.norm.sf(abs(_z)):.4f}") # -> unadjusted: iron 0.236 (0.126-0.443); KtV 0.392 (0.221-0.695) # -> fully: iron 0.482 (0.212-1.099); KtV 0.434 (0.189-0.995) # ---------------------------------------------------------------------- # 3e. Parsimonious adjusted model (age + sex + the two exposures only). # Addressed to EPV concerns: with 4 predictors and ~43 events the # events-per-variable ratio rises to ~10.7, trading confounding # control for estimator stability. Same standardisation (per pooled # SD) and same cluster-robust sandwich SEs as the primary analysis. # ---------------------------------------------------------------------- print("\n=== Parsimonious Adjusted Model (iron_sd + ktv_sd + age + male) ===") _preq = ["iron", "KtV", "sex", "age"] _d = df.dropna(subset=_preq) _b, _cm, _cc = fit_logit(_d, ["iron_sd", "ktv_sd", "age", "male"]) print(f"parsimonious: n={len(_d)} clusters={_d.patient_code.nunique()} " f"events={int(_d.y.sum())} EPV={_d.y.sum()/4:.1f}") for _i, _c in enumerate(["iron_sd", "ktv_sd", "age", "male"], start=1): _r = or_row(_b, _cc, _i) print(f" {_c:8s} cluster OR={_r['OR']:.3f} " f"({_r['lo']:.3f}-{_r['hi']:.3f}) P={_r['P']:.4f}") # -> parsimonious n=116, 44 events, EPV=11.0: # iron OR=0.244 (0.112-0.534) P=0.0004; KtV OR=0.329 (0.175-0.621) P=0.0006 # ---------------------------------------------------------------------- # 3f. Cluster-robust SE small-sample corrections (fully adjusted model). # The primary sandwich SE uses no finite-sample correction; with 43 # events and 102 clusters we compare two corrections: # (a) Stata-style multiplicative correction G/(G-1) * (N-1)/(N-K); # (b) CR3 leverage adjustment (each cluster score divided by 1 - h_g, # h_g = cluster leverage from W^{1/2} X (X'WX)^{-1} X' W^{1/2}). # ---------------------------------------------------------------------- print("\n=== Cluster-SE Small-Sample Corrections (fully adjusted) ===") _fareq = ["iron", "KtV", "sex", "age", "hemoglobin", "albumin", "log_crp"] _d = df.dropna(subset=_fareq) _xcols = ["iron_sd", "ktv_sd", "male", "age", "hemoglobin", "albumin", "log_crp"] _b, _cm, _cc = fit_logit(_d, _xcols) _X1 = np.column_stack([np.ones(len(_d))] + [_d[c].values for c in _xcols]) _y = _d["y"].values _z = np.clip(_X1 @ _b, -30, 30) _p = 1 / (1 + np.exp(-_z)) _W = _p * (1 - _p) _H = _X1.T @ (_W[:, None] * _X1) _Hinv = np.linalg.inv(_H) _G = _d.patient_code.nunique() _N, _K = _X1.shape _fc = _G / (_G - 1) * (_N - 1) / (_N - _K) # CR3: cluster leverage from hat matrix _A = np.sqrt(_W)[:, None] * _X1 @ _Hinv @ (_np_sqrtwX := (_X1 * np.sqrt(_W)[:, None]).T) # simpler: diag of hat = row-wise h_i _hhat = np.einsum("ij,jk,ik->i", _X1 * np.sqrt(_W)[:, None], _Hinv, _X1 * np.sqrt(_W)[:, None]) _meat_stata = np.zeros_like(_H) _meat_cr3 = np.zeros_like(_H) for g_ in np.unique(_d.patient_code.values): idx = np.where(_d.patient_code.values == g_)[0] sc = _X1[idx].T @ ((_y - _p)[idx]) _meat_stata += np.outer(sc, sc) _hg = _hhat[idx].sum() _meat_cr3 += np.outer(sc, sc) / (1 - _hg) ** 2 _cov_stata = _fc * _Hinv @ _meat_stata @ _Hinv _cov_cr3 = _Hinv @ _meat_cr3 @ _Hinv print(f"clusters={_G} n={_N} K={_K}; finite-correction factor={_fc:.4f}") for _lbl, _cv in [("primary (uncorrected)", _cc), ("Stata G/(G-1)*(N-1)/(N-K)", _cov_stata), ("CR3 leverage", _cov_cr3)]: for _i, _c in enumerate(["iron", "Kt/V"], start=1): _r = or_row(_b, _cv, _i) print(f" {_lbl:26s} {_c:5s} OR={_r['OR']:.3f} ({_r['lo']:.3f}-{_r['hi']:.3f}) P={_r['P']:.4f}") # -> corrections widen CIs modestly (CR3 Kt/V P = 0.048, borderline); conclusions qualitatively unchanged # ---------------------------------------------------------------------- # 3g. Pre-admission Kt/V restriction (data-gated; pre-specified) # Runs only if a boolean column `ktv_post_admission` (True = the case # episode's Kt/V value was obtained after that admission) is appended # to Data File S5. Case episodes are then restricted to pre-admission Kt/V # measurements and the fully adjusted model is refit with the primary # cluster-robust SEs; controls enter unchanged. # # Classification of cases with a Kt/V value but NO timing record: # they are NOT assumed to be pre-admission measurements (a naive # fillna(False) would silently fold them into the restricted sample # and bias it toward cases of unknown timing). They are EXCLUDED from # the restricted analysis and reported as a separate count, so the # restricted sample contains only episodes with a documented timing # record. (Alternative: mark them as a third category and run an # additional two-way sensitivity analysis.) # ---------------------------------------------------------------------- if "ktv_post_admission" in df.columns: _flag = df["ktv_post_admission"] _has_ktv = df["KtV"].notna() _no_time = (df["y"] == 1) & _has_ktv & _flag.isna() _d_pre = df[~(_flag.fillna(False).astype(bool)) & ~_no_time].dropna(subset=_fareq) print("\n=== 3g. Pre-admission Kt/V restriction ===") print(f"timing flags present for {int((_has_ktv & _flag.notna()).sum())} of " f"{int(_has_ktv.sum())} Kt/V values; {int(_no_time.sum())} case values " "without a timing record are EXCLUDED (not assumed pre-admission)") print(f"n={len(_d_pre)} cases={int(_d_pre['y'].sum())}") _bg, _cmg, _ccg = fit_logit(_d_pre, _xcols) for _i2, _c2 in enumerate(["iron", "Kt/V"], start=1): _r2 = or_row(_bg, _ccg, _i2) print(f" primary-SE {_c2:5s} OR={_r2['OR']:.3f} ({_r2['lo']:.3f}-{_r2['hi']:.3f}) P={_r2['P']:.4f}") else: print("\n=== 3g. Pre-admission Kt/V restriction: SKIPPED " "(ktv_post_admission column not present in Data File S5) ===") # ---------------------------------------------------------------------- # 4. Sensitivity: one observation per patient (first/earliest episode). # The source registry is ordered by admission date, so within each # patient the earliest observation is the first row (obs_id order). # Dedup is performed within the complete-case set requiring BOTH # exposures (n = 120 observations, 111 patients: 40 case patients, # 71 control patients); sequential models are then refit on this # independent-observation set with model-based (iid) SEs. Cluster # SEs are also printed and are numerically similar (each cluster # contains exactly one observation after dedup). # ---------------------------------------------------------------------- cc = df.dropna(subset=["iron", "KtV"]) fe = (cc.sort_values("obs_id") .groupby("patient_code", as_index=False).first()) print(f"\nfirst-episode set: n={len(fe)} " f"(case patients={int(fe.y.sum())}, " f"control patients={int((1 - fe.y).sum())})") for name, (xcols, req) in SPECS.items(): d = fe.dropna(subset=[c for c in req if c in fe.columns]) b, cov_m, cov_c = fit_logit(d, xcols, cluster_col="patient_code") for i, c in enumerate(xcols[:2], start=1): rm = or_row(b, cov_m, i) print(f" FE {name:11s} {c:8s} model OR={rm['OR']:.3f} " f"({rm['lo']:.3f}-{rm['hi']:.3f}) P={rm['P']:.4f}") # ---------------------------------------------------------------------- # 5. Multiple imputation (m=20, MICE with a linear main-effects Bayesian # ridge specification via sklearn IterativeImputer, # sample_posterior=True; Rubin's rules) # ---------------------------------------------------------------------- BASE = ["y", "iron", "KtV", "age", "male", "hemoglobin", "albumin", "log_crp"] D = df[BASE].copy() m_imp = 20 for label, xn in SPECS.items(): xcols = xn[0] qs, us = {v: [] for v in ["iron", "KtV"]}, {v: [] for v in ["iron", "KtV"]} for k in range(m_imp): # default estimator = BayesianRidge (chained equations with Bayesian # ridge regression); sample_posterior=True gives proper MI draws imp = IterativeImputer(max_iter=30, sample_posterior=True, random_state=2024 + k, min_value=0) Mi = pd.DataFrame(imp.fit_transform(D), columns=BASE) Mi["iron_sd"] = Mi["iron"] / SD_IRON Mi["ktv_sd"] = Mi["KtV"] / SD_KTV b, cov_model, _ = fit_logit(Mi, xcols, cluster_col=None) for idx, v in enumerate(["iron", "KtV"], start=1): qs[v].append(b[idx]) us[v].append(np.sqrt(cov_model[idx, idx])) for v in ["iron", "KtV"]: qbar = np.mean(qs[v]); ubar = np.mean(us[v]) B = np.var(qs[v], ddof=1) T = ubar + (1 + 1 / m_imp) * B se = np.sqrt(T) pv = 2 * stats.norm.sf(abs(qbar / se)) print(f"MI {label:11s} {v:5s} OR={np.exp(qbar):.2f} " f"({np.exp(qbar-1.96*se):.2f}-{np.exp(qbar+1.96*se):.2f}) " f"P={pv:.3f}") # ---------------------------------------------------------------------- # 5b. Correlation analysis (iron vs CRP) # ---------------------------------------------------------------------- print("\n=== Correlation Analysis (iron vs CRP) ===") # Full case set (iron + CRP both non-missing) cc_corr = df.dropna(subset=["iron", "CRP"]).copy() # log_crp is the CSV column ln(CRP+1), identical transform to the primary regression pearson_full, _ = stats.pearsonr(cc_corr["iron"], cc_corr["log_crp"]) spearman_full, _ = stats.spearmanr(cc_corr["iron"], cc_corr["CRP"]) print(f"Full case set (iron+CRP, n={len(cc_corr)}):") print(f" Pearson(iron vs log-CRP): {pearson_full:.3f}") print(f" Spearman(iron vs CRP): {spearman_full:.3f}") # Fully adjusted complete case set req_cols = ["iron", "KtV", "sex", "age", "hemoglobin", "albumin", "CRP"] cc_adj = df.dropna(subset=req_cols).copy() # log_crp is the CSV column ln(CRP+1), identical transform to the primary regression pearson_adj, _ = stats.pearsonr(cc_adj["iron"], cc_adj["log_crp"]) spearman_adj, _ = stats.spearmanr(cc_adj["iron"], cc_adj["CRP"]) print(f"\nFully adjusted complete case set (n={len(cc_adj)}):") print(f" Pearson(iron vs log-CRP): {pearson_adj:.3f}") print(f" Spearman(iron vs CRP): {spearman_adj:.3f}") # ---------------------------------------------------------------------- # 6. E-values (VanderWeele-Ding), OR-based, from primary cluster-robust CIs # In this case-control study the case proportion reflects sampling design, # not population incidence; OR-based E-values are therefore appropriate. # ---------------------------------------------------------------------- def evalue(or_, lo=None, hi=None): """VanderWeele-Ding E-value for OR (applicable to both risk and protective).""" def E(orv): orv = 1 / orv if orv < 1 else orv return orv + np.sqrt(orv * (orv - 1)) Ep = E(or_) Ec = (E(hi) if (or_ < 1 and hi is not None and hi < 1) else E(lo) if (or_ > 1 and lo is not None and lo > 1) else 1.0) return Ep, Ec # Calculate and output E-values for all models print("\n=== OR-based E-values ===") for name, (xcols, req) in SPECS.items(): d = df.dropna(subset=req) b, cov_model, cov_cl = fit_logit(d, xcols) print(f"\n{name}:") for i, c in enumerate(["iron", "KtV"], start=1): se_cl = np.sqrt(cov_cl[i, i]) or_val = np.exp(b[i]) lo = np.exp(b[i] - 1.96 * se_cl) hi = np.exp(b[i] + 1.96 * se_cl) ep, ec = evalue(or_val, lo, hi) print(f" {c}: OR={or_val:.4f} ({lo:.4f}-{hi:.4f}), E(point)={ep:.2f}, E(CI limit)={ec:.2f}") # ---------------------------------------------------------------------- # 7. Episode distribution (Supplementary Table S4), reproducible from S5 # ---------------------------------------------------------------------- _ep = df.groupby(["group", "patient_code"]).size() _case_ep = _ep.loc["case"]; _ctrl_ep = _ep.loc["ctrl"] print("\nepisode distribution (Supplementary Table S4):") print(f" PDAP episodes: {int((df.group=='case').sum())} from {_case_ep.size} " f"unique patients; {int((_case_ep>1).sum())} patients " f"({(_case_ep>1).mean()*100:.1f}%) contributed multiple episodes") for _k in sorted(int(x) for x in _case_ep.unique()): print(f" {_k} episode(s): {int((_case_ep==_k).sum())} patients") print(f" control observations: {int((df.group=='ctrl').sum())} from " f"{_ctrl_ep.size} unique patients ({int((_ctrl_ep>1).sum())} with two observations)") # ---------------------------------------------------------------------- # 8. Sensitivity: Robust Poisson regression (RR estimates) # ---------------------------------------------------------------------- print("\n=== Robust Poisson Regression (RR) ===") def robust_poisson(d, xcols, var_names): X = sm_glm.add_constant(d[xcols]) y = d["y"] model = sm_glm.GLM(y, X, family=sm_glm.families.Poisson()) result = model.fit(cov_type="HC0") for var, vn in zip(xcols, var_names): rr = np.exp(result.params[var]) ci = result.conf_int() lo = np.exp(ci.loc[var].iloc[0]) hi = np.exp(ci.loc[var].iloc[1]) p = result.pvalues[var] print(f" {vn}: RR={rr:.3f} ({lo:.3f}-{hi:.3f}) P={p:.4f}") # Unadjusted d = df.dropna(subset=["iron", "KtV"]).copy() print("Unadjusted:") robust_poisson(d, ["iron_sd", "ktv_sd"], ["iron", "KtV"]) # Adjusted d = df.dropna(subset=["iron", "KtV", "sex", "age", "hemoglobin", "albumin"]).copy() print("\nAdjusted (sex, age, Hb, Alb):") robust_poisson(d, ["iron_sd", "ktv_sd", "male", "age", "hemoglobin", "albumin"], ["iron", "KtV", "sex", "age", "hemoglobin", "albumin"]) # Fully adjusted d = df.dropna(subset=["iron", "KtV", "sex", "age", "hemoglobin", "albumin", "CRP"]).copy() # log_crp is the CSV column ln(CRP+1), identical transform to the primary regression print("\nFully adjusted (+log-CRP):") robust_poisson(d, ["iron_sd", "ktv_sd", "male", "age", "hemoglobin", "albumin", "log_crp"], ["iron", "KtV", "sex", "age", "hemoglobin", "albumin", "log-CRP"]) # ---------------------------------------------------------------------- # 8b. Tipping-point analysis: serum iron (missing in 30 cases, 7 controls) # ---------------------------------------------------------------------- print("\n=== Iron Tipping-Point Analysis ===") ctrl_iron = df.loc[df["group"] == "ctrl", "iron"].dropna() case_iron_missing = df.loc[(df["group"] == "case") & df["iron"].isna()] print(f"Iron missing: {len(case_iron_missing)} cases ({len(case_iron_missing)/(df.group=='case').sum()*100:.1f}%), " f"{(df.group=='ctrl').sum() - len(ctrl_iron)} controls") print(f"Control iron: min={ctrl_iron.min():.3f}, P5={np.percentile(ctrl_iron,5):.3f}, " f"mean={ctrl_iron.mean():.3f}, max={ctrl_iron.max():.3f}") # NOTE: P values printed below are cluster-robust (or_row on cov_cl), # consistent with the manuscript's primary SE specification. def iron_tipping(fill_value, label): d = df.copy() d.loc[(d.group == "case") & d.iron.isna(), "iron"] = fill_value d.loc[(d.group == "ctrl") & d.iron.isna(), "iron"] = fill_value d = d.dropna(subset=["iron"]) d["iron_sd"] = d["iron"] / SD_IRON b, cov_model, cov_cl = fit_logit(d, ["iron_sd"]) r = or_row(b, cov_cl, 1) print(f" {label}: fill={fill_value:.2f} -> OR={r['OR']:.3f} " f"({r['lo']:.3f}-{r['hi']:.3f}) P={r['P']:.4f}") return r iron_tipping(ctrl_iron.min(), "Control min") iron_tipping(np.percentile(ctrl_iron, 5), "Control P5") iron_tipping(np.percentile(ctrl_iron, 25), "Control P25") iron_tipping(ctrl_iron.mean(), "Control mean") iron_tipping(ctrl_iron.max(), "Control max") # Binary search for nullification point (CI lower bound crosses OR=1) from scipy.optimize import brentq def null_p_diff(fill_val): d = df.copy() d.loc[(d.group == "case") & d.iron.isna(), "iron"] = fill_val d.loc[(d.group == "ctrl") & d.iron.isna(), "iron"] = fill_val d = d.dropna(subset=["iron"]) d["iron_sd"] = d["iron"] / SD_IRON b, cov_model, _ = fit_logit(d, ["iron_sd"], cluster_col=None) se = np.sqrt(cov_model[1, 1]) return np.exp(b[1] - 1.96 * se) - 1.0 try: null_val = brentq(null_p_diff, 0.1, ctrl_iron.max() * 2, xtol=0.001) pct = (ctrl_iron < null_val).mean() * 100 print(f"\n Nullification point: fill = {null_val:.3f} umol/L " f"({pct:.1f}th percentile of control distribution)") except: print("\n Nullification point: above control max (association robust to all plausible fills)") # ---------------------------------------------------------------------- # 8c. Tipping-point analysis: total Kt/V (missing in 16 cases, 1 control) # ---------------------------------------------------------------------- print("\n=== Kt/V Tipping-Point Analysis ===") ctrl_ktv = df.loc[df["group"] == "ctrl", "KtV"].dropna() case_ktv_missing = df.loc[(df["group"] == "case") & df["KtV"].isna()] print(f"Kt/V missing: {len(case_ktv_missing)} cases ({len(case_ktv_missing)/(df.group=='case').sum()*100:.1f}%), " f"{(df.group=='ctrl').sum() - len(ctrl_ktv)} controls") print(f"Control Kt/V: min={ctrl_ktv.min():.3f}, P5={np.percentile(ctrl_ktv,5):.3f}, " f"mean={ctrl_ktv.mean():.3f}, max={ctrl_ktv.max():.3f}") def ktv_tipping(fill_value, label): d = df.copy() d.loc[(d.group == "case") & d.KtV.isna(), "KtV"] = fill_value d.loc[(d.group == "ctrl") & d.KtV.isna(), "KtV"] = fill_value d = d.dropna(subset=["KtV"]) d["ktv_sd"] = d["KtV"] / SD_KTV # Unadjusted logistic regression (case status on standardised Kt/V) b, cov_model, cov_cl = fit_logit(d, ["ktv_sd"]) r = or_row(b, cov_cl, 1) print(f" {label}: fill={fill_value:.2f} -> OR={r['OR']:.3f} " f"({r['lo']:.3f}-{r['hi']:.3f}) P={r['P']:.4f}") return r ktv_tipping(ctrl_ktv.min(), "Control min") ktv_tipping(np.percentile(ctrl_ktv, 5), "Control P5") ktv_tipping(np.percentile(ctrl_ktv, 25), "Control P25") ktv_tipping(ctrl_ktv.mean(), "Control mean") ktv_tipping(ctrl_ktv.max(), "Control max") # Binary search for point-estimate nullification (OR crosses 1) from scipy.optimize import brentq def ktv_null_point_est(fill_val): d = df.copy() d.loc[(d.group == "case") & d.KtV.isna(), "KtV"] = fill_val d.loc[(d.group == "ctrl") & d.KtV.isna(), "KtV"] = fill_val d = d.dropna(subset=["KtV"]) d["ktv_sd"] = d["KtV"] / SD_KTV b, cov_model, _ = fit_logit(d, ["ktv_sd"], cluster_col=None) se = np.sqrt(cov_model[1, 1]) return np.exp(b[1]) - 1.0 try: ktv_null_val = brentq(ktv_null_point_est, ctrl_ktv.min(), ctrl_ktv.max() * 2, xtol=0.001) ktv_pct = (ctrl_ktv < ktv_null_val).mean() * 100 print(f"\n Point-estimate nullification: fill = {ktv_null_val:.3f} " f"({ktv_pct:.1f}th percentile of control distribution)") except: print("\n Point-estimate nullification: above control max") # ---------------------------------------------------------------------- # 8d. Case characteristics by total Kt/V measurement availability # (Supplementary Table S3b). Cases with versus without a recent # total Kt/V measurement; Mann-Whitney U / Fisher's exact. # ---------------------------------------------------------------------- print("\n=== Case characteristics by Kt/V measurement availability (Table S3b) ===") cases = df[df.group == "case"].copy() ep_counts = cases.groupby("patient_code")["obs_id"].transform("size") cases["recurrent"] = (ep_counts > 1).astype(int) has_ktv = cases[cases.KtV.notna()] no_ktv = cases[cases.KtV.isna()] print(f"with recent Kt/V: n={len(has_ktv)}; without: n={len(no_ktv)}") for _v in ["age", "hemoglobin", "WBC", "CRP", "PCT", "IL6", "ferritin", "iron", "TIBC", "TSAT", "albumin", "creatinine", "pro_BNP"]: _a = has_ktv[_v].dropna(); _bq = no_ktv[_v].dropna() if len(_a) < 2 or len(_bq) < 2: continue _, _p = stats.mannwhitneyu(_a, _bq, alternative="two-sided") print(f" {_v:12s} with n={len(_a)} median={_a.median():.2f} | " f"without n={len(_bq)} median={_bq.median():.2f} | MWU P={_p:.4f}") _tab = np.array([[int(has_ktv.recurrent.sum()), int((1 - has_ktv.recurrent).sum())], [int(no_ktv.recurrent.sum()), int((1 - no_ktv.recurrent).sum())]]) print(f" recurrent: with {_tab[0,0]}/{len(has_ktv)} ({has_ktv.recurrent.mean()*100:.1f}%) | " f"without {_tab[1,0]}/{len(no_ktv)} ({no_ktv.recurrent.mean()*100:.1f}%) | " f"Fisher P={stats.fisher_exact(_tab)[1]:.4f}") _mh = has_ktv["male"].dropna(); _mn = no_ktv["male"].dropna() _tab2 = np.array([[int(_mh.sum()), int((1 - _mh).sum())], [int(_mn.sum()), int((1 - _mn).sum())]]) print(f" male: with {int(_mh.sum())}/{len(_mh)} | without {int(_mn.sum())}/{len(_mn)} | " f"Fisher P={stats.fisher_exact(_tab2)[1]:.4f}") # -> lower haemoglobin (99.3 vs 115.2 g/L, P=0.007) and TIBC (31.4 vs 39.9 # umol/L, P=0.016) among cases without a recent Kt/V measurement, with # numerically higher IL-6; other characteristics similar (Table S3b). # ---------------------------------------------------------------------- # 8e. Delta-adjustment (pattern-mixture style) sensitivity for total Kt/V. # Because S3b suggests missingness is informative (MNAR rather than # MAR), missing Kt/V values in CASE episodes are deterministically # filled at (observed case median + delta) and the single missing # CONTROL value at (observed control median - delta), shifting the # imputed values in the direction that WEAKENS the reported # association (i.e., assuming the missing cases were not truly # Kt/V-depleted). delta = 0 / 0.2 / 0.4 / 0.6 Kt/V units. The fully # adjusted specification is then refit on each completed dataset # with the same cluster-robust SEs. This brackets how large the # departure from MAR must be to erase the association. # ---------------------------------------------------------------------- print("\n=== Delta-Adjustment Sensitivity for Total Kt/V (fully adjusted) ===") print(" (cluster-robust SE without finite-sample correction, matching the primary analysis)") _fully_req = ["iron", "KtV", "sex", "age", "hemoglobin", "albumin", "log_crp"] _ktv_sd_div = 0.672 # pooled SD used for standardisation (ktv_sd = KtV / 0.672) _case_med = df.loc[(df.is_case == 1) & df.KtV.notna(), "KtV"].median() _ctrl_med = df.loc[(df.is_case == 0) & df.KtV.notna(), "KtV"].median() _mi_case = df.index[(df.is_case == 1) & df.KtV.isna()] _mi_ctrl = df.index[(df.is_case == 0) & df.KtV.isna()] print(f"case Kt/V missing n={len(_mi_case)} (observed case median {_case_med:.2f}); " f"control Kt/V missing n={len(_mi_ctrl)} (observed control median {_ctrl_med:.2f})") for _delta in [0.0, 0.2, 0.4, 0.6, 0.7, 0.8, 0.9]: _d = df.copy() _d.loc[_mi_case, "KtV"] = _case_med + _delta _d.loc[_mi_ctrl, "KtV"] = _ctrl_med - _delta _d["ktv_sd"] = _d["KtV"] / _ktv_sd_div _dd = _d.dropna(subset=_fully_req) _b, _cm, _cc = fit_logit(_dd, ["iron_sd", "ktv_sd", "male", "age", "hemoglobin", "albumin", "log_crp"]) _r = or_row(_b, _cc, 2) # ktv_sd is the 2nd covariate print(f" delta={_delta:.2f}: n={len(_dd)} Kt/V cluster OR={_r['OR']:.3f} " f"({_r['lo']:.3f}-{_r['hi']:.3f}) P={_r['P']:.4f}") # ---------------------------------------------------------------------- # 9. Sensitivity: Cluster bootstrap by patient # ---------------------------------------------------------------------- print("\n=== Cluster Bootstrap by Patient (500 iterations) ===") np.random.seed(2024) n_boot = 500 d = df.dropna(subset=["iron", "KtV", "sex", "age", "hemoglobin", "albumin", "CRP"]).copy() # log_crp is the CSV column ln(CRP+1), identical transform to the primary regression clusters = d["patient_code"].unique() boot_iron_rr = [] boot_ktv_rr = [] for i in range(n_boot): boot_clusters = np.random.choice(clusters, size=len(clusters), replace=True) boot_data = [] for c in boot_clusters: boot_data.append(d[d["patient_code"] == c]) boot_df = pd.concat(boot_data) X = sm_glm.add_constant(boot_df[["iron_sd", "ktv_sd", "male", "age", "hemoglobin", "albumin", "log_crp"]]) y = boot_df["y"] try: model = sm_glm.GLM(y, X, family=sm_glm.families.Poisson()) result = model.fit(cov_type="HC0", maxiter=100) boot_iron_rr.append(np.exp(result.params["iron_sd"])) boot_ktv_rr.append(np.exp(result.params["ktv_sd"])) except Exception as e: print(f"Bootstrap iteration failed: {e}") continue boot_iron_rr = np.array(boot_iron_rr) boot_ktv_rr = np.array(boot_ktv_rr) print(f"Iron: mean RR={np.mean(boot_iron_rr):.3f}, 95% CI=({np.percentile(boot_iron_rr, 2.5):.3f}-{np.percentile(boot_iron_rr, 97.5):.3f})") print(f"KtV: mean RR={np.mean(boot_ktv_rr):.3f}, 95% CI=({np.percentile(boot_ktv_rr, 2.5):.3f}-{np.percentile(boot_ktv_rr, 97.5):.3f})")