import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
from scipy.stats import shapiro, jarque_bera

# Load the diabetes dataset
diabetes = pd.read_csv('diabetes.csv')

# First look at our data
print("Dataset shape:", diabetes.shape)
print("\nFirst few rows:")
print(diabetes.head())

# Check for data quality issues
print("\nData types:")
print(diabetes.dtypes)

# Look for missing values
print(f"\nMissing values: {diabetes.isnull().sum().sum()}")

# Basic statistics
print("\nSummary statistics:")
print(diabetes.describe())

# These variables should NEVER be zero in living patients
impossible_zero_vars = ['Glucose', 'BloodPressure', 'SkinThickness', 'Insulin', 'BMI']

print("IMPOSSIBLE ZERO VALUES DETECTED:")
print("=" * 50)

total_impossible = 0
for var in impossible_zero_vars:
    zero_count = (diabetes[var] == 0).sum()
    zero_pct = round(100 * zero_count / len(diabetes), 1)
    
    # Severity assessment
    if zero_pct > 20:
        severity = "CRITICAL"
    elif zero_pct > 5:
        severity = "MODERATE" 
    else:
        severity = "MINOR"
    
    print(f"  {var:15}: {zero_count:3d} zeros ({zero_pct:4.1f}%) - {severity}")
    total_impossible += zero_count

print(f"\nTotal impossible values: {total_impossible:,}")
print(f"Percentage of affected measurements: {100*total_impossible/(len(diabetes)*len(impossible_zero_vars)):.1f}%")

# Load the massive cardiovascular dataset
cardio = pd.read_csv('cardio_train.csv', sep=';')
print(f"Full dataset: {cardio.shape[0]:,} patients")

# Create derived clinical variables
cardio['age_years'] = cardio['age'] / 365.25  # Convert days to years
cardio['bmi'] = cardio['weight'] / (cardio['height']/100)**2  # Calculate BMI
cardio['pulse_pressure'] = cardio['ap_hi'] - cardio['ap_lo']  # Pulse pressure

# Smart sampling for visualization (maintains statistical properties)
np.random.seed(123)  # Reproducible results
sample_size = 3000
cardio_sample = cardio.sample(sample_size)

print(f"Sample for analysis: {sample_size:,} patients")
print("Sample maintains population characteristics")

