Recipe 6.6: Python Implementation Example

pip install boto3 numpy pandas scikit-learn

import numpy as np
import pandas as pd
from sklearn.neighbors import NearestNeighbors
from sklearn.preprocessing import StandardScaler
from datetime import datetime, timezone
import json
import logging
import boto3
from decimal import Decimal

# Structured logging. In production, use JSON-formatted output for
# CloudWatch Logs Insights queries. Never log PHI field values.
logger = logging.getLogger(__name__)
logger.setLevel(logging.INFO)

# --- Feature Configuration ---
# These are the features we use for diabetes patient similarity.
# Each feature has a weight that encodes clinical importance.
# Higher weight = more influence on similarity score.
#
# These weights are NOT learned from data. They come from clinical SMEs
# who understand which factors drive diabetes outcomes. Treat this as
# a living document that gets refined as you validate the model.

FEATURE_WEIGHTS = {
    "age": 1.0,
    "bmi": 1.5,
    "a1c": 3.0,              # A1C is the primary diabetes marker
    "systolic_bp": 1.0,
    "years_since_diagnosis": 2.0,
    "medication_count": 1.0,
    "on_insulin": 2.0,       # insulin status changes trajectory significantly
    "on_metformin": 1.5,
    "on_glp1": 1.5,
    "heart_failure": 2.5,    # comorbidities that change the game
    "ckd": 2.5,
    "depression": 1.5,
    "hypertension": 1.0,
    "ed_visits_12mo": 1.5,
    "hospitalizations_12mo": 2.0,
}

# How many similar patients to retrieve per query.
K_NEIGHBORS = 20

# Maximum distance threshold. Neighbors beyond this are "not similar enough"
# and should not be presented to the care manager.
# This value depends on your feature scaling; tune it empirically.
MAX_DISTANCE_THRESHOLD = 3.0

# DynamoDB table for caching similarity results.
CACHE_TABLE_NAME = "patient-similarity-cache"

# S3 bucket and prefix for the feature store.
FEATURE_STORE_BUCKET = "my-health-system-features"
FEATURE_STORE_PREFIX = "patient-features/diabetes/v2026-03/"

# Cache TTL in seconds (24 hours). Results expire when the feature store refreshes.
CACHE_TTL_SECONDS = 86400

def generate_synthetic_cohort(n_patients: int = 500, seed: int = 42) -> pd.DataFrame:
    """
    Generate a synthetic patient cohort with realistic feature distributions.

    This creates fake-but-plausible patient data for demonstrating the
    similarity algorithm. The distributions are loosely based on published
    Type 2 diabetes population statistics, but they're synthetic.
    Do NOT use this for clinical validation.

    Args:
        n_patients: Number of synthetic patients to generate.
        seed: Random seed for reproducibility.

    Returns:
        DataFrame with one row per patient, columns matching FEATURE_WEIGHTS keys.
    """
    rng = np.random.default_rng(seed)

    patients = pd.DataFrame({
        "patient_id": [f"PAT-{i:05d}" for i in range(n_patients)],

        # Age: Type 2 diabetes population skews older. Mean ~58, std ~12.
        "age": rng.normal(58, 12, n_patients).clip(30, 90).astype(int),

        # BMI: overweight/obese population. Mean ~32, std ~6.
        "bmi": rng.normal(32, 6, n_patients).clip(18, 55).round(1),

        # A1C: the primary diabetes marker. Mean ~7.8, std ~1.5.
        # Values below 5.7 are non-diabetic; we clip to realistic range.
        "a1c": rng.normal(7.8, 1.5, n_patients).clip(5.7, 14.0).round(1),

        # Systolic blood pressure. Mean ~138, std ~18.
        "systolic_bp": rng.normal(138, 18, n_patients).clip(90, 200).astype(int),

        # Years since diabetes diagnosis. Exponential-ish distribution.
        "years_since_diagnosis": rng.exponential(5, n_patients).clip(0, 30).round(1),

        # Medication count (all medications, not just diabetes).
        "medication_count": rng.poisson(5, n_patients).clip(0, 20),

        # Binary flags for specific drug classes.
        "on_insulin": rng.binomial(1, 0.3, n_patients),
        "on_metformin": rng.binomial(1, 0.75, n_patients),
        "on_glp1": rng.binomial(1, 0.2, n_patients),

        # Comorbidity flags.
        "heart_failure": rng.binomial(1, 0.15, n_patients),
        "ckd": rng.binomial(1, 0.2, n_patients),
        "depression": rng.binomial(1, 0.25, n_patients),
        "hypertension": rng.binomial(1, 0.7, n_patients),

        # Utilization in the past 12 months.
        "ed_visits_12mo": rng.poisson(0.8, n_patients).clip(0, 10),
        "hospitalizations_12mo": rng.poisson(0.3, n_patients).clip(0, 5),
    })

    # Generate synthetic outcomes for the cohort.
    # Goal: A1C < 7.0 within 6 months. Probability depends on features.
    # This is a simplified outcome model for demonstration purposes.
    goal_probability = (
        0.5
        - 0.05 * (patients["a1c"] - 7.0)       # higher baseline A1C = harder
        + 0.1 * patients["on_metformin"]         # metformin helps
        + 0.15 * patients["on_glp1"]             # GLP-1 helps more
        - 0.1 * patients["heart_failure"]        # comorbidities make it harder
        - 0.08 * patients["ckd"]
        - 0.03 * patients["years_since_diagnosis"] / 10
    ).clip(0.1, 0.95)

    patients["goal_achieved"] = rng.binomial(1, goal_probability)

    # Time to goal (months) for those who achieved it. Faster for less severe cases.
    base_time = 3.0 + 0.5 * (patients["a1c"] - 6.5)
    patients["time_to_goal_months"] = np.where(
        patients["goal_achieved"] == 1,
        (base_time + rng.normal(0, 1, n_patients)).clip(1, 12).round(1),
        np.nan
    )

    # Interventions received (simplified: assign based on features).
    patients["intervention_metformin"] = patients["on_metformin"]
    patients["intervention_glp1"] = patients["on_glp1"]
    patients["intervention_lifestyle"] = rng.binomial(1, 0.6, n_patients)
    patients["intervention_insulin_titration"] = patients["on_insulin"]

    # Adverse events (simplified).
    patients["adverse_gi_intolerance"] = (
        patients["on_metformin"] * rng.binomial(1, 0.15, n_patients)
    )
    patients["adverse_hypoglycemia"] = (
        patients["on_insulin"] * rng.binomial(1, 0.1, n_patients)
    )

    return patients

