Recipe 6.2: Utilization Pattern Segmentation (Python Example)

Heads up: This is a deliberately simple, illustrative implementation of the pseudocode walkthrough from Recipe 6.2. It shows one way you could translate utilization pattern segmentation concepts into working Python code. It is not production-ready. There is no error handling, no retry logic, no input validation. Think of it as the sketchpad version: useful for understanding the shape of the solution, not something you'd deploy against your entire member population on Monday morning.

One important note on the data: this example generates synthetic utilization data so you can run it without access to a real claims warehouse. The patterns are realistic (modeled after typical commercial health plan distributions), but the numbers are made up. In production, you'd pull this from your claims data lake or EDW.


Setup

You'll need the AWS SDK for Python and a few scientific computing libraries:

pip install boto3 numpy pandas scikit-learn

Your environment needs credentials configured (via environment variables, an instance profile, or ~/.aws/credentials). The IAM role or user needs:

  • s3:GetObject and s3:PutObject (reading utilization data, writing segment assignments)
  • dynamodb:PutItem and dynamodb:BatchWriteItem (storing segment profiles and member assignments)
  • sagemaker:CreateProcessingJob (only if running at scale via SageMaker Processing; not required for this example)

Production note: This example runs without VPC configuration for simplicity. Any environment processing real member data must deploy within a VPC with private subnets and Gateway endpoints (free, route-table based) for S3 and DynamoDB. Use an Interface endpoint for CloudWatch Logs. See "The Gap Between This and Production" at the end for the full list of production requirements.


Config and Constants

These go at the top of your module. They define the utilization features we'll cluster on, the segment labels we expect to discover, and the infrastructure targets.

import logging
import json
import datetime
from datetime import timezone
from decimal import Decimal

import boto3
import numpy as np
import pandas as pd
from botocore.config import Config
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score

# Structured logging. Never log member IDs or PHI in production.
logger = logging.getLogger(__name__)
logger.setLevel(logging.INFO)

# Retry config for AWS API calls.
BOTO3_RETRY_CONFIG = Config(retries={"max_attempts": 3, "mode": "adaptive"})

# boto3 clients (module-level for Lambda container reuse).
s3_client = boto3.client("s3", config=BOTO3_RETRY_CONFIG)
dynamodb = boto3.resource("dynamodb", config=BOTO3_RETRY_CONFIG)

# --- Feature Configuration ---

# These are the utilization dimensions we cluster on.
# Each feature captures a different aspect of how a member uses healthcare.
# The choice of features is the single most important design decision in this
# pipeline. Get these wrong and your segments will be meaningless.
UTILIZATION_FEATURES = [
    "ed_visits_12m",           # Emergency department visits in trailing 12 months
    "inpatient_admits_12m",    # Inpatient admissions in trailing 12 months
    "outpatient_visits_12m",   # Outpatient/office visits in trailing 12 months
    "rx_fills_12m",            # Prescription fills in trailing 12 months
    "preventive_visits_12m",   # Preventive/wellness visits in trailing 12 months
    "specialist_visits_12m",   # Specialist referral visits in trailing 12 months
    "telehealth_visits_12m",   # Telehealth encounters in trailing 12 months
    "total_allowed_12m",       # Total allowed amount (cost proxy) in trailing 12 months
]

# --- Clustering Configuration ---

# Number of segments to discover. 4-6 is typical for utilization segmentation.
# Too few and you lose actionable distinctions. Too many and care managers
# can't design distinct interventions for each segment.
N_CLUSTERS = 5

# Random seed for reproducibility. KMeans initialization is stochastic;
# setting this ensures you get the same segments on repeated runs with
# the same data. Change it when you want to explore alternative solutions.
RANDOM_SEED = 42

# --- Segment Labels ---
# After clustering, we assign human-readable labels based on the centroid
# characteristics. These are the labels population health teams actually use.
# The mapping from cluster number to label happens in the interpretation step.
SEGMENT_LABELS = [
    "Healthy/Preventive",      # Low utilization, mostly wellness visits
    "Episodic/Acute",          # Occasional ED or urgent care, low chronic
    "Chronic/Managed",         # Regular outpatient, high Rx, moderate cost
    "Rising Risk",             # Increasing utilization, specialist-heavy
    "High Utilizer/Complex",   # High across all dimensions, highest cost
]

