Bionic Pancreas Analytics Pipeline

Introduction

30 Second Summary

Millions of people with diabetes rely on automated insulin delivery systems that make thousands of micro-decisions per day. Behind those devices, data teams analyze telemetry, clinical outcomes, and support records to catch safety signals and improve patient outcomes.

In this project, you will build a complete analytics pipeline for a simulated bionic pancreas company. You'll generate synthetic data, construct an ETL pipeline, perform safety signal detection and clinical outcomes analysis, and deliver insights through interactive dashboards.

What You'll Build

A portfolio-ready analytics system that fuses device telemetry, clinical outcomes, and CRM data into an analysis-ready warehouse, with statistical analysis and Tableau dashboards for three stakeholder audiences (safety, clinical, and commercial).

By the end of this project, you'll have:

• A synthetic data generator and ETL pipeline loading device, clinical, and CRM data into a SQLite analytical warehouse.
• Safety signal detection with Poisson confidence intervals and control charts, plus clinical outcomes analysis with pre/post HbA1c comparisons and time-in-range metrics.
• An interactive Tableau Public dashboard and stakeholder slide deck demonstrating integrated insights, with a Prefect-orchestrated pipeline for reproducible execution.
• Secret Mission: A regional early-warning composite risk score combining device alerts, support ticket velocity, and glycemic outcomes into a heatmap dashboard.

Are there any prerequisites?

This project assumes intermediate Python proficiency (comfortable with pandas, functions, and classes), basic SQL knowledge (SELECT, JOIN, GROUP BY), and familiarity with statistics (means, confidence intervals, hypothesis tests).

Scaffold the Repository and Data Dictionary

Your Python environment and all core packages are installed and ready. Now it's time to build the foundation that makes your project navigable and credible: a well-organized repository structure with clear documentation.

A hiring manager reviewing your portfolio will look at the repo layout and documentation before reading a single line of code. The data dictionary and README give that context upfront, proving you think about data governance and communication from day one.

In this step, get ready to:

• Initialize a Git repository and create the full directory structure.
• Write a data dictionary documenting all 4 synthetic datasets.
• Draft a README with project summary, methods, privacy language, and navigation.

Initialize the repository and create directories

A clean folder structure signals professionalism. Each directory has a dedicated purpose: raw data in, processed data out, source code separate from notebooks, and documentation alongside the code it describes.

• Initialize a Git repository and create all project directories by running these commands:

git init
mkdir -p data/raw data/processed src notebooks tests dashboards docs
touch README.md docs/data_dictionary.md

What do these commands do?

• git init creates a hidden .git folder that tracks all changes to your project.
• mkdir -p creates nested directories in one go. The -p flag means "create parent directories as needed."
• touch creates empty files that you'll fill with content in the next substeps.

• Confirm the directories and files exist by running:

ls

You should see README.md, data, dashboards, docs, notebooks, src, and tests listed in the output.

Don't see all the directories?

Make sure you ran the mkdir -p command with the full list of directory paths separated by spaces. If you see a "permission denied" error, check that you're in the correct project folder where your virtual environment lives.

Help me figure out why my directories weren't created.

Write the data dictionary

A data dictionary is the single source of truth for what each field means, where it comes from, and what values are valid. Anyone reviewing your project can look here to understand the data without digging through code.

• Open docs/data_dictionary.md in your text editor and paste the following content:

# Data Dictionary

## device_telemetry
| Field | Type | Range | Notes |
|---|---|---|---|
| user_id | str | — | Unique user identifier |
| timestamp | datetime | device_start_date to +180d | 5-minute intervals |
| glucose_mg_dl | float | 20-600 | CGM reading |
| insulin_delivered_u | float | 0-20 | Units delivered per interval |
| loop_mode | str | open/closed | AID mode |
| alert_type | str | none/low/high/urgent_low | Alert fired |
| battery_pct | int | 0-100 | Device battery |

## clinical_outcomes
| Field | Type | Range | Notes |
|---|---|---|---|
| user_id | str | — | |
| visit_date | date | — | Quarterly visit |
| hba1c_pct | float | 5.0-14.0 | Lab value |
| severe_hypo_event | int | 0/1 | Binary flag |
| dka_event | int | 0/1 | Binary flag, less than 1% rate |

## crm_data
| Field | Type | Range | Notes |
|---|---|---|---|
| user_id | str | — | |
| onboarding_date | date | — | |
| first_closed_loop_date | date | — | |
| account_status | str | active/inactive/churned | |
| region | str | — | Geographic region |
| sales_rep | str | — | |
| ticket_id | str | — | Support ticket |
| ticket_severity | str | low/medium/high | |
| ticket_resolution_days | float | — | |

## demographics_enrollment
| Field | Type | Range | Notes |
|---|---|---|---|
| user_id | str | — | |
| age | int | 18-75 | |
| diabetes_duration_yrs | float | 1-50 | |
| device_start_date | date | — | Staggered over 6 months |
| sex | str | M/F/Other | |

What does each dataset represent?

• device_telemetry holds continuous glucose monitor readings and insulin delivery data at 5-minute intervals. This is the highest-volume table.
• clinical_outcomes captures quarterly lab visits including HbA1c (a 3-month blood sugar average) and rare safety events like severe hypoglycemia and diabetic ketoacidosis (DKA).
• crm_data tracks the commercial relationship: onboarding progress, support tickets, and account status.
• demographics_enrollment links user identity to their age, diabetes duration, and the date they started using the device.

• Save docs/data_dictionary.md.
• Confirm the file has content by running:

cat docs/data_dictionary.md

You should see the full markdown table output in your terminal, starting with # Data Dictionary and ending with the demographics_enrollment table.

Seeing an empty file?

If cat shows nothing, you may have forgotten to save after pasting. Open the file again, paste the content, and save with your editor's save shortcut.

Why is my data dictionary file empty when I try to cat it?

✔️ Awesome, I've got everything!

Great. Make sure the file is saved before moving on.

ⓧ I'd like to double check the full code

# Data Dictionary

## device_telemetry
| Field | Type | Range | Notes |
|---|---|---|---|
| user_id | str | — | Unique user identifier |
| timestamp | datetime | device_start_date to +180d | 5-minute intervals |
| glucose_mg_dl | float | 20-600 | CGM reading |
| insulin_delivered_u | float | 0-20 | Units delivered per interval |
| loop_mode | str | open/closed | AID mode |
| alert_type | str | none/low/high/urgent_low | Alert fired |
| battery_pct | int | 0-100 | Device battery |

## clinical_outcomes
| Field | Type | Range | Notes |
|---|---|---|---|
| user_id | str | — | |
| visit_date | date | — | Quarterly visit |
| hba1c_pct | float | 5.0-14.0 | Lab value |
| severe_hypo_event | int | 0/1 | Binary flag |
| dka_event | int | 0/1 | Binary flag, less than 1% rate |

## crm_data
| Field | Type | Range | Notes |
|---|---|---|---|
| user_id | str | — | |
| onboarding_date | date | — | |
| first_closed_loop_date | date | — | |
| account_status | str | active/inactive/churned | |
| region | str | — | Geographic region |
| sales_rep | str | — | |
| ticket_id | str | — | Support ticket |
| ticket_severity | str | low/medium/high | |
| ticket_resolution_days | float | — | |

## demographics_enrollment
| Field | Type | Range | Notes |
|---|---|---|---|
| user_id | str | — | |
| age | int | 18-75 | |
| diabetes_duration_yrs | float | 1-50 | |
| device_start_date | date | — | Staggered over 6 months |
| sex | str | M/F/Other | |

Draft the README

The README is the front door of your portfolio project. It communicates what you built, how you built it, and where to find everything. It also includes HIPAA-aware language clarifying that all data is synthetic.

• Open README.md in your text editor and paste the following content:

# Integrated Clinical & Commercial Analytics for a Bionic Pancreas

## Project Summary
End-to-end analytics pipeline simulating the work of a data analyst at an automated insulin delivery (AID) device company. All data is fully synthetic and does not represent any real patient information (HIPAA-aware: no PHI).

## Methods Overview
- Synthetic data generation (500 users, ~180 days telemetry)
- ETL pipeline (Python + SQLite)
- Safety signal detection (Poisson CIs, control charts)
- Clinical outcomes analysis (pre/post HbA1c, time-in-range)
- CRM/commercial cohort analysis
- Pipeline orchestration with Prefect
- Interactive Tableau Public dashboard

## Privacy Considerations
All data is synthetically generated. No real patient data is used or represented.

## Repository Navigation
- `data/raw/` — synthetic CSV exports
- `data/processed/` — SQLite analytical warehouse
- `src/` — ETL and pipeline scripts
- `notebooks/` — analytical Jupyter notebooks
- `tests/` — pytest unit tests
- `dashboards/` — Tableau screenshots and exports
- `docs/` — data dictionary, validation report, demo script

Why include HIPAA-aware language?

Even though this project uses fully synthetic data, a reviewer at a healthcare or medical device company will look for evidence that you understand data privacy regulations. Stating "no PHI" (Protected Health Information) upfront shows awareness of the regulatory environment you'd work in.

This single sentence differentiates your portfolio from candidates who skip governance considerations entirely.

• Save README.md.
• Confirm the file has content by running:

cat README.md

You should see the full project summary, methods list, privacy section, and repository navigation printed in your terminal.

Seeing raw markdown instead of formatted text?

That's expected. The cat command shows raw file contents. The markdown formatting will render properly on GitHub once you push the repository later.

Help me understand why my README looks different in the terminal vs on GitHub.

✔️ Awesome, I've got everything!

Great. Double check that you've saved the file before continuing.

ⓧ I'd like to double check the full code

# Integrated Clinical & Commercial Analytics for a Bionic Pancreas

## Project Summary
End-to-end analytics pipeline simulating the work of a data analyst at an automated insulin delivery (AID) device company. All data is fully synthetic and does not represent any real patient information (HIPAA-aware: no PHI).

## Methods Overview
- Synthetic data generation (500 users, ~180 days telemetry)
- ETL pipeline (Python + SQLite)
- Safety signal detection (Poisson CIs, control charts)
- Clinical outcomes analysis (pre/post HbA1c, time-in-range)
- CRM/commercial cohort analysis
- Pipeline orchestration with Prefect
- Interactive Tableau Public dashboard

## Privacy Considerations
All data is synthetically generated. No real patient data is used or represented.

## Repository Navigation
- `data/raw/` — synthetic CSV exports
- `data/processed/` — SQLite analytical warehouse
- `src/` — ETL and pipeline scripts
- `notebooks/` — analytical Jupyter notebooks
- `tests/` — pytest unit tests
- `dashboards/` — Tableau screenshots and exports
- `docs/` — data dictionary, validation report, demo script

Your repository structure and documentation are in place. Next up, you'll build the synthetic data generator that populates these directories with realistic device telemetry, clinical outcomes, and CRM records for 500 simulated users.

Build the Synthetic Data Generator

Your project scaffolding and data dictionary are in place. Now you need the actual data to flow through the pipeline.

In this step, you will build a Python script that creates four interrelated synthetic datasets mimicking real automated insulin delivery (AID) device data. The script produces demographics, continuous glucose monitor readings, clinical visit outcomes, and customer relationship records for 500 simulated users.

In this step, get ready to:

• Build demographics and device telemetry generation functions.
• Create clinical outcomes and CRM data generators.
• Export all four datasets to CSV with a reproducible seed parameter.

Create the demographics and device telemetry generator

The generator script needs to produce realistic device data for 500 users. Each user gets a profile (age, diabetes duration, device start date) and approximately 180 days of 5-minute glucose readings. That is around 26 million telemetry rows total.

• Create a new file called src/generate_data.py in your project.
• Add the imports, configuration constants, argument parser, and demographics generation function by copying the following code:

"""Synthetic Data Generator for Bionic Pancreas Analytics Pipeline.
Generates demographics, device telemetry, clinical outcomes, and CRM data
for 500 simulated automated insulin delivery (AID) users.
"""

import argparse
import os
from datetime import datetime, timedelta

import numpy as np
import pandas as pd

# --- Configuration ---
NUM_USERS = 500
TELEMETRY_DAYS = 180
READINGS_PER_DAY = 288  # 5-minute intervals (24 * 60 / 5)
REGIONS = ["Northeast", "Southeast", "Midwest", "West", "Southwest"]
SALES_REPS = [f"REP_{i:03d}" for i in range(1, 21)]
ENROLLMENT_START = datetime(2023, 1, 1)
ENROLLMENT_WINDOW_DAYS = 180  # Stagger start dates over 6 months
OUTPUT_DIR = os.path.join(os.path.dirname(__file__), "..", "data", "raw")


def parse_args():
    """Parse command-line arguments."""
    parser = argparse.ArgumentParser(
        description="Generate synthetic bionic pancreas data."
    )
    parser.add_argument(
        "--seed", type=int, default=42,
        help="Random seed for reproducibility"
    )
    return parser.parse_args()


def generate_demographics(rng):
    """Generate 500 synthetic user demographics and enrollment records."""
    records = []
    for i in range(NUM_USERS):
        user_id = f"USR_{i:04d}"
        age = int(rng.integers(18, 76))
        sex = rng.choice(["M", "F", "Other"], p=[0.48, 0.48, 0.04])
        diabetes_duration_yrs = round(float(rng.uniform(1, 50)), 1)
        # Stagger device start dates across the enrollment window
        days_offset = int(rng.integers(0, ENROLLMENT_WINDOW_DAYS))
        device_start_date = (
            ENROLLMENT_START + timedelta(days=days_offset)
        ).strftime("%Y-%m-%d")
        region = rng.choice(REGIONS)
        records.append({
            "user_id": user_id,
            "age": age,
            "sex": sex,
            "diabetes_duration_yrs": diabetes_duration_yrs,
            "device_start_date": device_start_date,
            "region": region,
        })
    return pd.DataFrame(records)

What does this code do?

The configuration constants define the shape of the synthetic data: 500 users, 180 days of telemetry, 288 readings per day (one every 5 minutes), and 5 geographic regions.

The parse_args() function accepts a --seed command-line argument (defaulting to 42) that makes the entire output deterministic. Running the script twice with the same seed produces identical datasets.

The generate_demographics() function loops over 500 users, assigning each a random age (18-75), sex, diabetes duration (1-50 years), region, and a device start date staggered across a 6-month enrollment window starting January 2023.

Why use numpy.random.default_rng instead of random.seed?

NumPy's default_rng() creates a standalone Generator using the PCG64 algorithm. Unlike the legacy np.random.seed() approach, each Generator instance manages its own state. This means you can pass the rng object into functions without worrying about global state collisions.

It also produces higher-quality random numbers and is the recommended approach in modern NumPy code.

Now add the telemetry generation functions directly below the demographics function. Because telemetry is massive (approximately 26 million rows), the script writes one user at a time to avoid exhausting memory.

• Add the per-user telemetry helper and the batch writer by pasting this code below generate_demographics():

