VITALS-WATCH

Executable Code

#!/usr/bin/env python3
"""
Claw4S Skill: Apple Watch Vital Sign Flare Detection
Bayesian Change-Point Detection on HR/HRV/SpO2 Time Series

Detects autoimmune flare onset from wearable vital sign data using
Bayesian online change-point detection (Adams & MacKay 2007).

Author: Zamora-Tehozol EA (ORCID:0000-0002-7888-3961), DNAI
License: MIT
"""

import numpy as np
from scipy import stats
import warnings
warnings.filterwarnings('ignore')

# ══════════════════════════════════════════════════════════════════
# BAYESIAN ONLINE CHANGE-POINT DETECTION (BOCPD)
# Adams & MacKay, arXiv:0710.3742
# ══════════════════════════════════════════════════════════════════

class BOCPD:
    """Bayesian Online Change-Point Detection with Gaussian likelihood."""

    def __init__(self, hazard_rate=1/200, mu0=0.0, kappa0=1.0, alpha0=1.0, beta0=1.0):
        self.hazard = hazard_rate
        self.mu0 = mu0
        self.kappa0 = kappa0
        self.alpha0 = alpha0
        self.beta0 = beta0

    def detect(self, data):
        """Run BOCPD on 1D time series. Returns run-length posterior matrix."""
        T = len(data)
        R = np.zeros((T + 1, T + 1))
        R[0, 0] = 1.0

        mu = np.array([self.mu0])
        kappa = np.array([self.kappa0])
        alpha = np.array([self.alpha0])
        beta = np.array([self.beta0])

        change_points = []
        max_run_lengths = np.zeros(T)

        for t in range(T):
            x = data[t]
            # Predictive probability (Student-t)
            df = 2 * alpha
            scale = np.sqrt(beta * (kappa + 1) / (alpha * kappa))
            pred_probs = np.exp(stats.t.logpdf(x, df=df, loc=mu, scale=scale))

            # Growth probabilities
            R[1:t+2, t+1] = R[:t+1, t] * pred_probs * (1 - self.hazard)
            # Change-point probability
            R[0, t+1] = np.sum(R[:t+1, t] * pred_probs * self.hazard)
            # Normalize
            evidence = R[:t+2, t+1].sum()
            if evidence > 0:
                R[:t+2, t+1] /= evidence

            # Track max run length
            max_run_lengths[t] = np.argmax(R[:t+2, t+1])

            # Detect change-point: run length drops to near 0
            if t > 5 and max_run_lengths[t] < 3 and max_run_lengths[t-1] > 10:
                change_points.append(t)

            # Update sufficient statistics
            new_mu = np.append(self.mu0, (kappa * mu + x) / (kappa + 1))
            new_kappa = np.append(self.kappa0, kappa + 1)
            new_alpha = np.append(self.alpha0, alpha + 0.5)
            new_beta = np.append(self.beta0, beta + kappa * (x - mu)**2 / (2 * (kappa + 1)))
            mu, kappa, alpha, beta = new_mu, new_kappa, new_alpha, new_beta

        return change_points, max_run_lengths, R


# ══════════════════════════════════════════════════════════════════
# VITAL SIGN PROCESSOR
# ══════════════════════════════════════════════════════════════════