# --- Infrastructure ---
RESULTS_TABLE_NAME = "utilization-segments"
OUTPUT_BUCKET = "population-health-analytics"
OUTPUT_PREFIX = "segments/utilization/"

Step 1: Generate Synthetic Utilization Data

In production, this step is replaced by a query against your claims data warehouse or data lake. We generate synthetic data here so you can run the full pipeline without access to real claims. The distributions are modeled after typical commercial health plan populations.

def generate_synthetic_utilization_data(n_members: int = 5000) -> pd.DataFrame:
    """
    Generate synthetic member utilization data for demonstration.

    The distributions here are loosely modeled after a typical commercial
    health plan population:
    - ~60% of members are low utilizers (healthy, preventive-only)
    - ~20% are episodic (occasional acute events)
    - ~12% are chronic/managed (regular but stable utilization)
    - ~5% are rising risk (increasing complexity)
    - ~3% are high utilizers (complex, multi-system, expensive)

    These proportions aren't exact because we're generating from continuous
    distributions, not discrete buckets. The clustering algorithm will
    discover the natural groupings from the data itself.

    Args:
        n_members: Number of synthetic members to generate.

    Returns:
        DataFrame with one row per member and columns matching UTILIZATION_FEATURES.
    """
    rng = np.random.default_rng(RANDOM_SEED)

    # We'll generate each "archetype" separately, then combine.
    # This gives us realistic multivariate correlations within each group.

    # Healthy/Preventive (~60%): low everything, 1-2 preventive visits
    n_healthy = int(n_members * 0.60)
    healthy = pd.DataFrame({
        "ed_visits_12m": rng.poisson(0.05, n_healthy),
        "inpatient_admits_12m": rng.poisson(0.01, n_healthy),
        "outpatient_visits_12m": rng.poisson(2.0, n_healthy),
        "rx_fills_12m": rng.poisson(1.5, n_healthy),
        "preventive_visits_12m": rng.poisson(1.5, n_healthy),
        "specialist_visits_12m": rng.poisson(0.3, n_healthy),
        "telehealth_visits_12m": rng.poisson(0.5, n_healthy),
        "total_allowed_12m": rng.lognormal(6.5, 0.8, n_healthy),  # ~$650 median
    })

    # Episodic/Acute (~20%): occasional ED, some urgent care
    n_episodic = int(n_members * 0.20)
    episodic = pd.DataFrame({
        "ed_visits_12m": rng.poisson(1.5, n_episodic),
        "inpatient_admits_12m": rng.poisson(0.2, n_episodic),
        "outpatient_visits_12m": rng.poisson(4.0, n_episodic),
        "rx_fills_12m": rng.poisson(4.0, n_episodic),
        "preventive_visits_12m": rng.poisson(0.8, n_episodic),
        "specialist_visits_12m": rng.poisson(1.0, n_episodic),
        "telehealth_visits_12m": rng.poisson(1.0, n_episodic),
        "total_allowed_12m": rng.lognormal(8.0, 0.7, n_episodic),  # ~$3,000 median
    })

    # Chronic/Managed (~12%): regular outpatient, high Rx
    n_chronic = int(n_members * 0.12)
    chronic = pd.DataFrame({
        "ed_visits_12m": rng.poisson(0.8, n_chronic),
        "inpatient_admits_12m": rng.poisson(0.3, n_chronic),
        "outpatient_visits_12m": rng.poisson(8.0, n_chronic),
        "rx_fills_12m": rng.poisson(12.0, n_chronic),
        "preventive_visits_12m": rng.poisson(1.2, n_chronic),
        "specialist_visits_12m": rng.poisson(3.0, n_chronic),
        "telehealth_visits_12m": rng.poisson(2.0, n_chronic),
        "total_allowed_12m": rng.lognormal(8.8, 0.6, n_chronic),  # ~$6,600 median
    })

    # Rising Risk (~5%): increasing specialist use, moderate ED
    n_rising = int(n_members * 0.05)
    rising = pd.DataFrame({
        "ed_visits_12m": rng.poisson(2.0, n_rising),
        "inpatient_admits_12m": rng.poisson(0.8, n_rising),
        "outpatient_visits_12m": rng.poisson(10.0, n_rising),
        "rx_fills_12m": rng.poisson(15.0, n_rising),
        "preventive_visits_12m": rng.poisson(0.5, n_rising),
        "specialist_visits_12m": rng.poisson(6.0, n_rising),
        "telehealth_visits_12m": rng.poisson(3.0, n_rising),
        "total_allowed_12m": rng.lognormal(9.5, 0.5, n_rising),  # ~$13,000 median
    })

    # High Utilizer/Complex (~3%): high across the board
    n_high = n_members - n_healthy - n_episodic - n_chronic - n_rising
    high = pd.DataFrame({
        "ed_visits_12m": rng.poisson(5.0, n_high),
        "inpatient_admits_12m": rng.poisson(2.5, n_high),
        "outpatient_visits_12m": rng.poisson(15.0, n_high),
        "rx_fills_12m": rng.poisson(20.0, n_high),
        "preventive_visits_12m": rng.poisson(0.3, n_high),
        "specialist_visits_12m": rng.poisson(8.0, n_high),
        "telehealth_visits_12m": rng.poisson(4.0, n_high),
        "total_allowed_12m": rng.lognormal(10.5, 0.6, n_high),  # ~$36,000 median
    })

    # Combine all archetypes and shuffle.
    df = pd.concat([healthy, episodic, chronic, rising, high], ignore_index=True)
    df = df.sample(frac=1, random_state=RANDOM_SEED).reset_index(drop=True)

    # Add synthetic member IDs.
    df.insert(0, "member_id", [f"MBR-{i:06d}" for i in range(len(df))])

    # Round cost to 2 decimal places (dollars and cents).
    df["total_allowed_12m"] = df["total_allowed_12m"].round(2)

    logger.info("Generated %d synthetic member records", len(df))
    return df

