# Common Clinical Data Cleaning Patterns

Reference document for the clean-data skill. Covers recurring data quality issues
in electronic health records, registries, and research databases.

---

## 1. Missing Data Patterns

### Classification

- **MCAR (Missing Completely At Random)**: Missingness is unrelated to any variable.
  Example: random equipment failure during lab measurement.
  Test: Little's MCAR test (chi-square). If p > 0.05, MCAR is plausible.

- **MAR (Missing At Random)**: Missingness depends on observed variables but not the
  missing value itself. Example: younger patients less likely to have bone density measured.
  Cannot be directly tested; inferred from associations between missingness indicators
  and observed covariates.

- **MNAR (Missing Not At Random)**: Missingness depends on the unobserved value itself.
  Example: severely ill patients too sick to complete follow-up surveys.
  Cannot be tested from the data alone; requires domain knowledge.

### Heuristic Assessment

1. Compute missing percentage per variable.
2. Create missingness indicator (0/1) for each variable with >5% missing.
3. Correlate missingness indicators with observed variables (chi-square, t-test).
4. If strong correlations exist: likely MAR. If none: plausible MCAR. If clinical
   reasoning suggests the value itself drives missingness: suspect MNAR.

### When to Use Each Imputation Method

| Method | When appropriate | Caution |
|--------|-----------------|---------|
| Listwise deletion (complete case) | MCAR, low % missing (<5%), large sample | Biased if MAR/MNAR; reduces power |
| Mean/median imputation | Quick exploratory analysis only | Underestimates variance; distorts distributions |
| Last observation carried forward | Longitudinal data, slow-changing variables | Biased if trajectory is changing |
| Multiple imputation (MICE) | MAR, moderate missing (5-40%), multivariate | Requires careful model specification |
| Maximum likelihood (FIML) | MAR, SEM or regression contexts | Needs software support |
| Sensitivity analysis | Always for MNAR suspicion | Report results under multiple assumptions |

### Key Rule

Never impute the outcome variable in the primary analysis without explicit justification.
Report the missing data mechanism assumption in the methods section.

---

## 2. Outlier Detection

### Statistical Methods

**IQR Method (Tukey Fences)**:
- Lower fence: Q1 - 1.5 * IQR
- Upper fence: Q3 + 1.5 * IQR
- Robust to non-normal distributions
- Preferred for clinical data where normality is rarely guaranteed

**Z-Score Method**:
- Flag values with |z| > 3 (or |z| > 2.5 for smaller samples)
- Assumes approximate normality
- Sensitive to the outliers themselves (mean and SD are affected)

**Modified Z-Score (MAD-based)**:
- Uses median and Median Absolute Deviation instead of mean/SD
- More robust than standard Z-score
- Formula: M_i = 0.6745 * (x_i - median) / MAD

### Decision Framework

| Scenario | Recommended action |
|----------|--------------------|
| Data entry error (clearly impossible) | Correct if source available; else set to missing |
| Measurement error (instrument fault) | Set to missing; document in cleaning log |
| True extreme value (biologically plausible) | Keep in dataset; consider sensitivity analysis with/without |
| Ambiguous | Flag for domain expert review; do not remove without justification |

### Clinical Context Matters

A BMI of 50 is an outlier statistically but clinically real. An age of 200 is impossible.
A creatinine of 15 mg/dL is extreme but occurs in dialysis patients. Always consult the
codebook and clinical context before removing outliers.

For the **biologically-impossible** side of this distinction — hard compatible-with-life bounds
per organ system, plus cross-field logical-consistency rules (date ordering, derived-vs-source,
sex-specific) — see `implausible_value_rules.md`. That handles likely data-entry/unit errors
(correct-or-set-missing); this section handles statistically extreme but possible values (keep +
sensitivity).

---

## 3. Duplicate Detection

### Exact Duplicates

- Identical across ALL columns.
- Usually safe to remove (keep first occurrence).
- Common cause: accidental double-submission or ETL errors.

### Near-Duplicates

- Same patient identifier, different records.
- May be legitimate (multiple visits) or errors (same visit entered twice with typos).

### Detection Strategy

1. Check for exact row duplicates: `df.duplicated().sum()`
2. Check for duplicate patient IDs: `df['patient_id'].duplicated().sum()`
3. For near-duplicates: group by patient ID, sort by date, check for records within
   a suspiciously short time window (e.g., same day for what should be annual visits).
4. Fuzzy matching: consider Levenshtein distance on name fields if no unique ID exists.

### Resolution

- Exact duplicates: drop duplicates, log count.
- Same-patient near-duplicates: present to researcher for manual review.
- Never auto-merge patient records without explicit approval.

---

## 4. Date Handling

### Common Date Formats in Clinical Data

| Format | Example | Source |
|--------|---------|--------|
| YYYY-MM-DD | 2024-03-15 | ISO 8601, most databases |
| MM/DD/YYYY | 03/15/2024 | US clinical systems |
| DD/MM/YYYY | 15/03/2024 | European systems |
| YYYYMMDD | 20240315 | DICOM, HL7 |
| DD-Mon-YYYY | 15-Mar-2024 | Some EMR exports |
| Excel serial | 45366 | Excel numeric date |

### Common Issues

- **Ambiguous dates**: Is 03/04/2024 March 4th or April 3rd? Check the data source locale.
  Look for values >12 in the first or second position to disambiguate.
- **Impossible dates**: February 30, month 13, year 0001.
- **Future dates**: Dates after the data extraction date (except for scheduled appointments).
- **Timezone issues**: Rarely relevant for clinical research dates, but critical for timestamps
  in multi-site studies across time zones.
- **Two-digit years**: 24 could be 1924 or 2024. Use a pivot year (e.g., 30: <=30 means 2000s,
  >30 means 1900s) or infer from context.

### Standardization

1. Parse all date columns to datetime using `pd.to_datetime(col, format=..., errors='coerce')`.
2. Check for NaT (failed parses) and investigate.
3. Standardize to ISO 8601 (YYYY-MM-DD) for storage.
4. Calculate derived variables (age at event, follow-up duration) from standardized dates.

---

## 5. Category Harmonization

### Common Inconsistencies

