import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# Set visualization style for professional graphs
sns.set_theme(style="whitegrid")

# 1. Extraction: Load the 2021-2023 data directly from the NEW CDC URLs
base_url = "https://wwwn.cdc.gov/Nchs/Data/Nhanes/Public/2021/DataFiles/"

print("Downloading and reading 2021-2023 SAS XPT files from CDC...")
df_demo = pd.read_sas(base_url + 'DEMO_L.xpt') # Demographics
df_body = pd.read_sas(base_url + 'BMX_L.xpt')  # Body Measures
df_sleep = pd.read_sas(base_url + 'SLQ_L.xpt') # Sleep
df_activity = pd.read_sas(base_url + 'PAQ_L.xpt') # Physical Activity
df_depr = pd.read_sas(base_url + 'DPQ_L.xpt')  # Depression Screener (PHQ-9)
df_crp = pd.read_sas(base_url + 'HSCRP_L.xpt') # NEW: High-Sensitivity C-Reactive Protein (Inflammation)

# 2. Transform: Merge all 6 datasets on the unique sequence number (SEQN)
df = pd.merge(df_demo, df_body, on='SEQN', how='inner')
df = pd.merge(df, df_sleep, on='SEQN', how='inner')
df = pd.merge(df, df_activity, on='SEQN', how='inner')
df = pd.merge(df, df_depr, on='SEQN', how='inner')
df = pd.merge(df, df_crp, on='SEQN', how='inner') # Merging new lab data

# 3. Transform: Calculate PHQ-9 Depression Score
phq_cols = [f'DPQ0{i}0' for i in range(1, 10)]
df[phq_cols] = df[phq_cols].replace([7, 9], np.nan) 
df['PHQ9_Score'] = df[phq_cols].sum(axis=1)

# 4. Transform: Rename to readable columns and filter adults (20+)
col_mapping = {
    'SEQN': 'ID',
    'RIAGENDR': 'Gender',
    'RIDAGEYR': 'Age',
    'BMXBMI': 'BMI',
    'SLD012': 'Sleep_Hours',           
    'PAD820': 'Vigorous_Activity_Min', 
    'INDFMPIR': 'Income_Poverty_Ratio',
    'LBXHSCRP': 'CRP_mgL' # NEW: Clinical inflammation marker
}
df_clean = df.rename(columns=col_mapping)[['ID', 'Gender', 'Age', 'BMI', 'Sleep_Hours', 
                                           'Vigorous_Activity_Min', 'Income_Poverty_Ratio', 'CRP_mgL', 'PHQ9_Score']]
df_clean = df_clean[df_clean['Age'] >= 20]

# Map Gender
df_clean['Gender'] = df_clean['Gender'].map({1.0: 'Male', 2.0: 'Female'})

# Create Binary Target: Clinical cutoff for depression is a score >= 10
df_clean['Depressed'] = df_clean['PHQ9_Score'].apply(lambda x: 1 if x >= 10 else 0 if pd.notnull(x) else np.nan)

# 5. Handle NaNs & Set Dtypes
df_clean = df_clean.dropna(subset=['PHQ9_Score']) 

df_clean['Vigorous_Activity_Min'] = df_clean['Vigorous_Activity_Min'].fillna(0)

# Median imputation for remaining features including the new CRP clinical data
for col in ['BMI', 'Sleep_Hours', 'Income_Poverty_Ratio', 'CRP_mgL']:
    df_clean[col] = df_clean[col].fillna(df_clean[col].median())

# Set logical datatypes
df_clean = df_clean.astype({
    'ID': 'int64',
    'Age': 'int64',
    'Sleep_Hours': 'float64',
    'Vigorous_Activity_Min': 'float64',
    'CRP_mgL': 'float64',
    'PHQ9_Score': 'int64',
    'Depressed': 'int64'
})

print("\nETL Complete. Dataset Shape:", df_clean.shape)
display(df_clean.head())

Downloading and reading 2021-2023 SAS XPT files from CDC...

ETL Complete. Dataset Shape: (6064, 10)

# --- Summary Statistics (5 generated) ---
print("--- EDA Summary Statistics ---\n")

# 1. Base rate of depression
dep_rate = df_clean['Depressed'].mean() * 100
print(f"1. Proportion of cohort screening positive for depression (PHQ9 >= 10): {dep_rate:.2f}%")

# 2. Median Activity
med_activity = df_clean['Vigorous_Activity_Min'].median()
print(f"2. Median vigorous physical activity (minutes/day) across cohort: {med_activity:.1f}")

# 3. Median Income Ratio
med_income = df_clean['Income_Poverty_Ratio'].median()
print(f"3. Median Income-to-Poverty Ratio (higher = greater financial stability): {med_income:.2f}")

# 4 & 5. Grouped statistic: Average Activity by Depression Status
activity_by_dep = df_clean.groupby('Depressed')['Vigorous_Activity_Min'].mean()
print(f"4. Average vigorous activity for NON-depressed individuals: {activity_by_dep[0]:.1f} mins/day")
print(f"5. Average vigorous activity for DEPRESSED individuals: {activity_by_dep[1]:.1f} mins/day\n")


# --- Graphic: Correlation Heatmap ---
plt.figure(figsize=(10, 8))

# Select only the continuous numerical variables for the heatmap
numeric_cols = ['Age', 'BMI', 'Sleep_Hours', 'Vigorous_Activity_Min', 'Income_Poverty_Ratio', 'PHQ9_Score']
corr_matrix = df_clean[numeric_cols].corr()

# Generate a mask for the upper triangle to make it easier to read
mask = np.triu(np.ones_like(corr_matrix, dtype=bool))