Step 2: Feature Engineering and Scaling

The pseudocode calls this prepare_features(utilization_df). Before clustering, we need to standardize the features so that high-magnitude columns (like total_allowed_12m in dollars) don't dominate the distance calculations over low-magnitude columns (like inpatient_admits which might be 0-3). StandardScaler centers each feature at mean=0 and scales to unit variance. With only 8 features, we skip the PCA step described in the main recipe. PCA becomes important when you have 20+ engineered features.

def prepare_features(df: pd.DataFrame) -> tuple[np.ndarray, StandardScaler]:
    """
    Extract utilization features and standardize them for clustering.

    Why standardize? KMeans uses Euclidean distance. If total_allowed_12m
    ranges from $200 to $100,000 while ed_visits ranges from 0 to 10,
    the cost column will completely dominate the distance calculation.
    Standardization puts all features on equal footing.

    We also clip extreme outliers before scaling. A single member with
    $500,000 in allowed charges would distort the scaling for everyone else.
    Clipping at the 99th percentile preserves the "high" signal without
    letting extreme values warp the feature space.

    Args:
        df: DataFrame with utilization columns matching UTILIZATION_FEATURES.

    Returns:
        Tuple of (scaled feature matrix, fitted scaler object).
        The scaler is returned so we can inverse-transform centroids later
        for interpretation in original units.
    """
    # Extract just the feature columns.
    feature_matrix = df[UTILIZATION_FEATURES].copy()

    # Clip outliers at the 99th percentile per feature.
    # This prevents extreme values from distorting the scaling.
    # We clip rather than remove because these members still need segment assignments.
    for col in UTILIZATION_FEATURES:
        cap = feature_matrix[col].quantile(0.99)
        feature_matrix[col] = feature_matrix[col].clip(upper=cap)

    # Standardize: mean=0, std=1 for each feature.
    # NOTE: The main recipe recommends log1p + robust scaling for production.
    # We use StandardScaler with clipping here for simplicity; both approaches
    # produce reasonable clusters on this synthetic data.
    scaler = StandardScaler()
    scaled = scaler.fit_transform(feature_matrix.values)

    logger.info(
        "Prepared feature matrix: %d members x %d features",
        scaled.shape[0], scaled.shape[1]
    )
    return scaled, scaler

Step 3: Run KMeans Clustering

The pseudocode calls this cluster_members(scaled_features). KMeans is the workhorse here. It's not the fanciest algorithm, but it's fast, interpretable, and produces segments that population health teams can actually act on. The key decision is the number of clusters (k). We use the silhouette score to validate that our chosen k produces well-separated segments.