def build_similarity_index(
    cohort: pd.DataFrame,
    feature_weights: dict,
) -> tuple:
    """
    Build a nearest-neighbor index from the patient cohort.

    This applies feature weights, scales the features, and fits a
    NearestNeighbors model for fast similarity queries.

    Args:
        cohort: DataFrame with patient features (one row per patient).
        feature_weights: Dict mapping feature names to importance weights.

    Returns:
        Tuple of (fitted NearestNeighbors model, fitted StandardScaler, feature column list).
    """
    feature_cols = list(feature_weights.keys())

    # Extract the feature matrix. Each row is a patient, each column is a feature.
    X = cohort[feature_cols].values.astype(float)

    # Standardize features to zero mean, unit variance.
    # Without this, features with large ranges (age: 30-90) would dominate
    # features with small ranges (on_insulin: 0-1) in distance calculations.
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)

    # Apply clinical weights. Multiply each feature column by its weight.
    # This makes high-weight features (A1C: 3.0) contribute more to distance
    # than low-weight features (age: 1.0).
    weights = np.array([feature_weights[col] for col in feature_cols])
    X_weighted = X_scaled * weights

    # Fit the nearest-neighbor index.
    # algorithm="ball_tree" works well for moderate dimensions (<50 features).
    # For very high dimensions or very large cohorts, use "brute" or switch
    # to an ANN library like FAISS or Annoy.
    nn_model = NearestNeighbors(
        n_neighbors=K_NEIGHBORS + 1,  # +1 because the patient is their own neighbor
        metric="euclidean",
        algorithm="ball_tree",
    )
    nn_model.fit(X_weighted)

    logger.info(
        "Built similarity index: %d patients, %d features",
        X_weighted.shape[0],
        X_weighted.shape[1],
    )

    return nn_model, scaler, feature_cols