def analyze_clinical_normality(data, variable_name, remove_zeros=False, 
                                 clinical_context=None, save_plots=True):
    """
    Comprehensive normality analysis for clinical data
    
    Parameters:
    -----------
    data : pandas.Series
        Clinical data to analyze
    variable_name : str
        Name of the clinical variable
    remove_zeros : bool
        Remove physiologically impossible zero values
    clinical_context : str
        Clinical context for interpretation
    save_plots : bool
        Save generated plots
    
    Returns:
    --------
    dict : Comprehensive analysis results
    """
    
    import numpy as np
    import pandas as pd
    import matplotlib.pyplot as plt
    from scipy import stats
    from scipy.stats import shapiro, jarque_bera
    
    # Data preparation
    if remove_zeros:
        original_size = len(data)
        clean_data = data[data > 0]
        removed_zeros = original_size - len(clean_data)
        print(f"Removed {removed_zeros} impossible zero values")
    else:
        clean_data = data.dropna()
    
    # Create comprehensive visualization
    fig, axes = plt.subplots(2, 3, figsize=(16, 10))
    fig.suptitle(f'Clinical Normality Analysis: {variable_name}', 
                fontsize=16, fontweight='bold')
    
    # 1. Q-Q plot
    stats.probplot(clean_data, dist="norm", plot=axes[0, 0])
    axes[0, 0].set_title('Q-Q Plot vs Normal Distribution')
    axes[0, 0].grid(True, alpha=0.3)
    
    # 2. Histogram with normal overlay
    axes[0, 1].hist(clean_data, bins=30, density=True, alpha=0.7, 
                   color='lightblue', edgecolor='black')
    mu, sigma = clean_data.mean(), clean_data.std()
    x = np.linspace(clean_data.min(), clean_data.max(), 100)
    axes[0, 1].plot(x, stats.norm.pdf(x, mu, sigma), 'r-', lw=2, 
                   label=f'Normal (mean={mu:.2f}, std={sigma:.2f})')
    axes[0, 1].set_title('Histogram with Normal Overlay')
    axes[0, 1].legend()
    axes[0, 1].grid(True, alpha=0.3)
    
    # 3. Box plot
    axes[0, 2].boxplot(clean_data, patch_artist=True)
    axes[0, 2].set_title('Box Plot with Outliers')
    axes[0, 2].grid(True, alpha=0.3)
    
    # 4. Log-transformed Q-Q (if applicable)
    if (clean_data > 0).all():
        log_data = np.log(clean_data)
        stats.probplot(log_data, dist="norm", plot=axes[1, 0])
        axes[1, 0].set_title('Log-transformed Q-Q Plot')
        log_shapiro_p = shapiro(log_data)[1]
        axes[1, 0].text(0.05, 0.95, f'Shapiro p: {log_shapiro_p:.4f}',
                       transform=axes[1, 0].transAxes, 
                       bbox=dict(boxstyle='round', facecolor='lightgreen'))
    else:
        axes[1, 0].text(0.5, 0.5, 'Cannot log-transform\n(non-positive values)', 
                       ha='center', va='center')
        axes[1, 0].set_title('Log-transformed Q-Q Plot')
    
    # 5. Square root-transformed Q-Q
    sqrt_data = np.sqrt(np.abs(clean_data))
    stats.probplot(sqrt_data, dist="norm", plot=axes[1, 1])
    axes[1, 1].set_title('Square Root-transformed Q-Q Plot')
    axes[1, 1].grid(True, alpha=0.3)
    
    # 6. Statistical summary
    axes[1, 2].axis('off')
    
    # Calculate comprehensive statistics
    n = len(clean_data)
    mean_val = clean_data.mean()
    median_val = clean_data.median()
    std_val = clean_data.std()
    skewness = stats.skew(clean_data)
    kurtosis = stats.kurtosis(clean_data)
    
    # Normality tests
    shapiro_stat, shapiro_p = shapiro(clean_data) if n <= 5000 else (np.nan, np.nan)
    jb_stat, jb_p = jarque_bera(clean_data)
    
    # Clinical interpretation
    if shapiro_p > 0.05 and jb_p > 0.05:
        conclusion = "Approximately Normal"
        recommendation = "Parametric tests appropriate"
    else:
        conclusion = "Non-Normal"
        recommendation = "Consider transformations or non-parametric tests"
    
    # Summary text
    summary = f"""CLINICAL ANALYSIS SUMMARY
{'='*30}
Sample size: {n:,}
Mean: {mean_val:.2f}
Median: {median_val:.2f}
Std Dev: {std_val:.2f}
Skewness: {skewness:.2f}
Kurtosis: {kurtosis:.2f}

NORMALITY TESTS
{'='*15}
Shapiro-Wilk p: {shapiro_p:.6f}
Jarque-Bera p: {jb_p:.6f}

CONCLUSION: {conclusion}
RECOMMENDATION: {recommendation}"""
    
    if clinical_context:
        summary += f"\n\nCLINICAL CONTEXT:\n{clinical_context}"
    
    axes[1, 2].text(0.05, 0.95, summary, transform=axes[1, 2].transAxes,
                   fontsize=9, verticalalignment='top', fontfamily='monospace')
    
    plt.tight_layout()
    
    if save_plots:
        filename = f'clinical_analysis_{variable_name.lower().replace(" ", "_")}.png'
        plt.savefig(filename, dpi=300, bbox_inches='tight')
        print(f"Saved: {filename}")
    
    plt.show()
    
    return {
        'variable': variable_name,
        'n': n,
        'normality_conclusion': conclusion,
        'shapiro_p': shapiro_p,
        'jb_p': jb_p,
        'recommendation': recommendation
    }

# Load your clinical dataset
diabetes = pd.read_csv('diabetes.csv')

# Analyze different types of clinical variables
age_results = analyze_clinical_normality(
    diabetes['Age'], 
    'Age',
    clinical_context="Age in cross-sectional diabetes study"
)

insulin_results = analyze_clinical_normality(
    diabetes['Insulin'], 
    'Insulin',
    remove_zeros=True,  # Remove impossible zeros!
    clinical_context="Biomarker with high inter-individual variability"
)

bmi_results = analyze_clinical_normality(
    diabetes['BMI'], 
    'BMI',
    remove_zeros=True,
    clinical_context="Anthropometric measurement, right-skewed in populations"
)

# Quick comparison
variables = ['Age', 'Insulin', 'BMI']
results = [age_results, insulin_results, bmi_results]

print("\nCLINICAL NORMALITY SUMMARY")
print("=" * 50)
for var, result in zip(variables, results):
    print(f"{var:10}: {result['normality_conclusion']} (p={result['shapiro_p']:.4f})")