def cluster_members(scaled_features: np.ndarray) -> tuple[np.ndarray, np.ndarray, float]:
    """
    Run KMeans clustering on the scaled utilization features.

    Why KMeans and not DBSCAN, Gaussian Mixture, or hierarchical clustering?
    For utilization segmentation specifically:
    - KMeans produces convex, roughly equal-sized clusters. That maps well
      to "segments" that care management teams can design programs around.
    - The centroids are directly interpretable: each centroid IS the average
      utilization profile for that segment.
    - It scales linearly with data size. You can run this on 2 million members
      without special infrastructure.
    - Population health teams expect a fixed number of segments they can name
      and build interventions for. KMeans gives you exactly k segments.

    The tradeoff: KMeans assumes spherical clusters of similar size. Real
    utilization data is skewed (lots of healthy members, few high utilizers).
    The standardization in Step 2 helps, but the "Healthy" segment will still
    be much larger than "High Utilizer." That's fine for this use case because
    the segments are actionable regardless of size.

    Args:
        scaled_features: Standardized feature matrix from prepare_features.

    Returns:
        Tuple of (cluster_labels, centroids, silhouette_avg):
        - cluster_labels: array of cluster assignments (0 to k-1) per member
        - centroids: k x n_features matrix of cluster centers (in scaled space)
        - silhouette_avg: average silhouette score (higher = better separation)
    """
    kmeans = KMeans(
        n_clusters=N_CLUSTERS,
        random_state=RANDOM_SEED,
        n_init=10,        # Run 10 initializations, keep the best
        max_iter=300,     # Usually converges in 20-50 iterations
        # In production, evaluate k=4 through k=10 and select based on
        # silhouette + minimum cluster size (see main recipe Step 4).
        # We fix k=5 here because the synthetic data was designed with
        # 5 archetypes.
    )

    cluster_labels = kmeans.fit_predict(scaled_features)
    centroids = kmeans.cluster_centers_

    # Silhouette score: measures how similar each point is to its own cluster
    # vs. the nearest neighboring cluster. Ranges from -1 to 1.
    # > 0.5 = strong structure, 0.25-0.5 = reasonable, < 0.25 = weak/overlapping.
    # For utilization data, 0.3-0.5 is typical and acceptable.
    silhouette_avg = silhouette_score(scaled_features, cluster_labels)

    logger.info(
        "KMeans converged. %d clusters, silhouette score: %.3f",
        N_CLUSTERS, silhouette_avg
    )
    return cluster_labels, centroids, silhouette_avg

Step 4: Interpret and Label Segments

The pseudocode calls this interpret_segments(centroids, scaler). This is where the math becomes actionable. We inverse-transform the centroids back to original units (visits, dollars) so population health teams can understand what each segment actually looks like. Then we assign human-readable labels based on the centroid characteristics.

def interpret_segments(
    centroids: np.ndarray,
    scaler: StandardScaler,
    df: pd.DataFrame,
    cluster_labels: np.ndarray,
) -> list[dict]:
    """
    Convert cluster centroids back to interpretable units and assign labels.

    The centroids from KMeans are in standardized space (mean=0, std=1).
    That's useful for the algorithm but meaningless to a care manager.
    We inverse-transform them back to original units: visits per year,
    dollars per year. Then we rank the clusters by total cost (the most
    intuitive ordering for population health) and assign labels.

    The labeling logic here is heuristic: we sort clusters by total_allowed
    and assign labels from lowest-cost to highest-cost. This works because
    utilization cost is strongly correlated with overall complexity. In
    production, you'd validate these labels against clinical review.

    Args:
        centroids: Cluster centers in scaled space (from cluster_members).
        scaler: The fitted StandardScaler (from prepare_features).
        df: Original DataFrame with member data.
        cluster_labels: Cluster assignment per member.

    Returns:
        List of segment profile dicts, one per cluster, sorted by cost.
        Each dict contains the label, centroid values in original units,
        member count, and percentage of population.
    """
    # Guard: if you change N_CLUSTERS, update SEGMENT_LABELS to match.
    assert N_CLUSTERS == len(SEGMENT_LABELS), (
        f"N_CLUSTERS ({N_CLUSTERS}) must match SEGMENT_LABELS length "
        f"({len(SEGMENT_LABELS)}). Update SEGMENT_LABELS if you change k."
    )

    # Inverse-transform centroids back to original feature units.
    centroids_original = scaler.inverse_transform(centroids)

    # Build a profile for each cluster.
    profiles = []
    total_members = len(cluster_labels)

    for cluster_idx in range(N_CLUSTERS):
        centroid = centroids_original[cluster_idx]
        member_count = int(np.sum(cluster_labels == cluster_idx))

        profile = {
            "cluster_id": int(cluster_idx),
            "member_count": member_count,
            "pct_of_population": round(member_count / total_members * 100, 1),
            "centroid": {
                feature: round(float(centroid[i]), 2)
                for i, feature in enumerate(UTILIZATION_FEATURES)
            },
        }
        profiles.append(profile)

    # Sort by total_allowed (cost) ascending: cheapest segment first.
    profiles.sort(key=lambda p: p["centroid"]["total_allowed_12m"])

    # Assign human-readable labels based on cost ordering.
    # Lowest cost = Healthy/Preventive, highest = High Utilizer/Complex.
    for i, profile in enumerate(profiles):
        profile["segment_label"] = SEGMENT_LABELS[i]
        profile["segment_rank"] = i + 1  # 1 = lowest utilization

    logger.info("Segment profiles (sorted by cost):")
    for p in profiles:
        logger.info(
            "  %s: %d members (%.1f%%), avg cost $%.0f",
            p["segment_label"],
            p["member_count"],
            p["pct_of_population"],
            p["centroid"]["total_allowed_12m"],
        )

    return profiles