def find_similar_patients(
    query_patient_id: str,
    cohort: pd.DataFrame,
    nn_model: NearestNeighbors,
    scaler: StandardScaler,
    feature_cols: list,
    feature_weights: dict,
    k: int = K_NEIGHBORS,
    max_distance: float = MAX_DISTANCE_THRESHOLD,
) -> list:
    """
    Find the k most similar patients to a query patient.

    Args:
        query_patient_id: The patient ID to find neighbors for.
        cohort: Full patient cohort DataFrame.
        nn_model: Fitted NearestNeighbors model.
        scaler: Fitted StandardScaler (same one used during index build).
        feature_cols: List of feature column names in order.
        feature_weights: Dict of feature weights.
        k: Number of neighbors to return.
        max_distance: Maximum distance threshold; farther neighbors are excluded.

    Returns:
        List of dicts with patient_id, distance, and similarity score.
    """
    # Look up the query patient's features.
    query_row = cohort[cohort["patient_id"] == query_patient_id]
    if query_row.empty:
        logger.warning("Patient %s not found in cohort", query_patient_id)
        return []

    # Extract, scale, and weight the query vector (same transforms as index build).
    query_features = query_row[feature_cols].values.astype(float)
    query_scaled = scaler.transform(query_features)
    weights = np.array([feature_weights[col] for col in feature_cols])
    query_weighted = query_scaled * weights

    # Query the index. Returns distances and indices of nearest neighbors.
    distances, indices = nn_model.kneighbors(query_weighted, n_neighbors=k + 1)

    # Flatten (kneighbors returns 2D arrays even for single queries).
    distances = distances[0]
    indices = indices[0]

    # Build results, skipping the query patient and filtering by max distance.
    results = []
    for dist, idx in zip(distances, indices):
        neighbor_id = cohort.iloc[idx]["patient_id"]

        # Skip self-match.
        if neighbor_id == query_patient_id:
            continue

        # Stop if we've exceeded the distance threshold.
        # kneighbors() returns results sorted by distance ascending, so once
        # we exceed the threshold, all remaining neighbors are also too far.
        # Note: if using approximate NN libraries (FAISS, Annoy), verify sort
        # order before relying on this early-exit pattern.
        if dist > max_distance:
            break

        results.append({
            "patient_id": neighbor_id,
            "distance": round(float(dist), 4),
            "similarity": round(1.0 / (1.0 + float(dist)), 4),
        })

        if len(results) >= k:
            break

    logger.info(
        "Found %d similar patients for %s (closest distance: %.3f)",
        len(results),
        query_patient_id,
        results[0]["distance"] if results else float("inf"),
    )

    return results

def aggregate_outcomes(
    similar_patients: list,
    cohort: pd.DataFrame,
) -> dict:
    """
    Aggregate outcomes for a set of similar patients.

    Computes goal achievement rates, common interventions, adverse events,
    and a confidence indicator based on cohort size.

    Args:
        similar_patients: List of dicts from find_similar_patients.
        cohort: Full cohort DataFrame (includes outcome columns).

    Returns:
        Dict with aggregated outcome summary.
    """
    if not similar_patients:
        return {"cohort_size": 0, "confidence": "insufficient"}

    # Get the patient IDs of similar patients.
    neighbor_ids = [p["patient_id"] for p in similar_patients]

    # Filter the cohort to just the similar patients.
    neighbors_df = cohort[cohort["patient_id"].isin(neighbor_ids)]

    cohort_size = len(neighbors_df)

    # Goal achievement rate.
    goal_rate = neighbors_df["goal_achieved"].mean()

    # Median time to goal (only for those who achieved it).
    achievers = neighbors_df[neighbors_df["goal_achieved"] == 1]
    median_time = (
        float(achievers["time_to_goal_months"].median())
        if not achievers.empty
        else None
    )

    # Intervention frequency among successful patients.
    intervention_cols = [
        "intervention_metformin",
        "intervention_glp1",
        "intervention_lifestyle",
        "intervention_insulin_titration",
    ]
    if not achievers.empty:
        intervention_freq = {
            col.replace("intervention_", ""): round(float(achievers[col].mean()), 2)
            for col in intervention_cols
        }
    else:
        intervention_freq = {}

    # Sort interventions by frequency (most common first).
    intervention_freq = dict(
        sorted(intervention_freq.items(), key=lambda x: x[1], reverse=True)
    )

    # Adverse event rates.
    adverse_cols = ["adverse_gi_intolerance", "adverse_hypoglycemia"]
    adverse_rates = {
        col.replace("adverse_", ""): round(float(neighbors_df[col].mean()), 3)
        for col in adverse_cols
        if neighbors_df[col].sum() > 0
    }

    # Confidence based on cohort size.
    if cohort_size >= 20:
        confidence = "high"
    elif cohort_size >= 10:
        confidence = "moderate"
    elif cohort_size >= 5:
        confidence = "low"
    else:
        confidence = "insufficient"

    summary = {
        "cohort_size": cohort_size,
        "confidence": confidence,
        "goal_achievement_rate": round(float(goal_rate), 3),
        "median_time_to_goal_months": median_time,
        "intervention_frequency": intervention_freq,
        "adverse_event_rates": adverse_rates,
    }

    logger.info(
        "Outcome summary: %d patients, %.0f%% goal achievement, confidence=%s",
        cohort_size,
        goal_rate * 100,
        confidence,
    )

    return summary