| Raw values | Harmonized |
|-----------|-----------|
| "Male", "male", "M", "MALE", " Male " | "Male" |
| "Y", "Yes", "yes", "YES", "1", "True" | 1 or "Yes" |
| "Right", "Rt", "R", "right", "RT" | "Right" |
| "Non-small cell", "NSCLC", "non small cell" | "NSCLC" |

### Harmonization Steps

1. Strip whitespace: `series.str.strip()`
2. Normalize case: `series.str.lower()` or `series.str.title()`
3. Build a mapping dictionary for known synonyms.
4. Review unmapped values manually.
5. Apply mapping: `series.map(mapping_dict).fillna(series)`

### Encoding Standards

- **ICD-10**: Diagnosis codes. Watch for version differences (ICD-10-CM vs ICD-10-PCS).
- **SNOMED CT**: Clinical terminology. More granular than ICD-10.
- **LOINC**: Laboratory observations. Use for standardizing lab test names.
- **CPT/HCPCS**: Procedure codes.

When possible, map free-text categories to standard coding systems. Document the mapping
table and include it in supplementary materials.

---

## 6. Common Clinical Data Pitfalls

### Lab Values with Inequality Prefixes

Values like "<0.01", ">10000", "<=5" are common for lab results at detection limits.

**Handling options**:
- Replace with the limit value: "<0.01" -> 0.01 (conservative)
- Replace with half the limit: "<0.01" -> 0.005 (common in environmental studies)
- Replace with limit / sqrt(2): "<0.01" -> 0.00707 (EPA method)
- Keep as censored data and use appropriate statistical methods (Tobit regression)

Document the chosen method in the statistical analysis plan.

### Mixed Units

Common in multi-site studies or data merged from different systems.

| Analyte | Unit A | Unit B | Conversion |
|---------|--------|--------|-----------|
| Glucose | mg/dL | mmol/L | mg/dL = mmol/L * 18.018 |
| Creatinine | mg/dL | umol/L | mg/dL = umol/L / 88.4 |
| Hemoglobin | g/dL | g/L | g/dL = g/L / 10 |
| Calcium | mg/dL | mmol/L | mg/dL = mmol/L * 4.008 |

**Detection**: look for bimodal distributions in lab values -- one mode per unit system.

### Sentinel Values

Values used as placeholders for missing data in legacy systems:

| Sentinel | Meaning |
|----------|---------|
| 999, 9999, 99999 | Missing / not recorded |
| -1, -9, -99 | Missing / not applicable |
| 0 | Could be true zero OR missing -- context-dependent |
| 88, 77 | "Not applicable" or "Refused" in survey data |
| 8888 | "Not applicable" or "Missing" in health screening/institutional databases |
| 01/01/1900 | Default/missing date |

**Action**: Replace sentinel values with `NaN` BEFORE computing any statistics.
Document which values were treated as sentinel.

### Excel Date Corruption

Excel auto-converts certain strings to dates:
- Gene names: SEPT1 -> Sep-1, MARCH1 -> Mar-1, DEC1 -> Dec-1
- Sample IDs: 1-3 -> Jan-3, 2/4 -> Feb-4

**Prevention**: Open CSV in a text editor first to verify. Import with explicit dtypes
in pandas: `pd.read_csv(path, dtype={'gene': str})`.

**Detection**: Look for datetime values in columns that should contain gene names or
sample identifiers.

### Numeric Precision

- Floating point: 0.1 + 0.2 != 0.3. Use `np.isclose()` for comparisons.
- Rounding: Be consistent. Define rounding rules before analysis.
- Integer overflow: Rare in Python, but watch for 32-bit integer limits in R or
  database imports (max 2,147,483,647).

---

## 7. Recommended Workflow

The recommended end-to-end data cleaning workflow:

1. **Profile**: Run the profiling script. Understand what you have.
2. **Flag**: Identify potential issues. Categorize by type and severity.
3. **Review**: Present flags to the domain expert (you, the researcher).
4. **Approve**: Decide which flags to act on. Document rationale for each decision.
5. **Clean**: Generate and run cleaning code for approved actions only.
6. **Verify**: Compare before/after summaries. Check that cleaning did not introduce
   new problems.
7. **Document**: Save the cleaning log, mapping tables, and decision rationale.
   Include in supplementary materials or methods section.

### Documentation Checklist

- [ ] Number of rows before and after cleaning
- [ ] Number and percentage of missing values per variable (before/after)
- [ ] Outlier handling decisions with justification
- [ ] Duplicate removal count
- [ ] Category mapping tables
- [ ] Imputation method and variables imputed
- [ ] Any variables excluded from analysis and why

---

## 8. Key References

1. Van den Broeck J, Cunningham SA,"; R,"; AB. Data cleaning: detecting,
   diagnosing, and editing data abnormalities. *PLoS Med*. 2005;2(10):e267.
   DOI: 10.1371/journal.pmed.0020267

2. Kang H. The prevention and handling of the missing data.
   *Korean J Anesthesiol*. 2013;64(5):402-406.
   DOI: 10.4097/kjae.2013.64.5.402

3. Sterne JAC, White IR, Carlin JB, et al. Multiple imputation for missing data
   in epidemiological and clinical research: potential and pitfalls.
   *BMJ*. 2009;338:b2393.
   DOI: 10.1136/bmj.b2393

4. Altman DG, Bland JM. Missing data. *BMJ*. 2007;334(7590):424.
   DOI: 10.1136/bmj.38977.682025.2C

5. White IR, Royston P, Wood AM. Multiple imputation using chained equations:
   Issues and guidance for practice. *Stat Med*. 2011;30(4):377-399.
   DOI: 10.1002/sim.4067

6. Ziemann M, Eren Y, El-Osta A. Gene name errors are widespread in the
   scientific literature. *Genome Biol*. 2016;17(1):177.
   DOI: 10.1186/s13059-016-1044-7

---

*This reference is part of the clean-data skill for the medical-research-skills package.*

# Implausible-Value & Cross-Field Validity Rules

Reference for the clean-data skill, used in **Stage 2 (Flagging)** item 5 ("Implausible values").
It supplies domain-default **hard plausibility bounds** and **cross-field logical-consistency
rules** so the skill can flag likely data-entry errors *even when the codebook does not state a
valid range* — while never auto-correcting them.