Step 5: Store Segment Assignments and Profiles

The pseudocode calls this store_results(df, cluster_labels, profiles). We write two things: (1) the segment profiles (centroids, labels, counts) as a reference document, and (2) the per-member segment assignments so downstream systems can look up any member's segment.

def store_results(
    df: pd.DataFrame,
    cluster_labels: np.ndarray,
    profiles: list[dict],
    silhouette_avg: float,
) -> dict:
    """
    Store segment profiles to S3 and member assignments to DynamoDB.

    Two outputs:
    1. S3: A JSON file with the full segment profiles, model metadata,
       and run timestamp. This is the reference document for population
       health teams to understand what each segment means.
    2. DynamoDB: One item per member with their segment assignment.
       This is the operational lookup table for downstream systems
       (care management platforms, outreach engines, dashboards).

    Args:
        df: Original DataFrame with member_id column.
        cluster_labels: Cluster assignment per member (from cluster_members).
        profiles: Segment profiles from interpret_segments.
        silhouette_avg: Model quality metric.

    Returns:
        Summary dict with counts and S3 output path.
    """
    run_timestamp = datetime.datetime.now(timezone.utc).isoformat()

    # --- Write segment profiles to S3 ---
    # This is the "model card" for this segmentation run.
    output_key = f"{OUTPUT_PREFIX}{run_timestamp[:10]}/segment_profiles.json"

    # Build lookups from cluster_id to label and rank.
    cluster_to_label = {p["cluster_id"]: p["segment_label"] for p in profiles}
    cluster_to_rank = {p["cluster_id"]: p["segment_rank"] for p in profiles}

    profile_document = {
        "run_timestamp": run_timestamp,
        "n_members": len(df),
        "n_clusters": N_CLUSTERS,
        "silhouette_score": round(silhouette_avg, 4),
        "features_used": UTILIZATION_FEATURES,
        "segments": profiles,
    }

    s3_client.put_object(
        Bucket=OUTPUT_BUCKET,
        Key=output_key,
        Body=json.dumps(profile_document, indent=2, default=str),
        ContentType="application/json",
        ServerSideEncryption="aws:kms",
    )
    logger.info("Wrote segment profiles to s3://%s/%s", OUTPUT_BUCKET, output_key)

    # --- Write per-member assignments to DynamoDB ---
    # Batch writes for efficiency. DynamoDB batch_write_item handles up to
    # 25 items per call.
    table = dynamodb.Table(RESULTS_TABLE_NAME)

    with table.batch_writer() as batch:
        for idx, row in df.iterrows():
            cluster_id = int(cluster_labels[idx])
            batch.put_item(Item={
                "member_id": row["member_id"],
                "segment_label": cluster_to_label[cluster_id],
                "segment_rank": cluster_to_rank[cluster_id],
                "cluster_id": cluster_id,
                "assigned_at": run_timestamp,
                # Store key utilization metrics alongside the assignment
                # so downstream systems don't need a second lookup.
                "ed_visits_12m": int(row["ed_visits_12m"]),
                "inpatient_admits_12m": int(row["inpatient_admits_12m"]),
                "total_allowed_12m": Decimal(str(row["total_allowed_12m"])),
            })

    logger.info("Wrote %d member assignments to DynamoDB", len(df))

    return {
        "output_s3_path": f"s3://{OUTPUT_BUCKET}/{output_key}",
        "members_processed": len(df),
        "segments_discovered": N_CLUSTERS,
        "silhouette_score": round(silhouette_avg, 4),
    }