def cache_similarity_results(
    query_patient_id: str,
    similar_patients: list,
    outcome_summary: dict,
    feature_version: str = "v2026-03-15",
) -> dict:
    """
    Cache similarity results in DynamoDB with TTL-based expiration.

    Args:
        query_patient_id: The patient we queried for.
        similar_patients: List of similar patient matches.
        outcome_summary: Aggregated outcome data.
        feature_version: Version tag of the feature store used.

    Returns:
        The cached record (for confirmation).
    """
    dynamodb = boto3.resource("dynamodb")
    table = dynamodb.Table(CACHE_TABLE_NAME)

    now = int(datetime.now(timezone.utc).timestamp())

    # DynamoDB requires Decimal for numeric types, not float.
    # This helper converts floats in nested structures.
    record = {
        "patient_id": query_patient_id,
        "feature_version": feature_version,
        "computed_at": datetime.now(timezone.utc).isoformat(),
        "ttl": now + CACHE_TTL_SECONDS,
        "similar_patients": json.loads(
            json.dumps(similar_patients), parse_float=Decimal
        ),
        "outcome_summary": json.loads(
            json.dumps(outcome_summary), parse_float=Decimal
        ),
    }

    table.put_item(Item=record)

    logger.info(
        "Cached similarity results for %s (expires in %d seconds)",
        query_patient_id,
        CACHE_TTL_SECONDS,
    )

    return record

def check_cache(query_patient_id: str, feature_version: str) -> dict | None:
    """
    Check if we have cached similarity results for this patient.

    Returns the cached record if found and not expired, None otherwise.
    DynamoDB TTL handles expiration, but items may linger briefly after TTL.
    We check the feature_version to ensure results match current features.
    """
    dynamodb = boto3.resource("dynamodb")
    table = dynamodb.Table(CACHE_TABLE_NAME)

    response = table.get_item(
        Key={"patient_id": query_patient_id}
    )

    item = response.get("Item")
    if item and item.get("feature_version") == feature_version:
        logger.info("Cache hit for %s", query_patient_id)
        return item

    logger.info("Cache miss for %s", query_patient_id)
    return None

def run_patient_similarity_pipeline(
    query_patient_id: str,
    use_cache: bool = True,
    feature_version: str = "v2026-03-15",
) -> dict:
    """
    Run the full patient similarity pipeline for one patient.

    Generates a synthetic cohort (in production, this loads from S3),
    builds the similarity index, finds neighbors, aggregates outcomes,
    and returns a care-planning-ready summary.

    Args:
        query_patient_id: Patient ID to find similar patients for.
        use_cache: Whether to check/write DynamoDB cache.
        feature_version: Feature store version tag.

    Returns:
        Dict with similar patients and outcome summary.
    """

    # In production, check cache first to avoid redundant computation.
    # Commented out for local testing (requires DynamoDB table to exist).
    # if use_cache:
    #     cached = check_cache(query_patient_id, feature_version)
    #     if cached:
    #         return cached

    # Step 1: Load the patient cohort.
    # In production: load from S3 feature store (Parquet files).
    # Here: generate synthetic data for demonstration.
    print("Step 1: Loading patient cohort...")
    cohort = generate_synthetic_cohort(n_patients=500)
    print(f"  Loaded {len(cohort)} patients with {len(FEATURE_WEIGHTS)} features each")

    # Step 2: Build the similarity index.
    # In production: this is pre-built and loaded from a model artifact.
    # Here: we build it fresh each time (fast enough for 500 patients).
    print("Step 2: Building similarity index...")
    nn_model, scaler, feature_cols = build_similarity_index(cohort, FEATURE_WEIGHTS)
    print(f"  Index built with algorithm={nn_model.algorithm}")

    # Step 3: Find similar patients.
    print(f"Step 3: Finding patients similar to {query_patient_id}...")
    similar = find_similar_patients(
        query_patient_id=query_patient_id,
        cohort=cohort,
        nn_model=nn_model,
        scaler=scaler,
        feature_cols=feature_cols,
        feature_weights=FEATURE_WEIGHTS,
        k=K_NEIGHBORS,
        max_distance=MAX_DISTANCE_THRESHOLD,
    )
    print(f"  Found {len(similar)} similar patients")

    if not similar:
        print("  No similar patients found within distance threshold.")
        return {"query_patient_id": query_patient_id, "similar_patients": [], "outcome_summary": {}}

    # Print the top 5 matches for visibility.
    print("  Top 5 matches:")
    for match in similar[:5]:
        print(f"    {match['patient_id']}: similarity={match['similarity']}, distance={match['distance']}")

    # Step 4: Aggregate outcomes for the similar cohort.
    print("Step 4: Aggregating outcomes...")
    outcome_summary = aggregate_outcomes(similar, cohort)
    print(f"  Goal achievement rate: {outcome_summary['goal_achievement_rate']:.0%}")
    print(f"  Confidence: {outcome_summary['confidence']}")
    print(f"  Top interventions: {outcome_summary['intervention_frequency']}")

    # Step 5: Cache results (uncomment for production with DynamoDB).
    # if use_cache:
    #     cache_similarity_results(query_patient_id, similar, outcome_summary, feature_version)

    # Assemble the final result.
    # NOTE: similar_patients contains individual patient IDs. This is the
    # internal API response used for caching and audit logging. The care
    # planning UI should display only the aggregated outcome_summary by default.
    # If your organization permits drill-down into individual similar patient
    # records, gate that access behind break-the-glass authorization and log
    # each access as a HIPAA disclosure event.
    result = {
        "query_patient_id": query_patient_id,
        "feature_version": feature_version,
        "similar_patients": similar[:5],  # top 5 for display; outcome_summary uses all k neighbors
        "outcome_summary": outcome_summary,
    }

    print("\nDone. Result:")
    print(json.dumps(result, indent=2, default=str))

    return result

