import pandas as pd
from statsmodels.stats.weightstats import DescrStatsW
import numpy as np
# Weighted mean & SD
def w_mean_sd(series, weights):
s = DescrStatsW(series, weights=weights, ddof=0)
return s.mean, s.std
# Weighted % for categorical
def w_percent(df, var, target=None):
d = df.copy()
if target is not None:
d = d[d['IRAMDEYR']==target]
return (d.groupby(var)['ANALWT2_C'].sum() / d['ANALWT2_C'].sum() * 100).round(1)
# Weighted N (population estimate)
def weighted_n(df):
return round(df['ANALWT2_C'].sum() / 1e6, 1) # in million (approx population)
# Columns you need
cols = [
'IRAMDEYR','CATAG7','IRSEX','NEWRACE2',
'AKSSLR6WRST','WHODASDASC','UD5ILLANY','ANALWT2_C'
]
# Load dataset (tab-delimited)
df = pd.read_csv("../data/raw/NSDUH_2023_Tab.txt", sep='\t', usecols=cols, low_memory=False)
# Adults only (CATAG7 ≥ 4)
df = df[df['CATAG7'] >= 4]
# Drop missing target
df = df[df['IRAMDEYR'].notna()]
no_mde_raw = len(df[df.IRAMDEYR == 0])
mde_raw = len(df[df.IRAMDEYR == 1])
total_raw = no_mde_raw + mde_raw
# Percent share of total (based on sample)
no_mde_pct = round(no_mde_raw / total_raw * 100, 1)
mde_pct = round(mde_raw / total_raw * 100, 1)
# ---------- Build rows ----------
rows = []
# --- 1️⃣ Total participants row ---
rows.append([
"Total participants",
"Total respondents in the analytic sample",
f"{no_mde_raw:,} ({no_mde_pct:,} %)",
f"{mde_raw:,} ({mde_pct:,} %)"
])
# --- 2️⃣ Continuous variables (mean ± SD) ---
for var in ['AKSSLR6WRST','WHODASDASC']:
mean0, sd0 = w_mean_sd(df[df.IRAMDEYR==0][var], df[df.IRAMDEYR==0]['ANALWT2_C'])
mean1, sd1 = w_mean_sd(df[df.IRAMDEYR==1][var], df[df.IRAMDEYR==1]['ANALWT2_C'])
rows.append([var, "", f"{mean0:.2f} ± {sd0:.2f}", f"{mean1:.2f} ± {sd1:.2f}"])
# --- 3️⃣ Categorical variables (weighted %) ---
for var in ['CATAG7','IRSEX','NEWRACE2','UD5ILLANY']:
no_mde = w_percent(df, var, 0)
mde = w_percent(df, var, 1)
for level in sorted(no_mde.index):
rows.append([
f"{var}_{level}",
"",
f"{no_mde.get(level, np.nan)} %",
f"{mde.get(level, np.nan)} %"
])
# ---------- Convert to DataFrame ----------
table1 = pd.DataFrame(rows, columns=[
"Feature (code)",
"Code Meaning",
"No MDE (Weighted %) or Mean ± SD",
"MDE (Weighted %) or Mean ± SD"
])
# ---------- Add readable code meanings ----------
labels = {
'CATAG7_4':'Age 18–25 y','CATAG7_5':'Age 26–34 y',
'CATAG7_6':'Age 35–49 y','CATAG7_7':'Age 50+ y',
'IRSEX_1':'Male','IRSEX_2':'Female',
'NEWRACE2_1':'White','NEWRACE2_2':'Black or African American',
'NEWRACE2_3':'American Indian / Alaska Native',
'NEWRACE2_4':'Native Hawaiian / Pacific Islander',
'NEWRACE2_5':'Asian','NEWRACE2_6':'Two or more races','NEWRACE2_7':'Hispanic or Latino',
'UD5ILLANY_0':'No Drug Use Disorder','UD5ILLANY_1':'Has Drug Use Disorder',
'AKSSLR6WRST':'K6 Psychological Distress Score (0–24)',
'WHODASDASC':'WHODAS Disability Score (0–100)'
}
table1["Code Meaning"] = table1["Feature (code)"].map(labels).fillna(table1["Code Meaning"])
# ---------- Export clean final table ----------
table1.to_csv("../results/table1_Baseline_Characteristics.csv", index=False)
print("✅ Clean Table 1 saved as 'table1_Baseline_Characteristics.csv'")
print(table1.head(10))
✅ Clean Table 1 saved as 'table1_Baseline_Characteristics.csv'
Feature (code) Code Meaning \
0 Total participants Total respondents in the analytic sample
1 AKSSLR6WRST K6 Psychological Distress Score (0–24)
2 WHODASDASC WHODAS Disability Score (0–100)
3 CATAG7_4 Age 18–25 y
4 CATAG7_5 Age 26–34 y
5 CATAG7_6 Age 35–49 y
6 CATAG7_7 Age 50+ y
7 IRSEX_1 Male
8 IRSEX_2 Female
9 NEWRACE2_1 White
No MDE (Weighted %) or Mean ± SD MDE (Weighted %) or Mean ± SD
0 39,948 (88.5 %) 5,185 (11.5 %)
1 1.19 ± 2.96 9.01 ± 5.20
2 0.79 ± 1.78 4.93 ± 2.66
3 4.3 % 8.4 %
4 7.6 % 18.3 %
5 14.8 % 24.4 %
6 73.2 % 48.9 %
7 49.9 % 35.4 %
8 50.1 % 64.6 %
9 60.7 % 66.0 %