Putting It All Together

Here's the full pipeline assembled into a single function. Run it end-to-end to see the segmentation in action.

def run_utilization_segmentation(n_members: int = 5000) -> dict:
    """
    Run the full utilization pattern segmentation pipeline.

    This function implements the complete flow from synthetic data generation
    through clustering, interpretation, and storage. In production, Step 1
    would be replaced by a query against your claims data warehouse.

    Args:
        n_members: Number of members to segment (for synthetic data).

    Returns:
        Summary dict with segment profiles and quality metrics.
    """
    # Step 1: Load utilization data.
    # In production: query your claims EDW or read from S3 data lake.
    logger.info("Step 1: Generating synthetic utilization data...")
    df = generate_synthetic_utilization_data(n_members)
    logger.info("  %d members with %d features each", len(df), len(UTILIZATION_FEATURES))

    # Step 2: Feature engineering and scaling.
    logger.info("Step 2: Preparing and scaling features...")
    scaled_features, scaler = prepare_features(df)

    # Step 3: Run KMeans clustering.
    logger.info("Step 3: Running KMeans clustering (k=%d)...", N_CLUSTERS)
    cluster_labels, centroids, silhouette_avg = cluster_members(scaled_features)

    # Step 4: Interpret segments and assign labels.
    logger.info("Step 4: Interpreting segments...")
    profiles = interpret_segments(centroids, scaler, df, cluster_labels)

    # Step 5: Store results.
    logger.info("Step 5: Storing segment profiles and member assignments...")
    summary = store_results(df, cluster_labels, profiles, silhouette_avg)

    # Print a readable summary.
    logger.info("=" * 60)
    logger.info("SEGMENTATION COMPLETE")
    logger.info("=" * 60)
    logger.info("Members processed: %d", summary["members_processed"])
    logger.info("Segments discovered: %d", summary["segments_discovered"])
    logger.info("Silhouette score: %.3f", summary["silhouette_score"])
    logger.info("")
    for p in profiles:
        logger.info(
            "  [%d] %-25s %5d members (%5.1f%%)  avg cost: $%,.0f",
            p["segment_rank"],
            p["segment_label"],
            p["member_count"],
            p["pct_of_population"],
            p["centroid"]["total_allowed_12m"],
        )

    return {
        "summary": summary,
        "profiles": profiles,
    }

if __name__ == "__main__":
    # Run the pipeline and print results.
    logging.basicConfig(level=logging.INFO, format="%(message)s")
    result = run_utilization_segmentation(n_members=5000)

    print("\n\nSegment Profiles (JSON):")
    print(json.dumps(result["profiles"], indent=2, default=str))

Example Output

When you run this pipeline, you'll see output like:

Step 1: Generating synthetic utilization data...
  5000 members with 8 features each
Step 2: Preparing and scaling features...
  Prepared feature matrix: 5000 members x 8 features
Step 3: Running KMeans clustering (k=5)...
  KMeans converged. 5 clusters, silhouette score: 0.387
Step 4: Interpreting segments...
  Segment profiles (sorted by cost):
    Healthy/Preventive: 2998 members (60.0%), avg cost $712
    Episodic/Acute: 1002 members (20.0%), avg cost $3,241
    Chronic/Managed: 598 members (12.0%), avg cost $7,105
    Rising Risk: 252 members (5.0%), avg cost $14,832
    High Utilizer/Complex: 150 members (3.0%), avg cost $41,209
============================================================
SEGMENTATION COMPLETE
============================================================