class VitalSignFlareDetector:
    """Multi-channel vital sign analysis for autoimmune flare detection."""

    # Clinical thresholds based on literature
    THRESHOLDS = {
        'hr_rest_elevated': 90,      # Resting HR > 90 bpm
        'hrv_sdnn_low': 50,          # SDNN < 50 ms (reduced vagal tone)
        'spo2_low': 94,              # SpO2 < 94%
        'hr_variability_increase': 1.5,  # >1.5x baseline HR CV
    }

    def __init__(self, baseline_days=7):
        self.baseline_days = baseline_days
        self.bocpd_hr = BOCPD(hazard_rate=1/100, mu0=72, kappa0=1, alpha0=2, beta0=50)
        self.bocpd_hrv = BOCPD(hazard_rate=1/100, mu0=60, kappa0=1, alpha0=2, beta0=100)
        self.bocpd_spo2 = BOCPD(hazard_rate=1/100, mu0=97, kappa0=1, alpha0=2, beta0=2)

    def analyze(self, hr_series, hrv_series, spo2_series, timestamps_hours):
        """Analyze multi-channel vitals for flare signals."""
        n = len(hr_series)
        baseline_n = min(self.baseline_days * 24, n // 3)

        # Baseline statistics
        hr_baseline = {'mean': np.mean(hr_series[:baseline_n]),
                       'std': np.std(hr_series[:baseline_n])}
        hrv_baseline = {'mean': np.mean(hrv_series[:baseline_n]),
                        'std': np.std(hrv_series[:baseline_n])}
        spo2_baseline = {'mean': np.mean(spo2_series[:baseline_n]),
                         'std': np.std(spo2_series[:baseline_n])}

        # Run BOCPD on each channel
        hr_cp, hr_rl, _ = self.bocpd_hr.detect(hr_series)
        hrv_cp, hrv_rl, _ = self.bocpd_hrv.detect(hrv_series)
        spo2_cp, spo2_rl, _ = self.bocpd_spo2.detect(spo2_series)

        # Compute flare probability at each time point
        flare_scores = np.zeros(n)
        for t in range(baseline_n, n):
            score = 0.0
            window = slice(max(0, t-24), t+1)

            # HR elevation from baseline
            hr_delta = np.mean(hr_series[window]) - hr_baseline['mean']
            if hr_delta > 2 * hr_baseline['std']:
                score += 0.25
            if np.mean(hr_series[window]) > self.THRESHOLDS['hr_rest_elevated']:
                score += 0.15

            # HRV depression (autonomic dysfunction)
            hrv_delta = hrv_baseline['mean'] - np.mean(hrv_series[window])
            if hrv_delta > 2 * hrv_baseline['std']:
                score += 0.25
            if np.mean(hrv_series[window]) < self.THRESHOLDS['hrv_sdnn_low']:
                score += 0.15

            # SpO2 drops (pulmonary involvement)
            if np.mean(spo2_series[window]) < self.THRESHOLDS['spo2_low']:
                score += 0.20

            flare_scores[t] = min(score, 1.0)

        # Combine change-points across channels
        all_cp = set()
        for cp_list in [hr_cp, hrv_cp, spo2_cp]:
            all_cp.update(cp_list)

        # Multi-channel concordance: if 2+ channels have nearby CPs
        concordant_cp = []
        for cp in sorted(all_cp):
            channels_hit = 0
            if any(abs(cp - c) < 12 for c in hr_cp): channels_hit += 1
            if any(abs(cp - c) < 12 for c in hrv_cp): channels_hit += 1
            if any(abs(cp - c) < 12 for c in spo2_cp): channels_hit += 1
            if channels_hit >= 2:
                concordant_cp.append((cp, channels_hit))

        # Risk classification
        max_flare_score = np.max(flare_scores) if len(flare_scores) > 0 else 0
        if max_flare_score >= 0.6:
            risk_level = "HIGH"
            recommendation = ("Probable flare onset detected. Contact rheumatologist within 24h. "
                            "Consider increasing monitoring frequency to q1h.")
        elif max_flare_score >= 0.35:
            risk_level = "MODERATE"
            recommendation = ("Possible flare signal. Increase monitoring frequency. "
                            "Consider symptom diary and next-day clinical contact.")
        else:
            risk_level = "LOW"
            recommendation = "Vital signs within expected range. Continue routine monitoring."

        return {
            'baselines': {'hr': hr_baseline, 'hrv': hrv_baseline, 'spo2': spo2_baseline},
            'change_points': {'hr': hr_cp, 'hrv': hrv_cp, 'spo2': spo2_cp},
            'concordant_change_points': concordant_cp,
            'flare_scores': flare_scores,
            'max_flare_score': round(float(max_flare_score), 3),
            'risk_level': risk_level,
            'recommendation': recommendation,
            'n_observations': n,
            'baseline_hours': baseline_n,
        }


def generate_synthetic_vitals(n_hours=336, flare_onset_hour=200, seed=42):
    """Generate synthetic Apple Watch vital signs with embedded flare event."""
    rng = np.random.RandomState(seed)
    t = np.arange(n_hours)

    # Circadian rhythm component
    circadian_hr = 5 * np.sin(2 * np.pi * t / 24 - np.pi/2)
    circadian_hrv = 8 * np.sin(2 * np.pi * t / 24 + np.pi/4)

    # Baseline vitals
    hr = 72 + circadian_hr + rng.normal(0, 3, n_hours)
    hrv_sdnn = 65 + circadian_hrv + rng.normal(0, 5, n_hours)
    spo2 = 97.5 + rng.normal(0, 0.5, n_hours)

    # Inject flare: gradual HR increase, HRV decrease, SpO2 dip
    flare_mask = t >= flare_onset_hour
    flare_ramp = np.clip((t - flare_onset_hour) / 48, 0, 1)  # 48h ramp

    hr[flare_mask] += 15 * flare_ramp[flare_mask] + rng.normal(0, 2, flare_mask.sum())
    hrv_sdnn[flare_mask] -= 25 * flare_ramp[flare_mask] + rng.normal(0, 3, flare_mask.sum())
    spo2[flare_mask] -= 2.5 * flare_ramp[flare_mask] + rng.normal(0, 0.3, flare_mask.sum())

    # Clamp physiological ranges
    hr = np.clip(hr, 40, 180)
    hrv_sdnn = np.clip(hrv_sdnn, 10, 150)
    spo2 = np.clip(spo2, 85, 100)

    return hr, hrv_sdnn, spo2, t


# ══════════════════════════════════════════════════════════════════
# DEMO
# ══════════════════════════════════════════════════════════════════

if __name__ == "__main__":
    print("=" * 70)
    print("VITALS-WATCH: Apple Watch Flare Detection Skill")
    print("Bayesian Online Change-Point Detection (Adams & MacKay 2007)")
    print("Authors: Zamora-Tehozol EA (ORCID:0000-0002-7888-3961), DNAI")
    print("=" * 70)

    # Generate 14-day synthetic data with flare at day ~8
    hr, hrv, spo2, t = generate_synthetic_vitals(n_hours=336, flare_onset_hour=200)

    print(f"\n[DATA] Synthetic 14-day vitals: {len(t)} hourly observations")
    print(f"  HR range: {hr.min():.0f}-{hr.max():.0f} bpm")
    print(f"  HRV range: {hrv.min():.0f}-{hrv.max():.0f} ms (SDNN)")
    print(f"  SpO2 range: {spo2.min():.1f}-{spo2.max():.1f}%")
    print(f"  Flare injected at hour 200 (day 8.3)")

    detector = VitalSignFlareDetector(baseline_days=7)
    results = detector.analyze(hr, hrv, spo2, t)

    print(f"\n── BASELINE (first {results['baseline_hours']} hours) ──")
    for name, bl in results['baselines'].items():
        print(f"  {name.upper()}: mean={bl['mean']:.1f}, std={bl['std']:.1f}")

    print(f"\n── CHANGE-POINTS DETECTED ──")
    print(f"  HR:   {len(results['change_points']['hr'])} change-points")
    print(f"  HRV:  {len(results['change_points']['hrv'])} change-points")
    print(f"  SpO2: {len(results['change_points']['spo2'])} change-points")

    print(f"\n── CONCORDANT CHANGE-POINTS (2+ channels) ──")
    for cp, n_ch in results['concordant_change_points']:
        print(f"  Hour {cp} (day {cp/24:.1f}): {n_ch} channels concordant")

    print(f"\n── FLARE ASSESSMENT ──")
    print(f"  Max flare score: {results['max_flare_score']}")
    print(f"  Risk level: {results['risk_level']}")
    print(f"  Recommendation: {results['recommendation']}")

    # Scenario 2: No flare (stable vitals)
    print(f"\n{'='*70}")
    print("SCENARIO 2: Stable patient (no flare)")
    hr2, hrv2, spo2_2, t2 = generate_synthetic_vitals(n_hours=336, flare_onset_hour=9999)
    results2 = detector.analyze(hr2, hrv2, spo2_2, t2)
    print(f"  Max flare score: {results2['max_flare_score']}")
    print(f"  Risk level: {results2['risk_level']}")

    print(f"\n── LIMITATIONS ──")
    print("  • Synthetic data only — not validated on real Apple Watch exports")
    print("  • BOCPD assumes Gaussian likelihoods (real vital signs may be non-Gaussian)")
    print("  • Circadian rhythm modeling is simplified (single harmonic)")
    print("  • Does not account for medication effects on HR/HRV")
    print("  • Not FDA-cleared; decision-support only, not diagnostic")
    print("  • Requires ≥7 days baseline for reliable detection")
    print(f"\n{'='*70}")
    print("END — Vitals-Watch Skill v1.0")

Demo Output

Max flare score: 0.65, Risk level: HIGH
Scenario 2 (stable): Max flare score: 0.0, Risk level: LOW