"""Phase 5 - Predictive modeling of normalized billing.
This module models the row-level target ``log1p(NBA)`` for provider-service
observations. NBA is Agrawal-normalized billing:
NBA_h = avg_submitted_charge_h * US_mean_payment
/ (US_mean_billing * provider_payment_h)
The target is ``log1p(max(NBA, 0))``. It is not a residual target and it is
not a state-level NB statistic. The charge-to-payment ratio is intentionally
excluded because it is mechanically adjacent to NBA.
Both models are evaluated with provider-grouped held-out validation. Fold
assignments, out-of-fold predictions, and their metadata are written as audit
artifacts under DATA_PROC. Full-data fits are used only for descriptive
coefficients and random-forest MDI; they are not used to claim held-out
performance.
"""
import json
from collections import OrderedDict
import numpy as np
import pandas as pd
from sklearn.base import clone
from sklearn.compose import ColumnTransformer
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score
from sklearn.model_selection import GroupKFold
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import OneHotEncoder, StandardScaler
from .config import DATA_PROC, DEFAULT_SEED, MIN_PROVIDERS_PER_SERVICE
from .analysis_insights import save_json
# ---------------------------------------------------------------------------
# Modeling contract
# ---------------------------------------------------------------------------
TARGET_COLUMN = "log_NBA"
GROUP_COLUMN = "provider_ccn"
TARGET_FORMULA = (
"NBA_h = avg_submitted_charge_h * US_mean_payment / "
"(US_mean_billing * provider_payment_h); "
"target_h = log1p(max(NBA_h, 0))"
)
NORMALIZATION_REFERENCE = (
"Fixed service-level means computed from the complete FY2023 inpatient or "
"CY2023 outpatient analysis frame; this is a cross-sectional benchmark, "
"not a future-period deployment simulation."
)
VALIDATION_ESTIMAND = (
"Mean fold-level row-wise R2 for log1p(NBA) under provider-held-out "
"GroupKFold; each provider_ccn is assigned to one test fold and never "
"appears in that fold's training data."
)
VALIDATION_N_SPLITS = 5
PERMUTATION_REPEATS = 2
PORTABLE_AUDIT_PREFIX = "site/data/audit"
# census_region is intentionally absent. State is retained because it is the
# finest geography available in this modeling frame and is explicitly coded.
CATEGORICAL_FEATS = ["urban_rural", "service_family", "state"]
NUMERIC_FEATS = ["log_volume"]
MODEL_FEATURES = CATEGORICAL_FEATS + NUMERIC_FEATS
PROVIDER_LEVEL_FEATURES = {"state", "urban_rural"}
# These are semantic preferred references. If a preferred level is absent
# from a dataset, the first observed level in deterministic lexical order is
# used and recorded in predictors.json.
REFERENCE_LEVEL_PREFERENCES = {
"urban_rural": ("Urban",),
"service_family": ("Other", "Other / Unclassified", "Unknown"),
"state": ("AL",),
}
# ---------------------------------------------------------------------------
# Build modeling frame
# ---------------------------------------------------------------------------
def build_modeling_frame(df: pd.DataFrame, service_col: str,
billing_col: str, payment_col: str,
n_volume_col: str, family_col: str,
dataset_name: str) -> pd.DataFrame:
"""Build provider-service rows with target ``log1p(NBA)``.
NBA is computed from service-level US means before the service-size
threshold is applied. The returned frame keeps the excluded
``charge_to_payment_ratio`` only for audit visibility; it is never part of
``MODEL_FEATURES``. Cleaning and clipping counts are attached to
``result.attrs['predictors_metadata']`` for the orchestration layer.
"""
required = [
"provider_ccn", "state", "urban_rural", service_col, family_col,
billing_col, payment_col, n_volume_col,
]
missing_columns = [column for column in required if column not in df.columns]
if missing_columns:
raise KeyError(f"missing predictor input columns: {missing_columns}")
input_rows = int(len(df))
work = df.copy()
work["_billing_value"] = pd.to_numeric(work[billing_col], errors="coerce")
work["_payment_value"] = pd.to_numeric(work[payment_col], errors="coerce")
work["_volume_value"] = pd.to_numeric(work[n_volume_col], errors="coerce")
missing_service = work[service_col].isna()
missing_provider = work[GROUP_COLUMN].isna()
rows_excluded_missing_service = int(missing_service.sum())
rows_excluded_missing_provider = int(missing_provider.sum())
work = work.loc[~missing_service & ~missing_provider].copy()
mapping_quarantined = work.get(
"mapping_quarantined", pd.Series(False, index=work.index)
).fillna(False).astype(bool)
rows_excluded_mapping_quarantine = int(mapping_quarantined.sum())
work = work.loc[~mapping_quarantined].copy()
# Compute service-level means and the provider count on the eligible input
# rows. The threshold is based on unique providers, not merely row count.
grouped = work.groupby(service_col, dropna=False, sort=False)
work["n_providers_per_service"] = grouped[GROUP_COLUMN].transform("nunique")
work["_us_p_mean"] = grouped["_payment_value"].transform("mean")
work["_us_b_mean"] = grouped["_billing_value"].transform("mean")
eligible_service = work["n_providers_per_service"] >= MIN_PROVIDERS_PER_SERVICE
rows_excluded_service_threshold = int((~eligible_service).sum())
services_below_threshold = int(
work.loc[~eligible_service, service_col].nunique(dropna=False)
)
work = work.loc[eligible_service].copy()
payment_valid = work["_payment_value"] > 0
rows_excluded_payment = int((~payment_valid).sum())
work = work.loc[payment_valid].copy()
with np.errstate(divide="ignore", invalid="ignore", over="ignore"):
work["NBA"] = (
work["_billing_value"] * work["_us_p_mean"]
/ (work["_us_b_mean"] * work["_payment_value"])
)
work["NBA"] = work["NBA"].replace([np.inf, -np.inf], np.nan)
finite_nba = work["NBA"].notna()
rows_excluded_nonfinite_nba = int((~finite_nba).sum())
work = work.loc[finite_nba].copy()
negative_nba = work["NBA"] < 0
rows_clipped_nba = int(negative_nba.sum())
work[TARGET_COLUMN] = np.log1p(work["NBA"].clip(lower=0))
finite_volume = np.isfinite(work["_volume_value"])
rows_excluded_nonfinite_volume = int((~finite_volume).sum())
work = work.loc[finite_volume].copy()
rows_clipped_volume = int((work["_volume_value"] < 0).sum())
work["log_volume"] = np.log1p(work["_volume_value"].clip(lower=0))
# Preserve missing categorical geography/classification as an explicit
# level rather than silently dropping rows. The replacement is recorded.
missing_category_counts = OrderedDict()
for column in ["state", "urban_rural", family_col]:
missing_category_counts[column] = int(work[column].isna().sum())
work[column] = work[column].where(work[column].notna(), "Unknown").astype(str)
# The ratio is retained for audit/debugging but explicitly excluded from X.
if "charge_to_payment_ratio" not in work.columns:
work["charge_to_payment_ratio"] = np.nan
result = work[
[
"provider_ccn", "state", "urban_rural", service_col, family_col,
"log_volume", "charge_to_payment_ratio", "NBA", TARGET_COLUMN,
"n_providers_per_service",
]
].rename(columns={
service_col: "service_code",
family_col: "service_family",
})
result["dataset"] = dataset_name
result = result.reset_index(drop=True)
metadata = OrderedDict(
dataset=dataset_name,
target_column=TARGET_COLUMN,
target_label="log1p(NBA)",
target_formula=TARGET_FORMULA,
service_provider_threshold=OrderedDict(
minimum_unique_providers=int(MIN_PROVIDERS_PER_SERVICE),
group_column=GROUP_COLUMN,
),
exclusion_counts=OrderedDict(
input_rows=input_rows,
missing_service_code=rows_excluded_missing_service,
missing_provider_group=rows_excluded_missing_provider,
mapping_quarantined=rows_excluded_mapping_quarantine,
below_minimum_providers_per_service=rows_excluded_service_threshold,
nonpositive_or_missing_payment=rows_excluded_payment,
nonfinite_NBA=rows_excluded_nonfinite_nba,
nonfinite_volume=rows_excluded_nonfinite_volume,
modeled_rows=int(len(result)),
),
clipping_counts=OrderedDict(
NBA_lower_bound_zero=rows_clipped_nba,
volume_lower_bound_zero=rows_clipped_volume,
),
services_below_threshold=services_below_threshold,
missing_categorical_values_replaced_with_Unknown=missing_category_counts,
excluded_features=["charge_to_payment_ratio", "census_region"],
excluded_geography=OrderedDict(
feature="census_region",
retained_geography=["state"],
reason="State is retained; census_region is nested geography and removed to avoid redundant coding.",
),
payment_column=payment_col,
)
result.attrs["predictors_metadata"] = metadata
return result
# ---------------------------------------------------------------------------
# Preprocessing and reference coding
# ---------------------------------------------------------------------------
def _prepare_model_inputs(df: pd.DataFrame) -> pd.DataFrame:
"""Return model columns with categorical missing values made explicit."""
missing = [column for column in MODEL_FEATURES if column not in df.columns]
if missing:
raise KeyError(f"missing model feature columns: {missing}")
X = df[MODEL_FEATURES].copy()
for column in CATEGORICAL_FEATS:
X[column] = X[column].where(X[column].notna(), "Unknown").astype(str)
X["log_volume"] = pd.to_numeric(X["log_volume"], errors="coerce")
if not np.isfinite(X["log_volume"]).all():
raise ValueError("log_volume contains non-finite values")
return X
def _reference_spec(X: pd.DataFrame):
"""Return deterministic category levels and explicit reference levels."""
category_levels = OrderedDict()
reference_levels = OrderedDict()
for feature in CATEGORICAL_FEATS:
levels = sorted(X[feature].astype(str).unique().tolist())
if not levels:
raise ValueError(f"no observed levels for categorical feature {feature}")
category_levels[feature] = levels
preferred = REFERENCE_LEVEL_PREFERENCES.get(feature, ())
reference_levels[feature] = next(
(level for level in preferred if level in levels), levels[0]
)
return category_levels, reference_levels
def _pipeline(model, category_levels=None, reference_levels=None):
"""Build a pipeline with one explicit reference category per categorical.
``min_frequency`` is deliberately not supplied: no rare-category grouping
is performed. The main validation path supplies observed category levels
and explicit drops; the fallback keeps the helper usable in isolation.
"""
if category_levels is None:
encoder = OneHotEncoder(handle_unknown="ignore", drop="first")
else:
if reference_levels is None:
raise ValueError("reference_levels required with category_levels")
encoder = OneHotEncoder(
categories=[category_levels[feature] for feature in CATEGORICAL_FEATS],
drop=[reference_levels[feature] for feature in CATEGORICAL_FEATS],
handle_unknown="ignore",
)
pre = ColumnTransformer(
transformers=[
("cat", encoder, CATEGORICAL_FEATS),
("num", StandardScaler(), NUMERIC_FEATS),
],
remainder="drop",
)
return Pipeline([("pre", pre), ("model", model)])
# ---------------------------------------------------------------------------
# Grouped validation and held-out importance
# ---------------------------------------------------------------------------
def _grouped_splits(df: pd.DataFrame, n_splits: int = VALIDATION_N_SPLITS):
"""Create and validate provider-disjoint GroupKFold splits."""
if GROUP_COLUMN not in df.columns:
raise KeyError(f"missing validation group column: {GROUP_COLUMN}")
groups = df[GROUP_COLUMN]
if groups.isna().any():
raise ValueError("provider groups must be non-missing for GroupKFold")
n_groups = int(groups.nunique())
if n_groups < n_splits:
raise ValueError(
f"GroupKFold requires at least {n_splits} provider groups; found {n_groups}"
)
splitter = GroupKFold(n_splits=n_splits)
X = _prepare_model_inputs(df)
y = df[TARGET_COLUMN]
splits = list(splitter.split(X, y, groups))
for train_idx, test_idx in splits:
train_groups = set(groups.iloc[train_idx])
test_groups = set(groups.iloc[test_idx])
if train_groups.intersection(test_groups):
raise RuntimeError("provider groups overlap between train and test folds")
return splits
def _safe_r2(y_true, prediction):
if len(y_true) < 2:
return None
score = r2_score(y_true, prediction, force_finite=True)
return float(score) if np.isfinite(score) else None
def _aggregate_permutation_importance(fold_importances, raw_features):
"""Aggregate fold-level held-out permutation means for site compatibility."""
if not fold_importances:
return []
values = np.asarray([
[row[feature] for feature in raw_features]
for row in fold_importances
], dtype=float)
aggregate = []
for position, feature in enumerate(raw_features):
fold_values = values[:, position]
aggregate.append(OrderedDict(
feature=feature,
importance=float(np.mean(fold_values)),
std=float(np.std(fold_values, ddof=1)) if len(fold_values) > 1 else 0.0,
))
aggregate.sort(key=lambda row: row["importance"], reverse=True)
return aggregate
def _provider_group_permutation_importance(
model, X_test: pd.DataFrame, y_test: pd.Series, groups_test: pd.Series,
raw_features: list, n_repeats: int, random_state: int,
) -> dict:
"""Permute features without moving observations between providers.
Provider-level geography is shuffled as a whole provider block. Service
family and log volume are shuffled within each provider, preserving the
held-out provider composition while still perturbing row-level signals.
"""
X_test = X_test.reset_index(drop=True).copy()
y_test = pd.Series(y_test).reset_index(drop=True)
groups_test = pd.Series(groups_test).reset_index(drop=True)
if len(X_test) != len(y_test) or len(X_test) != len(groups_test):
raise ValueError("held-out feature, target, and group lengths differ")
baseline = float(model.score(X_test, y_test))
if not np.isfinite(baseline):
baseline = 0.0
rng = np.random.default_rng(random_state)
group_values = groups_test.drop_duplicates().to_numpy()
group_array = groups_test.to_numpy()
importances = {feature: [] for feature in raw_features}
for _ in range(n_repeats):
for feature in raw_features:
permuted = X_test.copy()
column_position = permuted.columns.get_loc(feature)
if feature in PROVIDER_LEVEL_FEATURES:
donor_groups = rng.permutation(group_values)
donor_values = {}
for donor in group_values:
donor_positions = np.flatnonzero(group_array == donor)
donor_values[donor] = X_test.iloc[donor_positions][feature].iloc[0]
for target, donor in zip(group_values, donor_groups):
target_positions = np.flatnonzero(group_array == target)
permuted.iloc[target_positions, column_position] = donor_values[donor]
else:
for group in group_values:
positions = np.flatnonzero(group_array == group)
values = X_test.iloc[positions][feature].to_numpy(copy=True)
rng.shuffle(values)
permuted.iloc[positions, column_position] = values
score = float(model.score(permuted, y_test))
if not np.isfinite(score):
score = 0.0
importances[feature].append(baseline - score)
return {
"importances_mean": np.asarray(
[np.mean(importances[feature]) for feature in raw_features], dtype=float
),
"importances_std": np.asarray(
[
np.std(importances[feature], ddof=1) if n_repeats > 1 else 0.0
for feature in raw_features
],
dtype=float,
),
"baseline_score": baseline,
}
def cv_score(df: pd.DataFrame, model, n_splits: int = VALIDATION_N_SPLITS,
splits=None, category_levels=None, reference_levels=None,
n_repeats: int = PERMUTATION_REPEATS) -> dict:
"""Run provider-grouped held-out CV and fold-by-fold raw-feature importance.
The returned full-data pipeline is explicitly descriptive. R2 and
permutation importance are calculated only from fold-specific pipelines on
their held-out rows.
"""
X = _prepare_model_inputs(df)
y = pd.to_numeric(df[TARGET_COLUMN], errors="coerce")
if not np.isfinite(y).all():
raise ValueError(f"{TARGET_COLUMN} contains non-finite values")
if category_levels is None or reference_levels is None:
category_levels, reference_levels = _reference_spec(X)
if splits is None:
splits = _grouped_splits(df, n_splits=n_splits)
if len(splits) != n_splits:
raise ValueError("provided splits do not match n_splits")
if n_repeats < 1:
raise ValueError("n_repeats must be at least 1")
groups = df[GROUP_COLUMN]
oof_predictions = np.full(len(df), np.nan, dtype=float)
scores = []
fold_summaries = []
fold_importance_values = []
fold_importance_details = []
for fold_number, (train_idx, test_idx) in enumerate(splits, start=1):
train_groups = set(groups.iloc[train_idx])
test_groups = set(groups.iloc[test_idx])
overlap = train_groups.intersection(test_groups)
if overlap:
raise RuntimeError("provider groups overlap between train and test folds")
fold_pipe = _pipeline(
clone(model), category_levels=category_levels,
reference_levels=reference_levels,
)
fold_pipe.fit(X.iloc[train_idx], y.iloc[train_idx])
predictions = fold_pipe.predict(X.iloc[test_idx])
oof_predictions[test_idx] = predictions
fold_r2 = _safe_r2(y.iloc[test_idx], predictions)
if fold_r2 is not None:
scores.append(fold_r2)
# Evaluate on the held-out partition and preserve provider groups
# during each feature perturbation.
perm = _provider_group_permutation_importance(
fold_pipe,
X.iloc[test_idx],
y.iloc[test_idx],
groups.iloc[test_idx],
MODEL_FEATURES,
n_repeats,
DEFAULT_SEED + fold_number,
)
fold_means = {
feature: float(importance)
for feature, importance in zip(
MODEL_FEATURES, perm["importances_mean"]
)
}
fold_importance_values.append(fold_means)
fold_importance_details.append(OrderedDict(
fold=fold_number,
n_test_rows=int(len(test_idx)),
test_r2=fold_r2,
importance=[
OrderedDict(
feature=feature,
importance=float(perm["importances_mean"][position]),
std=float(perm["importances_std"][position]),
)
for position, feature in enumerate(MODEL_FEATURES)
],
))
fold_summaries.append(OrderedDict(
fold=fold_number,
n_train_rows=int(len(train_idx)),
n_test_rows=int(len(test_idx)),
n_train_groups=int(len(train_groups)),
n_test_groups=int(len(test_groups)),
group_overlap_count=int(len(overlap)),
r2=fold_r2,
))
if np.isnan(oof_predictions).any():
raise RuntimeError("GroupKFold did not produce a complete OOF prediction")
full_pipe = _pipeline(
clone(model), category_levels=category_levels,
reference_levels=reference_levels,
)
full_pipe.fit(X, y)
mean_r2 = float(np.mean(scores)) if scores else None
oof_r2 = _safe_r2(y, oof_predictions)
importance_evaluation = OrderedDict(
method="held-out provider-group-preserving permutation importance",
scoring="r2",
evaluated_on="test partition within each provider-grouped fold",
fitted_on="training partition within the same fold",
aggregation="mean of fold-level importances; std is across fold means",
full_data_permutation=False,
n_repeats=int(n_repeats),
features=MODEL_FEATURES,
provider_level_features=sorted(PROVIDER_LEVEL_FEATURES),
within_provider_features=[
feature for feature in MODEL_FEATURES
if feature not in PROVIDER_LEVEL_FEATURES
],
)
return {
"r2_mean": mean_r2,
"r2_std": float(np.std(scores, ddof=1)) if len(scores) > 1 else 0.0,
"r2_folds": [float(score) for score in scores],
"r2_oof": oof_r2,
"n_samples": int(len(df)),
"n_groups": int(groups.nunique()),
"splitter": "GroupKFold",
"n_splits": int(n_splits),
"group_column": GROUP_COLUMN,
"fold_summaries": fold_summaries,
"permutation_importance": _aggregate_permutation_importance(
fold_importance_values, MODEL_FEATURES
),
"permutation_importance_by_fold": fold_importance_details,
"importance_evaluation": importance_evaluation,
"oof_predictions": oof_predictions,
"fitted_pipeline": full_pipe,
"full_fit_scope": "descriptive_full_sample",
"category_levels": category_levels,
"reference_levels": reference_levels,
}
# ---------------------------------------------------------------------------
# Descriptive full-fit summaries
# ---------------------------------------------------------------------------
def linear_coefficients(df: pd.DataFrame, fitted_pipe,
reference_levels: dict = None) -> list:
"""Extract dummy coefficients with explicit reference metadata."""
del df # retained in the signature for compatibility with earlier callers
pre = fitted_pipe.named_steps["pre"]
model = fitted_pipe.named_steps["model"]
feature_names_out = pre.get_feature_names_out()
coefs = np.asarray(model.coef_).ravel()
out = []
for fname, coef in zip(feature_names_out, coefs):
token = fname.split("__", 1)[-1]
raw_feature = next(
(feature for feature in MODEL_FEATURES
if token == feature or token.startswith(f"{feature}_")),
token,
)
coefficient = float(coef)
reference = None
if reference_levels and raw_feature in reference_levels:
reference = f"{raw_feature}={reference_levels[raw_feature]}"
out.append(OrderedDict(
feature=token,
reference=reference,
coefficient=float(round(coefficient, 4)),
coefficient_log_difference=float(round(coefficient, 4)),
multiplicative_effect=float(round(np.exp(coefficient), 4)),
percent_effect=float(round((np.exp(coefficient) - 1) * 100, 2)),
))
out.sort(key=lambda row: abs(row["coefficient"]), reverse=True)
return out
def _mdi_aggregated(fitted_pipe) -> list:
pre = fitted_pipe.named_steps["pre"]
model = fitted_pipe.named_steps["model"]
raw_names = pre.get_feature_names_out()
mdi = pd.Series(model.feature_importances_, index=raw_names)
mdi_agg = {feature: 0.0 for feature in MODEL_FEATURES}
for encoded_name, value in mdi.items():
token = encoded_name.split("__", 1)[-1]
for feature in MODEL_FEATURES:
if token == feature or token.startswith(f"{feature}_"):
mdi_agg[feature] += float(value)
break
return [
OrderedDict(feature=feature, importance=float(round(value, 4)))
for feature, value in sorted(
mdi_agg.items(), key=lambda item: item[1], reverse=True
)
]
def rf_importance(df: pd.DataFrame, fitted_pipe, n_repeats: int = 10,
heldout_permutation=None,
heldout_permutation_by_fold=None,
importance_evaluation=None) -> dict:
"""Combine descriptive full-fit MDI with supplied held-out importance.
``heldout_permutation`` must come from ``cv_score``. This function no
longer computes permutation importance on the full fitted data.
``n_repeats`` remains accepted for caller compatibility but is not used to
perform another evaluation here.
"""
del df, n_repeats
return {
"permutation_importance": heldout_permutation or [],
"permutation_importance_by_fold": heldout_permutation_by_fold or [],
"mdi_aggregated": _mdi_aggregated(fitted_pipe),
"mdi_fit_scope": "descriptive_full_sample",
"importance_evaluation": importance_evaluation or OrderedDict(
method="held-out provider-group-preserving permutation importance",
full_data_permutation=False,
),
}
# ---------------------------------------------------------------------------
# Audit artifacts
# ---------------------------------------------------------------------------
def _combined_fold_summaries(lr_res, rf_res):
summaries = []
for lr_fold, rf_fold in zip(
lr_res["fold_summaries"], rf_res["fold_summaries"]
):
summaries.append(OrderedDict(
fold=lr_fold["fold"],
n_train_rows=lr_fold["n_train_rows"],
n_test_rows=lr_fold["n_test_rows"],
n_train_groups=lr_fold["n_train_groups"],
n_test_groups=lr_fold["n_test_groups"],
group_overlap_count=lr_fold["group_overlap_count"],
models=OrderedDict(
linear_regression=OrderedDict(r2=lr_fold["r2"]),
random_forest=OrderedDict(r2=rf_fold["r2"]),
),
))
return summaries
def _save_validation_artifacts(dataset_name: str, frame: pd.DataFrame,
splits, lr_res: dict, rf_res: dict,
reference_levels: dict, frame_metadata: dict,
fold_summaries: list) -> dict:
"""Write row-level fold assignments and OOF predictions for audit use."""
fold_assignments = np.full(len(frame), -1, dtype=int)
for fold_number, (_, test_idx) in enumerate(splits, start=1):
if (fold_assignments[test_idx] != -1).any():
raise RuntimeError("a row was assigned to multiple audit folds")
fold_assignments[test_idx] = fold_number
if (fold_assignments == -1).any():
raise RuntimeError("some rows have no audit fold assignment")
audit = pd.DataFrame({
"dataset": dataset_name,
"row_id": np.arange(len(frame), dtype=int),
"provider_ccn": frame[GROUP_COLUMN].to_numpy(),
"fold": fold_assignments,
"state": frame["state"].to_numpy(),
"service_code": frame["service_code"].to_numpy(),
"service_family": frame["service_family"].to_numpy(),
"NBA": frame["NBA"].to_numpy(),
"target_log1p_NBA": frame[TARGET_COLUMN].to_numpy(),
"linear_regression_oof_prediction": lr_res["oof_predictions"],
"random_forest_oof_prediction": rf_res["oof_predictions"],
})
DATA_PROC.mkdir(parents=True, exist_ok=True)
audit_path = DATA_PROC / f"predictors_validation_{dataset_name}.parquet"
metadata_path = DATA_PROC / f"predictors_validation_{dataset_name}.json"
audit.to_parquet(audit_path, index=False)
audit_metadata = OrderedDict(
artifact_type="provider_grouped_fold_assignments_and_oof_predictions",
dataset=dataset_name,
format="parquet",
path=f"{PORTABLE_AUDIT_PREFIX}/{audit_path.name}",
row_count=int(len(audit)),
columns=list(audit.columns),
target_column="target_log1p_NBA",
target_formula=TARGET_FORMULA,
group_column=GROUP_COLUMN,
n_groups=int(frame[GROUP_COLUMN].nunique()),
splitter="GroupKFold",
n_splits=int(len(splits)),
group_disjoint=True,
fold_column="fold",
oof_prediction_columns=[
"linear_regression_oof_prediction",
"random_forest_oof_prediction",
],
reference_levels=reference_levels,
clipping_counts=frame_metadata.get("clipping_counts", {}),
exclusion_counts=frame_metadata.get("exclusion_counts", {}),
fold_summaries=fold_summaries,
)
with open(metadata_path, "w") as handle:
json.dump(audit_metadata, handle, indent=2, default=_json_default)
return {
"oof_predictions": OrderedDict(
format="parquet",
location="SITE_DATA_AUDIT",
path=f"{PORTABLE_AUDIT_PREFIX}/{audit_path.name}",
rows=int(len(audit)),
),
"metadata": OrderedDict(
format="json",
location="SITE_DATA_AUDIT",
path=f"{PORTABLE_AUDIT_PREFIX}/{metadata_path.name}",
),
}
# ---------------------------------------------------------------------------
# Orchestration
# ---------------------------------------------------------------------------
def run() -> dict:
print("[Phase 5] Loading feature parquets ...")
inp = pd.read_parquet(DATA_PROC / "inpatient_features.parquet")
out = pd.read_parquet(DATA_PROC / "outpatient_features.parquet")
print("[Phase 5] Building modeling frames (target log1p(NBA)) ...")
inp_m = build_modeling_frame(
inp, "drg_code", "avg_submitted_charge", "avg_total_payment",
"total_discharges", "drg_mdc", "inpatient"
)
out_m = build_modeling_frame(
out, "apc_code", "avg_submitted_charge", "avg_allowed_amount",
"apc_services", "apc_family", "outpatient"
)
print(f" inpatient modeling rows: {len(inp_m):,}")
print(f" outpatient modeling rows: {len(out_m):,}")
artifacts = {}
for name, frame in [("inpatient", inp_m), ("outpatient", out_m)]:
print(f"\n[Phase 5] --- {name} ---")
X = _prepare_model_inputs(frame)
category_levels, reference_levels = _reference_spec(X)
splits = _grouped_splits(frame, n_splits=VALIDATION_N_SPLITS)
n_groups = int(frame[GROUP_COLUMN].nunique())
print("[Phase 5] Linear regression with provider-grouped CV ...")
lr = LinearRegression()
lr_res = cv_score(
frame, lr, n_splits=VALIDATION_N_SPLITS, splits=splits,
category_levels=category_levels, reference_levels=reference_levels,
n_repeats=PERMUTATION_REPEATS,
)
print(
f" held-out R2 (GroupKFold) = {lr_res['r2_mean']} "
f"+/- {lr_res['r2_std']} on N={lr_res['n_samples']:,}"
)
lr_coefs = linear_coefficients(
frame, lr_res["fitted_pipeline"], reference_levels
)
print("[Phase 5] Random forest (150 trees) with provider-grouped CV ...")
rf = RandomForestRegressor(
n_estimators=150,
max_depth=10,
n_jobs=1,
random_state=DEFAULT_SEED,
min_samples_leaf=50,
)
rf_res = cv_score(
frame, rf, n_splits=VALIDATION_N_SPLITS, splits=splits,
category_levels=category_levels, reference_levels=reference_levels,
n_repeats=PERMUTATION_REPEATS,
)
print(
f" held-out R2 (GroupKFold) = {rf_res['r2_mean']} "
f"+/- {rf_res['r2_std']}"
)
rf_imp = rf_importance(
frame,
rf_res["fitted_pipeline"],
heldout_permutation=rf_res["permutation_importance"],
heldout_permutation_by_fold=rf_res["permutation_importance_by_fold"],
importance_evaluation=rf_res["importance_evaluation"],
)
fold_summaries = _combined_fold_summaries(lr_res, rf_res)
frame_metadata = frame.attrs.get("predictors_metadata", {})
audit_artifacts = _save_validation_artifacts(
name, frame, splits, lr_res, rf_res, reference_levels,
frame_metadata, fold_summaries,
)
coefs_sorted = sorted(lr_coefs, key=lambda row: row["coefficient"])
validation = OrderedDict(
estimand=VALIDATION_ESTIMAND,
splitter="GroupKFold",
splitter_parameters=OrderedDict(
n_splits=VALIDATION_N_SPLITS,
shuffle=False,
),
group_column=GROUP_COLUMN,
n_groups=n_groups,
groups_disjoint=True,
fold_summaries=fold_summaries,
oof_r2=OrderedDict(
linear_regression=lr_res["r2_oof"],
random_forest=rf_res["r2_oof"],
),
)
artifacts[name] = OrderedDict(
n_samples=int(len(frame)),
target="log1p(NBA) - Agrawal-normalized billing",
target_column=TARGET_COLUMN,
target_formula=TARGET_FORMULA,
payment_column=frame_metadata.get("payment_column"),
payment_measure={
"avg_total_payment": "average total payment",
"avg_medicare_payment": "average Medicare payment",
"avg_allowed_amount": "average Medicare allowed amount",
}.get(frame_metadata.get("payment_column"), frame_metadata.get("payment_column")),
normalization_reference=NORMALIZATION_REFERENCE,
features=MODEL_FEATURES,
excluded_features=["charge_to_payment_ratio", "census_region"],
feature_exclusions=OrderedDict(
charge_to_payment_ratio=OrderedDict(
excluded=True,
reason="Mechanically derived from billing/payment, which enter NBA; including it leaks the target.",
),
census_region=OrderedDict(
excluded=True,
reason="Removed because state is retained and census_region is nested geography.",
),
),
excluded_geography=frame_metadata.get("excluded_geography", {}),
reference_levels=reference_levels,
category_levels=category_levels,
categorical_encoding=OrderedDict(
encoder="OneHotEncoder",
handle_unknown="ignore",
explicit_reference_coding=True,
min_frequency=None,
rare_category_grouping=False,
),
clipping_counts=frame_metadata.get("clipping_counts", {}),
exclusion_counts=frame_metadata.get("exclusion_counts", {}),
validation_estimand=VALIDATION_ESTIMAND,
splitter="GroupKFold",
group_column=GROUP_COLUMN,
n_groups=n_groups,
fold_summaries=fold_summaries,
validation=validation,
audit_artifacts=audit_artifacts,
importance_evaluation=OrderedDict(
linear_regression=lr_res["importance_evaluation"],
random_forest=rf_res["importance_evaluation"],
),
linear_regression=OrderedDict(
cv_r2_mean=lr_res["r2_mean"],
cv_r2_std=lr_res["r2_std"],
cv_r2_folds=lr_res["r2_folds"],
cv_r2_oof=lr_res["r2_oof"],
fold_summaries=lr_res["fold_summaries"],
top_positive_coefs=[
row for row in coefs_sorted if row["coefficient"] > 0
][-10:][::-1],
top_negative_coefs=coefs_sorted[:10],
all_coefs_count=len(lr_coefs),
fit_scope="descriptive_full_sample",
coefficients_fit_scope="descriptive_full_sample",
held_out_performance_source="provider-grouped CV above; coefficients are not held-out estimates",
permutation_importance=lr_res["permutation_importance"],
permutation_importance_by_fold=lr_res["permutation_importance_by_fold"],
),
random_forest=OrderedDict(
cv_r2_mean=rf_res["r2_mean"],
cv_r2_std=rf_res["r2_std"],
cv_r2_folds=rf_res["r2_folds"],
cv_r2_oof=rf_res["r2_oof"],
fold_summaries=rf_res["fold_summaries"],
permutation_importance=rf_imp["permutation_importance"],
permutation_importance_by_fold=rf_imp["permutation_importance_by_fold"],
mdi_aggregated=rf_imp["mdi_aggregated"],
fit_scope="descriptive_full_sample",
mdi_fit_scope=rf_imp["mdi_fit_scope"],
held_out_performance_source="provider-grouped CV above; MDI is descriptive full-sample importance",
),
)
# Keep the site's existing limitations shape, while describing the actual
# target and held-out estimand rather than a residual or state-level target.
artifacts["limitations"] = {
"data_constraints": [
"CMS PUFs omit length-of-stay (LOS), so we cannot adjust for case severity the way OIG did.",
"No case-mix index (CMI) at the row level - variations in patient acuity within a DRG are unobservable.",
"No patient-outcome data, so we cannot validate whether high-charge providers deliver correspondingly better care.",
"charge_to_payment_ratio is partly a function of geography (including wage-index effects) and is not a cost or markup measure.",
"Geography is at state granularity - within-state hospital-level cost variation is unobservable.",
],
"what_would_help": [
"CMS Hospital Cost Reports with cost-to-charge ratios (CCR) per facility - actual OIG-level analysis.",
"Hospital Compare / Care Compare outcome metrics - could correlate outliers with mortality / readmission.",
"Hospital Referral Region (HRR) granularity (Mulani showed it) - sub-state geographic precision.",
"Additional annual source files would be needed before estimating a temporal trend or future-period performance.",
],
"interpretive_framing": [
"Reported R2 is provider-grouped held-out performance for log1p(NBA), not a causal estimate.",
"Coefficients and random-forest MDI are descriptive summaries from full-data fits; they must not be read as held-out performance.",
"Cross-dataset outlier co-occurrence (Phase 4) is a separate screening hypothesis only; it is selected, cross-period, and not an outcome validation.",
"charge_to_payment_ratio was deliberately EXCLUDED from features because it is mechanically derived from billing/payment inputs to NBA.",
],
"open_questions": [
"Which cross-dataset co-occurrence logic could payer/plan integrity teams adopt as a screening tool?",
"Could this be linked to a charge-audit workflow as a high-charge-authorization flag?",
"How does provider behavior evolve over multiple years?",
],
}
model_validation_manifest = OrderedDict(
artifact_type="provider_grouped_predictor_validation",
target_column=TARGET_COLUMN,
target_formula=TARGET_FORMULA,
normalization_reference=NORMALIZATION_REFERENCE,
validation_estimand=VALIDATION_ESTIMAND,
splitter="GroupKFold",
group_column=GROUP_COLUMN,
n_splits=VALIDATION_N_SPLITS,
datasets=OrderedDict(
(name, OrderedDict(
n_samples=artifacts[name]["n_samples"],
n_groups=artifacts[name]["n_groups"],
payment_column=artifacts[name]["payment_column"],
payment_measure=artifacts[name]["payment_measure"],
reference_levels=artifacts[name]["reference_levels"],
fold_summaries=artifacts[name]["fold_summaries"],
audit_artifacts=artifacts[name]["audit_artifacts"],
importance_evaluation=artifacts[name]["importance_evaluation"],
))
for name in ("inpatient", "outpatient")
),
)
artifacts["model_validation_manifest"] = model_validation_manifest
# Persist the site-compatible summary. The audit Parquets and metadata
# JSONs are written separately under DATA_PROC by the loop above.
save_json(artifacts, "predictors.json")
save_json(model_validation_manifest, "model_validation_manifest.json")
print("\n[Phase 5] DONE.")
return artifacts
def _json_default(o):
if isinstance(o, (np.integer,)):
return int(o)
if isinstance(o, (np.floating,)):
value = float(o)
return None if pd.isna(value) else value
if isinstance(o, (np.ndarray,)):
return o.tolist()
if pd.isna(o):
return None
raise TypeError(f"not serializable: {type(o)}")
if __name__ == "__main__": # pragma: no cover
run()