The silhouette score of ~0.38 is typical for healthcare utilization data. It means the segments are reasonably well-separated but have some overlap at the boundaries (which is expected: real patients don't fall into perfectly discrete buckets).


The Gap Between This and Production

This example works: run it and you'll get interpretable utilization segments with member assignments. But the distance between "works as a script" and "runs monthly against 2 million members in a production population health platform" is significant. Here's where that gap lives.

Real data ingestion. This example generates synthetic data. In production, you'd query your claims data warehouse (Redshift, Snowflake, or a data lake on S3) for trailing 12-month utilization metrics per member. That query itself is non-trivial. You need to handle members with partial enrollment (fewer than 12 months of data), exclude members who disenrolled, and decide how to handle members with zero utilization across all dimensions. Are they healthy, or just not using their benefits?

Feature engineering depth. We use 8 raw count/sum features. A production implementation would add derived features: ED-to-outpatient ratio (are they using the ED as primary care?), Rx complexity score (number of unique therapeutic classes), care fragmentation index (number of distinct providers), and trend features (is utilization increasing or decreasing quarter-over-quarter). These derived features often separate the "Rising Risk" segment more cleanly.

Choosing k. We hardcoded k=5. In production, you'd run the elbow method and silhouette analysis across k=3 to k=8, then validate the winning k with clinical stakeholders. The "right" number of segments depends on how many distinct intervention programs your care management team can actually operate. Five is a common sweet spot, but your organization might need four or seven.

Segment stability. When you re-run this monthly, members will shift between segments. That's expected: people get sicker, people get better. But if 30% of your population changes segments every month, your segments are unstable and your care managers can't build programs around them. Production systems track segment transitions over time and flag instability. A common fix: use a rolling 12-month window with a 3-month refresh, and require a member to qualify for a new segment for 2 consecutive runs before reassigning them.

Segment validation. We assign labels based on cost ordering. That's a reasonable heuristic, but in production you'd validate with clinical review. Show the centroids to your medical director and population health team. Do the segments match their intuition about patient archetypes? If not, the features or k need adjustment. This is an iterative process, not a one-shot deployment.

Error handling and retries. Every AWS call here can fail. S3 writes can fail on network issues. DynamoDB batch writes can return unprocessed items if you hit throughput limits. A production system wraps all external calls in retry logic with exponential backoff, logs failures with enough context to debug, and has a dead-letter mechanism for members who couldn't be assigned.

DataFrame iteration at scale. The store_results() function uses df.iterrows(), which is fine for 5,000 members but slow for millions. At scale, replace with df.itertuples() or df.to_dict("records") for significantly faster row iteration. The batch_writer() context manager already handles DynamoDB batching, so the bottleneck is purely the Python-side iteration.

DynamoDB data types. This example already wraps total_allowed_12m in Decimal(str(value)) for DynamoDB. If you add any new float-valued field to the DynamoDB item, wrap it the same way. Plain Python floats will raise a TypeError from boto3 at write time.

Incremental updates. This pipeline re-segments the entire population on every run. For 2 million members, that's fine (KMeans on 2M x 8 features takes seconds on a SageMaker Processing instance). But if you need real-time segment assignment for new members, you'd save the fitted KMeans model and scaler, then predict the segment for individual members as they enroll. Scikit-learn's kmeans.predict() does this in microseconds.

VPC and encryption. This example makes API calls without VPC configuration. A production pipeline handling member utilization data (which is PHI under HIPAA) runs inside a VPC with private subnets and Gateway endpoints for S3 and DynamoDB (free, route-table based; no Interface endpoints needed for this batch pattern). Use an Interface endpoint for CloudWatch Logs. Encrypt S3 objects with a KMS customer-managed key (specify SSEKMSKeyId in the put_object call for key rotation control and cross-account access policies). DynamoDB encryption at rest is enabled by default. All API calls travel over TLS.

Bias and equity. Utilization-based segmentation can encode existing access disparities. Members in underserved areas may appear "healthy" (low utilization) when they actually cannot access care. Production systems cross-reference segments with social determinant data and flag populations where low utilization might indicate access barriers rather than good health. This is not just an ethical concern; it's a clinical accuracy concern.

Access control. The DynamoDB table stores member IDs alongside utilization metrics, which constitutes PHI. In production, scope IAM policies to the minimum set of roles that need segment lookups. Consider using opaque surrogate keys rather than direct member identifiers if downstream consumers don't need the original ID. Classify the table as PHI-containing in your data catalog and apply appropriate retention policies.

Testing. There are no tests here. A production pipeline has unit tests for feature engineering (does clipping work correctly at boundaries?), integration tests for the full pipeline with known synthetic data (do you get the expected number of clusters?), and regression tests that verify segment stability across code changes.


Part of the Healthcare AI/ML Cookbook. See Recipe 6.2: Utilization Pattern Segmentation for the full architectural walkthrough, pseudocode, and honest take on where this gets hard.