# Draw the heatmap
sns.heatmap(corr_matrix, mask=mask, annot=True, fmt=".2f", cmap='coolwarm', 
            vmax=0.3, vmin=-0.3, square=True, linewidths=.5, cbar_kws={"shrink": .8})

plt.title('Correlation Heatmap: Lifestyle, Socioeconomic Factors, and Depression (2021-2023)', fontsize=14, pad=20)
plt.show()

--- EDA Summary Statistics ---

1. Proportion of cohort screening positive for depression (PHQ9 >= 10): 11.05%
2. Median vigorous physical activity (minutes/day) across cohort: 0.0
3. Median Income-to-Poverty Ratio (higher = greater financial stability): 2.82
4. Average vigorous activity for NON-depressed individuals: 37.8 mins/day
5. Average vigorous activity for DEPRESSED individuals: 66.8 mins/day

# Set custom aesthetic palette
custom_palette = ['#232D4B', '#E57200']
sns.set_palette(custom_palette)

fig, axes = plt.subplots(2, 2, figsize=(16, 12))
fig.suptitle('Exploratory Data Analysis: Distributions and Relationships', fontsize=16, fontweight='bold', y=1.02)

# Graph 2: Distribution of PHQ-9 Scores
sns.histplot(data=df_clean, x='PHQ9_Score', bins=27, kde=True, color='#232D4B', ax=axes[0, 0])
axes[0, 0].set_title('Graph 2: Distribution of PHQ-9 Depression Scores', fontsize=12)
axes[0, 0].set_xlabel('PHQ-9 Score')
axes[0, 0].set_ylabel('Count of Individuals')

# Graph 3: Depression Base Rate by Gender
sns.barplot(data=df_clean, x='Gender', y='Depressed', errorbar=None, palette=custom_palette, hue='Gender', legend=False, ax=axes[0, 1])
axes[0, 1].set_title('Graph 3: Clinical Depression Rate by Gender', fontsize=12)
axes[0, 1].set_ylabel('Proportion Depressed (PHQ9 >= 10)')

# Graph 4: BMI vs Depression Status (Boxplot)
sns.boxplot(data=df_clean, x='Depressed', y='BMI', palette=custom_palette, hue='Depressed', legend=False, ax=axes[1, 0])
axes[1, 0].set_title('Graph 4: Clinical BMI Distribution by Depression Status', fontsize=12)
axes[1, 0].set_xticks([0, 1])
axes[1, 0].set_xticklabels(['Not Depressed (0)', 'Depressed (1)'])
axes[1, 0].set_xlabel('Depression Status') # Added to override the x-axis label
axes[1, 0].set_ylabel('Body Mass Index (BMI)')

# Graph 5: Inflammation (CRP) vs Depression Status (Violin Plot)
# Capping CRP at 15 mg/L for visualization purposes to cut extreme outliers
df_vis = df_clean[df_clean['CRP_mgL'] <= 15]
sns.violinplot(data=df_vis, x='Depressed', y='CRP_mgL', palette=custom_palette, hue='Depressed', legend=False, inner="quartile", ax=axes[1, 1])
axes[1, 1].set_title('Graph 5: Systemic Inflammation (CRP) by Depression Status', fontsize=12)
axes[1, 1].set_xticks([0, 1])
axes[1, 1].set_xticklabels(['Not Depressed (0)', 'Depressed (1)'])
axes[1, 1].set_xlabel('Depression Status') # Added to override the x-axis label
axes[1, 1].set_ylabel('C-Reactive Protein (mg/L)')

plt.tight_layout()
plt.show()

from scipy import stats

# 1. Isolate the CRP data for the two groups
crp_not_depressed = df_clean[df_clean['Depressed'] == 0]['CRP_mgL'].dropna()
crp_depressed = df_clean[df_clean['Depressed'] == 1]['CRP_mgL'].dropna()

# 2. Run the Mann-Whitney U test
# We use 'greater' because our hypothesis is that the depressed group has HIGHER inflammation
stat, p_value = stats.mannwhitneyu(crp_depressed, crp_not_depressed, alternative='greater')

# 3. Print the results cleanly
print(f"Mann-Whitney U Statistic: {stat}")
print(f"P-Value: {p_value}")

if p_value < 0.05:
    print("\nResult: STATISTICALLY SIGNIFICANT (p < 0.05).")
    print("The depressed cohort has significantly higher systemic inflammation.")
else:
    print("\nResult: NOT STATISTICALLY SIGNIFICANT.")
    print("We cannot confidently say the inflammation levels are different.")

Mann-Whitney U Statistic: 2009218.0
P-Value: 1.1068184701508483e-06

Result: STATISTICALLY SIGNIFICANT (p < 0.05).
The depressed cohort has significantly higher systemic inflammation.

	ID	Gender	Age	BMI	Sleep_Hours	Vigorous_Activity_Min	Income_Poverty_Ratio	CRP_mgL	PHQ9_Score
0	130378	Male	43	27.0	9.5	45.0	5.00	1.78	0
1	130379	Male	66	33.5	9.0	45.0	5.00	2.03	1
2	130380	Female	44	29.7	8.0	0.0	1.41	5.62	2
3	130386	Male	34	30.2	7.5	30.0	1.33	1.05	1
4	130387	Female	68	42.6	3.0	0.0	1.32	3.96	0

Predictive Modeling of Mental Health and Wellness: Synergies of Lifestyle and Socioeconomic Factors (2021-2023 Cohort)¶

Project Plan¶

Extraction, Transform, and Load (ETL)¶

Exploratory Data Analysis (EDA)¶

Discussion of EDA Findings¶

Proposed Machine Learning Models (Milestone 2 Outline)¶

Model 1: Binary Classification of Clinical Depression (Logistic Regression / Random Forest)¶

Model 2: Continuous Prediction of Mental Wellness (Multiple Linear Regression)¶