def generate_telemetry_for_user(rng, user_id, device_start_date):
    """Generate ~180 days of 5-minute glucose readings for one user."""
    start = datetime.strptime(device_start_date, "%Y-%m-%d")
    # User-specific baseline glucose and insulin parameters
    base_glucose = float(rng.uniform(140, 220))
    insulin_base = float(rng.uniform(0.5, 2.0))
    # Day when user transitions from open-loop to closed-loop
    loop_transition_day = int(rng.integers(7, 45))

    records = []
    for day in range(TELEMETRY_DAYS):
        in_closed_loop = day >= loop_transition_day
        # Glucose center improves over time once in closed-loop mode
        improvement = (
            min(day - loop_transition_day, 90) * 0.3
            if in_closed_loop else 0
        )
        day_center = max(100, base_glucose - improvement)
        current_date = start + timedelta(days=day)
        # Battery resets each day (simulating daily charge)
        battery_start = int(rng.integers(90, 101))

        for reading in range(READINGS_PER_DAY):
            timestamp = current_date + timedelta(minutes=reading * 5)
            # Normally distributed glucose with clipping to physiological range
            glucose = float(np.clip(rng.normal(day_center, 30), 40, 400))
            insulin = round(
                float(np.clip(rng.normal(insulin_base, 0.4), 0, 20)), 2
            )
            loop_mode = "closed" if in_closed_loop else "open"
            # Alert logic based on glucose thresholds
            if glucose < 54:
                alert_type = "urgent_low"
            elif glucose < 70:
                alert_type = "low"
            elif glucose > 250:
                alert_type = "high"
            else:
                alert_type = "none"
            # Battery drains through the day, resets next day
            battery_pct = max(
                5,
                battery_start - (reading * 85 // READINGS_PER_DAY)
                + int(rng.integers(-2, 3))
            )
            records.append({
                "user_id": user_id,
                "timestamp": timestamp.strftime("%Y-%m-%d %H:%M:%S"),
                "glucose_mg_dl": round(glucose, 1),
                "insulin_delivered_u": insulin,
                "loop_mode": loop_mode,
                "alert_type": alert_type,
                "battery_pct": battery_pct,
            })
    return records


def generate_telemetry(rng, demographics_df, output_path):
    """Generate device telemetry for all users, writing per-user to manage memory."""
    header_written = False
    for idx, user in demographics_df.iterrows():
        records = generate_telemetry_for_user(
            rng, user["user_id"], user["device_start_date"]
        )
        chunk_df = pd.DataFrame(records)
        chunk_df.to_csv(
            output_path, mode="a", header=not header_written, index=False
        )
        header_written = True
        if (idx + 1) % 50 == 0:
            print(f"  Telemetry: {idx + 1}/{NUM_USERS} users complete")
    print(f"  Telemetry: {NUM_USERS}/{NUM_USERS} users complete")

What does this code do?

The generate_telemetry_for_user() function creates 180 days of 5-minute readings for a single user. Each user starts in open-loop mode and transitions to closed-loop mode after 7-45 days. Once in closed-loop, glucose levels gradually improve (the day_center drops over time), simulating realistic AID device behavior.

Alert types follow clinical thresholds: below 54 mg/dL triggers urgent_low, below 70 triggers low, and above 250 triggers high.

The generate_telemetry() wrapper writes each user's data to the CSV file in append mode. This keeps memory usage manageable despite the total dataset being approximately 26 million rows.

• Save src/generate_data.py.

You will not run the script yet. The next substep adds the remaining two generators before the main execution block.

Seeing an import error for numpy or pandas?

Make sure your virtual environment is active before opening the file. If you deactivated it since Step 1, reactivate it from your project directory.

help me fix my import error for numpy or pandas in generate_data.py

Generate clinical outcomes and CRM records

Clinical outcomes track quarterly HbA1c lab visits and rare adverse events. CRM data captures the commercial side: onboarding timelines, support tickets, and account retention. Both datasets link back to the same 500 users through user_id.

• Add the clinical outcomes generator below the telemetry functions in src/generate_data.py:

def generate_clinical_outcomes(rng, demographics_df):
    """Generate quarterly clinical visit data with improvement trajectory."""
    records = []
    for _, user in demographics_df.iterrows():
        user_id = user["user_id"]
        start = datetime.strptime(user["device_start_date"], "%Y-%m-%d")
        # Baseline HbA1c correlates loosely with diabetes duration
        baseline_hba1c = float(rng.uniform(7.5, 10.5))

        for quarter in range(4):
            visit_date = (
                start + timedelta(days=90 * (quarter + 1))
            ).strftime("%Y-%m-%d")
            # HbA1c improves over quarters (typical AID improvement)
            improvement = quarter * float(rng.uniform(0.2, 0.5))
            hba1c = round(max(5.5, baseline_hba1c - improvement), 1)
            # Severe hypoglycemia at ~3.5% rate per visit
            severe_hypo_event = 1 if rng.random() < 0.035 else 0
            # DKA at ~0.8% rate per visit
            dka_event = 1 if rng.random() < 0.008 else 0
            records.append({
                "user_id": user_id,
                "visit_date": visit_date,
                "hba1c_pct": hba1c,
                "severe_hypo_event": severe_hypo_event,
                "dka_event": dka_event,
            })
    return pd.DataFrame(records)

What does this code do?

Each of the 500 users gets 4 quarterly visits (2000 total rows). The baseline HbA1c starts between 7.5% and 10.5%, then drops by 0.2-0.5% each quarter as the AID device helps the user manage glucose better.

Severe hypoglycemia events fire at approximately 3.5% probability per visit, and DKA events at approximately 0.8%. These rates match the ranges specified in your data dictionary (2-5% for severe hypo, less than 1% for DKA).

• Add the CRM data generator directly below the clinical outcomes function:

def generate_crm_data(rng, demographics_df):
    """Generate CRM records including onboarding, account status, and support tickets."""
    records = []
    for _, user in demographics_df.iterrows():
        user_id = user["user_id"]
        region = user["region"]
        start = datetime.strptime(user["device_start_date"], "%Y-%m-%d")
        # Onboarding happens within first 2 weeks
        onboarding_days = int(rng.integers(1, 15))
        onboarding_date = (
            start + timedelta(days=onboarding_days)
        ).strftime("%Y-%m-%d")
        # First closed-loop session 1-6 weeks after device start
        first_cl_days = int(rng.integers(7, 43))
        first_closed_loop_date = (
            start + timedelta(days=first_cl_days)
        ).strftime("%Y-%m-%d")
        # Account status weighted toward active
        account_status = rng.choice(
            ["active", "inactive", "churned"], p=[0.75, 0.15, 0.10]
        )
        sales_rep = rng.choice(SALES_REPS)
        # Generate 1-7 support tickets per user
        num_tickets = int(rng.integers(1, 8))
        for t in range(num_tickets):
            ticket_id = f"TKT_{user_id}_{t:02d}"
            ticket_severity = rng.choice(
                ["low", "medium", "high"], p=[0.5, 0.35, 0.15]
            )
            # Resolution time depends on severity
            if ticket_severity == "low":
                resolution_days = round(float(rng.uniform(0.5, 5)), 1)
            elif ticket_severity == "medium":
                resolution_days = round(float(rng.uniform(2, 14)), 1)
            else:
                resolution_days = round(float(rng.uniform(5, 30)), 1)
            records.append({
                "user_id": user_id,
                "onboarding_date": onboarding_date,
                "first_closed_loop_date": first_closed_loop_date,
                "account_status": account_status,
                "region": region,
                "sales_rep": sales_rep,
                "ticket_id": ticket_id,
                "ticket_severity": ticket_severity,
                "ticket_resolution_days": resolution_days,
            })
    return pd.DataFrame(records)

What does this code do?

Each user generates 1-7 support tickets (producing approximately 1750 total rows). Ticket severity follows a realistic distribution: 50% low, 35% medium, 15% high. Resolution time scales with severity.

The CRM record also captures onboarding speed (days until initial setup) and first closed-loop date (days until the user enables the automated algorithm). These metrics become critical in Step 7 when you analyze whether slow onboarding predicts churn.

• Save src/generate_data.py.

Seeing an IndentationError?

Make sure each new function starts at column 0 (no indentation before def). The functions should be at the module level, not nested inside another function or class.

help me fix my IndentationError in generate_data.py

Export all datasets and run the generator

The final piece is the main execution block that ties all four generators together, handles the seed, and writes the CSV files to data/raw/.

• Add the main block at the very bottom of src/generate_data.py:

if __name__ == "__main__":
    args = parse_args()
    rng = np.random.default_rng(args.seed)
    print(f"Generating synthetic data with seed={args.seed}...")

    # Ensure output directory exists
    os.makedirs(OUTPUT_DIR, exist_ok=True)

    # Step 1: Demographics
    print("Generating demographics...")
    demographics_df = generate_demographics(rng)
    demographics_df.to_csv(
        os.path.join(OUTPUT_DIR, "demographics_enrollment.csv"), index=False
    )
    print(f"  Saved {len(demographics_df)} user records.")

    # Step 2: Device telemetry (written per-user to manage memory)
    print("Generating device telemetry (this may take several minutes)...")
    telemetry_path = os.path.join(OUTPUT_DIR, "device_telemetry.csv")
    # Remove existing file to avoid appending to stale data
    if os.path.exists(telemetry_path):
        os.remove(telemetry_path)
    generate_telemetry(rng, demographics_df, telemetry_path)

    # Step 3: Clinical outcomes
    print("Generating clinical outcomes...")
    clinical_df = generate_clinical_outcomes(rng, demographics_df)
    clinical_df.to_csv(
        os.path.join(OUTPUT_DIR, "clinical_outcomes.csv"), index=False
    )
    print(f"  Saved {len(clinical_df)} visit records.")

    # Step 4: CRM data
    print("Generating CRM data...")
    crm_df = generate_crm_data(rng, demographics_df)
    crm_df.to_csv(
        os.path.join(OUTPUT_DIR, "crm_data.csv"), index=False
    )
    print(f"  Saved {len(crm_df)} CRM records.")

    print("Done! All datasets saved to data/raw/")

What does this code do?

The if __name__ == "__main__" block initializes a single rng Generator seeded from the command-line argument. All four generators share this same rng instance, so the entire output is deterministic for a given seed value.

The telemetry writer removes any existing CSV before starting to prevent appending duplicate rows if you re-run the script. Demographics, clinical outcomes, and CRM data are small enough to hold in memory and write in one pass.

• Save src/generate_data.py.
• Run the generator by executing this command from your project root directory:

python src/generate_data.py --seed 42

What should I see?

The script prints progress messages as it generates each dataset. Telemetry generation takes several minutes because it is producing approximately 26 million rows. You should see output similar to:

Generating synthetic data with seed=42...
Generating demographics...
  Saved 500 user records.
Generating device telemetry (this may take several minutes)...
  Telemetry: 50/500 users complete
  Telemetry: 100/500 users complete
  ...
  Telemetry: 500/500 users complete
Generating clinical outcomes...
  Saved 2000 visit records.
Generating CRM data...
  Saved 1750 CRM records.
Done! All datasets saved to data/raw/

• Confirm the files were created by running:

ls data/raw/

You should see clinical_outcomes.csv, crm_data.csv, demographics_enrollment.csv, and device_telemetry.csv listed.

Script fails with FileNotFoundError?

Make sure you are running the command from your project root directory (the folder containing src/ and data/). The script uses relative paths from its own location to find data/raw/.

help me fix FileNotFoundError when running generate_data.py

Verify reproducibility

To confirm the seed parameter works, run the script again with the same --seed 42 flag. The output files should be byte-for-byte identical. You can verify with md5sum data/raw/demographics_enrollment.csv on macOS/Linux. Both runs should produce the same hash.

✔️ Awesome, I've got everything!

Great. Double check you have saved src/generate_data.py and that all four CSVs exist in data/raw/.

ⓧ I'd like to double check the full code

Here is the complete src/generate_data.py file for comparison:

"""Synthetic Data Generator for Bionic Pancreas Analytics Pipeline.
Generates demographics, device telemetry, clinical outcomes, and CRM data
for 500 simulated automated insulin delivery (AID) users.
"""

import argparse
import os
from datetime import datetime, timedelta

import numpy as np
import pandas as pd

# --- Configuration ---
NUM_USERS = 500
TELEMETRY_DAYS = 180
READINGS_PER_DAY = 288  # 5-minute intervals (24 * 60 / 5)
REGIONS = ["Northeast", "Southeast", "Midwest", "West", "Southwest"]
SALES_REPS = [f"REP_{i:03d}" for i in range(1, 21)]
ENROLLMENT_START = datetime(2023, 1, 1)
ENROLLMENT_WINDOW_DAYS = 180  # Stagger start dates over 6 months
OUTPUT_DIR = os.path.join(os.path.dirname(__file__), "..", "data", "raw")


def parse_args():
    """Parse command-line arguments."""
    parser = argparse.ArgumentParser(
        description="Generate synthetic bionic pancreas data."
    )
    parser.add_argument(
        "--seed", type=int, default=42,
        help="Random seed for reproducibility"
    )
    return parser.parse_args()


def generate_demographics(rng):
    """Generate 500 synthetic user demographics and enrollment records."""
    records = []
    for i in range(NUM_USERS):
        user_id = f"USR_{i:04d}"
        age = int(rng.integers(18, 76))
        sex = rng.choice(["M", "F", "Other"], p=[0.48, 0.48, 0.04])
        diabetes_duration_yrs = round(float(rng.uniform(1, 50)), 1)
        # Stagger device start dates across the enrollment window
        days_offset = int(rng.integers(0, ENROLLMENT_WINDOW_DAYS))
        device_start_date = (
            ENROLLMENT_START + timedelta(days=days_offset)
        ).strftime("%Y-%m-%d")
        region = rng.choice(REGIONS)
        records.append({
            "user_id": user_id,
            "age": age,
            "sex": sex,
            "diabetes_duration_yrs": diabetes_duration_yrs,
            "device_start_date": device_start_date,
            "region": region,
        })
    return pd.DataFrame(records)


def generate_telemetry_for_user(rng, user_id, device_start_date):
    """Generate ~180 days of 5-minute glucose readings for one user."""
    start = datetime.strptime(device_start_date, "%Y-%m-%d")
    # User-specific baseline glucose and insulin parameters
    base_glucose = float(rng.uniform(140, 220))
    insulin_base = float(rng.uniform(0.5, 2.0))
    # Day when user transitions from open-loop to closed-loop
    loop_transition_day = int(rng.integers(7, 45))

    records = []
    for day in range(TELEMETRY_DAYS):
        in_closed_loop = day >= loop_transition_day
        # Glucose center improves over time once in closed-loop mode
        improvement = (
            min(day - loop_transition_day, 90) * 0.3
            if in_closed_loop else 0
        )
        day_center = max(100, base_glucose - improvement)
        current_date = start + timedelta(days=day)
        # Battery resets each day (simulating daily charge)
        battery_start = int(rng.integers(90, 101))

        for reading in range(READINGS_PER_DAY):
            timestamp = current_date + timedelta(minutes=reading * 5)
            # Normally distributed glucose with clipping to physiological range
            glucose = float(np.clip(rng.normal(day_center, 30), 40, 400))
            insulin = round(
                float(np.clip(rng.normal(insulin_base, 0.4), 0, 20)), 2
            )
            loop_mode = "closed" if in_closed_loop else "open"
            # Alert logic based on glucose thresholds
            if glucose < 54:
                alert_type = "urgent_low"
            elif glucose < 70:
                alert_type = "low"
            elif glucose > 250:
                alert_type = "high"
            else:
                alert_type = "none"
            # Battery drains through the day, resets next day
            battery_pct = max(
                5,
                battery_start - (reading * 85 // READINGS_PER_DAY)
                + int(rng.integers(-2, 3))
            )
            records.append({
                "user_id": user_id,
                "timestamp": timestamp.strftime("%Y-%m-%d %H:%M:%S"),
                "glucose_mg_dl": round(glucose, 1),
                "insulin_delivered_u": insulin,
                "loop_mode": loop_mode,
                "alert_type": alert_type,
                "battery_pct": battery_pct,
            })
    return records


def generate_telemetry(rng, demographics_df, output_path):
    """Generate device telemetry for all users, writing per-user to manage memory."""
    header_written = False
    for idx, user in demographics_df.iterrows():
        records = generate_telemetry_for_user(
            rng, user["user_id"], user["device_start_date"]
        )
        chunk_df = pd.DataFrame(records)
        chunk_df.to_csv(
            output_path, mode="a", header=not header_written, index=False
        )
        header_written = True
        if (idx + 1) % 50 == 0:
            print(f"  Telemetry: {idx + 1}/{NUM_USERS} users complete")
    print(f"  Telemetry: {NUM_USERS}/{NUM_USERS} users complete")


def generate_clinical_outcomes(rng, demographics_df):
    """Generate quarterly clinical visit data with improvement trajectory."""
    records = []
    for _, user in demographics_df.iterrows():
        user_id = user["user_id"]
        start = datetime.strptime(user["device_start_date"], "%Y-%m-%d")
        # Baseline HbA1c correlates loosely with diabetes duration
        baseline_hba1c = float(rng.uniform(7.5, 10.5))

        for quarter in range(4):
            visit_date = (
                start + timedelta(days=90 * (quarter + 1))
            ).strftime("%Y-%m-%d")
            # HbA1c improves over quarters (typical AID improvement)
            improvement = quarter * float(rng.uniform(0.2, 0.5))
            hba1c = round(max(5.5, baseline_hba1c - improvement), 1)
            # Severe hypoglycemia at ~3.5% rate per visit
            severe_hypo_event = 1 if rng.random() < 0.035 else 0
            # DKA at ~0.8% rate per visit
            dka_event = 1 if rng.random() < 0.008 else 0
            records.append({
                "user_id": user_id,
                "visit_date": visit_date,
                "hba1c_pct": hba1c,
                "severe_hypo_event": severe_hypo_event,
                "dka_event": dka_event,
            })
    return pd.DataFrame(records)


def generate_crm_data(rng, demographics_df):
    """Generate CRM records including onboarding, account status, and support tickets."""
    records = []
    for _, user in demographics_df.iterrows():
        user_id = user["user_id"]
        region = user["region"]
        start = datetime.strptime(user["device_start_date"], "%Y-%m-%d")
        # Onboarding happens within first 2 weeks
        onboarding_days = int(rng.integers(1, 15))
        onboarding_date = (
            start + timedelta(days=onboarding_days)
        ).strftime("%Y-%m-%d")
        # First closed-loop session 1-6 weeks after device start
        first_cl_days = int(rng.integers(7, 43))
        first_closed_loop_date = (
            start + timedelta(days=first_cl_days)
        ).strftime("%Y-%m-%d")
        # Account status weighted toward active
        account_status = rng.choice(
            ["active", "inactive", "churned"], p=[0.75, 0.15, 0.10]
        )
        sales_rep = rng.choice(SALES_REPS)
        # Generate 1-7 support tickets per user
        num_tickets = int(rng.integers(1, 8))
        for t in range(num_tickets):
            ticket_id = f"TKT_{user_id}_{t:02d}"
            ticket_severity = rng.choice(
                ["low", "medium", "high"], p=[0.5, 0.35, 0.15]
            )
            # Resolution time depends on severity
            if ticket_severity == "low":
                resolution_days = round(float(rng.uniform(0.5, 5)), 1)
            elif ticket_severity == "medium":
                resolution_days = round(float(rng.uniform(2, 14)), 1)
            else:
                resolution_days = round(float(rng.uniform(5, 30)), 1)
            records.append({
                "user_id": user_id,
                "onboarding_date": onboarding_date,
                "first_closed_loop_date": first_closed_loop_date,
                "account_status": account_status,
                "region": region,
                "sales_rep": sales_rep,
                "ticket_id": ticket_id,
                "ticket_severity": ticket_severity,
                "ticket_resolution_days": resolution_days,
            })
    return pd.DataFrame(records)


if __name__ == "__main__":
    args = parse_args()
    rng = np.random.default_rng(args.seed)
    print(f"Generating synthetic data with seed={args.seed}...")

    # Ensure output directory exists
    os.makedirs(OUTPUT_DIR, exist_ok=True)

    # Step 1: Demographics
    print("Generating demographics...")
    demographics_df = generate_demographics(rng)
    demographics_df.to_csv(
        os.path.join(OUTPUT_DIR, "demographics_enrollment.csv"), index=False
    )
    print(f"  Saved {len(demographics_df)} user records.")

    # Step 2: Device telemetry (written per-user to manage memory)
    print("Generating device telemetry (this may take several minutes)...")
    telemetry_path = os.path.join(OUTPUT_DIR, "device_telemetry.csv")
    # Remove existing file to avoid appending to stale data
    if os.path.exists(telemetry_path):
        os.remove(telemetry_path)
    generate_telemetry(rng, demographics_df, telemetry_path)

    # Step 3: Clinical outcomes
    print("Generating clinical outcomes...")
    clinical_df = generate_clinical_outcomes(rng, demographics_df)
    clinical_df.to_csv(
        os.path.join(OUTPUT_DIR, "clinical_outcomes.csv"), index=False
    )
    print(f"  Saved {len(clinical_df)} visit records.")

    # Step 4: CRM data
    print("Generating CRM data...")
    crm_df = generate_crm_data(rng, demographics_df)
    crm_df.to_csv(
        os.path.join(OUTPUT_DIR, "crm_data.csv"), index=False
    )
    print(f"  Saved {len(crm_df)} CRM records.")

    print("Done! All datasets saved to data/raw/")

Your synthetic data foundation is ready. Next up, you will build the ETL pipeline that cleans these raw CSVs and loads them into a SQLite analytical warehouse.

Build the ETL Pipeline

Your synthetic data generator is producing four realistic CSVs. Raw data alone is not analysis-ready.

In this step, you'll build an ETL pipeline that extracts data from those CSVs, transforms it (type casting, null removal, range validation), and loads it into a SQLite database. You'll also create a user_summary table that joins all four sources into one analysis-ready layer.

In this step, get ready to:

• Build data cleaning functions that validate and standardize each CSV.
• Create a user-level summary table joining device metrics, clinical outcomes, and CRM data.
• Write example SQL queries demonstrating cohort extraction and aggregation.

Write the data cleaning functions

Each of the four raw CSVs needs type casting, date parsing, and null handling before it can be loaded into the database. The telemetry data also needs range validation to exclude impossible glucose readings.

• Create a new file called src/etl.py by running:

touch src/etl.py

• Open src/etl.py in your editor and add the imports, path constants, and first two cleaning functions:

import pandas as pd
import sqlite3
from pathlib import Path

# Database and raw data paths
DB_PATH = 'data/processed/bionic_pancreas.db'
RAW_PATH = Path('data/raw')


def load_and_clean_demographics():
    """Load and clean demographics_enrollment.csv."""
    df = pd.read_csv(RAW_PATH / 'demographics_enrollment.csv')

    # Parse dates and validate types
    df['device_start_date'] = pd.to_datetime(df['device_start_date'])
    df['age'] = df['age'].astype(int)
    df['diabetes_duration_yrs'] = df['diabetes_duration_yrs'].astype(float)

    # Drop rows with null user_id
    df = df.dropna(subset=['user_id'])

    return df


def load_and_clean_telemetry():
    """Load and clean device_telemetry.csv."""
    df = pd.read_csv(RAW_PATH / 'device_telemetry.csv')

    # Parse timestamps
    df['timestamp'] = pd.to_datetime(df['timestamp'])

    # Validate glucose range (20-600 mg/dL)
    df = df[(df['glucose_mg_dl'] >= 20) & (df['glucose_mg_dl'] <= 600)]

    # Cast numeric columns
    df['insulin_delivered_u'] = df['insulin_delivered_u'].astype(float)
    df['battery_pct'] = df['battery_pct'].astype(int)

    # Drop rows with null user_id
    df = df.dropna(subset=['user_id'])

    return df

What does this code do?

DB_PATH and RAW_PATH centralize the file locations so every function can reference them without hardcoding paths.

load_and_clean_demographics() parses device_start_date into a proper datetime, casts age to integer, and removes any rows missing a user_id.

load_and_clean_telemetry() applies range validation that discards glucose readings outside 20-600 mg/dL. These bounds match physiological reality. Any reading outside this range is sensor noise or data corruption.

• Below the load_and_clean_telemetry() function, add the clinical and CRM cleaning functions:

def load_and_clean_clinical():
    """Load and clean clinical_outcomes.csv."""
    df = pd.read_csv(RAW_PATH / 'clinical_outcomes.csv')

    # Parse dates
    df['visit_date'] = pd.to_datetime(df['visit_date'])

    # Cast event flags as integers
    df['severe_hypo_event'] = df['severe_hypo_event'].astype(int)
    df['dka_event'] = df['dka_event'].astype(int)

    # Drop rows with null user_id
    df = df.dropna(subset=['user_id'])

    return df


def load_and_clean_crm():
    """Load and clean crm_data.csv."""
    df = pd.read_csv(RAW_PATH / 'crm_data.csv')

    # Parse dates
    df['onboarding_date'] = pd.to_datetime(df['onboarding_date'])
    df['first_closed_loop_date'] = pd.to_datetime(df['first_closed_loop_date'])

    # Cast resolution days as float, coerce errors to NaN
    df['ticket_resolution_days'] = pd.to_numeric(
        df['ticket_resolution_days'], errors='coerce'
    )

    # Drop rows with null user_id
    df = df.dropna(subset=['user_id'])

    return df

What does this code do?

load_and_clean_clinical() ensures severe_hypo_event and dka_event are proper integers (0 or 1) for downstream aggregation.

load_and_clean_crm() uses pd.to_numeric(..., errors='coerce') for ticket_resolution_days. This converts unparseable strings to NaN instead of crashing the pipeline.

You won't see a change yet. The next chunk adds the orchestration that actually runs these functions.

Create the user summary and run the pipeline

The analysis-ready layer collapses millions of telemetry rows into one row per user, joining demographics, device metrics, and CRM status into a single table.

• Below load_and_clean_crm(), add the create_user_summary() function:

def create_user_summary(conn):
    """Create user_summary table joining all sources."""
    query = """
    CREATE TABLE IF NOT EXISTS user_summary AS
    SELECT
        d.user_id,
        d.age,
        d.sex,
        d.diabetes_duration_yrs,
        d.device_start_date,
        d.region,
        t.mean_glucose,
        t.time_in_range_pct,
        t.pct_closed_loop,
        t.mean_insulin,
        c.account_status,
        c.onboarding_date,
        c.first_closed_loop_date,
        c.days_to_closed_loop
    FROM demographics_enrollment d
    LEFT JOIN (
        SELECT
            user_id,
            AVG(glucose_mg_dl) AS mean_glucose,
            100.0 * SUM(CASE WHEN glucose_mg_dl BETWEEN 70 AND 180 THEN 1 ELSE 0 END)
                / COUNT(*) AS time_in_range_pct,
            100.0 * SUM(CASE WHEN loop_mode = 'closed' THEN 1 ELSE 0 END)
                / COUNT(*) AS pct_closed_loop,
            AVG(insulin_delivered_u) AS mean_insulin
        FROM device_telemetry
        GROUP BY user_id
    ) t ON d.user_id = t.user_id
    LEFT JOIN (
        SELECT DISTINCT
            user_id,
            account_status,
            onboarding_date,
            first_closed_loop_date,
            JULIANDAY(first_closed_loop_date) - JULIANDAY(onboarding_date)
                AS days_to_closed_loop
        FROM crm_data
    ) c ON d.user_id = c.user_id
    """
    conn.execute(query)
    conn.commit()

What does this SQL query do?

The query creates a user_summary table with one row per user. It starts from demographics_enrollment as the base and LEFT JOINs two subqueries.

The telemetry subquery computes four device metrics: mean_glucose, time_in_range_pct (70-180 mg/dL), pct_closed_loop, and mean_insulin.

The CRM subquery calculates days_to_closed_loop using SQLite's JULIANDAY() function, which measures how many days each user took from onboarding to their first closed-loop session.

• Below create_user_summary(), add the run_etl() orchestrator and the main entry point:

def run_etl():
    """Run the full ETL pipeline."""
    # Ensure output directory exists
    Path('data/processed').mkdir(parents=True, exist_ok=True)

    # Connect to SQLite database
    conn = sqlite3.connect(DB_PATH)

    # Load, clean, and write each dataset
    print("Loading demographics...")
    demographics = load_and_clean_demographics()
    demographics.to_sql('demographics_enrollment', conn, if_exists='replace', index=False)

    print("Loading telemetry (this may take a moment)...")
    telemetry = load_and_clean_telemetry()
    telemetry.to_sql('device_telemetry', conn, if_exists='replace', index=False)

    print("Loading clinical outcomes...")
    clinical = load_and_clean_clinical()
    clinical.to_sql('clinical_outcomes', conn, if_exists='replace', index=False)

    print("Loading CRM data...")
    crm = load_and_clean_crm()
    crm.to_sql('crm_data', conn, if_exists='replace', index=False)

    # Create the analysis-ready summary table
    print("Creating user_summary table...")
    create_user_summary(conn)

    print("ETL complete. Database saved to:", DB_PATH)
    conn.close()


if __name__ == '__main__':
    run_etl()

What does this code do?

run_etl() orchestrates the full pipeline in sequence. It creates the output directory, opens a single SQLite connection, loads each cleaned DataFrame into its own table using to_sql(), then builds the summary table on top.

The if_exists='replace' parameter means running the script multiple times drops and recreates each table. This makes the pipeline idempotent.

• Save src/etl.py.
• Run the ETL pipeline by executing:

python src/etl.py

You should see progress messages printed for each stage, ending with the line ETL complete. Database saved to: data/processed/bionic_pancreas.db.

Why does telemetry take longer?

The telemetry CSV contains approximately 26 million rows (500 users x 180 days x 288 readings per day). Loading this into SQLite takes noticeably longer than the other tables. On most machines, expect 30-90 seconds for this step.

• Verify the database was created by running:

python -c "import sqlite3; conn = sqlite3.connect('data/processed/bionic_pancreas.db'); print([r[0] for r in conn.execute('SELECT name FROM sqlite_master WHERE type=\"table\"').fetchall()]); conn.close()"

You should see a list containing demographics_enrollment, device_telemetry, clinical_outcomes, crm_data, and user_summary.

Not seeing all five tables?

Make sure you ran the script from your project root directory (the folder containing src/ and data/). The paths are relative.

If you see a FileNotFoundError, confirm your raw CSVs exist in data/raw/ by running ls data/raw/.

help me debug my ETL pipeline errors

✔️ Awesome, I've got everything!

Great. Double check you've saved src/etl.py and that your database has all five tables.

ⓧ I'd like to double check the full code

Here is the complete src/etl.py file:

import pandas as pd
import sqlite3
from pathlib import Path

# Database and raw data paths
DB_PATH = 'data/processed/bionic_pancreas.db'
RAW_PATH = Path('data/raw')


def load_and_clean_demographics():
    """Load and clean demographics_enrollment.csv."""
    df = pd.read_csv(RAW_PATH / 'demographics_enrollment.csv')

    # Parse dates and validate types
    df['device_start_date'] = pd.to_datetime(df['device_start_date'])
    df['age'] = df['age'].astype(int)
    df['diabetes_duration_yrs'] = df['diabetes_duration_yrs'].astype(float)

    # Drop rows with null user_id
    df = df.dropna(subset=['user_id'])

    return df


def load_and_clean_telemetry():
    """Load and clean device_telemetry.csv."""
    df = pd.read_csv(RAW_PATH / 'device_telemetry.csv')

    # Parse timestamps
    df['timestamp'] = pd.to_datetime(df['timestamp'])

    # Validate glucose range (20-600 mg/dL)
    df = df[(df['glucose_mg_dl'] >= 20) & (df['glucose_mg_dl'] <= 600)]

    # Cast numeric columns
    df['insulin_delivered_u'] = df['insulin_delivered_u'].astype(float)
    df['battery_pct'] = df['battery_pct'].astype(int)

    # Drop rows with null user_id
    df = df.dropna(subset=['user_id'])

    return df


def load_and_clean_clinical():
    """Load and clean clinical_outcomes.csv."""
    df = pd.read_csv(RAW_PATH / 'clinical_outcomes.csv')

    # Parse dates
    df['visit_date'] = pd.to_datetime(df['visit_date'])

    # Cast event flags as integers
    df['severe_hypo_event'] = df['severe_hypo_event'].astype(int)
    df['dka_event'] = df['dka_event'].astype(int)

    # Drop rows with null user_id
    df = df.dropna(subset=['user_id'])

    return df


def load_and_clean_crm():
    """Load and clean crm_data.csv."""
    df = pd.read_csv(RAW_PATH / 'crm_data.csv')

    # Parse dates
    df['onboarding_date'] = pd.to_datetime(df['onboarding_date'])
    df['first_closed_loop_date'] = pd.to_datetime(df['first_closed_loop_date'])

    # Cast resolution days as float, coerce errors to NaN
    df['ticket_resolution_days'] = pd.to_numeric(
        df['ticket_resolution_days'], errors='coerce'
    )

    # Drop rows with null user_id
    df = df.dropna(subset=['user_id'])

    return df


def create_user_summary(conn):
    """Create user_summary table joining all sources."""
    query = """
    CREATE TABLE IF NOT EXISTS user_summary AS
    SELECT
        d.user_id,
        d.age,
        d.sex,
        d.diabetes_duration_yrs,
        d.device_start_date,
        d.region,
        t.mean_glucose,
        t.time_in_range_pct,
        t.pct_closed_loop,
        t.mean_insulin,
        c.account_status,
        c.onboarding_date,
        c.first_closed_loop_date,
        c.days_to_closed_loop
    FROM demographics_enrollment d
    LEFT JOIN (
        SELECT
            user_id,
            AVG(glucose_mg_dl) AS mean_glucose,
            100.0 * SUM(CASE WHEN glucose_mg_dl BETWEEN 70 AND 180 THEN 1 ELSE 0 END)
                / COUNT(*) AS time_in_range_pct,
            100.0 * SUM(CASE WHEN loop_mode = 'closed' THEN 1 ELSE 0 END)
                / COUNT(*) AS pct_closed_loop,
            AVG(insulin_delivered_u) AS mean_insulin
        FROM device_telemetry
        GROUP BY user_id
    ) t ON d.user_id = t.user_id
    LEFT JOIN (
        SELECT DISTINCT
            user_id,
            account_status,
            onboarding_date,
            first_closed_loop_date,
            JULIANDAY(first_closed_loop_date) - JULIANDAY(onboarding_date)
                AS days_to_closed_loop
        FROM crm_data
    ) c ON d.user_id = c.user_id
    """
    conn.execute(query)
    conn.commit()


def run_etl():
    """Run the full ETL pipeline."""
    # Ensure output directory exists
    Path('data/processed').mkdir(parents=True, exist_ok=True)

    # Connect to SQLite database
    conn = sqlite3.connect(DB_PATH)

    # Load, clean, and write each dataset
    print("Loading demographics...")
    demographics = load_and_clean_demographics()
    demographics.to_sql('demographics_enrollment', conn, if_exists='replace', index=False)

    print("Loading telemetry (this may take a moment)...")
    telemetry = load_and_clean_telemetry()
    telemetry.to_sql('device_telemetry', conn, if_exists='replace', index=False)

    print("Loading clinical outcomes...")
    clinical = load_and_clean_clinical()
    clinical.to_sql('clinical_outcomes', conn, if_exists='replace', index=False)

    print("Loading CRM data...")
    crm = load_and_clean_crm()
    crm.to_sql('crm_data', conn, if_exists='replace', index=False)

    # Create the analysis-ready summary table
    print("Creating user_summary table...")
    create_user_summary(conn)

    print("ETL complete. Database saved to:", DB_PATH)
    conn.close()


if __name__ == '__main__':
    run_etl()

Write example SQL queries

A reviewer cloning your repository should be able to run queries directly against the database without reading all your Python code. A standalone SQL file demonstrates your SQL proficiency and makes the analysis reproducible.

• Create notebooks/sql_queries.sql by running:

touch notebooks/sql_queries.sql

• Open notebooks/sql_queries.sql and add the following queries:

-- Example SQL Queries for Bionic Pancreas Analytics
-- Run these against data/processed/bionic_pancreas.db

-- 1. Preview the user summary table
SELECT * FROM user_summary LIMIT 10;

-- 2. Closed-loop vs open-loop cohort comparison
SELECT
    CASE WHEN pct_closed_loop >= 50 THEN 'closed_loop' ELSE 'open_loop' END AS cohort,
    COUNT(*) AS n_users,
    AVG(mean_glucose) AS avg_glucose,
    AVG(time_in_range_pct) AS avg_tir
FROM user_summary
GROUP BY cohort;

-- 3. Safety events by region
SELECT
    d.region,
    COUNT(*) AS total_visits,
    SUM(c.severe_hypo_event) AS severe_hypo_count,
    SUM(c.dka_event) AS dka_count,
    ROUND(100.0 * SUM(c.severe_hypo_event) / COUNT(*), 2) AS severe_hypo_rate_pct
FROM demographics_enrollment d
JOIN clinical_outcomes c ON d.user_id = c.user_id
GROUP BY d.region
ORDER BY severe_hypo_rate_pct DESC;

-- 4. Support ticket burden by severity
SELECT
    ticket_severity,
    COUNT(*) AS ticket_count,
    AVG(ticket_resolution_days) AS avg_resolution_days,
    MAX(ticket_resolution_days) AS max_resolution_days
FROM crm_data
GROUP BY ticket_severity
ORDER BY ticket_count DESC;

What do these queries demonstrate?

Query 1 is a simple preview. Query 2 splits users into cohorts based on whether they spent 50%+ of time in closed-loop mode, then compares their glucose outcomes.

Query 3 joins demographics_enrollment with clinical_outcomes to compute safety event rates per region. Query 4 aggregates support ticket resolution times by severity.

• Save notebooks/sql_queries.sql.
• Verify one of the queries works by running:

python -c "import sqlite3; conn = sqlite3.connect('data/processed/bionic_pancreas.db'); print(conn.execute('SELECT * FROM user_summary LIMIT 3').fetchall()); conn.close()"

You should see three rows of user-level summary data with columns for age, region, mean glucose, time-in-range percentage, and more.

Seeing empty results or an error?

If user_summary returns no rows, re-run python src/etl.py and check that the previous tables loaded correctly by querying SELECT COUNT(*) FROM demographics_enrollment. You should see 500.

help me debug my SQL queries against the bionic pancreas database

Your analytical warehouse is ready. Next up, you'll use these tables to detect safety signals with Poisson confidence intervals and control charts.

Perform Safety Signal Detection

Your ETL pipeline is running and the analytical warehouse is ready. Now it is time to answer the first critical question any medical device company must track: are there safety signals in the data?

Safety signal detection looks for unexpected spikes in adverse events (severe hypoglycemia, DKA, device alerts) across different patient cohorts. In this step, you will compute event rates with confidence intervals, build a control chart to catch week-over-week spikes, and estimate relative risk between device operating modes.

In this step, get ready to:

• Compute safety event rates per 1,000 device-days by cohort with Poisson confidence intervals.
• Build a c-chart for weekly safety event counts with 2-sigma and 3-sigma thresholds.
• Estimate relative risk ratios comparing closed-loop vs. open-loop users with bootstrapped 95% confidence intervals.

Compute safety event rates with Poisson confidence intervals

Before you can detect a signal, you need a baseline. Event rates express how often adverse events occur relative to total exposure time. Expressing them per 1,000 device-days makes rates comparable across cohorts of different sizes.

• Open Jupyter by running this command in your terminal (with your virtual environment active):

jupyter notebook

• Create a new Python 3 notebook.
• Save it as notebooks/01_safety_signal_detection.ipynb.

• Set up the notebook by pasting this imports cell:

import sqlite3
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import chi2

# Connect to the analytical warehouse
DB_PATH = 'data/processed/bionic_pancreas.db'
conn = sqlite3.connect(DB_PATH)

# Set plot style
sns.set_style('whitegrid')
print('Connected to database successfully.')

What does this code do?

This cell imports the core analysis libraries and opens a connection to your SQLite warehouse. The key import here is scipy.stats.chi2, which provides the chi-squared distribution quantiles needed for exact Poisson confidence intervals.

• Run the cell. You should see Connected to database successfully. printed below.

• In a new cell, document your inclusion/exclusion criteria and load the safety events data:

# =============================================================
# INCLUSION / EXCLUSION CRITERIA
# =============================================================
# Inclusion:
#   - All enrolled users with >= 30 days of device telemetry
#   - Safety events: severe_hypo_event, dka_event (clinical_outcomes),
#     and device alerts classified as adverse (alert_type != None)
# Exclusion:
#   - Users with < 30 days of telemetry (insufficient exposure)
#   - Battery/informational alerts (only clinically significant alerts count)
# =============================================================

# Load demographics and compute device-days per user
demographics = pd.read_sql('SELECT * FROM demographics_enrollment', conn)
telemetry_summary = pd.read_sql("""
    SELECT user_id,
           COUNT(DISTINCT DATE(timestamp)) as device_days,
           SUM(CASE WHEN alert_type IS NOT NULL THEN 1 ELSE 0 END) as device_alerts
    FROM device_telemetry
    GROUP BY user_id
""", conn)

# Load clinical safety events
clinical = pd.read_sql('SELECT * FROM clinical_outcomes', conn)

# Apply inclusion criteria: >= 30 device-days
eligible_users = telemetry_summary[telemetry_summary['device_days'] >= 30].copy()
print(f"Eligible users (>= 30 device-days): {len(eligible_users)} / {len(telemetry_summary)}")

# Aggregate safety events per user
user_safety = clinical.groupby('user_id').agg(
    severe_hypo_events=('severe_hypo_event', 'sum'),
    dka_events=('dka_event', 'sum')
).reset_index()

# Merge all sources
safety_df = eligible_users.merge(user_safety, on='user_id', how='left')
safety_df = safety_df.merge(demographics[['user_id', 'age', 'region']], on='user_id', how='left')

# Add device mode from user_summary
user_summary = pd.read_sql('SELECT user_id, pct_closed_loop FROM user_summary', conn)
safety_df = safety_df.merge(user_summary, on='user_id', how='left')

# Fill NaN events with 0
safety_df[['severe_hypo_events', 'dka_events', 'device_alerts']] = (
    safety_df[['severe_hypo_events', 'dka_events', 'device_alerts']].fillna(0)
)

# Create cohort columns
safety_df['age_group'] = pd.cut(safety_df['age'], bins=[17, 30, 50, 75],
                                labels=['18-30', '31-50', '51-75'])
safety_df['device_mode'] = np.where(safety_df['pct_closed_loop'] >= 50,
                                    'closed_loop', 'open_loop')

print(f"Safety dataframe shape: {safety_df.shape}")
safety_df.head()

What does this code do?

This cell documents the study's inclusion/exclusion criteria (a regulatory best practice) and builds a unified safety dataframe. It merges device-days exposure, clinical adverse events, and device alerts for each eligible user.

Users are classified into cohorts by age group, device mode (closed-loop if pct_closed_loop >= 50), and region. These cohorts form the basis for rate comparisons.

• In a new cell, compute event rates per 1,000 device-days with exact Poisson confidence intervals using the Garwood method:

def poisson_ci(k, alpha=0.05):
    """Compute exact Poisson CI using chi-squared quantiles (Garwood 1936).

    Parameters:
        k: observed count of events
        alpha: significance level (default 0.05 for 95% CI)

    Returns:
        (lower, upper) bounds on the Poisson mean
    """
    if k == 0:
        lower = 0.0
    else:
        lower = chi2.ppf(alpha / 2, 2 * k) / 2
    upper = chi2.ppf(1 - alpha / 2, 2 * (k + 1)) / 2
    return lower, upper

def compute_cohort_rates(df, cohort_col, event_col, exposure_col='device_days', per=1000):
    """Compute event rates per `per` device-days by cohort with Poisson CIs."""
    results = []
    for cohort, group in df.groupby(cohort_col):
        total_events = int(group[event_col].sum())
        total_exposure = group[exposure_col].sum()
        rate = (total_events / total_exposure) * per

        # Exact Poisson CI on the count, then convert to rate
        ci_low, ci_high = poisson_ci(total_events)
        rate_low = (ci_low / total_exposure) * per
        rate_high = (ci_high / total_exposure) * per

        results.append({
            'cohort': cohort,
            'events': total_events,
            'device_days': total_exposure,
            'rate_per_1000': rate,
            'ci_lower': rate_low,
            'ci_upper': rate_high
        })
    return pd.DataFrame(results)

# Compute rates by age group
print("=== Severe Hypoglycemia Rates by Age Group ===")
rates_age_hypo = compute_cohort_rates(safety_df, 'age_group', 'severe_hypo_events')
print(rates_age_hypo.to_string(index=False))

print("\n=== DKA Rates by Device Mode ===")
rates_mode_dka = compute_cohort_rates(safety_df, 'device_mode', 'dka_events')
print(rates_mode_dka.to_string(index=False))

print("\n=== Device Alert Rates by Region ===")
rates_region_alerts = compute_cohort_rates(safety_df, 'region', 'device_alerts')
print(rates_region_alerts.to_string(index=False))

How does the Garwood Poisson CI work?

The Garwood (1936) method is the standard exact confidence interval for Poisson counts. For k observed events, the 95% CI on the true mean is:

Lower = chi2.ppf(0.025, 2k) / 2 and Upper = chi2.ppf(0.975, 2(k+1)) / 2. Dividing by total exposure and multiplying by 1,000 converts from a count CI to a rate CI per 1,000 device-days. This approach is preferred over normal approximations because safety events are rare.

• Run both cells. You should see three tables printed showing event rates with confidence intervals broken down by age group, device mode, and region.

Seeing empty results or all zeros?

Check that your ETL ran successfully and that clinical_outcomes contains rows with severe_hypo_event = 1 or dka_event = 1. Run clinical['severe_hypo_event'].sum() to verify the raw event count is non-zero.

Also confirm the database path is correct relative to where Jupyter launched. If you opened Jupyter from a different folder, the relative path to data/processed/bionic_pancreas.db may not resolve.

Help me debug my safety event rate computation.

Build a c-chart for weekly safety event counts

A c-chart is a statistical process control tool that detects whether event counts deviate from their expected range. You will aggregate all safety events into weekly bins and plot them against control limits to spot weeks with unusually high counts.

• In a new cell, aggregate events by week and build the control chart:

# Aggregate safety events by week using clinical visit dates
clinical_dated = pd.read_sql("""
    SELECT visit_date,
           SUM(severe_hypo_event) + SUM(dka_event) as total_safety_events
    FROM clinical_outcomes
    GROUP BY visit_date
    ORDER BY visit_date
""", conn)
clinical_dated['visit_date'] = pd.to_datetime(clinical_dated['visit_date'])

# Resample to weekly counts
clinical_dated = clinical_dated.set_index('visit_date')
weekly_events = clinical_dated['total_safety_events'].resample('W').sum().reset_index()
weekly_events.columns = ['week', 'event_count']

# C-chart parameters
c_bar = weekly_events['event_count'].mean()
sigma = np.sqrt(c_bar)
ucl_3 = c_bar + 3 * sigma  # 3-sigma upper control limit
lcl_3 = max(0, c_bar - 3 * sigma)  # 3-sigma lower control limit
ucl_2 = c_bar + 2 * sigma  # 2-sigma warning limit
lcl_2 = max(0, c_bar - 2 * sigma)  # 2-sigma warning limit

print(f"Center line (c-bar): {c_bar:.2f}")
print(f"2-sigma limits: [{lcl_2:.2f}, {ucl_2:.2f}]")
print(f"3-sigma limits: [{lcl_3:.2f}, {ucl_3:.2f}]")
print(f"Weeks exceeding 3-sigma UCL: {(weekly_events['event_count'] > ucl_3).sum()}")

What is a c-chart?

A c-chart (counts control chart) monitors the number of events per equal time period. Because safety events follow a Poisson distribution, the standard deviation equals the square root of the mean. The center line is c-bar (mean weekly count).

The 2-sigma limits (c-bar +/- 2 * sqrt(c-bar)) serve as warning thresholds. The 3-sigma limits (c-bar +/- 3 * sqrt(c-bar)) are action thresholds. Points outside 3-sigma signal a statistically unusual week that warrants investigation.

• In the next cell, plot the c-chart:

# Plot the c-chart
fig, ax = plt.subplots(figsize=(14, 5))
ax.plot(weekly_events['week'], weekly_events['event_count'],
        marker='o', linewidth=1, markersize=4, label='Weekly events')

# Control limits
ax.axhline(c_bar, color='green', linestyle='-', linewidth=1.5, label=f'Center (c-bar={c_bar:.1f})')
ax.axhline(ucl_3, color='red', linestyle='--', linewidth=1.5, label=f'UCL 3-sigma={ucl_3:.1f}')
ax.axhline(lcl_3, color='red', linestyle='--', linewidth=1.5, label=f'LCL 3-sigma={lcl_3:.1f}')
ax.axhline(ucl_2, color='orange', linestyle=':', linewidth=1, label=f'UCL 2-sigma={ucl_2:.1f}')
ax.axhline(lcl_2, color='orange', linestyle=':', linewidth=1, label=f'LCL 2-sigma={lcl_2:.1f}')

# Highlight points exceeding 3-sigma
exceeding = weekly_events[weekly_events['event_count'] > ucl_3]
ax.scatter(exceeding['week'], exceeding['event_count'],
           color='red', s=80, zorder=5, label='Exceeds 3-sigma')

ax.set_xlabel('Week')
ax.set_ylabel('Safety Event Count')
ax.set_title('C-Chart: Weekly Safety Events with Control Limits')
ax.legend(loc='upper right', fontsize=8)
plt.tight_layout()
plt.savefig('notebooks/c_chart_safety_events.png', dpi=150)
plt.show()

print("C-chart saved to notebooks/c_chart_safety_events.png")

• Run both cells. You should see a line chart with green center line, orange 2-sigma warning limits, and red 3-sigma action limits. Any red highlighted dots represent weeks exceeding the upper control limit.

Chart is empty or shows a flat line?

If all events fall on the same date, the weekly resample will only produce one data point. Check your synthetic data generator produced visit dates spread across the full ~180-day window. Run clinical_dated.index.min() and clinical_dated.index.max() to confirm the date range spans multiple weeks.

Help me debug my c-chart visualization.

Estimate relative risk ratios with bootstrapped confidence intervals

Rates by cohort are useful, but stakeholders want to know: does closed-loop mode reduce the risk of adverse events compared to open-loop? A relative risk ratio with a confidence interval answers this directly.

• In a new cell, compute the relative risk and bootstrap a 95% confidence interval:

def compute_relative_risk(df, exposure_group, reference_group, event_col, exposure_col='device_days'):
    """Compute relative risk: rate(exposed) / rate(reference)."""
    exp = df[df['device_mode'] == exposure_group]
    ref = df[df['device_mode'] == reference_group]

    rate_exp = exp[event_col].sum() / exp[exposure_col].sum()
    rate_ref = ref[event_col].sum() / ref[exposure_col].sum()

    if rate_ref == 0:
        return np.inf
    return rate_exp / rate_ref

def bootstrap_rr_ci(df, exposure_group, reference_group, event_col,
                    exposure_col='device_days', n_boot=2000, alpha=0.05, seed=42):
    """Bootstrap 95% CI for relative risk by resampling users."""
    rng = np.random.default_rng(seed)
    rr_samples = []

    for _ in range(n_boot):
        # Resample users with replacement within each group
        exp_mask = df['device_mode'] == exposure_group
        ref_mask = df['device_mode'] == reference_group

        boot_exp = df[exp_mask].sample(n=exp_mask.sum(), replace=True, random_state=rng.integers(1e9))
        boot_ref = df[ref_mask].sample(n=ref_mask.sum(), replace=True, random_state=rng.integers(1e9))

        rate_exp = boot_exp[event_col].sum() / boot_exp[exposure_col].sum()
        rate_ref = boot_ref[event_col].sum() / boot_ref[exposure_col].sum()

        if rate_ref > 0:
            rr_samples.append(rate_exp / rate_ref)

    rr_samples = np.array(rr_samples)
    ci_lower = np.percentile(rr_samples, 100 * (alpha / 2))
    ci_upper = np.percentile(rr_samples, 100 * (1 - alpha / 2))
    return ci_lower, ci_upper

# Compute RR for severe hypoglycemia: closed-loop vs open-loop
rr_hypo = compute_relative_risk(safety_df, 'closed_loop', 'open_loop', 'severe_hypo_events')
ci_hypo = bootstrap_rr_ci(safety_df, 'closed_loop', 'open_loop', 'severe_hypo_events')

print("=== Relative Risk: Closed-Loop vs Open-Loop ===")
print(f"Severe Hypoglycemia RR: {rr_hypo:.3f} (95% CI: [{ci_hypo[0]:.3f}, {ci_hypo[1]:.3f}])")

# Compute RR for DKA
rr_dka = compute_relative_risk(safety_df, 'closed_loop', 'open_loop', 'dka_events')
ci_dka = bootstrap_rr_ci(safety_df, 'closed_loop', 'open_loop', 'dka_events')
print(f"DKA RR: {rr_dka:.3f} (95% CI: [{ci_dka[0]:.3f}, {ci_dka[1]:.3f}])")

# Compute RR for device alerts
rr_alerts = compute_relative_risk(safety_df, 'closed_loop', 'open_loop', 'device_alerts')
ci_alerts = bootstrap_rr_ci(safety_df, 'closed_loop', 'open_loop', 'device_alerts')
print(f"Device Alerts RR: {rr_alerts:.3f} (95% CI: [{ci_alerts[0]:.3f}, {ci_alerts[1]:.3f}])")

# Interpretation guide
print("\nInterpretation: RR < 1 means closed-loop has LOWER risk than open-loop.")
print("If the CI does not cross 1.0, the difference is statistically significant at alpha=0.05.")

What does this code do?

The relative risk (RR) is the ratio of event rates between two groups. An RR of 0.6 means the exposed group (closed-loop) has 40% lower risk than the reference group (open-loop).

The bootstrap resamples users (not individual events) 2,000 times to build a distribution of plausible RR values. The 2.5th and 97.5th percentiles of this distribution form the 95% CI. If the interval excludes 1.0, the difference between groups is statistically significant.

• Run the cell. You should see RR values for each event type with their confidence intervals. An RR below 1.0 indicates closed-loop mode has a lower event rate.
• Save the notebook with Ctrl+S (Windows) or Cmd+S (macOS).

Getting inf or NaN for relative risk?

An infinite RR means the reference group (open-loop) had zero events, making division impossible. Check that your cohort split using pct_closed_loop >= 50 produces enough users in both groups. Run safety_df['device_mode'].value_counts() to verify both groups have adequate sample sizes.

Help me fix my relative risk calculation returning inf or NaN.

✔️ Awesome, I've got everything!

Great. Double check you have saved your notebook.

ⓧ I'd like to double check the full code

Here is the full notebook code for notebooks/01_safety_signal_detection.ipynb:

import sqlite3
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import chi2

# Connect to the analytical warehouse
DB_PATH = 'data/processed/bionic_pancreas.db'
conn = sqlite3.connect(DB_PATH)

# Set plot style
sns.set_style('whitegrid')
print('Connected to database successfully.')

# =============================================================
# INCLUSION / EXCLUSION CRITERIA
# =============================================================
# Inclusion:
#   - All enrolled users with >= 30 days of device telemetry
#   - Safety events: severe_hypo_event, dka_event (clinical_outcomes),
#     and device alerts classified as adverse (alert_type != None)
# Exclusion:
#   - Users with < 30 days of telemetry (insufficient exposure)
#   - Battery/informational alerts (only clinically significant alerts count)
# =============================================================

# Load demographics and compute device-days per user
demographics = pd.read_sql('SELECT * FROM demographics_enrollment', conn)
telemetry_summary = pd.read_sql("""
    SELECT user_id,
           COUNT(DISTINCT DATE(timestamp)) as device_days,
           SUM(CASE WHEN alert_type IS NOT NULL THEN 1 ELSE 0 END) as device_alerts
    FROM device_telemetry
    GROUP BY user_id
""", conn)

# Load clinical safety events
clinical = pd.read_sql('SELECT * FROM clinical_outcomes', conn)

# Apply inclusion criteria: >= 30 device-days
eligible_users = telemetry_summary[telemetry_summary['device_days'] >= 30].copy()
print(f"Eligible users (>= 30 device-days): {len(eligible_users)} / {len(telemetry_summary)}")

# Aggregate safety events per user
user_safety = clinical.groupby('user_id').agg(
    severe_hypo_events=('severe_hypo_event', 'sum'),
    dka_events=('dka_event', 'sum')
).reset_index()

# Merge all sources
safety_df = eligible_users.merge(user_safety, on='user_id', how='left')
safety_df = safety_df.merge(demographics[['user_id', 'age', 'region']], on='user_id', how='left')

# Add device mode from user_summary
user_summary = pd.read_sql('SELECT user_id, pct_closed_loop FROM user_summary', conn)
safety_df = safety_df.merge(user_summary, on='user_id', how='left')

# Fill NaN events with 0
safety_df[['severe_hypo_events', 'dka_events', 'device_alerts']] = (
    safety_df[['severe_hypo_events', 'dka_events', 'device_alerts']].fillna(0)
)

# Create cohort columns
safety_df['age_group'] = pd.cut(safety_df['age'], bins=[17, 30, 50, 75],
                                labels=['18-30', '31-50', '51-75'])
safety_df['device_mode'] = np.where(safety_df['pct_closed_loop'] >= 50,
                                    'closed_loop', 'open_loop')

print(f"Safety dataframe shape: {safety_df.shape}")
safety_df.head()

# --- Poisson CI and cohort rates ---

def poisson_ci(k, alpha=0.05):
    """Compute exact Poisson CI using chi-squared quantiles (Garwood 1936).

    Parameters:
        k: observed count of events
        alpha: significance level (default 0.05 for 95% CI)

    Returns:
        (lower, upper) bounds on the Poisson mean
    """
    if k == 0:
        lower = 0.0
    else:
        lower = chi2.ppf(alpha / 2, 2 * k) / 2
    upper = chi2.ppf(1 - alpha / 2, 2 * (k + 1)) / 2
    return lower, upper

def compute_cohort_rates(df, cohort_col, event_col, exposure_col='device_days', per=1000):
    """Compute event rates per `per` device-days by cohort with Poisson CIs."""
    results = []
    for cohort, group in df.groupby(cohort_col):
        total_events = int(group[event_col].sum())
        total_exposure = group[exposure_col].sum()
        rate = (total_events / total_exposure) * per

        # Exact Poisson CI on the count, then convert to rate
        ci_low, ci_high = poisson_ci(total_events)
        rate_low = (ci_low / total_exposure) * per
        rate_high = (ci_high / total_exposure) * per

        results.append({
            'cohort': cohort,
            'events': total_events,
            'device_days': total_exposure,
            'rate_per_1000': rate,
            'ci_lower': rate_low,
            'ci_upper': rate_high
        })
    return pd.DataFrame(results)

# Compute rates by age group
print("=== Severe Hypoglycemia Rates by Age Group ===")
rates_age_hypo = compute_cohort_rates(safety_df, 'age_group', 'severe_hypo_events')
print(rates_age_hypo.to_string(index=False))

print("\n=== DKA Rates by Device Mode ===")
rates_mode_dka = compute_cohort_rates(safety_df, 'device_mode', 'dka_events')
print(rates_mode_dka.to_string(index=False))

print("\n=== Device Alert Rates by Region ===")
rates_region_alerts = compute_cohort_rates(safety_df, 'region', 'device_alerts')
print(rates_region_alerts.to_string(index=False))

# --- C-Chart ---

# Aggregate safety events by week using clinical visit dates
clinical_dated = pd.read_sql("""
    SELECT visit_date,
           SUM(severe_hypo_event) + SUM(dka_event) as total_safety_events
    FROM clinical_outcomes
    GROUP BY visit_date
    ORDER BY visit_date
""", conn)
clinical_dated['visit_date'] = pd.to_datetime(clinical_dated['visit_date'])

# Resample to weekly counts
clinical_dated = clinical_dated.set_index('visit_date')
weekly_events = clinical_dated['total_safety_events'].resample('W').sum().reset_index()
weekly_events.columns = ['week', 'event_count']

# C-chart parameters
c_bar = weekly_events['event_count'].mean()
sigma = np.sqrt(c_bar)
ucl_3 = c_bar + 3 * sigma  # 3-sigma upper control limit
lcl_3 = max(0, c_bar - 3 * sigma)  # 3-sigma lower control limit
ucl_2 = c_bar + 2 * sigma  # 2-sigma warning limit
lcl_2 = max(0, c_bar - 2 * sigma)  # 2-sigma warning limit

print(f"Center line (c-bar): {c_bar:.2f}")
print(f"2-sigma limits: [{lcl_2:.2f}, {ucl_2:.2f}]")
print(f"3-sigma limits: [{lcl_3:.2f}, {ucl_3:.2f}]")
print(f"Weeks exceeding 3-sigma UCL: {(weekly_events['event_count'] > ucl_3).sum()}")

# Plot the c-chart
fig, ax = plt.subplots(figsize=(14, 5))
ax.plot(weekly_events['week'], weekly_events['event_count'],
        marker='o', linewidth=1, markersize=4, label='Weekly events')

# Control limits
ax.axhline(c_bar, color='green', linestyle='-', linewidth=1.5, label=f'Center (c-bar={c_bar:.1f})')
ax.axhline(ucl_3, color='red', linestyle='--', linewidth=1.5, label=f'UCL 3-sigma={ucl_3:.1f}')
ax.axhline(lcl_3, color='red', linestyle='--', linewidth=1.5, label=f'LCL 3-sigma={lcl_3:.1f}')
ax.axhline(ucl_2, color='orange', linestyle=':', linewidth=1, label=f'UCL 2-sigma={ucl_2:.1f}')
ax.axhline(lcl_2, color='orange', linestyle=':', linewidth=1, label=f'LCL 2-sigma={lcl_2:.1f}')

# Highlight points exceeding 3-sigma
exceeding = weekly_events[weekly_events['event_count'] > ucl_3]
ax.scatter(exceeding['week'], exceeding['event_count'],
           color='red', s=80, zorder=5, label='Exceeds 3-sigma')

ax.set_xlabel('Week')
ax.set_ylabel('Safety Event Count')
ax.set_title('C-Chart: Weekly Safety Events with Control Limits')
ax.legend(loc='upper right', fontsize=8)
plt.tight_layout()
plt.savefig('notebooks/c_chart_safety_events.png', dpi=150)
plt.show()

print("C-chart saved to notebooks/c_chart_safety_events.png")

# --- Relative Risk with Bootstrap CI ---

def compute_relative_risk(df, exposure_group, reference_group, event_col, exposure_col='device_days'):
    """Compute relative risk: rate(exposed) / rate(reference)."""
    exp = df[df['device_mode'] == exposure_group]
    ref = df[df['device_mode'] == reference_group]

    rate_exp = exp[event_col].sum() / exp[exposure_col].sum()
    rate_ref = ref[event_col].sum() / ref[exposure_col].sum()

    if rate_ref == 0:
        return np.inf
    return rate_exp / rate_ref

def bootstrap_rr_ci(df, exposure_group, reference_group, event_col,
                    exposure_col='device_days', n_boot=2000, alpha=0.05, seed=42):
    """Bootstrap 95% CI for relative risk by resampling users."""
    rng = np.random.default_rng(seed)
    rr_samples = []

    for _ in range(n_boot):
        # Resample users with replacement within each group
        exp_mask = df['device_mode'] == exposure_group
        ref_mask = df['device_mode'] == reference_group

        boot_exp = df[exp_mask].sample(n=exp_mask.sum(), replace=True, random_state=rng.integers(1e9))
        boot_ref = df[ref_mask].sample(n=ref_mask.sum(), replace=True, random_state=rng.integers(1e9))

        rate_exp = boot_exp[event_col].sum() / boot_exp[exposure_col].sum()
        rate_ref = boot_ref[event_col].sum() / boot_ref[exposure_col].sum()

        if rate_ref > 0:
            rr_samples.append(rate_exp / rate_ref)

    rr_samples = np.array(rr_samples)
    ci_lower = np.percentile(rr_samples, 100 * (alpha / 2))
    ci_upper = np.percentile(rr_samples, 100 * (1 - alpha / 2))
    return ci_lower, ci_upper

# Compute RR for severe hypoglycemia: closed-loop vs open-loop
rr_hypo = compute_relative_risk(safety_df, 'closed_loop', 'open_loop', 'severe_hypo_events')
ci_hypo = bootstrap_rr_ci(safety_df, 'closed_loop', 'open_loop', 'severe_hypo_events')

print("=== Relative Risk: Closed-Loop vs Open-Loop ===")
print(f"Severe Hypoglycemia RR: {rr_hypo:.3f} (95% CI: [{ci_hypo[0]:.3f}, {ci_hypo[1]:.3f}])")

# Compute RR for DKA
rr_dka = compute_relative_risk(safety_df, 'closed_loop', 'open_loop', 'dka_events')
ci_dka = bootstrap_rr_ci(safety_df, 'closed_loop', 'open_loop', 'dka_events')
print(f"DKA RR: {rr_dka:.3f} (95% CI: [{ci_dka[0]:.3f}, {ci_dka[1]:.3f}])")

# Compute RR for device alerts
rr_alerts = compute_relative_risk(safety_df, 'closed_loop', 'open_loop', 'device_alerts')
ci_alerts = bootstrap_rr_ci(safety_df, 'closed_loop', 'open_loop', 'device_alerts')
print(f"Device Alerts RR: {rr_alerts:.3f} (95% CI: [{ci_alerts[0]:.3f}, {ci_alerts[1]:.3f}])")

# Interpretation guide
print("\nInterpretation: RR < 1 means closed-loop has LOWER risk than open-loop.")
print("If the CI does not cross 1.0, the difference is statistically significant at alpha=0.05.")

You now have a complete safety signal detection notebook. Your c-chart can catch unusual weeks in real time, and the relative risk analysis quantifies whether closed-loop mode provides a safety benefit. Next up, you will use similar statistical rigor to evaluate clinical outcomes with pre/post HbA1c analysis.

Analyze Clinical Outcomes

In the previous step, you built a safety signal detection system that flags abnormal event rates across patient cohorts. Now you'll answer the core clinical question: does the bionic pancreas actually improve glycemic outcomes after patients start using it?

Clinical outcomes analysis compares each user's health metrics before and after they started the device. You'll define pre/post cohorts, run a formal statistical test, calculate an effect size, and compute time-in-range metrics that clinicians care about.

In this step, get ready to:

• Define pre/post cohorts and compute paired HbA1c changes with a Wilcoxon signed-rank test and Cohen's d effect size.
• Calculate time-in-range and time-below-range metrics with a difference-in-differences approach.
• Perform sensitivity analyses varying lookback windows, documenting missingness, and stratifying by demographics.

Define pre/post cohorts and compute paired HbA1c change

To measure whether the device improves outcomes, you need to compare each user's HbA1c before they started the device versus after they've been using it for a meaningful period. This requires strict inclusion criteria so the comparison is valid.

• Create a new notebook file called 02_clinical_outcomes.ipynb inside your notebooks/ folder.
• Add the following imports and data loading cell to connect to your analytical warehouse:

import pandas as pd
import numpy as np
import sqlite3
from scipy import stats
import matplotlib.pyplot as plt
import seaborn as sns

# Connect to the analytical warehouse
conn = sqlite3.connect("../data/processed/bionic_pancreas.db")

# Load required tables
demographics = pd.read_sql("SELECT * FROM demographics_enrollment", conn)
clinical = pd.read_sql("SELECT * FROM clinical_outcomes", conn)
telemetry = pd.read_sql("SELECT * FROM device_telemetry", conn)

# Convert date columns to datetime
demographics["device_start_date"] = pd.to_datetime(demographics["device_start_date"])
clinical["visit_date"] = pd.to_datetime(clinical["visit_date"])
telemetry["timestamp"] = pd.to_datetime(telemetry["timestamp"])

print(f"Loaded {len(demographics)} users, {len(clinical)} clinical visits, {len(telemetry)} telemetry readings")

What does this code do?

This cell imports all the libraries you'll need and pulls data from three tables in your SQLite warehouse.

You're loading demographics (for device start dates), clinical outcomes (for HbA1c measurements), and device telemetry (for glucose readings) so you can compare pre-device vs. post-device metrics for each user.

• Run the cell. You should see a message confirming 500 users, approximately 2000 clinical visits, and millions of telemetry readings loaded.

Next, you'll define the cohort logic. The pre-period is the 90-day window before each user's device start date. The post-period is 90 or more days after they started using the device.

• Add a new cell with the cohort definition and inclusion criteria:

# Define pre/post cohort parameters
LOOKBACK_DAYS = 90
POST_MIN_DAYS = 90

# Inclusion criteria:
# 1. At least one HbA1c measurement in the 90-day pre-period
# 2. At least one HbA1c measurement 90+ days post device start
# 3. Valid device_start_date

# Merge clinical visits with device start dates
clinical_merged = clinical.merge(
    demographics[["user_id", "device_start_date"]], on="user_id"
)

# Calculate days relative to device start for each visit
clinical_merged["days_from_start"] = (
    clinical_merged["visit_date"] - clinical_merged["device_start_date"]
).dt.days

# Pre-period: within LOOKBACK_DAYS before device_start_date
pre_visits = clinical_merged[
    (clinical_merged["days_from_start"] >= -LOOKBACK_DAYS) &
    (clinical_merged["days_from_start"] < 0)
]

# Post-period: 90+ days after device_start_date
post_visits = clinical_merged[
    clinical_merged["days_from_start"] >= POST_MIN_DAYS
]

# Get last pre-HbA1c and last post-HbA1c per user
pre_hba1c = pre_visits.groupby("user_id")["hba1c_pct"].last().reset_index()
pre_hba1c.columns = ["user_id", "pre_hba1c"]

post_hba1c = post_visits.groupby("user_id")["hba1c_pct"].last().reset_index()
post_hba1c.columns = ["user_id", "post_hba1c"]

# Apply inclusion criteria: users must have both pre and post measurements
paired_cohort = pre_hba1c.merge(post_hba1c, on="user_id")
paired_cohort["hba1c_change"] = paired_cohort["post_hba1c"] - paired_cohort["pre_hba1c"]

print(f"Paired cohort: {len(paired_cohort)} users with both pre and post HbA1c")
print(f"Excluded: {len(demographics) - len(paired_cohort)} users (missing pre or post measurement)")

What are the inclusion criteria doing?

The 90-day lookback window captures each user's baseline HbA1c before they started the device. The 90-day minimum post-period ensures users have had enough time on the device for measurable glycemic improvement.

Users missing either measurement are excluded because you can't compute a valid paired difference without both data points. This is a standard approach in clinical research to avoid bias from incomplete data.

• Run the cell. You should see the number of users included in your paired cohort and how many were excluded.

Now you'll run the formal statistical test. The Wilcoxon signed-rank test is a non-parametric alternative to the paired t-test that doesn't assume the differences are normally distributed.

• Add a new cell for the statistical test and effect size calculation:

# Wilcoxon signed-rank test (non-parametric paired test)
# Tests whether the distribution of pre-post differences is symmetric about zero
stat, p_value = stats.wilcoxon(paired_cohort["pre_hba1c"], paired_cohort["post_hba1c"])

# Cohen's d effect size for paired samples
# d = mean(differences) / std(differences)
mean_change = paired_cohort["hba1c_change"].mean()
std_change = paired_cohort["hba1c_change"].std()
cohens_d = mean_change / std_change

print(f"=== Pre/Post HbA1c Analysis ===")
print(f"Mean HbA1c change: {mean_change:.3f}%")
print(f"Wilcoxon signed-rank statistic: {stat:.1f}")
print(f"P-value: {p_value:.6f}")
print(f"Cohen's d: {cohens_d:.3f}")
print(f"\nInterpretation: {'Statistically significant' if p_value < 0.05 else 'Not significant'} (alpha=0.05)")
print(f"Effect size: {'Small' if abs(cohens_d) < 0.5 else 'Medium' if abs(cohens_d) < 0.8 else 'Large'} (|d|={abs(cohens_d):.3f})")

Why Wilcoxon instead of a paired t-test?

The Wilcoxon signed-rank test makes fewer assumptions about your data's distribution. In clinical data, HbA1c changes are often skewed, so a non-parametric test is more appropriate than assuming normality.

Cohen's d tells you the practical magnitude of the effect. A p-value alone only tells you the change is non-zero. Cohen's d tells you whether it's clinically meaningful (0.2 = small, 0.5 = medium, 0.8 = large).

• Run the cell. You should see the mean HbA1c change (likely negative, indicating improvement), a p-value, and a Cohen's d value.

Seeing a ValueError from wilcoxon?

If you get an error about sample sizes, check that your paired cohort has more than zero users. Run len(paired_cohort) to confirm it's not empty.

If the cohort is empty, verify your date parsing produced valid datetime objects by checking demographics["device_start_date"].dtype is datetime64. Help me debug why my paired cohort is empty.

Calculate time-in-range and time-below-range metrics

HbA1c is a quarterly snapshot. Continuous glucose data tells a richer story through time-in-range (percentage of readings between 70-180 mg/dL) and time-below-range (percentage below 54 mg/dL, indicating dangerous hypoglycemia).

• Add a new cell to compute glycemic metrics from the telemetry data using the same pre/post periods:

# Time-in-range: % of glucose readings between 70-180 mg/dL
# Time-below-range: % of glucose readings below 54 mg/dL

# Merge telemetry with device start dates
telemetry_merged = telemetry.merge(
    demographics[["user_id", "device_start_date"]], on="user_id"
)
telemetry_merged["days_from_start"] = (
    telemetry_merged["timestamp"] - telemetry_merged["device_start_date"]
).dt.days

# Define pre and post periods for telemetry
pre_telemetry = telemetry_merged[
    (telemetry_merged["days_from_start"] >= -LOOKBACK_DAYS) &
    (telemetry_merged["days_from_start"] < 0)
]
post_telemetry = telemetry_merged[
    telemetry_merged["days_from_start"] >= POST_MIN_DAYS
]

def calculate_glycemic_metrics(df):
    """Calculate time-in-range and time-below-range per user."""
    metrics = df.groupby("user_id")["glucose_mg_dl"].agg(
        total_readings="count",
        tir_count=lambda x: ((x >= 70) & (x <= 180)).sum(),
        tbr_count=lambda x: (x < 54).sum()
    ).reset_index()
    metrics["time_in_range_pct"] = (metrics["tir_count"] / metrics["total_readings"]) * 100
    metrics["time_below_range_pct"] = (metrics["tbr_count"] / metrics["total_readings"]) * 100
    return metrics

# Calculate metrics for both periods
pre_metrics = calculate_glycemic_metrics(pre_telemetry)
post_metrics = calculate_glycemic_metrics(post_telemetry)

pre_metrics.columns = ["user_id", "pre_readings", "pre_tir_count", "pre_tbr_count", "pre_tir_pct", "pre_tbr_pct"]
post_metrics.columns = ["user_id", "post_readings", "post_tir_count", "post_tbr_count", "post_tir_pct", "post_tbr_pct"]

print(f"Pre-period telemetry: {len(pre_metrics)} users with glucose data")
print(f"Post-period telemetry: {len(post_metrics)} users with glucose data")

What are time-in-range and time-below-range?

Time-in-range (70-180 mg/dL) is the international consensus target for automated insulin delivery devices. A higher percentage means better glucose control.

Time-below-range (below 54 mg/dL) measures dangerous hypoglycemia. This metric should decrease with device use. Even a small reduction in time-below-range is clinically significant because severe hypoglycemia can be life-threatening.

• Run the cell. You should see two user counts confirming telemetry data exists for both periods.
• Add the difference-in-differences cell to combine the pre/post metrics:

# Difference-in-differences approach for users with sufficient data
tir_paired = pre_metrics.merge(post_metrics, on="user_id")
tir_paired["tir_change"] = tir_paired["post_tir_pct"] - tir_paired["pre_tir_pct"]
tir_paired["tbr_change"] = tir_paired["post_tbr_pct"] - tir_paired["pre_tbr_pct"]

print(f"Users with both pre and post telemetry: {len(tir_paired)}")
print(f"\n{'='*50}")
print(f"Time-in-Range (70-180 mg/dL):")
print(f"  Pre:  {tir_paired['pre_tir_pct'].mean():.1f}% +/- {tir_paired['pre_tir_pct'].std():.1f}%")
print(f"  Post: {tir_paired['post_tir_pct'].mean():.1f}% +/- {tir_paired['post_tir_pct'].std():.1f}%")
print(f"  Mean change: {tir_paired['tir_change'].mean():.1f}%")
print(f"\nTime-Below-Range (<54 mg/dL):")
print(f"  Pre:  {tir_paired['pre_tbr_pct'].mean():.2f}% +/- {tir_paired['pre_tbr_pct'].std():.2f}%")
print(f"  Post: {tir_paired['post_tbr_pct'].mean():.2f}% +/- {tir_paired['post_tbr_pct'].std():.2f}%")
print(f"  Mean change: {tir_paired['tbr_change'].mean():.2f}%")

# Statistical tests on TIR and TBR changes
tir_stat, tir_p = stats.wilcoxon(tir_paired["pre_tir_pct"], tir_paired["post_tir_pct"])
tbr_stat, tbr_p = stats.wilcoxon(tir_paired["pre_tbr_pct"], tir_paired["post_tbr_pct"])

print(f"\nWilcoxon test - TIR change: p={tir_p:.6f}")
print(f"Wilcoxon test - TBR change: p={tbr_p:.6f}")

What is difference-in-differences?

Difference-in-differences compares each user's own pre-period metrics to their post-period metrics, then aggregates across the cohort. This controls for individual baseline variation because each user serves as their own control.

You expect time-in-range to increase (glucose better controlled) and time-below-range to decrease (fewer dangerous lows) after the device is activated.

• Run the cell. You should see pre and post percentages for both metrics, the mean change, and Wilcoxon p-values for each.

Getting all zeros for time-below-range?

Severe hypoglycemia (below 54 mg/dL) is intentionally rare in the synthetic data (2-5% of users). If time-below-range is very close to zero in both periods, that's expected. The Wilcoxon test may return a warning about identical values.

Help me understand why my TBR values are all near zero.

Perform sensitivity analyses

A single analysis with one set of parameters isn't enough for a credible portfolio piece. Sensitivity analyses show that your conclusions hold under different assumptions about the lookback window, and that you've investigated potential biases from missing data.

• Add a cell to test different lookback windows (60, 90, and 120 days):

# Sensitivity analysis: vary lookback windows
windows = [60, 90, 120]
sensitivity_results = []

for window in windows:
    # Redefine pre-period with varied window
    pre_w = clinical_merged[
        (clinical_merged["days_from_start"] >= -window) &
        (clinical_merged["days_from_start"] < 0)
    ]
    pre_w_hba1c = pre_w.groupby("user_id")["hba1c_pct"].last().reset_index()
    pre_w_hba1c.columns = ["user_id", "pre_hba1c"]

    # Post stays the same (90+ days after)
    paired_w = pre_w_hba1c.merge(post_hba1c, on="user_id")
    paired_w["hba1c_change"] = paired_w["post_hba1c"] - paired_w["pre_hba1c"]

    if len(paired_w) > 1:
        stat_w, p_w = stats.wilcoxon(paired_w["pre_hba1c"], paired_w["post_hba1c"])
        d_w = paired_w["hba1c_change"].mean() / paired_w["hba1c_change"].std()
    else:
        stat_w, p_w, d_w = np.nan, np.nan, np.nan

    sensitivity_results.append({
        "lookback_days": window,
        "n_users": len(paired_w),
        "mean_change": paired_w["hba1c_change"].mean() if len(paired_w) > 0 else np.nan,
        "p_value": p_w,
        "cohens_d": d_w
    })

sensitivity_df = pd.DataFrame(sensitivity_results)
print("Sensitivity Analysis: Lookback Window Variation")
print("=" * 55)
print(sensitivity_df.to_string(index=False))

Why vary the lookback window?

A 90-day window is one reasonable choice, but results should be robust to that choice. If the effect disappears when you use 60 or 120 days instead, your conclusion is fragile and depends too much on an arbitrary parameter.

Consistent results across windows strengthen your credibility with reviewers and hiring managers.

• Run the cell. You should see a table with three rows showing how the cohort size, mean change, p-value, and Cohen's d vary across 60, 90, and 120-day windows.
• Add a cell to document missingness patterns:

# Document missingness patterns
users_with_pre = set(pre_hba1c["user_id"])
users_with_post = set(post_hba1c["user_id"])
all_users = set(demographics["user_id"])

missing_pre = all_users - users_with_pre
missing_post = all_users - users_with_post
missing_both = missing_pre & missing_post

print("Missingness Report")
print("=" * 40)
print(f"Total users: {len(all_users)}")
print(f"Users with pre-period HbA1c: {len(users_with_pre)} ({len(users_with_pre)/len(all_users)*100:.1f}%)")
print(f"Users with post-period HbA1c: {len(users_with_post)} ({len(users_with_post)/len(all_users)*100:.1f}%)")
print(f"Users missing pre-period: {len(missing_pre)} ({len(missing_pre)/len(all_users)*100:.1f}%)")
print(f"Users missing post-period: {len(missing_post)} ({len(missing_post)/len(all_users)*100:.1f}%)")
print(f"Users missing both: {len(missing_both)} ({len(missing_both)/len(all_users)*100:.1f}%)")
print(f"\nFinal paired cohort: {len(paired_cohort)} users ({len(paired_cohort)/len(all_users)*100:.1f}%)")

# Check if missingness correlates with demographics
missing_users = demographics[demographics["user_id"].isin(missing_pre | missing_post)]
included_users = demographics[demographics["user_id"].isin(paired_cohort["user_id"])]

print(f"\nDemographic comparison (missing vs. included):")
print(f"  Mean age - Missing: {missing_users['age'].mean():.1f}, Included: {included_users['age'].mean():.1f}")
print(f"  Mean diabetes duration - Missing: {missing_users['diabetes_duration_yrs'].mean():.1f}, Included: {included_users['diabetes_duration_yrs'].mean():.1f}")

Why document missingness?

If the users excluded from your analysis are systematically different from those included (older, sicker, newer to the device), your results may not generalize. Documenting this is standard practice in clinical research and shows reviewers you've considered selection bias.

• Run the cell. You should see a complete missingness report with demographic comparisons between missing and included groups.
• Add a final cell to stratify results by diabetes duration and age:

# Stratification by diabetes duration and age
paired_with_demo = paired_cohort.merge(
    demographics[["user_id", "age", "diabetes_duration_yrs"]], on="user_id"
)

# Diabetes duration strata
paired_with_demo["duration_group"] = pd.cut(
    paired_with_demo["diabetes_duration_yrs"],
    bins=[0, 5, 15, 50],
    labels=["Short (0-5 yrs)", "Medium (5-15 yrs)", "Long (15+ yrs)"]
)

# Age strata
paired_with_demo["age_group"] = pd.cut(
    paired_with_demo["age"],
    bins=[0, 30, 50, 100],
    labels=["Young (18-30)", "Middle (31-50)", "Older (51+)"]
)

print("Stratified HbA1c Change by Diabetes Duration")
print("=" * 55)
for group, data in paired_with_demo.groupby("duration_group"):
    if len(data) > 1:
        stat_s, p_s = stats.wilcoxon(data["pre_hba1c"], data["post_hba1c"])
        d_s = data["hba1c_change"].mean() / data["hba1c_change"].std()
        print(f"  {group}: n={len(data)}, mean change={data['hba1c_change'].mean():.3f}, p={p_s:.4f}, d={d_s:.3f}")

print(f"\nStratified HbA1c Change by Age Group")
print("=" * 55)
for group, data in paired_with_demo.groupby("age_group"):
    if len(data) > 1:
        stat_s, p_s = stats.wilcoxon(data["pre_hba1c"], data["post_hba1c"])
        d_s = data["hba1c_change"].mean() / data["hba1c_change"].std()
        print(f"  {group}: n={len(data)}, mean change={data['hba1c_change'].mean():.3f}, p={p_s:.4f}, d={d_s:.3f}")

# Close the database connection
conn.close()
print("\nAnalysis complete. Connection closed.")

What does stratification reveal?

Stratification breaks the cohort into subgroups to see whether the device's effect is consistent across different populations. Patients with longer diabetes duration might respond differently than newly diagnosed patients.

Similarly, age-based stratification can reveal whether younger or older users benefit more from automated insulin delivery. These insights are exactly what a medical device company would use to target specific patient populations.

• Run the cell. You should see stratified results for three diabetes duration groups and three age groups, each with sample size, mean change, p-value, and Cohen's d.
• Save your notebook.

Wilcoxon warning about too few observations in a stratum?

Some strata may have very few users. The Wilcoxon test needs at least a handful of observations to produce meaningful results. If a stratum has fewer than 10 users, the p-value is unreliable. Note this as a limitation in your analysis rather than treating it as a code error.

Help me handle small sample sizes in my stratified analysis.

✔️ Awesome, I've got everything!

Great. Make sure your notebook is saved and all cells have been executed in order.

ⓧ I'd like to double check the full code

Here is the complete notebook code for notebooks/02_clinical_outcomes.ipynb:

import pandas as pd
import numpy as np
import sqlite3
from scipy import stats
import matplotlib.pyplot as plt
import seaborn as sns

# Connect to the analytical warehouse
conn = sqlite3.connect("../data/processed/bionic_pancreas.db")

# Load required tables
demographics = pd.read_sql("SELECT * FROM demographics_enrollment", conn)
clinical = pd.read_sql("SELECT * FROM clinical_outcomes", conn)
telemetry = pd.read_sql("SELECT * FROM device_telemetry", conn)

# Convert date columns to datetime
demographics["device_start_date"] = pd.to_datetime(demographics["device_start_date"])
clinical["visit_date"] = pd.to_datetime(clinical["visit_date"])
telemetry["timestamp"] = pd.to_datetime(telemetry["timestamp"])

print(f"Loaded {len(demographics)} users, {len(clinical)} clinical visits, {len(telemetry)} telemetry readings")

# --- Cell 2: Define pre/post cohorts ---

# Define pre/post cohort parameters
LOOKBACK_DAYS = 90
POST_MIN_DAYS = 90

# Inclusion criteria:
# 1. At least one HbA1c measurement in the 90-day pre-period
# 2. At least one HbA1c measurement 90+ days post device start
# 3. Valid device_start_date

# Merge clinical visits with device start dates
clinical_merged = clinical.merge(
    demographics[["user_id", "device_start_date"]], on="user_id"
)

# Calculate days relative to device start for each visit
clinical_merged["days_from_start"] = (
    clinical_merged["visit_date"] - clinical_merged["device_start_date"]
).dt.days

# Pre-period: within LOOKBACK_DAYS before device_start_date
pre_visits = clinical_merged[
    (clinical_merged["days_from_start"] >= -LOOKBACK_DAYS) &
    (clinical_merged["days_from_start"] < 0)
]

# Post-period: 90+ days after device_start_date
post_visits = clinical_merged[
    clinical_merged["days_from_start"] >= POST_MIN_DAYS
]

# Get last pre-HbA1c and last post-HbA1c per user
pre_hba1c = pre_visits.groupby("user_id")["hba1c_pct"].last().reset_index()
pre_hba1c.columns = ["user_id", "pre_hba1c"]

post_hba1c = post_visits.groupby("user_id")["hba1c_pct"].last().reset_index()
post_hba1c.columns = ["user_id", "post_hba1c"]

# Apply inclusion criteria: users must have both pre and post measurements
paired_cohort = pre_hba1c.merge(post_hba1c, on="user_id")
paired_cohort["hba1c_change"] = paired_cohort["post_hba1c"] - paired_cohort["pre_hba1c"]

print(f"Paired cohort: {len(paired_cohort)} users with both pre and post HbA1c")
print(f"Excluded: {len(demographics) - len(paired_cohort)} users (missing pre or post measurement)")

# --- Cell 3: Wilcoxon test and Cohen's d ---

# Wilcoxon signed-rank test (non-parametric paired test)
# Tests whether the distribution of pre-post differences is symmetric about zero
stat, p_value = stats.wilcoxon(paired_cohort["pre_hba1c"], paired_cohort["post_hba1c"])

# Cohen's d effect size for paired samples
# d = mean(differences) / std(differences)
mean_change = paired_cohort["hba1c_change"].mean()
std_change = paired_cohort["hba1c_change"].std()
cohens_d = mean_change / std_change

print(f"=== Pre/Post HbA1c Analysis ===")
print(f"Mean HbA1c change: {mean_change:.3f}%")
print(f"Wilcoxon signed-rank statistic: {stat:.1f}")
print(f"P-value: {p_value:.6f}")
print(f"Cohen's d: {cohens_d:.3f}")
print(f"\nInterpretation: {'Statistically significant' if p_value < 0.05 else 'Not significant'} (alpha=0.05)")
print(f"Effect size: {'Small' if abs(cohens_d) < 0.5 else 'Medium' if abs(cohens_d) < 0.8 else 'Large'} (|d|={abs(cohens_d):.3f})")

# --- Cell 4: Time-in-range and time-below-range ---

# Time-in-range: % of glucose readings between 70-180 mg/dL
# Time-below-range: % of glucose readings below 54 mg/dL

# Merge telemetry with device start dates
telemetry_merged = telemetry.merge(
    demographics[["user_id", "device_start_date"]], on="user_id"
)
telemetry_merged["days_from_start"] = (
    telemetry_merged["timestamp"] - telemetry_merged["device_start_date"]
).dt.days

# Define pre and post periods for telemetry
pre_telemetry = telemetry_merged[
    (telemetry_merged["days_from_start"] >= -LOOKBACK_DAYS) &
    (telemetry_merged["days_from_start"] < 0)
]
post_telemetry = telemetry_merged[
    telemetry_merged["days_from_start"] >= POST_MIN_DAYS
]

def calculate_glycemic_metrics(df):
    """Calculate time-in-range and time-below-range per user."""
    metrics = df.groupby("user_id")["glucose_mg_dl"].agg(
        total_readings="count",
        tir_count=lambda x: ((x >= 70) & (x <= 180)).sum(),
        tbr_count=lambda x: (x < 54).sum()
    ).reset_index()
    metrics["time_in_range_pct"] = (metrics["tir_count"] / metrics["total_readings"]) * 100
    metrics["time_below_range_pct"] = (metrics["tbr_count"] / metrics["total_readings"]) * 100
    return metrics

# Calculate metrics for both periods
pre_metrics = calculate_glycemic_metrics(pre_telemetry)
post_metrics = calculate_glycemic_metrics(post_telemetry)

pre_metrics.columns = ["user_id", "pre_readings", "pre_tir_count", "pre_tbr_count", "pre_tir_pct", "pre_tbr_pct"]
post_metrics.columns = ["user_id", "post_readings", "post_tir_count", "post_tbr_count", "post_tir_pct", "post_tbr_pct"]

print(f"Pre-period telemetry: {len(pre_metrics)} users with glucose data")
print(f"Post-period telemetry: {len(post_metrics)} users with glucose data")

# --- Cell 5: Difference-in-differences ---

# Difference-in-differences approach for users with sufficient data
tir_paired = pre_metrics.merge(post_metrics, on="user_id")
tir_paired["tir_change"] = tir_paired["post_tir_pct"] - tir_paired["pre_tir_pct"]
tir_paired["tbr_change"] = tir_paired["post_tbr_pct"] - tir_paired["pre_tbr_pct"]

print(f"Users with both pre and post telemetry: {len(tir_paired)}")
print(f"\n{'='*50}")
print(f"Time-in-Range (70-180 mg/dL):")
print(f"  Pre:  {tir_paired['pre_tir_pct'].mean():.1f}% +/- {tir_paired['pre_tir_pct'].std():.1f}%")
print(f"  Post: {tir_paired['post_tir_pct'].mean():.1f}% +/- {tir_paired['post_tir_pct'].std():.1f}%")
print(f"  Mean change: {tir_paired['tir_change'].mean():.1f}%")
print(f"\nTime-Below-Range (<54 mg/dL):")
print(f"  Pre:  {tir_paired['pre_tbr_pct'].mean():.2f}% +/- {tir_paired['pre_tbr_pct'].std():.2f}%")
print(f"  Post: {tir_paired['post_tbr_pct'].mean():.2f}% +/- {tir_paired['post_tbr_pct'].std():.2f}%")
print(f"  Mean change: {tir_paired['tbr_change'].mean():.2f}%")

# Statistical tests on TIR and TBR changes
tir_stat, tir_p = stats.wilcoxon(tir_paired["pre_tir_pct"], tir_paired["post_tir_pct"])
tbr_stat, tbr_p = stats.wilcoxon(tir_paired["pre_tbr_pct"], tir_paired["post_tbr_pct"])

print(f"\nWilcoxon test - TIR change: p={tir_p:.6f}")
print(f"Wilcoxon test - TBR change: p={tbr_p:.6f}")

# --- Cell 6: Sensitivity - vary lookback windows ---

windows = [60, 90, 120]
sensitivity_results = []

for window in windows:
    # Redefine pre-period with varied window
    pre_w = clinical_merged[
        (clinical_merged["days_from_start"] >= -window) &
        (clinical_merged["days_from_start"] < 0)
    ]
    pre_w_hba1c = pre_w.groupby("user_id")["hba1c_pct"].last().reset_index()
    pre_w_hba1c.columns = ["user_id", "pre_hba1c"]

    # Post stays the same (90+ days after)
    paired_w = pre_w_hba1c.merge(post_hba1c, on="user_id")
    paired_w["hba1c_change"] = paired_w["post_hba1c"] - paired_w["pre_hba1c"]

    if len(paired_w) > 1:
        stat_w, p_w = stats.wilcoxon(paired_w["pre_hba1c"], paired_w["post_hba1c"])
        d_w = paired_w["hba1c_change"].mean() / paired_w["hba1c_change"].std()
    else:
        stat_w, p_w, d_w = np.nan, np.nan, np.nan

    sensitivity_results.append({
        "lookback_days": window,
        "n_users": len(paired_w),
        "mean_change": paired_w["hba1c_change"].mean() if len(paired_w) > 0 else np.nan,
        "p_value": p_w,
        "cohens_d": d_w
    })

sensitivity_df = pd.DataFrame(sensitivity_results)
print("Sensitivity Analysis: Lookback Window Variation")
print("=" * 55)
print(sensitivity_df.to_string(index=False))

# --- Cell 7: Document missingness ---

users_with_pre = set(pre_hba1c["user_id"])
users_with_post = set(post_hba1c["user_id"])
all_users = set(demographics["user_id"])

missing_pre = all_users - users_with_pre
missing_post = all_users - users_with_post
missing_both = missing_pre & missing_post

print("Missingness Report")
print("=" * 40)
print(f"Total users: {len(all_users)}")
print(f"Users with pre-period HbA1c: {len(users_with_pre)} ({len(users_with_pre)/len(all_users)*100:.1f}%)")
print(f"Users with post-period HbA1c: {len(users_with_post)} ({len(users_with_post)/len(all_users)*100:.1f}%)")
print(f"Users missing pre-period: {len(missing_pre)} ({len(missing_pre)/len(all_users)*100:.1f}%)")
print(f"Users missing post-period: {len(missing_post)} ({len(missing_post)/len(all_users)*100:.1f}%)")
print(f"Users missing both: {len(missing_both)} ({len(missing_both)/len(all_users)*100:.1f}%)")
print(f"\nFinal paired cohort: {len(paired_cohort)} users ({len(paired_cohort)/len(all_users)*100:.1f}%)")

# Check if missingness correlates with demographics
missing_users = demographics[demographics["user_id"].isin(missing_pre | missing_post)]
included_users = demographics[demographics["user_id"].isin(paired_cohort["user_id"])]

print(f"\nDemographic comparison (missing vs. included):")
print(f"  Mean age - Missing: {missing_users['age'].mean():.1f}, Included: {included_users['age'].mean():.1f}")
print(f"  Mean diabetes duration - Missing: {missing_users['diabetes_duration_yrs'].mean():.1f}, Included: {included_users['diabetes_duration_yrs'].mean():.1f}")

# --- Cell 8: Stratification ---

# Stratification by diabetes duration and age
paired_with_demo = paired_cohort.merge(
    demographics[["user_id", "age", "diabetes_duration_yrs"]], on="user_id"
)

# Diabetes duration strata
paired_with_demo["duration_group"] = pd.cut(
    paired_with_demo["diabetes_duration_yrs"],
    bins=[0, 5, 15, 50],
    labels=["Short (0-5 yrs)", "Medium (5-15 yrs)", "Long (15+ yrs)"]
)

# Age strata
paired_with_demo["age_group"] = pd.cut(
    paired_with_demo["age"],
    bins=[0, 30, 50, 100],
    labels=["Young (18-30)", "Middle (31-50)", "Older (51+)"]
)

print("Stratified HbA1c Change by Diabetes Duration")
print("=" * 55)
for group, data in paired_with_demo.groupby("duration_group"):
    if len(data) > 1:
        stat_s, p_s = stats.wilcoxon(data["pre_hba1c"], data["post_hba1c"])
        d_s = data["hba1c_change"].mean() / data["hba1c_change"].std()
        print(f"  {group}: n={len(data)}, mean change={data['hba1c_change'].mean():.3f}, p={p_s:.4f}, d={d_s:.3f}")

print(f"\nStratified HbA1c Change by Age Group")
print("=" * 55)
for group, data in paired_with_demo.groupby("age_group"):
    if len(data) > 1:
        stat_s, p_s = stats.wilcoxon(data["pre_hba1c"], data["post_hba1c"])
        d_s = data["hba1c_change"].mean() / data["hba1c_change"].std()
        print(f"  {group}: n={len(data)}, mean change={data['hba1c_change'].mean():.3f}, p={p_s:.4f}, d={d_s:.3f}")

# Close the database connection
conn.close()
print("\nAnalysis complete. Connection closed.")

Your clinical outcomes notebook now contains a complete pre/post analysis with formal statistical testing, glycemic metrics, and robustness checks. Next, you'll integrate the CRM data to connect clinical outcomes with commercial insights like onboarding speed and support ticket burden.

Integrate CRM and Commercial Analysis

Your safety and clinical analyses revealed how the device performs medically. Now it's time to bring in the business dimension.

In this step, you'll connect CRM data (onboarding timelines, support tickets, account status) to clinical outcomes. The goal is to answer questions a commercial team cares about: does faster onboarding improve retention? Do high-ticket users have worse glycemic control? Which regions need intervention?

In this step, get ready to:

• Build a logistic regression model linking onboarding speed to 90-day retention.
• Analyze support ticket burden by region and severity, correlating tickets with glycemic outcomes.
• Compare fast-onboarding regions vs. slow-onboarding regions on glycemic outcomes and ticket rates.

Build a logistic regression for onboarding speed vs. retention

A device company's commercial team wants to know whether getting users into closed-loop mode faster leads to better retention. Logistic regression is the right tool here because the outcome (retained or not at 90 days) is binary.

• Create a new Jupyter notebook at notebooks/03_crm_commercial_analysis.ipynb by running:

jupyter notebook notebooks/03_crm_commercial_analysis.ipynb

• In the first cell, add the imports and data loading code:

import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
import seaborn as sns
import sqlite3

# Connect to the analytical warehouse
conn = sqlite3.connect("../data/processed/bionic_pancreas.db")

# Load user summary with device metrics and outcomes
user_summary = pd.read_sql("SELECT * FROM user_summary", conn)

# Load CRM data for onboarding and ticket details
crm_data = pd.read_sql("SELECT * FROM crm_data", conn)

print(f"User summary: {len(user_summary)} rows")
print(f"CRM data: {len(crm_data)} rows")

What does this code do?

This cell connects to the same SQLite warehouse your ETL built in Step 4. It loads two tables: user_summary (the joined analytical layer with device metrics and clinical outcomes per user) and crm_data (raw CRM records with onboarding dates, tickets, and account status).

• Run the cell. You should see row counts printed for both DataFrames (500 for user_summary, approximately 1750 for crm_data).

• In the next cell, compute the onboarding speed metric and fit the logistic regression:

# Aggregate CRM data to one row per user
onboarding = crm_data.groupby("user_id").agg(
    first_closed_loop_date=("first_closed_loop_date", "first"),
    account_status=("account_status", "first"),
    region=("region", "first")
).reset_index()

# Merge with user_summary to get device_start_date and time_in_range
onboarding = onboarding.merge(
    user_summary[["user_id", "device_start_date", "time_in_range_pct"]],
    on="user_id",
    how="left"
)

# Calculate days from device start to first closed-loop session
onboarding["device_start_date"] = pd.to_datetime(onboarding["device_start_date"])
onboarding["first_closed_loop_date"] = pd.to_datetime(onboarding["first_closed_loop_date"])
onboarding["days_to_closed_loop"] = (
    onboarding["first_closed_loop_date"] - onboarding["device_start_date"]
).dt.days

# Define 90-day retention as binary outcome
onboarding["retained_90d"] = (onboarding["account_status"] == "active").astype(int)

print(f"Median days to closed loop: {onboarding['days_to_closed_loop'].median():.1f}")
print(f"Retention rate: {onboarding['retained_90d'].mean():.1%}")

What does this code do?

The days_to_closed_loop column measures how quickly each user transitioned from receiving the device to running in fully automated mode. This is the predictor variable for our regression.

The retained_90d column is the binary outcome: 1 if the user's account is still active, 0 otherwise.

• Run the cell. You should see the median days to closed loop and overall retention rate printed.

• In the next cell, fit the logistic regression model:

# Fit logistic regression: does onboarding speed predict retention?
model_data = onboarding[["days_to_closed_loop", "retained_90d"]].dropna()

# Add constant for intercept term
X = sm.add_constant(model_data["days_to_closed_loop"])
y = model_data["retained_90d"]

# Fit the Logit model using maximum likelihood estimation
logit_model = sm.Logit(y, X).fit()
print(logit_model.summary())

What does this code do?

The sm.Logit() function from statsmodels fits a logistic regression using maximum likelihood estimation. The summary() output shows the coefficient, standard error, z-statistic, p-value, and 95% confidence interval for days_to_closed_loop.

A negative coefficient means longer onboarding time reduces the probability of retention. The p-value tells you whether this relationship is statistically significant.

• Run the cell. You should see a Logit Regression Results table with the model summary, including the coefficient for days_to_closed_loop, its p-value, and pseudo R-squared.

Seeing a convergence warning?

If the model prints a convergence warning, it usually means there isn't enough variation in the outcome variable or the predictor range is very narrow. Verify that retained_90d has both 0s and 1s by checking model_data["retained_90d"].value_counts().

If fewer than 10 users churned, the synthetic data may need regeneration with a different seed. Help me debug my logistic regression convergence issue

Analyze support ticket burden by region and severity

High support ticket volume can indicate onboarding friction or device usability problems. If users who file more tickets also have worse glycemic control, that suggests support isn't resolving the underlying clinical issue.

• In a new cell, compute the ticket burden by region and severity, then correlate ticket count with time-in-range:

# Count tickets per user and capture their region
ticket_counts = crm_data.groupby("user_id").agg(
    ticket_count=("ticket_id", "count"),
    region=("region", "first"),
    mean_severity=("ticket_severity", "mean")
).reset_index()

# Merge with glycemic outcomes from user_summary
ticket_outcomes = ticket_counts.merge(
    user_summary[["user_id", "time_in_range_pct"]],
    on="user_id",
    how="left"
)

# Correlation between ticket count and time-in-range percentage
correlation = ticket_outcomes["ticket_count"].corr(ticket_outcomes["time_in_range_pct"])
print(f"Correlation (ticket count vs. time-in-range): {correlation:.3f}")

What does this correlation tell you?

A negative correlation means users who file more support tickets tend to spend less time in the target glucose range (70-180 mg/dL). This is a signal that support volume may be a proxy for clinical struggle, not just administrative friction.

A near-zero correlation would suggest tickets and glycemic outcomes are independent. Either result is useful for stakeholder recommendations.

• Run the cell. You should see the Pearson correlation coefficient printed.

• In the next cell, break down ticket burden by region and severity for a more granular view:

# Ticket burden by region and severity
region_severity = crm_data.groupby(["region", "ticket_severity"]).agg(
    ticket_count=("ticket_id", "count"),
    mean_resolution_days=("ticket_resolution_days", "mean")
).reset_index()

# Pivot to see severity distribution across regions
print("Ticket count by region and severity:")
print(region_severity.pivot_table(
    index="region", columns="ticket_severity",
    values="ticket_count", aggfunc="sum"
).fillna(0))

# Visualize the relationship between tickets and glycemic outcomes
fig, ax = plt.subplots(figsize=(8, 5))
sns.scatterplot(
    data=ticket_outcomes, x="ticket_count", y="time_in_range_pct",
    hue="region", alpha=0.6, ax=ax
)
ax.set_xlabel("Support Ticket Count")
ax.set_ylabel("Time in Range (%)")
ax.set_title("Support Tickets vs. Glycemic Outcomes by Region")
plt.tight_layout()
plt.show()

What does this code do?

The pivot table shows how ticket volume distributes across severity levels in each region. Regions with disproportionately high severe tickets may need targeted intervention.

The scatter plot visualizes whether the correlation you computed is driven by a few outliers or is a consistent pattern across users and regions.

• Run the cell. You should see a pivot table of ticket counts by region and severity, followed by a scatter plot with colored dots for each region.

Scatter plot looks like a single blob?

If the plot shows all points clustered together, try adjusting the alpha (transparency) or reducing point size with s=20 in the scatterplot() call. Synthetic data sometimes produces tight distributions.

Help me improve my scatter plot visualization for ticket count vs. time-in-range

Build a cohort comparison of fast vs. slow onboarding regions

The final analysis segments regions into two groups: those where users reach closed-loop mode faster than the median, and those where onboarding is slower. This directly answers whether regional onboarding differences translate into clinical or support outcomes.

• In a new cell, define the fast and slow cohorts and compare their outcomes:

# Calculate median days_to_closed_loop per region
region_onboarding = onboarding.groupby("region")["days_to_closed_loop"].median().reset_index()
region_onboarding.columns = ["region", "median_days_to_closed_loop"]

# Overall median as the threshold for fast vs. slow
overall_median = onboarding["days_to_closed_loop"].median()
print(f"Overall median days to closed loop: {overall_median:.1f}")

# Classify regions as fast or slow onboarding
region_onboarding["onboarding_speed"] = np.where(
    region_onboarding["median_days_to_closed_loop"] < overall_median,
    "fast", "slow"
)
print(f"\nRegion classifications:")
print(region_onboarding.to_string(index=False))

Why use the overall median as threshold?

The median splits regions into roughly equal groups, making the comparison balanced. A region is "fast" if its typical user reaches closed-loop mode quicker than the fleet-wide median. This relative approach avoids arbitrary cutoffs and adapts to the data's natural distribution.

• Run the cell. You should see the overall median printed, followed by a table showing each region's median onboarding days and its fast/slow classification.

• In the final cell, merge the cohort labels back to user-level data and compare outcomes:

# Merge cohort classification to user-level data
cohort_data = onboarding.merge(
    region_onboarding[["region", "onboarding_speed"]],
    on="region",
    how="left"
)

# Merge ticket counts for each user
cohort_data = cohort_data.merge(
    ticket_counts[["user_id", "ticket_count"]],
    on="user_id",
    how="left"
)

# Compare cohorts on key metrics
cohort_comparison = cohort_data.groupby("onboarding_speed").agg(
    mean_time_in_range=("time_in_range_pct", "mean"),
    mean_ticket_count=("ticket_count", "mean"),
    retention_rate=("retained_90d", "mean"),
    n_users=("user_id", "count")
).reset_index()

print("\nCohort comparison (fast vs. slow onboarding regions):")
print(cohort_comparison.to_string(index=False))

# Close the database connection
conn.close()

How to interpret the cohort comparison

The table shows four metrics side by side for fast vs. slow regions: mean time-in-range (higher is better), mean ticket count (lower suggests fewer support issues), retention rate (proportion of users still active), and sample size.

If fast-onboarding regions show better glycemic outcomes AND fewer tickets, that's strong evidence for investing in onboarding acceleration programs.

• Run the cell. You should see a comparison table with columns for onboarding_speed, mean_time_in_range, mean_ticket_count, retention_rate, and n_users.

• Save the notebook with Cmd+S (macOS) or Ctrl+S (Windows).

Getting a KeyError on a column name?

Make sure the column names in your groupby() and merge() calls match the exact column names in your database. Run crm_data.columns.tolist() and user_summary.columns.tolist() to verify.

Help me debug a KeyError in my CRM commercial analysis notebook

✔️ Awesome, I've got everything!

Great. Double check that you have saved notebooks/03_crm_commercial_analysis.ipynb and all cells ran without errors.

ⓧ I'd like to double check the full code

Here is the complete notebook code across all cells:

# Cell 1: Imports and data loading
import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
import seaborn as sns
import sqlite3

# Connect to the analytical warehouse
conn = sqlite3.connect("../data/processed/bionic_pancreas.db")

# Load user summary with device metrics and outcomes
user_summary = pd.read_sql("SELECT * FROM user_summary", conn)

# Load CRM data for onboarding and ticket details
crm_data = pd.read_sql("SELECT * FROM crm_data", conn)

print(f"User summary: {len(user_summary)} rows")
print(f"CRM data: {len(crm_data)} rows")

# Cell 2: Compute onboarding speed metric
# Aggregate CRM data to one row per user
onboarding = crm_data.groupby("user_id").agg(
    first_closed_loop_date=("first_closed_loop_date", "first"),
    account_status=("account_status", "first"),
    region=("region", "first")
).reset_index()

# Merge with user_summary to get device_start_date and time_in_range
onboarding = onboarding.merge(
    user_summary[["user_id", "device_start_date", "time_in_range_pct"]],
    on="user_id",
    how="left"
)

# Calculate days from device start to first closed-loop session
onboarding["device_start_date"] = pd.to_datetime(onboarding["device_start_date"])
onboarding["first_closed_loop_date"] = pd.to_datetime(onboarding["first_closed_loop_date"])
onboarding["days_to_closed_loop"] = (
    onboarding["first_closed_loop_date"] - onboarding["device_start_date"]
).dt.days

# Define 90-day retention as binary outcome
onboarding["retained_90d"] = (onboarding["account_status"] == "active").astype(int)

print(f"Median days to closed loop: {onboarding['days_to_closed_loop'].median():.1f}")
print(f"Retention rate: {onboarding['retained_90d'].mean():.1%}")

# Cell 3: Logistic regression - onboarding speed predicting retention
model_data = onboarding[["days_to_closed_loop", "retained_90d"]].dropna()

# Add constant for intercept term
X = sm.add_constant(model_data["days_to_closed_loop"])
y = model_data["retained_90d"]

# Fit the Logit model using maximum likelihood estimation
logit_model = sm.Logit(y, X).fit()
print(logit_model.summary())

# Cell 4: Support ticket burden - count per user and correlate with outcomes
ticket_counts = crm_data.groupby("user_id").agg(
    ticket_count=("ticket_id", "count"),
    region=("region", "first"),
    mean_severity=("ticket_severity", "mean")
).reset_index()

# Merge with glycemic outcomes from user_summary
ticket_outcomes = ticket_counts.merge(
    user_summary[["user_id", "time_in_range_pct"]],
    on="user_id",
    how="left"
)

# Correlation between ticket count and time-in-range percentage
correlation = ticket_outcomes["ticket_count"].corr(ticket_outcomes["time_in_range_pct"])
print(f"Correlation (ticket count vs. time-in-range): {correlation:.3f}")

# Cell 5: Ticket burden by region and severity with visualization
region_severity = crm_data.groupby(["region", "ticket_severity"]).agg(
    ticket_count=("ticket_id", "count"),
    mean_resolution_days=("ticket_resolution_days", "mean")
).reset_index()

# Pivot to see severity distribution across regions
print("Ticket count by region and severity:")
print(region_severity.pivot_table(
    index="region", columns="ticket_severity",
    values="ticket_count", aggfunc="sum"
).fillna(0))

# Visualize the relationship between tickets and glycemic outcomes
fig, ax = plt.subplots(figsize=(8, 5))
sns.scatterplot(
    data=ticket_outcomes, x="ticket_count", y="time_in_range_pct",
    hue="region", alpha=0.6, ax=ax
)
ax.set_xlabel("Support Ticket Count")
ax.set_ylabel("Time in Range (%)")
ax.set_title("Support Tickets vs. Glycemic Outcomes by Region")
plt.tight_layout()
plt.show()

# Cell 6: Define fast vs. slow onboarding regions
region_onboarding = onboarding.groupby("region")["days_to_closed_loop"].median().reset_index()
region_onboarding.columns = ["region", "median_days_to_closed_loop"]

# Overall median as the threshold for fast vs. slow
overall_median = onboarding["days_to_closed_loop"].median()
print(f"Overall median days to closed loop: {overall_median:.1f}")

# Classify regions as fast or slow onboarding
region_onboarding["onboarding_speed"] = np.where(
    region_onboarding["median_days_to_closed_loop"] < overall_median,
    "fast", "slow"
)
print(f"\nRegion classifications:")
print(region_onboarding.to_string(index=False))

# Cell 7: Cohort comparison - fast vs. slow onboarding regions
# Merge cohort classification to user-level data
cohort_data = onboarding.merge(
    region_onboarding[["region", "onboarding_speed"]],
    on="region",
    how="left"
)

# Merge ticket counts for each user
cohort_data = cohort_data.merge(
    ticket_counts[["user_id", "ticket_count"]],
    on="user_id",
    how="left"
)

# Compare cohorts on key metrics
cohort_comparison = cohort_data.groupby("onboarding_speed").agg(
    mean_time_in_range=("time_in_range_pct", "mean"),
    mean_ticket_count=("ticket_count", "mean"),
    retention_rate=("retained_90d", "mean"),
    n_users=("user_id", "count")
).reset_index()

print("\nCohort comparison (fast vs. slow onboarding regions):")
print(cohort_comparison.to_string(index=False))

# Close the database connection
conn.close()

You now have a complete commercial picture linking onboarding speed, support burden, and clinical outcomes across regions. Next up, you'll wire all of these analyses into an automated Prefect pipeline so the whole thing runs on a schedule.

Orchestrate the Pipeline with Prefect

You have a working data generator, ETL scripts, and three analytical notebooks. Right now, running the full pipeline means executing each piece manually in the right order.

In this step, you'll wrap everything into a single Prefect flow so that one command regenerates data, loads it into the warehouse, runs all analyses, and exports results.

In this step, get ready to:

• Convert the ETL and analysis pipeline into a Prefect flow with @task and @flow decorators.
• Add a weekly cron schedule to simulate automated safety reporting.
• Run the pipeline end-to-end and verify it in the Prefect UI dashboard.

Create the Prefect pipeline file

Prefect uses pure Python decorators to turn ordinary functions into observable, recoverable pipeline stages. The @task decorator marks individual units of work. The @flow decorator marks the orchestrator that calls them in order.

• Create a new file called src/pipeline.py by running:

touch src/pipeline.py

• Confirm the file exists by running:

ls src/pipeline.py

You should see src/pipeline.py listed in the output.

• Open src/pipeline.py in your editor and add the imports and the ingest task:

"""Prefect pipeline orchestrating end-to-end analytics for bionic pancreas data."""
import subprocess
import sqlite3
from pathlib import Path

import pandas as pd
from prefect import flow, task

# Project paths
PROJECT_ROOT = Path(__file__).resolve().parent.parent
DATA_RAW = PROJECT_ROOT / "data" / "raw"
DATA_PROCESSED = PROJECT_ROOT / "data" / "processed"
DB_PATH = DATA_PROCESSED / "bionic_pancreas.db"
NOTEBOOKS = PROJECT_ROOT / "notebooks"


@task(log_prints=True)
def ingest_task():
    """Generate synthetic data by running the data generator."""
    print("Generating synthetic data...")
    subprocess.run(
        ["python", str(PROJECT_ROOT / "src" / "generate_data.py"), "--seed", "42"],
        check=True
    )
    print("Synthetic data generated successfully.")

What does this code do?

The imports bring in Prefect's flow and task decorators, plus standard libraries for running subprocesses and reading the database.

The ingest_task calls generate_data.py with the --seed 42 flag to produce reproducible synthetic datasets. Prefect tracks this as a discrete pipeline stage with its own logs and state.

• Save src/pipeline.py.

Seeing an import error on prefect?

Make sure your virtual environment is active. Run pip install prefect if the import fails.

help me fix my Prefect import error

• Below ingest_task, add the ETL and analysis tasks:

@task(log_prints=True)
def clean_and_load_task():
    """Run the ETL pipeline to clean data and load into SQLite."""
    print("Running ETL pipeline...")
    # Import and run the ETL module
    import sys
    sys.path.insert(0, str(PROJECT_ROOT / "src"))
    from etl import run_etl
    run_etl()
    print("ETL complete. Database updated.")


@task(log_prints=True)
def analyze_safety_task():
    """Execute the safety signal detection notebook."""
    print("Running safety signal detection analysis...")
    notebook = NOTEBOOKS / "01_safety_signal_detection.ipynb"
    subprocess.run(
        ["jupyter", "nbconvert", "--to", "notebook", "--execute", str(notebook)],
        check=True
    )
    print("Safety analysis complete.")


@task(log_prints=True)
def analyze_clinical_task():
    """Execute the clinical outcomes analysis notebook."""
    print("Running clinical outcomes analysis...")
    notebook = NOTEBOOKS / "02_clinical_outcomes.ipynb"
    subprocess.run(
        ["jupyter", "nbconvert", "--to", "notebook", "--execute", str(notebook)],
        check=True
    )
    print("Clinical analysis complete.")


@task(log_prints=True)
def analyze_crm_task():
    """Execute the CRM commercial analysis notebook."""
    print("Running CRM commercial analysis...")
    notebook = NOTEBOOKS / "03_crm_commercial_analysis.ipynb"
    subprocess.run(
        ["jupyter", "nbconvert", "--to", "notebook", "--execute", str(notebook)],
        check=True
    )
    print("CRM analysis complete.")

What does this code do?

The clean_and_load_task imports your existing run_etl() function from src/etl.py and executes it. This reads the raw CSVs, cleans them, and loads everything into bionic_pancreas.db.

Each analyze_* task uses jupyter nbconvert --to notebook --execute to programmatically run a notebook. This is the standard way to execute Jupyter notebooks in an automated pipeline without opening a browser.

• Save src/pipeline.py.

• Add the export task and the main flow below the analysis tasks:

@task(log_prints=True)
def export_task():
    """Export the user_summary table to CSV for dashboarding."""
    print("Exporting user_summary to CSV...")
    conn = sqlite3.connect(str(DB_PATH))
    df = pd.read_sql_query("SELECT * FROM user_summary", conn)
    conn.close()
    output_path = DATA_PROCESSED / "user_summary_export.csv"
    df.to_csv(output_path, index=False)
    print(f"Exported {len(df)} rows to {output_path}")


@flow(name="bionic-pancreas-pipeline", log_prints=True)
def main_flow():
    """Orchestrate the full analytics pipeline."""
    # Stage 1: Generate synthetic data
    ingest_task()

    # Stage 2: Clean and load into warehouse
    clean_and_load_task()

    # Stage 3: Run analyses
    analyze_safety_task()
    analyze_clinical_task()
    analyze_crm_task()

    # Stage 4: Export results
    export_task()

    print("Pipeline complete. All stages finished successfully.")


if __name__ == "__main__":
    main_flow()

What does this code do?

The export_task reads the user_summary table from your SQLite database and writes it to data/processed/user_summary_export.csv. This CSV feeds Tableau in the next step.

The main_flow function decorated with @flow is the orchestrator. It calls each task in sequence. Prefect builds a dependency graph and tracks the state of every task run, so if one fails you can see exactly which stage broke.

Why Prefect over Airflow?

Prefect uses pure Python decorators and requires no separate scheduler daemon, metadata database, or webserver. You can run the pipeline locally with a single python pipeline.py command.

Airflow requires standing up multiple processes and a configuration layer before you can run anything. For a portfolio project that a reviewer should be able to clone and run immediately, Prefect demonstrates the same orchestration concepts with far less infrastructure overhead.

• Save src/pipeline.py.
• Run the pipeline end-to-end by executing:

python src/pipeline.py

You should see Prefect log output for each task stage, ending with a message confirming the pipeline completed successfully.

Pipeline failing on notebook execution?

If you see a timeout error, notebooks default to a 30-second cell timeout. Add --ExecutePreprocessor.timeout=600 to the subprocess.run arguments in each analysis task to allow longer-running cells.

If you see an import error from etl.py, make sure you are running the command from your project root directory (the folder containing src/).

help me debug my Prefect pipeline execution errors

• Confirm the export file was created by running:

ls data/processed/user_summary_export.csv

You should see data/processed/user_summary_export.csv listed in the output.

Add a weekly schedule

In a production environment, a safety reporting pipeline would run automatically on a recurring schedule. Prefect's .serve() method lets you attach a cron schedule to any flow without deploying to a remote server.

• In src/pipeline.py, find the if __name__ == "__main__": block at the bottom of the file.
• Find this line:

if __name__ == "__main__":
    main_flow()

• Replace main_flow() with the .serve() method, so the block now looks like this:

if __name__ == "__main__":
    main_flow.serve(name="weekly-safety-pipeline", cron="0 8 * * 1")

What does this schedule mean?

The cron expression 0 8 * * 1 means "run at 8:00 AM every Monday." This simulates a weekly automated safety report cadence typical in medical device analytics.

When you run the script, Prefect starts a lightweight process that listens for scheduled triggers. Press Ctrl+C to stop it when you are done.

• Save src/pipeline.py.

Run the pipeline and verify in the Prefect UI

Prefect includes a local dashboard where you can see flow runs, task states, and logs. You will start the Prefect server in one terminal, then trigger your pipeline from another.

• Open a new terminal window (keep it separate from your main terminal).
• Start the Prefect local server by running:

prefect server start

You should see output confirming the server is running, with a URL (typically http://127.0.0.1:4200) where you can access the dashboard.

• Open http://127.0.0.1:4200 in your browser. You should see the Prefect dashboard.
• Switch back to your original terminal.
• Run the pipeline with the serve command to register the deployment and trigger an initial run:

python src/pipeline.py

The serve process starts and registers your weekly schedule. Prefect triggers the first run immediately.

• Switch to your browser and refresh the Prefect dashboard. You should see your bionic-pancreas-pipeline flow run appear with all six task stages listed with green completed states.
• Press Ctrl+C in your original terminal to stop the serve process.

Don't see your flow run in the dashboard?

Make sure prefect server start is running in a separate terminal before you execute python src/pipeline.py. The pipeline needs the server to be running so it can report its state.

If the dashboard shows the flow as "Late" rather than running, the scheduled time has not arrived yet. You can trigger a one-off execution by temporarily changing the if __name__ block back to main_flow() and running the script directly.

help me troubleshoot my Prefect server and dashboard connection

✔️ Awesome, I've got everything!

Great. Double check you have saved src/pipeline.py and confirmed the pipeline ran successfully in the Prefect dashboard.

ⓧ I'd like to double check the full code

Here is the complete src/pipeline.py file:

"""Prefect pipeline orchestrating end-to-end analytics for bionic pancreas data."""
import subprocess
import sqlite3
from pathlib import Path

import pandas as pd
from prefect import flow, task

# Project paths
PROJECT_ROOT = Path(__file__).resolve().parent.parent
DATA_RAW = PROJECT_ROOT / "data" / "raw"
DATA_PROCESSED = PROJECT_ROOT / "data" / "processed"
DB_PATH = DATA_PROCESSED / "bionic_pancreas.db"
NOTEBOOKS = PROJECT_ROOT / "notebooks"


@task(log_prints=True)
def ingest_task():
    """Generate synthetic data by running the data generator."""
    print("Generating synthetic data...")
    subprocess.run(
        ["python", str(PROJECT_ROOT / "src" / "generate_data.py"), "--seed", "42"],
        check=True
    )
    print("Synthetic data generated successfully.")


@task(log_prints=True)
def clean_and_load_task():
    """Run the ETL pipeline to clean data and load into SQLite."""
    print("Running ETL pipeline...")
    # Import and run the ETL module
    import sys
    sys.path.insert(0, str(PROJECT_ROOT / "src"))
    from etl import run_etl
    run_etl()
    print("ETL complete. Database updated.")


@task(log_prints=True)
def analyze_safety_task():
    """Execute the safety signal detection notebook."""
    print("Running safety signal detection analysis...")
    notebook = NOTEBOOKS / "01_safety_signal_detection.ipynb"
    subprocess.run(
        ["jupyter", "nbconvert", "--to", "notebook", "--execute", str(notebook)],
        check=True
    )
    print("Safety analysis complete.")


@task(log_prints=True)
def analyze_clinical_task():
    """Execute the clinical outcomes analysis notebook."""
    print("Running clinical outcomes analysis...")
    notebook = NOTEBOOKS / "02_clinical_outcomes.ipynb"
    subprocess.run(
        ["jupyter", "nbconvert", "--to", "notebook", "--execute", str(notebook)],
        check=True
    )
    print("Clinical analysis complete.")


@task(log_prints=True)
def analyze_crm_task():
    """Execute the CRM commercial analysis notebook."""
    print("Running CRM commercial analysis...")
    notebook = NOTEBOOKS / "03_crm_commercial_analysis.ipynb"
    subprocess.run(
        ["jupyter", "nbconvert", "--to", "notebook", "--execute", str(notebook)],
        check=True
    )
    print("CRM analysis complete.")


@task(log_prints=True)
def export_task():
    """Export the user_summary table to CSV for dashboarding."""
    print("Exporting user_summary to CSV...")
    conn = sqlite3.connect(str(DB_PATH))
    df = pd.read_sql_query("SELECT * FROM user_summary", conn)
    conn.close()
    output_path = DATA_PROCESSED / "user_summary_export.csv"
    df.to_csv(output_path, index=False)
    print(f"Exported {len(df)} rows to {output_path}")


@flow(name="bionic-pancreas-pipeline", log_prints=True)
def main_flow():
    """Orchestrate the full analytics pipeline."""
    # Stage 1: Generate synthetic data
    ingest_task()

    # Stage 2: Clean and load into warehouse
    clean_and_load_task()

    # Stage 3: Run analyses
    analyze_safety_task()
    analyze_clinical_task()
    analyze_crm_task()

    # Stage 4: Export results
    export_task()

    print("Pipeline complete. All stages finished successfully.")


if __name__ == "__main__":
    main_flow.serve(name="weekly-safety-pipeline", cron="0 8 * * 1")

Your pipeline is now fully orchestrated. In the next step, you will write tests to validate your data transforms and statistical functions.

Write Tests and Validation Report

Your pipeline runs end-to-end and produces analytical outputs. But how do you know the underlying data is correct?

In this step, you will write automated tests that validate your ETL output and statistical functions. You will also create a validation report that documents data quality checks and known limitations for anyone reviewing your repository.

In this step, get ready to:

• Write pytest unit tests for ETL data quality (row counts, null checks, value ranges).
• Write pytest unit tests for statistical functions (Poisson CI, rate computation, cohort assignment).
• Create a validation report documenting test results, quality checks, and limitations.

Write ETL data quality tests

Automated tests catch regressions early. If someone re-runs the pipeline with new parameters and the user count changes or null values creep in, these tests surface the problem immediately.

• Create the file tests/test_etl.py inside your tests/ folder by pasting the following code:

import sqlite3
import os
import pandas as pd

# Path to the processed database
DB_PATH = os.path.join("data", "processed", "bionic_pancreas.db")


def get_connection():
    """Return a connection to the analytics database."""
    return sqlite3.connect(DB_PATH)


def test_user_summary_row_count():
    """Verify user_summary has exactly 500 rows (one per user)."""
    conn = get_connection()
    df = pd.read_sql("SELECT * FROM user_summary", conn)
    conn.close()
    assert len(df) == 500, f"Expected 500 rows, got {len(df)}"


def test_no_null_user_ids():
    """Verify no null user_ids in the user_summary table."""
    conn = get_connection()
    df = pd.read_sql(
        "SELECT user_id FROM user_summary WHERE user_id IS NULL", conn
    )
    conn.close()
    assert len(df) == 0, f"Found {len(df)} null user_ids"


def test_glucose_range():
    """Verify all glucose values fall within physiological range (20-600 mg/dL)."""
    conn = get_connection()
    df = pd.read_sql(
        "SELECT glucose_mg_dl FROM device_telemetry "
        "WHERE glucose_mg_dl < 20 OR glucose_mg_dl > 600",
        conn
    )
    conn.close()
    assert len(df) == 0, f"Found {len(df)} out-of-range glucose values"


def test_date_range():
    """Verify device_start_date is not null for any user."""
    conn = get_connection()
    df = pd.read_sql(
        "SELECT user_id FROM demographics_enrollment "
        "WHERE device_start_date IS NULL",
        conn
    )
    conn.close()
    assert len(df) == 0, f"Found {len(df)} users with null device_start_date"

What does this code do?

Each function prefixed with test_ is automatically discovered by pytest. The tests query the SQLite database directly and assert data quality expectations.

• test_user_summary_row_count confirms exactly 500 users exist after all joins.
• test_no_null_user_ids catches any broken joins that would produce null keys.
• test_glucose_range ensures glucose readings stay within the physiological range of 20-600 mg/dL.
• test_date_range ensures every enrolled user has a valid device start date.

• Save tests/test_etl.py.
• Run the ETL tests by running this command from the project root:

pytest tests/test_etl.py -v

• You should see 4 tests pass, each showing PASSED next to the test function name.

Seeing test failures?

• If test_user_summary_row_count fails, confirm you ran the full pipeline at least once so that user_summary is populated in the database.
• If test_glucose_range fails, check your data generator's glucose bounds in src/generate_data.py.
• Make sure you are running pytest from the project root folder so the relative DB_PATH resolves correctly.

Help me debug my failing ETL tests.

Write analysis function tests

Testing statistical functions ensures your Poisson confidence interval formula and rate calculations produce correct results. These tests verify the math independently from the data.

• Create the file tests/test_analysis.py inside your tests/ folder by pasting the following code:

from scipy.stats import chi2


def poisson_ci(events, exposure, alpha=0.05):
    """Calculate Poisson confidence interval for an event rate.

    Uses the chi-squared distribution to derive exact bounds.
    """
    rate = events / exposure
    # Lower bound using chi-squared inverse CDF
    lo = chi2.ppf(alpha / 2, 2 * events) / (2 * exposure)
    # Upper bound using chi-squared inverse CDF
    hi = chi2.ppf(1 - alpha / 2, 2 * (events + 1)) / (2 * exposure)
    return rate, lo, hi


def test_poisson_ci_known_value():
    """Verify Poisson CI calculation against a known value."""
    rate, lo, hi = poisson_ci(10, 1000)
    assert rate == 0.01, f"Expected rate 0.01, got {rate}"
    assert lo < rate, f"Lower bound {lo} should be less than rate {rate}"
    assert hi > rate, f"Upper bound {hi} should be greater than rate {rate}"


def test_rate_computation():
    """Verify rate per 1,000 device-days computation."""
    events = 5
    device_days = 500
    rate_per_1000 = (events / device_days) * 1000
    assert rate_per_1000 == 10.0, f"Expected 10.0, got {rate_per_1000}"


def test_cohort_assignment():
    """Verify cohort assignment logic based on pct_closed_loop."""
    pct_closed_loop = 50
    cohort = "closed_loop" if pct_closed_loop >= 50 else "open_loop"
    assert cohort == "closed_loop", f"Expected closed_loop, got {cohort}"

What does this code do?

The poisson_ci helper reimplements the Poisson confidence interval formula using scipy.stats.chi2. This is the same Garwood interval method you used in your safety signal notebook.

• test_poisson_ci_known_value passes 10 events over 1,000 exposure units and confirms the rate is 0.01 with the lower bound below and upper bound above.
• test_rate_computation verifies the "events per 1,000 device-days" arithmetic directly.
• test_cohort_assignment confirms the threshold logic that maps pct_closed_loop >= 50 to the closed_loop cohort.

• Save tests/test_analysis.py.
• Run all 7 tests across both files by running:

pytest tests/ -v

• You should see 7 tests pass (4 from test_etl.py and 3 from test_analysis.py), each showing PASSED.

Seeing an import error for scipy?

Make sure your virtual environment is activated before running pytest. The scipy package must be installed in the active environment.

Help me fix my pytest import errors.

Create the validation report

A validation report gives reviewers confidence that you checked your own work. It documents what passed, what the data looks like in aggregate, and where the limitations are.

• Create the file docs/validation_report.md inside your docs/ folder by pasting the following:

# Validation Report

## Test Results

All 7 tests pass across two test modules:

### test_etl.py (4 tests)
- `test_user_summary_row_count`: Confirms 500 rows in user_summary (one per synthetic user).
- `test_no_null_user_ids`: Confirms zero null user_id values in the final table.
- `test_glucose_range`: Confirms all glucose readings fall within 20-600 mg/dL.
- `test_date_range`: Confirms no null device_start_date values.

### test_analysis.py (3 tests)
- `test_poisson_ci_known_value`: Verifies Poisson CI math against known input (10 events, 1000 exposure).
- `test_rate_computation`: Verifies rate-per-1000 calculation (5 events / 500 days * 1000 = 10.0).
- `test_cohort_assignment`: Verifies pct_closed_loop >= 50 maps to closed_loop cohort.

## Data Quality Checks

| Metric | Expected | Observed |
|--------|----------|----------|
| Total users | 500 | 500 |
| Null user_ids | 0 | 0 |
| Glucose range | 20-600 mg/dL | Confirmed |
| HbA1c range | 5.0-14.0% | Confirmed |
| Severe hypo rate | ~2-5% | ~3.5% |
| DKA rate | <1% | ~0.8% |

## Known Limitations

1. **Synthetic data**: All data is generated. Patterns are illustrative, not derived from real patients.
2. **Simulated pre-device HbA1c**: Pre-device baseline values are generated independently, not from longitudinal records.
3. **Logistic regression assumptions**: The onboarding-retention model assumes a linear log-odds relationship without accounting for all potential confounders.
4. **Illustrative findings only**: Statistical results demonstrate methodology, not real clinical evidence.

Why document known limitations?

In medical device analytics, transparency about data provenance and methodological assumptions is critical. Listing limitations upfront shows a reviewer that you understand the boundaries of your analysis.

A hiring manager reading this report sees that you can self-audit your work. That is a skill that separates junior analysts from senior ones.

• Save docs/validation_report.md.
• Confirm the file exists by running:

ls docs/validation_report.md

• You should see docs/validation_report.md in the output.

✔️ Awesome, I've got everything!

Great. Double check that all three new files are saved.

ⓧ I'd like to double check the full code

Here is the complete tests/test_etl.py:

import sqlite3
import os
import pandas as pd

# Path to the processed database
DB_PATH = os.path.join("data", "processed", "bionic_pancreas.db")


def get_connection():
    """Return a connection to the analytics database."""
    return sqlite3.connect(DB_PATH)


def test_user_summary_row_count():
    """Verify user_summary has exactly 500 rows (one per user)."""
    conn = get_connection()
    df = pd.read_sql("SELECT * FROM user_summary", conn)
    conn.close()
    assert len(df) == 500, f"Expected 500 rows, got {len(df)}"


def test_no_null_user_ids():
    """Verify no null user_ids in the user_summary table."""
    conn = get_connection()
    df = pd.read_sql(
        "SELECT user_id FROM user_summary WHERE user_id IS NULL", conn
    )
    conn.close()
    assert len(df) == 0, f"Found {len(df)} null user_ids"


def test_glucose_range():
    """Verify all glucose values fall within physiological range (20-600 mg/dL)."""
    conn = get_connection()
    df = pd.read_sql(
        "SELECT glucose_mg_dl FROM device_telemetry "
        "WHERE glucose_mg_dl < 20 OR glucose_mg_dl > 600",
        conn
    )
    conn.close()
    assert len(df) == 0, f"Found {len(df)} out-of-range glucose values"


def test_date_range():
    """Verify device_start_date is not null for any user."""
    conn = get_connection()
    df = pd.read_sql(
        "SELECT user_id FROM demographics_enrollment "
        "WHERE device_start_date IS NULL",
        conn
    )
    conn.close()
    assert len(df) == 0, f"Found {len(df)} users with null device_start_date"

Here is the complete tests/test_analysis.py:

from scipy.stats import chi2


def poisson_ci(events, exposure, alpha=0.05):
    """Calculate Poisson confidence interval for an event rate.

    Uses the chi-squared distribution to derive exact bounds.
    """
    rate = events / exposure
    # Lower bound using chi-squared inverse CDF
    lo = chi2.ppf(alpha / 2, 2 * events) / (2 * exposure)
    # Upper bound using chi-squared inverse CDF
    hi = chi2.ppf(1 - alpha / 2, 2 * (events + 1)) / (2 * exposure)
    return rate, lo, hi


def test_poisson_ci_known_value():
    """Verify Poisson CI calculation against a known value."""
    rate, lo, hi = poisson_ci(10, 1000)
    assert rate == 0.01, f"Expected rate 0.01, got {rate}"
    assert lo < rate, f"Lower bound {lo} should be less than rate {rate}"
    assert hi > rate, f"Upper bound {hi} should be greater than rate {rate}"


def test_rate_computation():
    """Verify rate per 1,000 device-days computation."""
    events = 5
    device_days = 500
    rate_per_1000 = (events / device_days) * 1000
    assert rate_per_1000 == 10.0, f"Expected 10.0, got {rate_per_1000}"


def test_cohort_assignment():
    """Verify cohort assignment logic based on pct_closed_loop."""
    pct_closed_loop = 50
    cohort = "closed_loop" if pct_closed_loop >= 50 else "open_loop"
    assert cohort == "closed_loop", f"Expected closed_loop, got {cohort}"

Here is the complete docs/validation_report.md:

# Validation Report

## Test Results

All 7 tests pass across two test modules:

### test_etl.py (4 tests)
- `test_user_summary_row_count`: Confirms 500 rows in user_summary (one per synthetic user).
- `test_no_null_user_ids`: Confirms zero null user_id values in the final table.
- `test_glucose_range`: Confirms all glucose readings fall within 20-600 mg/dL.
- `test_date_range`: Confirms no null device_start_date values.

### test_analysis.py (3 tests)
- `test_poisson_ci_known_value`: Verifies Poisson CI math against known input (10 events, 1000 exposure).
- `test_rate_computation`: Verifies rate-per-1000 calculation (5 events / 500 days * 1000 = 10.0).
- `test_cohort_assignment`: Verifies pct_closed_loop >= 50 maps to closed_loop cohort.

## Data Quality Checks

| Metric | Expected | Observed |
|--------|----------|----------|
| Total users | 500 | 500 |
| Null user_ids | 0 | 0 |
| Glucose range | 20-600 mg/dL | Confirmed |
| HbA1c range | 5.0-14.0% | Confirmed |
| Severe hypo rate | ~2-5% | ~3.5% |
| DKA rate | <1% | ~0.8% |

## Known Limitations

1. **Synthetic data**: All data is generated. Patterns are illustrative, not derived from real patients.
2. **Simulated pre-device HbA1c**: Pre-device baseline values are generated independently, not from longitudinal records.
3. **Logistic regression assumptions**: The onboarding-retention model assumes a linear log-odds relationship without accounting for all potential confounders.
4. **Illustrative findings only**: Statistical results demonstrate methodology, not real clinical evidence.

Your pipeline is now backed by automated tests and documented quality checks. Next up, you will build an interactive Tableau dashboard to visualize your findings for stakeholders.

Build the Tableau Dashboard and Finalize Deliverables

Your analysis notebooks are complete, your tests pass, and your pipeline runs end-to-end. Now it is time to translate all of that work into assets that a hiring manager can actually interact with.

In this step, you will build an interactive Tableau dashboard, publish it to the web, and create the stakeholder deliverables (one-pager, demo script, finalized README) that turn your repository into a polished portfolio piece.

In this step, get ready to:

• Build three Tableau dashboard sheets (Safety, Outcomes, Commercial KPIs) and publish them to Tableau Public.
• Create a prioritization backlog and export dashboard screenshots.
• Write stakeholder deliverables and finalize the project README.

Connect to your data and build the dashboard sheets

Each of your three audiences (safety, clinical, commercial) needs a tailored view of the data. You will create one worksheet per audience inside a single Tableau workbook, then combine them into a dashboard.

• Open Tableau Desktop Public Edition.
• On the Start page, under Connect, click Text File.
• Navigate to your data/processed/ folder and select user_summary_export.csv.
• Click Open.
• On the Data Source page, verify that Tableau correctly parsed your columns (user_id, age, region, mean_glucose, pct_time_in_range, hba1c_pre, hba1c_post, etc.).
• Click the Sheet 1 tab at the bottom to begin building your first view.

Why connect to the CSV export instead of SQLite directly?

Tableau Desktop Public Edition supports text files (CSV) as a data source. Using the pre-built user_summary_export.csv produced by your Prefect pipeline gives you a single flat file with all the joined metrics already computed.

When you publish to Tableau Public, the data is extracted into Tableau's format automatically. No live connection is maintained.

Now build the three dashboard sheets. Each sheet focuses on one audience.

For the Safety dashboard sheet:

• Create a bar chart showing safety event rates (severe hypoglycemia, DKA) per 1,000 device-days broken down by cohort (age group, device mode, region).
• Add a line chart for the weekly safety event control chart. Include reference lines for the 2-sigma and 3-sigma thresholds.
• Create a dot plot or bar chart as a forest-plot mockup showing relative risk ratios with confidence interval error bars for each cohort comparison.

For the Outcomes dashboard sheet:

• Build a waterfall chart showing the pre-to-post HbA1c change distribution across users.
• Create a grouped bar chart of time-in-range percentage by cohort (age group, diabetes duration).
• Add a line chart showing the improvement trajectory (mean glucose or time-in-range) over the first 180 days post-device start.

For the Commercial KPIs sheet:

• Build a bar chart of 90-day retention rate by region.
• Create a funnel visualization for the onboarding pipeline (enrolled, first closed-loop session, 90-day active).
• Add a heatmap of support ticket counts by region and severity level.

Use your notebooks as a blueprint

You already built each of these visualizations in matplotlib/seaborn inside your analysis notebooks. Use the same axis labels, cohort definitions, and color coding so your Tableau dashboard stays consistent with the notebook outputs.

• Rename each sheet tab at the bottom (double-click the tab name) to Safety, Outcomes, and Commercial KPIs.

Fields showing as strings instead of numbers?

Tableau sometimes misclassifies numeric columns imported from CSV. Right-click the field in the Data pane and select Change Data Type to set it to Number (decimal) or Number (whole). Repeat for any date fields that should be Date type.

Still stuck on building a chart?

Publish to Tableau Public and create the prioritization backlog

With all three sheets complete, you will now publish your workbook so a hiring manager can interact with it from a browser link.

• In Tableau Desktop Public Edition, select Server > Tableau Public > Save to Tableau Public.
• Sign in with your Tableau Public account. If you do not have one, click the link to create a free account.
• Enter a workbook name (e.g. Bionic Pancreas Analytics) and click Save.
• After publishing completes, your browser opens to the live dashboard on Tableau Public. Copy the URL from your browser address bar. You will add this to your README shortly.

What happens when you publish?

Tableau Public creates a data extract from your CSV. The dashboard becomes publicly accessible, and anyone with the link can view and interact with it. This is exactly what you want for a portfolio piece.

Tableau Public is limited to 15 million rows per workbook. Your user summary CSV is well within that limit.

Next, capture screenshots and create the backlog.

• Take a screenshot of each dashboard sheet (Safety, Outcomes, Commercial KPIs) and save them to dashboards/ as safety_dashboard.png, outcomes_dashboard.png, and commercial_kpis_dashboard.png.
• Create the prioritization backlog file by running:

touch docs/backlog.csv

• Open docs/backlog.csv and add the following JIRA-style intake content:

id,title,priority,owner,eta,status
BP-001,Expand safety signal detection to sub-populations,High,Data Analyst,2025-Q3,In Progress
BP-002,Automate weekly HbA1c improvement report,Medium,Analytics Engineer,2025-Q3,Backlog
BP-003,Add NPS survey data to CRM integration,Medium,Data Analyst,2025-Q4,Backlog
BP-004,Build regional early-warning composite score,High,Data Scientist,2025-Q3,In Progress
BP-005,Validate control chart thresholds with clinical team,Low,Data Analyst,2025-Q4,Backlog

What is this backlog for?

In real analytics teams, work arrives as intake requests with priority, owner, and estimated delivery. Including a sample backlog in your portfolio demonstrates that you understand how analytical work is triaged and managed in a team setting.

• Save docs/backlog.csv.
• Confirm your dashboards/ folder now contains the three PNG screenshots by running:

ls dashboards/

You should see safety_dashboard.png, outcomes_dashboard.png, and commercial_kpis_dashboard.png listed.

Publish option greyed out?

The Save to Tableau Public option only appears when you have at least one field on a worksheet. Make sure you are on a sheet tab (not the Data Source tab) and that at least one visualization has been built.

Help me troubleshoot publishing to Tableau Public

Write stakeholder deliverables and finalize the README

A strong portfolio project is not just code. It includes the communication artifacts that prove you can translate analysis into action. You will now create a one-pager, a demo script, and a finalized README.

• Create a one-pager slide deck (using Google Slides, PowerPoint, or any tool you prefer) and export it as docs/stakeholder_onepager.pdf. Include the following sections:

• 2-3 key findings from your analysis (e.g. "Closed-loop users show 0.8% mean HbA1c reduction vs. open-loop", "Region West has 2x the support ticket rate").
• Recommended actions tied to each finding.
• Confidence levels (high/medium/low) for each recommendation based on sample size and statistical significance.
• Limitations section noting synthetic data, sample size constraints, and any assumptions.

• Create the demo script by running:

touch docs/demo_script.md

• Open docs/demo_script.md and write a 10-minute interview walkthrough. Structure it with the following time blocks:

# Demo Script: 10-Minute Interview Walkthrough

## Project Overview (1 min)
- Introduce the problem: analytics for an automated insulin delivery company
- Mention the three data domains (device telemetry, clinical outcomes, CRM)
- State the goal: fuse data to surface safety, efficacy, and business insights

## Data Generation and ETL (2 min)
- Show generate_data.py: 500 users, ~180 days telemetry, reproducible with seed
- Walk through etl.py: CSV ingestion, cleaning, SQLite warehouse
- Highlight the user_summary table joining all domains

## Safety Signal Detection (2 min)
- Show Poisson CI calculations for event rates
- Highlight the control chart with 2-sigma/3-sigma thresholds
- Discuss relative risk ratios between closed-loop and open-loop cohorts

## Clinical Outcomes (2 min)
- Explain pre/post cohort definition (90-day windows)
- Show HbA1c paired change with Wilcoxon test and effect size
- Mention sensitivity analyses and stratification

## CRM/Commercial Analysis (1 min)
- Logistic regression: onboarding speed predicts retention
- Support ticket burden correlates with worse glycemic outcomes
- Regional comparison of fast vs. slow onboarding

## Tableau Dashboard (1 min)
- Open live Tableau Public link
- Walk through Safety, Outcomes, and Commercial KPIs sheets
- Point out interactive filters

## Pipeline and Testing (1 min)
- Show Prefect flow with @task and @flow decorators
- Mention pytest suite (7 tests covering ETL and analysis)
- Highlight reproducibility: clone, install, run one command

What makes a good demo script?

A 10-minute walkthrough should tell a story, not just list features. Start with the problem, show how you solved it, and end with impact. Practice delivering it aloud at least twice before an interview.

• Save docs/demo_script.md.
• Now finalize README.md by replacing the current contents with the following:

# Integrated Clinical & Commercial Analytics for a Bionic Pancreas

![Python](https://img.shields.io/badge/Python-3.11%2B-blue)
![SQLite](https://img.shields.io/badge/Database-SQLite-lightgrey)
![Prefect](https://img.shields.io/badge/Orchestration-Prefect-blue)
![Tableau](https://img.shields.io/badge/Dashboard-Tableau%20Public-orange)
![Tests](https://img.shields.io/badge/Tests-7%20passing-brightgreen)

## Project Summary
End-to-end analytics pipeline simulating the work of a data analyst at an automated insulin delivery (AID) device company. Fuses synthetic device telemetry, clinical outcomes, and CRM data to surface safety signals, measure clinical efficacy, and drive commercial decisions.

## Quick Start
```bash
# Clone and set up
git clone <your-repo-url>
cd bionic-pancreas-analytics
python -m venv venv && source venv/bin/activate
pip install pandas==2.2.3 numpy==2.1.3 scipy==1.14.1 statsmodels==0.14.6 matplotlib==3.10.9 seaborn==0.13.2 jupyter pytest prefect

# Run the full pipeline
python src/pipeline.py
```

## Live Dashboard
[View the interactive Tableau Public dashboard](YOUR_TABLEAU_PUBLIC_URL)

## Methods Overview
- **Synthetic data generation**: 500 users, ~180 days of 5-minute telemetry, clinical visits, CRM records
- **ETL pipeline**: Python + SQLite with cleaning, validation, and joined user summary table
- **Safety signal detection**: Event rates with Poisson CIs, c-chart control charts, relative risk ratios
- **Clinical outcomes**: Pre/post HbA1c (Wilcoxon signed-rank, Cohen's d), time-in-range analysis
- **CRM/commercial**: Onboarding speed vs. retention (logistic regression), support ticket analysis
- **Orchestration**: Prefect @flow/@task with weekly cron schedule
- **Testing**: 7 pytest tests covering ETL integrity and statistical functions

## Privacy Considerations
All data is synthetically generated using `src/generate_data.py` with a fixed random seed. No real patient data is used, referenced, or represented. This project is HIPAA-aware by design.

## Repository Navigation
| Directory | Contents |
|-----------|----------|
| `data/raw/` | Synthetic CSV exports (4 datasets) |
| `data/processed/` | SQLite analytical warehouse + CSV export |
| `src/` | Data generation, ETL, and Prefect pipeline |
| `notebooks/` | Analytical Jupyter notebooks (safety, outcomes, CRM) |
| `tests/` | pytest unit tests for ETL and analysis |
| `dashboards/` | Tableau dashboard screenshots |
| `docs/` | Data dictionary, validation report, demo script, backlog |

What does the finalized README include?

The README now includes shield badges for visual scanning, a one-command quick-start, a direct link to your live Tableau Public dashboard, a methods summary, privacy statement, and a navigation table. These elements let a hiring manager assess the project in under 60 seconds.

• Replace YOUR_TABLEAU_PUBLIC_URL in the README with the actual URL you copied from Tableau Public earlier.
• Save README.md.
• Confirm all deliverables exist by running:

ls docs/

You should see backlog.csv, data_dictionary.md, demo_script.md, stakeholder_onepager.pdf, and validation_report.md listed.

Missing the stakeholder one-pager?

The one-pager is created outside your code editor (Google Slides, PowerPoint, etc.) and exported as a PDF. Make sure you export it to docs/stakeholder_onepager.pdf before checking your file list.

Help me create a one-pager slide for my analytics project

Your project is now complete. You have a published interactive dashboard, a stakeholder one-pager, a rehearsal-ready demo script, and a README that gets to the point in under a minute. This is exactly the kind of deliverable set that demonstrates analytics ownership in an interview.

Secret mission

Add a Regional Early-Warning System

What if a specific region is silently failing? Users might not complain directly, but rising device alerts, increased support ticket volume, and worsening glucose readings could signal a problem before it reaches a crisis. In this extension, you'll build a composite early-warning system that flags at-risk regions automatically.

Clean Up Your Resources

Clean Up Your Resources

This project runs entirely locally with free tools, so there are no ongoing costs. Decide whether to keep your resources running, pause them to come back later, or delete them entirely.

Resources you used:

• Local project folder (all source code, SQLite database, raw CSVs, notebooks, test files, and dashboard screenshots).
• Python virtual environment (inside your project folder).
• Tableau Public published dashboard (hosted online at public.tableau.com).
• Tableau Desktop Public Edition (installed application).

Keep everything running

No action needed. Choose this if you are still actively building on this project or plan to use it as a live portfolio piece.

• Your local project folder remains on your Desktop.
• Your Tableau Public dashboard stays live at its published URL for hiring managers to view.
• Tableau Desktop Public Edition stays installed for future visualization work.

Pause - I'll come back to this later

Shut down any running processes to free up memory. Your files and published dashboard stay intact.

• If you have a Prefect server running, stop it by pressing Ctrl+C in the terminal where you ran the server.
• If you have Jupyter Notebook running, stop the kernel and close the browser tab.
• Deactivate your Python virtual environment by running:

deactivate

• Your project files and Tableau Public dashboard remain available when you return.

Delete - I don't want to use this again

Remove all project resources and start fresh. Follow each step below for a complete cleanup.

Delete your local project folder

• Open your terminal and navigate to the location where your project folder lives (your Desktop, if you followed the setup steps).
• Delete the entire project folder by running:

macOS

cd ~/Desktop
rm -rf bionic-pancreas-analytics

Windows

cd ~\Desktop
Remove-Item -Recurse -Force bionic-pancreas-analytics

Replace bionic-pancreas-analytics with the actual name of your project folder if you named it differently.

Delete your published Tableau Public dashboard

• Sign in to Tableau Public at public.tableau.com.
• Go to your profile page.
• Select the More Actions icon on the viz you want to delete.
• Select Delete This Viz.
• Repeat for each dashboard you published during this project (Safety, Outcomes, Commercial KPIs, Regional Alert).

Uninstall Tableau Desktop Public Edition

macOS

• Open Finder and go to your Applications folder.
• Delete the Tableau Public application by dragging it to the Trash.
• Empty your Trash.

Windows

• Open Control Panel and go to Programs, then Programs and Features.
• Right-click Tableau Public and select Uninstall/Change.
• From the Modify Setup menu, select Uninstall.

Nice Work!

Nice Work!

You just built a complete, portfolio-ready analytics pipeline for a simulated medical device company. This project demonstrates end-to-end ownership of the analytics lifecycle, from synthetic data generation and governance through statistical analysis to executive-facing delivery.

Here's what you accomplished:

• Built an end-to-end ETL pipeline that generates synthetic device telemetry, clinical outcomes, and CRM data for 500 users, then cleans and loads it into a SQLite analytical warehouse with analysis-ready summary tables.
• Performed safety signal detection with Poisson confidence intervals and control charts, clinical outcomes analysis with pre/post HbA1c comparisons and time-in-range metrics, and commercial cohort analysis linking onboarding speed to retention.
• Orchestrated the full pipeline with Prefect for reproducible scheduled execution, built interactive Tableau dashboards for three stakeholder audiences, and delivered unit tests, a validation report, and a stakeholder one-pager.
• Secret Mission: Created a regional early-warning composite risk score combining device alert frequency, support ticket velocity, and glycemic deterioration into a 0-100 index surfaced as a Tableau alert view.

Ready to quiz yourself?