This complements `cleaning_patterns.md` §2 (statistical outlier detection): that section finds
values that are *statistically extreme*; this one finds values that are *biologically impossible
or internally contradictory*. The two are different and must not be conflated (§ "Implausible vs
outlier" below).

---

## Implausible vs outlier (the distinction that drives the action)

- **Implausible (compatible-with-life violation / impossible)** — outside the range a living
  human can take, or an impossible code (age 200, SBP 0 in a live patient, BMI 4, % > 100,
  negative count, future birth date). Almost always a **data-entry / unit / sentinel error**.
  Action: correct from source if available, else set to missing and log — do **not** keep it in
  the analytic distribution. (A frank sentinel like 999/9999 is handled in `cleaning_patterns.md`
  §6; treat it as missing *before* range-checking.)
- **Outlier (statistically extreme but biologically possible)** — a BMI of 55, creatinine 12
  mg/dL in dialysis, WBC 300 in leukemia. Action: **keep**; consider a with/without sensitivity
  analysis (`cleaning_patterns.md` §2). Removing a real extreme value as if it were an error
  biases the result.
- The same variable has **both** a hard implausibility bound and a soft outlier fence; apply the
  hard bound first (error removal), then assess outliers among the survivors.

## Population & codebook first

These are **screening bounds to catch errors, not reference/normal ranges, and not analysis
cut-points.** Always reconcile against the data dictionary (dictionary-first) and the **study
population**: neonatal/pediatric, pregnancy, dialysis, ICU, and post-surgical cohorts legitimately
exceed several adult bounds. When the codebook states a valid range, the codebook wins; use this
table only to fill silence, and present every hit as a **flag for review**, never an auto-fix.

## 1. Single-field hard plausibility bounds (adult human, screening)

> Units matter — a bound violation is often a **unit error** (e.g., height in cm stored as m,
> temperature in °F stored as °C, glucose mg/dL vs mmol/L). Check units before treating as an
> error (`cleaning_patterns.md` §6 Mixed Units).

| Domain | Variable (unit) | Implausible if outside | Notes / common error |
|---|---|---|---|
| Demographics | Age (years) | 0–120 | >120 entry error; negative impossible; check DOB-derived age |
| Demographics | Body weight (kg) | 1–400 | <1 = unit error (g?); adult typically 30–250 |
| Demographics | Height (cm) | 30–250 | values 1.5–2.5 = metres stored as cm; >250 entry error |
| Demographics | BMI (kg/m²) | 8–80 | derive and cross-check vs height/weight (§2) |
| Vitals | Systolic BP (mmHg) | 40–300 | 0 in a live patient = error; >300 entry error |
| Vitals | Diastolic BP (mmHg) | 20–200 | must be < systolic (§2) |
| Vitals | Heart rate (bpm) | 20–300 | 0 in a live patient = error |
| Vitals | Respiratory rate (/min) | 4–80 | |
| Vitals | Temperature (°C) | 30–45 | 95–105 = °F stored as °C; outside 30–45 incompatible with life |
| Vitals | SpO₂ (%) | 50–100 | >100 impossible; <50 rarely survivable, suspect error |
| Labs | Blood pH | 6.6–7.8 | outside is incompatible with life |
| Labs | Glucose (mg/dL) | 10–2000 | ÷18 for mmol/L unit error; <10/>2000 suspect |
| Labs | Sodium (mmol/L) | 100–190 | |
| Labs | Potassium (mmol/L) | 1.5–9.5 | |
| Labs | Creatinine (mg/dL) | 0.1–25 | high real in renal failure (outlier, keep) |
| Labs | Hemoglobin (g/dL) | 2–25 | g/L stored as g/dL = ×10 unit error |
| Labs | WBC (×10⁹/L) | 0–500 | high real in leukemia (outlier) |
| Labs | Platelets (×10⁹/L) | 0–3000 | |
| General | Any percentage / proportion | 0–100 (or 0–1) | >100% or >1 impossible; mixed 0–1 vs 0–100 scales |
| General | Any count / duration / dose | ≥ 0 | negative impossible (except intentionally signed change scores) |

(Extend per dataset from the codebook; keep the table version-controlled with the project.)

## 2. Cross-field logical-consistency rules

Single-field bounds miss contradictions *between* fields. Flag a record when:

- **Temporal ordering**: birth date ≤ index/event date ≤ death/last-seen date; admission ≤
  discharge; surgery date ≤ complication date; first ≤ last visit. Any future date beyond the
  extraction date (except scheduled appointments) is suspect (`cleaning_patterns.md` §4).
- **Derived-vs-source**: recomputed BMI from height/weight matches the stored BMI (±rounding);
  age matches (event date − DOB); eGFR consistent with creatinine/age/sex if both stored;
  a "total" equals the sum of its parts; a subset count ≤ its superset count.
- **Sex-specific / state-specific**: pregnancy, parity, or gestational fields present for
  `sex == male`; PSA recorded for females; menopause age < current age; pediatric-only or
  adult-only fields outside the eligible age.
- **Mutually exclusive / dependency**: `deceased == no` with a death date populated;
  `smoking_status == never` with `pack_years > 0` (vs the structural-zero case where it is NULL —
  see SKILL Stage 2 item 7); treatment-stop date before treatment-start date.
- **Range pair**: a "min" field exceeding its paired "max"; diastolic ≥ systolic.

## 3. Flagging discipline (no auto-fix)

- Emit each violation as a **flag** in the Stage 2 table (`Issue Type = Implausible value` or
  `Cross-field inconsistency`), with the rule that fired, the count, and a **suggested** action —
  never an automatic edit (clean-data philosophy: researcher approves each action).
- Severity: a compatible-with-life violation or a hard cross-field contradiction is **High**; a
  bound exceeded only under a special population (possible unit issue) is **Medium** pending unit
  confirmation.
- Prefer **"unit-error" correction over deletion** when a unit mistake is the likely cause (a
  height of 1.75 in a cm column → ×100), but only after confirming the unit; otherwise set to
  missing and log. Record every decision in the cleaning log (`cleaning_patterns.md` §7).

## 4. Key reference

- Van den Broeck J, et al. Data cleaning: detecting, diagnosing, and editing data abnormalities.
  *PLoS Med* 2005;2(10):e267. DOI: 10.1371/journal.pmed.0020267 (the diagnose-vs-edit framework
  this section operationalizes).

---

*This reference is part of the clean-data skill for the medical-research-skills package. The
bounds are error-screening defaults, not clinical reference ranges — verify against the codebook
and study population.*

#!/usr/bin/env python3
"""
Data Profiling Template for Clinical Research Datasets
======================================================

Generates a structured profile of a CSV or Excel dataset.
Outputs a summary table to the console and saves it as CSV.

Usage:
    python profiling_template.py <file_path> [--output <output_dir>]

Requirements:
    - pandas
    - numpy
    - matplotlib, seaborn (optional, for plots)

This script does NOT modify the input data. It is read-only.
"""

import argparse
import os
import sys
import random
from pathlib import Path

import numpy as np
import pandas as pd

# Reproducibility
np.random.seed(42)
random.seed(42)


# ---------------------------------------------------------------------------
# 1. Data Loading
# ---------------------------------------------------------------------------

def load_data(file_path: str) -> pd.DataFrame:
    """Auto-detect CSV vs Excel and load into a DataFrame."""
    path = Path(file_path)
    ext = path.suffix.lower()

    if ext in (".csv", ".tsv"):
        sep = "\t" if ext == ".tsv" else ","
        df = pd.read_csv(path, sep=sep, low_memory=False)
    elif ext in (".xls", ".xlsx", ".xlsm"):
        df = pd.read_excel(path, engine="openpyxl")
    else:
        raise ValueError(f"Unsupported file format: {ext}. Use CSV, TSV, or Excel.")

    print(f"Loaded {len(df)} rows x {len(df.columns)} columns from {path.name}")
    return df


# ---------------------------------------------------------------------------
# 2. Variable Summary
# ---------------------------------------------------------------------------

def build_variable_summary(df: pd.DataFrame) -> pd.DataFrame:
    """Build a per-variable summary with type, missingness, and descriptive stats."""
    records = []

    for col in df.columns:
        series = df[col]
        n_missing = int(series.isna().sum())
        pct_missing = round(100 * n_missing / len(df), 2) if len(df) > 0 else 0.0
        n_unique = int(series.nunique(dropna=True))
        inferred_type = _infer_variable_type(series)

        rec = {
            "variable": col,
            "dtype": str(series.dtype),
            "inferred_type": inferred_type,
            "n_total": len(df),
            "n_missing": n_missing,
            "pct_missing": pct_missing,
            "n_unique": n_unique,
            "min": None,
            "max": None,
            "mean": None,
            "median": None,
            "sd": None,
        }

        # Numeric descriptive statistics
        if inferred_type == "numeric":
            numeric = pd.to_numeric(series, errors="coerce")
            rec["min"] = round(float(numeric.min()), 4) if numeric.notna().any() else None
            rec["max"] = round(float(numeric.max()), 4) if numeric.notna().any() else None
            rec["mean"] = round(float(numeric.mean()), 4) if numeric.notna().any() else None
            rec["median"] = round(float(numeric.median()), 4) if numeric.notna().any() else None
            rec["sd"] = round(float(numeric.std()), 4) if numeric.notna().any() else None

        records.append(rec)

    summary = pd.DataFrame(records)
    return summary


def _infer_variable_type(series: pd.Series) -> str:
    """Heuristic type inference: numeric, categorical, datetime, or text."""
    if pd.api.types.is_numeric_dtype(series):
        return "numeric"
    if pd.api.types.is_datetime64_any_dtype(series):
        return "datetime"

    # Try to parse as numeric (catches numeric-stored-as-string)
    coerced = pd.to_numeric(series.dropna(), errors="coerce")
    if coerced.notna().sum() > 0.8 * series.dropna().shape[0]:
        return "numeric"

    # Try to parse as datetime
    try:
        parsed = pd.to_datetime(series.dropna(), infer_datetime_format=True, errors="coerce")
        if parsed.notna().sum() > 0.8 * series.dropna().shape[0]:
            return "datetime"
    except Exception:
        pass

    # Categorical vs free text heuristic
    n_unique = series.nunique(dropna=True)
    n_rows = len(series.dropna())
    if n_rows > 0 and n_unique / n_rows < 0.05:
        return "categorical"
    if n_unique <= 20:
        return "categorical"

    return "text"


# ---------------------------------------------------------------------------
# 3. Flag Detection
# ---------------------------------------------------------------------------

def flag_missing(summary: pd.DataFrame, threshold: float = 5.0) -> pd.DataFrame:
    """Flag variables with missing percentage above the threshold."""
    flagged = summary[summary["pct_missing"] > threshold].copy()
    flagged["issue"] = "Missing > " + str(threshold) + "%"
    flagged["severity"] = flagged["pct_missing"].apply(
        lambda x: "High" if x > 30 else ("Medium" if x > 10 else "Low")
    )
    return flagged[["variable", "issue", "n_missing", "pct_missing", "severity"]]


def flag_outliers_iqr(df: pd.DataFrame, summary: pd.DataFrame) -> pd.DataFrame:
    """Flag numeric variables with outliers using the IQR method."""
    results = []
    numeric_vars = summary[summary["inferred_type"] == "numeric"]["variable"].tolist()

    for col in numeric_vars:
        numeric = pd.to_numeric(df[col], errors="coerce").dropna()
        if len(numeric) < 10:
            continue
        q1 = numeric.quantile(0.25)
        q3 = numeric.quantile(0.75)
        iqr = q3 - q1
        if iqr == 0:
            continue
        lower = q1 - 1.5 * iqr
        upper = q3 + 1.5 * iqr
        outliers = numeric[(numeric < lower) | (numeric > upper)]
        if len(outliers) > 0:
            results.append({
                "variable": col,
                "issue": f"Outlier (IQR): {len(outliers)} values outside [{lower:.2f}, {upper:.2f}]",
                "count": len(outliers),
                "severity": "Medium",
            })

    return pd.DataFrame(results) if results else pd.DataFrame(
        columns=["variable", "issue", "count", "severity"]
    )


# ---------------------------------------------------------------------------
# 4. Distribution Plots (optional)
# ---------------------------------------------------------------------------

def plot_distributions(df: pd.DataFrame, summary: pd.DataFrame, output_dir: str):
    """Generate histograms for numeric and bar charts for categorical variables."""
    try:
        import matplotlib
        matplotlib.use("Agg")
        import matplotlib.pyplot as plt
        import seaborn as sns
    except ImportError:
        print("[INFO] matplotlib/seaborn not installed. Skipping distribution plots.")
        return

    plot_dir = os.path.join(output_dir, "profile_plots")
    os.makedirs(plot_dir, exist_ok=True)

    # Numeric histograms
    numeric_vars = summary[summary["inferred_type"] == "numeric"]["variable"].tolist()
    for col in numeric_vars[:20]:  # Limit to first 20 to avoid excessive plots
        fig, ax = plt.subplots(figsize=(6, 4))
        numeric = pd.to_numeric(df[col], errors="coerce").dropna()
        if len(numeric) == 0:
            plt.close(fig)
            continue
        ax.hist(numeric, bins=30, edgecolor="black", alpha=0.7)
        ax.set_title(f"Distribution: {col}")
        ax.set_xlabel(col)
        ax.set_ylabel("Frequency")
        fig.tight_layout()
        fig.savefig(os.path.join(plot_dir, f"hist_{col}.png"), dpi=100)
        plt.close(fig)

    # Categorical bar charts
    cat_vars = summary[summary["inferred_type"] == "categorical"]["variable"].tolist()
    for col in cat_vars[:20]:
        fig, ax = plt.subplots(figsize=(6, 4))
        counts = df[col].value_counts().head(15)
        counts.plot(kind="barh", ax=ax, color="steelblue", edgecolor="black")
        ax.set_title(f"Categories: {col}")
        ax.set_xlabel("Count")
        fig.tight_layout()
        fig.savefig(os.path.join(plot_dir, f"bar_{col}.png"), dpi=100)
        plt.close(fig)

    print(f"[INFO] Distribution plots saved to {plot_dir}/")


# ---------------------------------------------------------------------------
# 5. Report Output
# ---------------------------------------------------------------------------

def print_summary(summary: pd.DataFrame):
    """Print a formatted summary table to the console."""
    display_cols = [
        "variable", "inferred_type", "n_missing", "pct_missing",
        "n_unique", "min", "max", "mean", "median", "sd"
    ]
    print("\n" + "=" * 80)
    print("VARIABLE SUMMARY")
    print("=" * 80)
    print(summary[display_cols].to_string(index=False))
    print("=" * 80)


def save_outputs(summary: pd.DataFrame, flags_missing: pd.DataFrame,
                 flags_outlier: pd.DataFrame, output_dir: str):
    """Save profiling results as CSV files."""
    os.makedirs(output_dir, exist_ok=True)

    summary_path = os.path.join(output_dir, "variable_summary.csv")
    summary.to_csv(summary_path, index=False)
    print(f"[SAVED] Variable summary -> {summary_path}")

    if len(flags_missing) > 0:
        missing_path = os.path.join(output_dir, "flags_missing.csv")
        flags_missing.to_csv(missing_path, index=False)
        print(f"[SAVED] Missing flags -> {missing_path}")

    if len(flags_outlier) > 0:
        outlier_path = os.path.join(output_dir, "flags_outliers.csv")
        flags_outlier.to_csv(outlier_path, index=False)
        print(f"[SAVED] Outlier flags -> {outlier_path}")


# ---------------------------------------------------------------------------
# Main
# ---------------------------------------------------------------------------

def main():
    parser = argparse.ArgumentParser(description="Profile a clinical research dataset.")
    parser.add_argument("file_path", help="Path to CSV or Excel file")
    parser.add_argument("--output", default="./data_profile",
                        help="Output directory for profiling results (default: ./data_profile)")
    parser.add_argument("--no-plots", action="store_true",
                        help="Skip distribution plots")
    args = parser.parse_args()

    # Load
    df = load_data(args.file_path)

    # Profile
    summary = build_variable_summary(df)
    print_summary(summary)

    # Flag
    flags_missing = flag_missing(summary, threshold=5.0)
    flags_outlier = flag_outliers_iqr(df, summary)

    if len(flags_missing) > 0:
        print("\n[FLAGS] Variables with >5% missing:")
        print(flags_missing.to_string(index=False))

    if len(flags_outlier) > 0:
        print("\n[FLAGS] Variables with IQR outliers:")
        print(flags_outlier.to_string(index=False))

    # Save
    save_outputs(summary, flags_missing, flags_outlier, args.output)

    # Plot (optional)
    if not args.no_plots:
        plot_distributions(df, summary, args.output)

    print("\n[DONE] Profiling complete. Review outputs before proceeding to cleaning.")


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""Reverse-coded-item detector for multi-item scale reliability (clean-data Stage 2).

A multi-item scale (Likert Trust/Satisfaction/Burden, etc.) mixes positively- and
negatively-worded items. A negatively-worded ("reverse") item must be recoded
`(min+max) - x` before the scale total or Cronbach's alpha is computed. When that
recoding is skipped, the reverse item correlates negatively with the rest of the
scale and Cronbach's alpha collapses — often turning *negative*. A negative alpha
is almost never a real measurement phenomenon: it is a reverse-coding bug. Authors
who instead defend it as "multidimensional structure" lose a review round (the
motivating incident: a Trust scale shipped at alpha = -0.57 until one item was
recoded, after which alpha = 0.58).

This script flags reverse-code suspects *before* alpha is reported, so they are
recoded at cleaning time rather than mis-explained at revision time.

INPUTS
  --data        CSV with one row per respondent.
  --items       scale item columns (repeatable / space-separated). >= 2 required.
  --min         lowest point of the response scale (default 1).
  --max         highest point of the response scale (default: inferred per item
                from the observed maximum across all items).
  --threshold   item-rest correlation at or below which an item is a reverse-code
                suspect (default 0.0 — i.e. a negative item-rest correlation).

OUTPUT  (per scale)
  - alpha_raw            Cronbach's alpha on the items AS-GIVEN (un-recoded)
  - per item: item_rest_r (corrected item-total correlation)
  - suspects             items with item_rest_r <= threshold
  - verdict:
      REVERSE_CODING_LIKELY   alpha_raw < 0           (recode then re-run)
      REVERSE_CODING_SUSPECT  suspects present, alpha_raw >= 0
      OK                      no suspects, alpha_raw >= 0
  Exit 1 (with --strict) if verdict != OK.

Stdlib-only (csv / json / argparse / math / statistics). Exit codes: 0 clean
(or report-only), 1 reverse-coding flagged (with --strict), 2 input/usage error.
"""

from __future__ import annotations

import argparse
import csv
import json
import math
import sys
from pathlib import Path

NULLISH = {"", "na", "n/a", "nan", "null", "none", ".", "missing"}


def _to_float(raw: str):
    if raw is None:
        return None
    s = raw.strip()
    if s.lower() in NULLISH:
        return None
    try:
        return float(s)
    except ValueError:
        return None


def _pearson(xs: list[float], ys: list[float]):
    """Pearson r over paired complete observations. None if undefined."""
    n = len(xs)
    if n < 2:
        return None
    mx = sum(xs) / n
    my = sum(ys) / n
    sxx = sum((x - mx) ** 2 for x in xs)
    syy = sum((y - my) ** 2 for y in ys)
    if sxx <= 0 or syy <= 0:        # a constant column has no correlation
        return None
    sxy = sum((x - mx) * (y - my) for x, y in zip(xs, ys))
    return sxy / math.sqrt(sxx * syy)


def _variance(vals: list[float]):
    n = len(vals)
    if n < 2:
        return 0.0
    m = sum(vals) / n
    return sum((v - m) ** 2 for v in vals) / (n - 1)


def cronbach_alpha(rows: list[list[float]]):
    """rows = list of complete item vectors (one per respondent). None if undefined."""
    if not rows:
        return None
    k = len(rows[0])
    if k < 2:
        return None
    item_vars = [_variance([r[j] for r in rows]) for j in range(k)]
    totals = [sum(r) for r in rows]
    total_var = _variance(totals)
    if total_var <= 0:
        return None
    return (k / (k - 1.0)) * (1.0 - sum(item_vars) / total_var)


def analyze(data_path: Path, items: list[str], scale_min: float,
            scale_max, threshold: float) -> dict:
    with data_path.open(newline="", encoding="utf-8-sig") as fh:
        reader = csv.DictReader(fh)
        header = reader.fieldnames or []
        missing = [c for c in items if c not in header]
        if missing:
            raise SystemExit(f"USAGE-ERR: item columns not in CSV: {missing}")
        raw_rows = list(reader)

    # Listwise-complete matrix (drop any respondent missing >=1 item) — matches
    # the default complete-case alpha so the flag mirrors the reported alpha.
    complete: list[list[float]] = []
    n_total = len(raw_rows)
    for row in raw_rows:
        vec = [_to_float(row.get(c)) for c in items]
        if all(v is not None for v in vec):
            complete.append(vec)  # type: ignore[arg-type]

    n_used = len(complete)
    alpha_raw = cronbach_alpha(complete) if n_used >= 2 else None

    # Per-item corrected item-total (item-rest) correlation.
    per_item = []
    suspects = []
    k = len(items)
    for j, name in enumerate(items):
        item_rest_r = None
        if n_used >= 2 and k >= 2:
            xs = [r[j] for r in complete]
            rest = [sum(r[t] for t in range(k) if t != j) for r in complete]
            item_rest_r = _pearson(xs, rest)
        flagged = item_rest_r is not None and item_rest_r <= threshold
        if flagged:
            suspects.append(name)
        per_item.append({
            "item": name,
            "item_rest_r": None if item_rest_r is None else round(item_rest_r, 3),
            "reverse_suspect": bool(flagged),
        })

    if alpha_raw is not None and alpha_raw < 0:
        verdict = "REVERSE_CODING_LIKELY"
    elif suspects:
        verdict = "REVERSE_CODING_SUSPECT"
    else:
        verdict = "OK"

    return {
        "items": items,
        "n_total": n_total,
        "n_complete": n_used,
        "scale_min": scale_min,
        "scale_max": scale_max,  # echoed only; flag is correlation-based
        "threshold": threshold,
        "alpha_raw": None if alpha_raw is None else round(alpha_raw, 3),
        "per_item": per_item,
        "suspects": suspects,
        "verdict": verdict,
        "recode_hint": (
            None if verdict == "OK"
            else "Recode suspect items as (min+max)-x, then recompute alpha "
                 "(see ~/.claude/rules/survey-scale-reliability.md)."
        ),
    }


def main() -> int:
    ap = argparse.ArgumentParser(description="Reverse-coded-item / negative-alpha detector.")
    ap.add_argument("--data", required=True)
    ap.add_argument("--items", nargs="+", required=True)
    ap.add_argument("--min", type=float, default=1.0, dest="scale_min")
    ap.add_argument("--max", type=float, default=None, dest="scale_max")
    ap.add_argument("--threshold", type=float, default=0.0)
    ap.add_argument("--out", default=None, help="write JSON report to this path")
    ap.add_argument("--strict", action="store_true",
                    help="exit 1 if any reverse-coding is flagged")
    args = ap.parse_args()

    if len(args.items) < 2:
        print("USAGE-ERR: need >= 2 scale items", file=sys.stderr)
        return 2

    data_path = Path(args.data)
    if not data_path.is_file():
        print(f"USAGE-ERR: data file not found: {data_path}", file=sys.stderr)
        return 2

    report = analyze(data_path, args.items, args.scale_min, args.scale_max, args.threshold)

    if args.out:
        Path(args.out).write_text(json.dumps({"detector": "check_reverse_coding", **report}, indent=2), encoding="utf-8")

    a = report["alpha_raw"]
    print(f"Scale items: {', '.join(report['items'])}  (n={report['n_complete']} complete)")
    print(f"Cronbach's alpha (as-given): {a}")
    for it in report["per_item"]:
        mark = "  <-- reverse-code suspect" if it["reverse_suspect"] else ""
        print(f"  item-rest r  {it['item']:<16} {it['item_rest_r']}{mark}")
    print(f"VERDICT: {report['verdict']}")
    if report["recode_hint"]:
        print(f"  {report['recode_hint']}")

    if args.strict and report["verdict"] != "OK":
        return 1
    return 0


if __name__ == "__main__":
    sys.exit(main())

#!/usr/bin/env python3
"""Categorical-implied-zero (structural-zero) detector for clean-data Stage 2.

A dose/duration variable anchored to a categorical exposure has a known zero at
the reference level: a never-smoker's pack-years is 0 *by definition*, not
missing. When that implied zero is instead stored as NULL/blank, two downstream
failures follow — a complete-case model silently drops the whole unexposed
stratum, or MICE imputes a non-zero dose for people who have no exposure — both
of which corrupt the exposure contrast. This script finds those contradiction
rows so they are fixed at cleaning time (set 0 at the reference level; impute
only the residual missingness among the exposed).

INPUTS
  --data            CSV with one row per subject.
  --category-col    categorical exposure column (e.g. smoking_status).
  --reference-level value of the category that implies a zero dose (e.g. never).
  --dose-col        dose/duration column that should be 0 at the reference level
                    (e.g. pack_years). Repeatable for several dose columns.

OUTPUT
  Per (reference-level, dose) pair: counts of
    - implied_zero_missing : category == reference AND dose is NULL/blank   (FIX)
    - implied_zero_nonzero : category == reference AND dose > 0  (mislabeled — review)
    - implied_zero_ok      : category == reference AND dose == 0
  Exit 1 (with --strict) if any implied_zero_missing rows exist.

Stdlib-only (csv / json / argparse). Exit codes: 0 clean (or report-only),
1 contradiction rows found (with --strict), 2 input/usage error.
"""

from __future__ import annotations

import argparse
import csv
import json
import sys
from pathlib import Path

NULLISH = {"", "na", "n/a", "nan", "null", "none", ".", "missing"}


def _is_null(v: str) -> bool:
    return v.strip().lower() in NULLISH


def _to_float(v: str):
    try:
        return float(v.strip())
    except (ValueError, AttributeError):
        return None


def _norm(s: str) -> str:
    return s.strip().lower()


def analyze(data: str, category_col: str, reference: str, dose_cols: list[str]) -> dict:
    p = Path(data)
    if not p.is_file():
        sys.stderr.write(f"ERROR: data file not found: {data}\n")
        sys.exit(2)
    with p.open(encoding="utf-8-sig", newline="") as f:
        reader = csv.DictReader(f)
        if reader.fieldnames is None:
            sys.stderr.write("ERROR: empty CSV\n")
            sys.exit(2)
        fields = {_norm(c): c for c in reader.fieldnames}
        cat_key = fields.get(_norm(category_col))
        if cat_key is None:
            sys.stderr.write(f"ERROR: category column '{category_col}' not in {reader.fieldnames}\n")
            sys.exit(2)
        dose_keys = {}
        for d in dose_cols:
            k = fields.get(_norm(d))
            if k is None:
                sys.stderr.write(f"ERROR: dose column '{d}' not in {reader.fieldnames}\n")
                sys.exit(2)
            dose_keys[d] = k
        rows = list(reader)

    ref = _norm(reference)
    results = []
    total_missing = 0
    for d, dkey in dose_keys.items():
        miss = nonzero = ok = ref_n = 0
        for r in rows:
            if _norm(r.get(cat_key, "")) != ref:
                continue
            ref_n += 1
            raw = r.get(dkey, "")
            if _is_null(raw):
                miss += 1
            else:
                val = _to_float(raw)
                if val is None:
                    continue            # non-numeric, not our concern here
                if val > 0:
                    nonzero += 1
                else:
                    ok += 1
        total_missing += miss
        results.append({
            "dose_col": d,
            "reference_level": reference,
            "reference_n": ref_n,
            "implied_zero_missing": miss,
            "implied_zero_nonzero": nonzero,
            "implied_zero_ok": ok,
            "verdict": "FIX_STRUCTURAL_ZERO" if miss else ("REVIEW_MISLABEL" if nonzero else "OK"),
        })

    return {
        "data": str(p),
        "category_col": category_col,
        "n_rows": len(rows),
        "results": results,
        "total_implied_zero_missing": total_missing,
        "suggested_fix": (
            f"Set dose = 0 where {category_col} == '{reference}' (structural zero), "
            "then impute only the residual missingness among the exposed."
        ) if total_missing else None,
    }


def render(result: dict) -> str:
    lines = [
        "| Dose col | ref n | implied-zero MISSING | nonzero (mislabel) | zero (ok) | Verdict |",
        "|---|---|---|---|---|---|",
    ]
    for r in result["results"]:
        mark = {"FIX_STRUCTURAL_ZERO": "✗ Fix", "REVIEW_MISLABEL": "△ Review", "OK": "✓"}[r["verdict"]]
        lines.append(
            f"| {r['dose_col']} | {r['reference_n']} | {r['implied_zero_missing']} | "
            f"{r['implied_zero_nonzero']} | {r['implied_zero_ok']} | {mark} |"
        )
    return "\n".join(lines)


def main() -> int:
    ap = argparse.ArgumentParser(description="Categorical-implied-zero (structural-zero) detector.")
    ap.add_argument("--data", required=True, help="subject-level CSV")
    ap.add_argument("--category-col", required=True, help="categorical exposure column")
    ap.add_argument("--reference-level", required=True, help="category value implying a zero dose")
    ap.add_argument("--dose-col", required=True, action="append",
                    help="dose/duration column (repeatable)")
    ap.add_argument("--out", help="write JSON artifact to this path")
    ap.add_argument("--strict", action="store_true", help="exit 1 if any implied-zero-missing rows")
    args = ap.parse_args()

    result = analyze(args.data, args.category_col, args.reference_level, args.dose_col)

    print("=" * 41)
    print(" Categorical-Implied-Zero (structural zero)")
    print("=" * 41)
    print(f"category: {args.category_col} == '{args.reference_level}'")
    print(render(result))
    print()
    if result["total_implied_zero_missing"]:
        print(f"FIX: {result['total_implied_zero_missing']} reference-level row(s) store the implied "
              f"zero as missing.")
        print(result["suggested_fix"])
    else:
        print("OK: no categorical-implied-zero stored as missing.")

    if args.out:
        Path(args.out).parent.mkdir(parents=True, exist_ok=True)
        Path(args.out).write_text(json.dumps({"detector": "check_structural_zero", **result}, indent=2), encoding="utf-8")
        print(f"\nwrote {args.out}")

    return 1 if (args.strict and result["total_implied_zero_missing"]) else 0


if __name__ == "__main__":
    sys.exit(main())

#!/usr/bin/env bash
# Regression test for the reverse-coded-item / negative-alpha detector.
# Synthetic fixture: a 3-item scale where E3 is reverse-worded (stored 6-x) so the
# raw alpha is negative and E3 has a negative item-rest correlation; plus a clean
# 3-item scale (G1-G3) that should pass. Stdlib-only (python3).
set -u

HERE="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
SCRIPT="$HERE/../scripts/check_reverse_coding.py"
FIXTURE="$HERE/fixtures/scale_reverse.csv"
OUT="$(mktemp -t rc_XXXX).json"
trap 'rm -f "$OUT"' EXIT

fail=0
check() { local label="$1"; shift
    if "$@" >/dev/null 2>&1; then printf '  PASS  %s\n' "$label"
    else printf '  FAIL  %s\n' "$label"; fail=$((fail+1)); fi
}

[[ -f "$SCRIPT" ]]  || { echo "ENV-ERR: script missing"  >&2; exit 2; }
[[ -f "$FIXTURE" ]] || { echo "ENV-ERR: fixture missing" >&2; exit 2; }

# Reverse-coded scale: must flag under --strict (exit 1) and write the report.
python3 "$SCRIPT" --data "$FIXTURE" --items E1 E2 E3 --out "$OUT" --strict >/dev/null 2>&1
check "exit 1 under --strict (reverse coding present)" test "$?" -eq 1
check "JSON artifact written" test -s "$OUT"

check "verdict REVERSE_CODING_LIKELY" python3 -c "
import json; d=json.load(open('$OUT'))
assert d['verdict']=='REVERSE_CODING_LIKELY', d['verdict']"
check "alpha_raw is negative" python3 -c "
import json; d=json.load(open('$OUT'))
assert d['alpha_raw'] is not None and d['alpha_raw']<0, d['alpha_raw']"
check "E3 is the sole suspect" python3 -c "
import json; d=json.load(open('$OUT'))
assert d['suspects']==['E3'], d['suspects']"
check "E3 item-rest r is negative" python3 -c "
import json; d=json.load(open('$OUT'))
r={it['item']:it['item_rest_r'] for it in d['per_item']}
assert r['E3']<0 and r['E1']>0 and r['E2']>0, r"
check "n_complete == 8" python3 -c "
import json; d=json.load(open('$OUT'))
assert d['n_complete']==8, d['n_complete']"

# Clean scale: no reverse item -> verdict OK, exit 0 under --strict.
python3 "$SCRIPT" --data "$FIXTURE" --items G1 G2 G3 --out "$OUT" --strict >/dev/null 2>&1
check "exit 0 on a clean (aligned) scale" test "$?" -eq 0
check "verdict OK on clean scale" python3 -c "
import json; d=json.load(open('$OUT'))
assert d['verdict']=='OK' and not d['suspects'], d"

# Usage guard: a single item is not a scale.
python3 "$SCRIPT" --data "$FIXTURE" --items E1 >/dev/null 2>&1
check "exit 2 on <2 items (usage error)" test "$?" -eq 2

echo "fail=$fail"; [[ "$fail" -eq 0 ]] && echo "ALL PASS" || echo "FAILURES: $fail"
exit "$fail"

#!/usr/bin/env bash
# Regression test for the categorical-implied-zero (structural-zero) detector.
# Synthetic fixture: never-smokers with NULL pack-years (the bug), one with an
# explicit 0 (ok), one mislabeled with a positive dose. Stdlib-only (python3).
set -u

HERE="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
SCRIPT="$HERE/../scripts/check_structural_zero.py"
FIXTURE="$HERE/fixtures/smoking.csv"
OUT="$(mktemp -t sz_XXXX).json"
trap 'rm -f "$OUT"' EXIT

fail=0
check() { local label="$1"; shift
    if "$@" >/dev/null 2>&1; then printf '  PASS  %s\n' "$label"
    else printf '  FAIL  %s\n' "$label"; fail=$((fail+1)); fi
}

[[ -f "$SCRIPT" ]]  || { echo "ENV-ERR: script missing"  >&2; exit 2; }
[[ -f "$FIXTURE" ]] || { echo "ENV-ERR: fixture missing" >&2; exit 2; }

python3 "$SCRIPT" --data "$FIXTURE" --category-col smoking_status \
    --reference-level never --dose-col pack_years --out "$OUT" --strict >/dev/null 2>&1
check "exit 1 under --strict (implied-zero missing present)" test "$?" -eq 1
check "JSON artifact written" test -s "$OUT"

check "implied_zero_missing == 2" python3 -c "
import json; d=json.load(open('$OUT'))
assert d['total_implied_zero_missing']==2, d['total_implied_zero_missing']"
check "implied_zero_nonzero (mislabel) == 1" python3 -c "
import json; d=json.load(open('$OUT'))
assert d['results'][0]['implied_zero_nonzero']==1"
check "implied_zero_ok == 1" python3 -c "
import json; d=json.load(open('$OUT'))
assert d['results'][0]['implied_zero_ok']==1"
check "reference_n == 4 (never-smokers)" python3 -c "
import json; d=json.load(open('$OUT'))
assert d['results'][0]['reference_n']==4"
check "verdict FIX_STRUCTURAL_ZERO" python3 -c "
import json; d=json.load(open('$OUT'))
assert d['results'][0]['verdict']=='FIX_STRUCTURAL_ZERO'"

# Clean case: a reference level with no missing dose -> exit 0 under --strict.
python3 "$SCRIPT" --data "$FIXTURE" --category-col smoking_status \
    --reference-level former --dose-col pack_years --strict >/dev/null 2>&1
check "exit 0 when reference level has no missing dose" test "$?" -eq 0

echo "fail=$fail"; [[ "$fail" -eq 0 ]] && echo "ALL PASS" || echo "FAILURES: $fail"
exit "$fail"