# --- Run it ---
if __name__ == "__main__":
    # Pick a patient from our synthetic cohort to query.
    # PAT-00042 is arbitrary; any valid patient_id works.
    result = run_patient_similarity_pipeline(query_patient_id="PAT-00042")

def explain_similarity(
    query_patient_id: str,
    neighbor_patient_id: str,
    cohort: pd.DataFrame,
    feature_cols: list,
    feature_weights: dict,
) -> list:
    """
    Explain which features make two patients similar (or different).

    Returns a ranked list of features by their contribution to the
    distance between the two patients. Features with small differences
    (high agreement) are what makes them similar. Features with large
    differences are where they diverge.

    Args:
        query_patient_id: The query patient.
        neighbor_patient_id: A similar patient to explain.
        cohort: Full cohort DataFrame.
        feature_cols: Feature column names.
        feature_weights: Feature importance weights.

    Returns:
        List of dicts with feature name, query value, neighbor value,
        and weighted contribution to distance.
    """
    query_row = cohort[cohort["patient_id"] == query_patient_id][feature_cols].iloc[0]
    neighbor_row = cohort[cohort["patient_id"] == neighbor_patient_id][feature_cols].iloc[0]

    explanations = []
    for col in feature_cols:
        q_val = float(query_row[col])
        n_val = float(neighbor_row[col])
        weight = feature_weights[col]
        # Weighted absolute difference on RAW (unscaled) values.
        # Note: the actual similarity model uses StandardScaler-transformed
        # values before applying weights, so this explanation is an
        # approximation. The ranking of contributions here may differ from
        # the true distance decomposition. For precise explanations, pass
        # the fitted scaler and compute contributions on scaled values.
        contribution = abs(q_val - n_val) * weight

        explanations.append({
            "feature": col,
            "query_value": q_val,
            "neighbor_value": n_val,
            "weight": weight,
            "distance_contribution": round(contribution, 3),
        })

    # Sort by contribution: smallest first (most similar features at top).
    explanations.sort(key=lambda x: x["distance_contribution"])

    return explanations

# After running the pipeline, explain why the top match is similar.
if __name__ == "__main__":
    cohort = generate_synthetic_cohort(n_patients=500)
    _, _, feature_cols = build_similarity_index(cohort, FEATURE_WEIGHTS)

    explanation = explain_similarity(
        query_patient_id="PAT-00042",
        neighbor_patient_id="PAT-00187",  # replace with actual top match
        cohort=cohort,
        feature_cols=feature_cols,
        feature_weights=FEATURE_WEIGHTS,
    )

    print("\nSimilarity explanation (most similar features first):")
    for item in explanation[:5]:
        print(f"  {item['feature']}: query={item['query_value']}, "
              f"neighbor={item['neighbor_value']} (contribution={item['distance_contribution']})")

    print("\nBiggest differences:")
    for item in explanation[-3:]:
        print(f"  {item['feature']}: query={item['query_value']}, "
              f"neighbor={item['neighbor_value']} (contribution={item['distance_contribution']})")