Data Loading & Initial Cleaning


# Import Libraries
import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.impute import SimpleImputer
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score, f1_score
import matplotlib.pyplot as plt
import seaborn as sns

# Import NSUDH dataset
df = pd.read_csv("../data/NSDUH_2023_Tab.txt", 
                 sep="\t", 
                 low_memory=False)

# Checking Data
df.head

<bound method NDFrame.head of        QUESTID2    FILEDATE     ANALWT2_C  VESTR_C  VEREP  PDEN10  COUTYP4  \
0      10000053  03/25/2025   3276.469874    40031      2       2        2   
1      10000679  03/25/2025  15630.082955    40021      2       2        3   
2      10001208  03/25/2025   4018.172390    40043      1       2        2   
3      10001260  03/25/2025  10711.709540    40030      2       2        2   
4      10001588  03/25/2025   8195.104779    40023      2       2        2   
...         ...         ...           ...      ...    ...     ...      ...   
56700  50556120  03/25/2025    417.630119    40018      2       2        2   
56701  50557151  03/25/2025   7625.717934    40017      1       2        2   
56702  50558694  03/25/2025  23556.083908    40042      1       1        1   
56703  50558696  03/25/2025   5193.882625    40040      2       1        1   
56704  50558785  03/25/2025   2676.234296    40027      1       2        2   

       MAIIN102  AIIND102  AGE3  ...  COSUTELE2  COSUAPTDL2  COSURXDL2  \
0             2         2    10  ...        3.0         3.0        3.0   
1             2         2     9  ...        3.0         3.0        3.0   
2             2         2     9  ...        3.0         3.0        3.0   
3             2         2     1  ...        2.0         2.0        2.0   
4             2         2    10  ...        3.0         3.0        3.0   
...         ...       ...   ...  ...        ...         ...        ...   
56700         2         2     9  ...        NaN         NaN        NaN   
56701         2         2    11  ...        3.0         3.0        3.0   
56702         2         2    10  ...        2.0         2.0        2.0   
56703         2         2    11  ...        3.0         3.0        3.0   
56704         2         2     3  ...        2.0         2.0        2.0   

       COSUSVHLT2  COHCTELE2  COHCAPTDL2  COHCRXDL2  COHCSVHLT2  LANGVER  \
0             3.0        3.0         3.0        3.0         3.0        1   
1             3.0        3.0         3.0        3.0         3.0        1   
2             3.0        3.0         3.0        3.0         3.0        1   
3             2.0        2.0         2.0        2.0         2.0        1   
4             3.0        1.0         3.0        3.0         3.0        1   
...           ...        ...         ...        ...         ...      ...   
56700         NaN        NaN         NaN        NaN         NaN        2   
56701         3.0        2.0         2.0        2.0         2.0        1   
56702         2.0        2.0         2.0        2.0         2.0        1   
56703         3.0        2.0         2.0        2.0         2.0        1   
56704         2.0        2.0         2.0        2.0         2.0        1   

       GQTYPE2  
0          NaN  
1          NaN  
2          NaN  
3          NaN  
4          NaN  
...        ...  
56700      NaN  
56701      NaN  
56702      NaN  
56703      NaN  
56704      NaN  

[56705 rows x 2636 columns]>

missing_summary = df.isna().mean() * 100
missing_summary.sort_values(ascending=False)

SRCCLFRSED     99.871264
SRCFRSEDNM     99.871264
GQTYPE2        99.850101
SRCSEDNM2      99.728419
SRCFRTRQNM     99.370426
                 ...    
SEXIDENT22      0.000000
SPEAKENGL       0.000000
LVLDIFSEE2      0.000000
LVLDIFHEAR2     0.000000
LVLDIFWALK2     0.000000
Length: 2636, dtype: float64

# Checking for duplicate rows
df.duplicated().sum()

np.int64(0)

# Filtering dataset for relevant age group (18 and above)
'''I used AGE3 instead of CATAG7 because AGE3 is the final edited age variable in NSDUH. It’s the most accurate age measure because it incorporates multiple consistency checks—birthdate, roster age, screener age, and internal corrections. CATAG7 is just a categorical recode derived from AGE3. For filtering out respondents under 18, it’s better to filter at the source (AGE3) and let the model work with the true final age before any recoding.

AGE3 Len : 2 RECODE - FINAL EDITED AGE
Freq Pct
1 = Respondent is 12 or 13 years old
2 = Respondent is 14 or 15 years old
3 = Respondent is 16 or 17 years old
4 = Respondent is between 18 and 20 years old
5 = Respondent is between 21 and 23 years old
6 = Respondent is 24 or 25 years old
7 = Respondent is between 26 and 29 years old
8 = Respondent is between 30 and 34 years old
9 = Respondent is between 35 and 49 years old
10 = Respondent is between 50 and 64 years old
11 = Respondent is 65 years old or older'''

df_no_age = df[df['AGE3'] >= 4]   
df_no_age.shape

(45133, 2636)

# Checking how many error codes are present in the dataset
'''Code	Meaning	Should Convert to NaN?
94 / 994 / 9994	Don't know
97 / 997 / 9997	Refused
98 / 998 / 9998	Blank / Not answered
85 / 985 / 9985	Bad / inconsistent data
99 / 999 / 9999	Legitimate skip	
89 / 989 / 9989	Legitimate skip - logically assigned'''

nan_codes = [
    94, 97, 98, 85,
    994, 997, 998, 985,
    9994, 9997, 9998, 9985,
    99, 999, 9999,
    89, 989, 9989
]

for code in nan_codes:
    count = (df_no_age == code).sum().sum()
    print(code, count)

94 72183
97 44145
98 1898217
85 9487
994 7015
997 3857
998 33083
985 5817
9994 4675
9997 1899
9998 20242
9985 1774
99 17222835
999 1267661
9999 1453080
89 4185
989 154
9989 1332

# Converting the following error codes to NaN

nan_codes = [
    94, 97, 98, 85,
    994, 997, 998, 985,
    9994, 9997, 9998, 9985,
    99, 999, 9999,
    89, 989, 9989
]

df_ErrorCode_NaN = df_no_age.replace(nan_codes, np.nan)

nan_codes = [
    94, 97, 98, 85,
    994, 997, 998, 985,
    9994, 9997, 9998, 9985,
    99, 999, 9999,
    89, 989, 9989
]

for code in nan_codes:
    count = (df_ErrorCode_NaN == code).sum().sum()
    print(code, count)

# Checking for any NaNs in the Target Variable (IRAMDEYR)
df_ErrorCode_NaN['IRAMDEYR'].isna().sum()

np.int64(0)

# 2. DROP rows where the target is NaN
df_dropped_targetNaN = df_ErrorCode_NaN.dropna(subset=['IRAMDEYR'])

before = df.shape[0]
after = df_dropped_targetNaN.shape[0]

print("Rows before:", before)
print("Rows after :", after)
print("Rows removed:", before - after)

Rows before: 56705
Rows after : 45133
Rows removed: 11572

df_dropped_targetNaN.describe()

	QUESTID2	ANALWT2_C	VESTR_C	VEREP	PDEN10	COUTYP4	MAIIN102	AIIND102	AGE3	SERVICE	...	COMHSVHLT2	COSUTELE2	COSUAPTDL2	COSURXDL2	COSUSVHLT2	COHCTELE2	COHCAPTDL2	COHCRXDL2	COHCSVHLT2	LANGVER
count	4.513300e+04	45133.000000	45133.000000	45133.000000	45133.000000	45133.000000	45133.000000	45133.000000	45133.000000	45108.000000	...	43333.000000	43354.000000	43364.000000	43362.000000	43347.000000	43378.000000	43375.000000	43380.000000	43347.00000	45133.000000
mean	3.027962e+07	5706.297779	40025.494849	1.429774	1.628631	1.713602	1.984579	1.984357	7.816697	1.948967	...	2.410772	2.533884	2.537450	2.563004	2.578518	1.968809	2.036749	2.196358	2.26424	1.042519
std	1.172008e+07	8511.336263	14.499402	0.495049	0.586271	0.730958	0.123222	0.124090	2.194038	0.220068	...	0.592895	0.618880	0.603669	0.563047	0.542894	0.761610	0.720678	0.611436	0.56382	0.201772
min	1.000005e+07	1.516955	40001.000000	1.000000	1.000000	1.000000	1.000000	1.000000	4.000000	1.000000	...	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.00000	1.000000
25%	2.013831e+07	973.740425	40013.000000	1.000000	1.000000	1.000000	2.000000	2.000000	6.000000	2.000000	...	2.000000	2.000000	2.000000	2.000000	2.000000	1.000000	2.000000	2.000000	2.00000	1.000000
50%	3.027958e+07	2697.651023	40025.000000	1.000000	2.000000	2.000000	2.000000	2.000000	8.000000	2.000000	...	2.000000	3.000000	3.000000	3.000000	3.000000	2.000000	2.000000	2.000000	2.00000	1.000000
75%	4.039094e+07	6788.874141	40038.000000	2.000000	2.000000	2.000000	2.000000	2.000000	9.000000	2.000000	...	3.000000	3.000000	3.000000	3.000000	3.000000	3.000000	3.000000	3.000000	3.00000	1.000000
max	5.055870e+07	118941.430150	40050.000000	2.000000	3.000000	3.000000	2.000000	2.000000	11.000000	2.000000	...	3.000000	3.000000	3.000000	3.000000	3.000000	3.000000	3.000000	3.000000	3.00000	2.000000

8 rows × 2634 columns

df_dropped_targetNaN.isna().sum().sort_values(ascending=False)

YOLOSEV         45133
YODSCEV         45133
YODPREV         45133
YOWRCHR         45133
YOWRDST         45133
                ...  
FLVVAPMON           0
FLVVAPYR            0
IREDUHIGHST2        0
IIMARIT             0
MAIIN102            0
Length: 2636, dtype: int64

# Drop ID, weight, design, and metadata variables
cols_to_drop = []

# Unique identifiers
cols_to_drop += [c for c in df_dropped_targetNaN.columns if 'CASEID' in c or 'QUESTID' in c]

# Weight variables
cols_to_drop += [c for c in df_dropped_targetNaN.columns if 'WT' in c or 'WGT' in c or 'ANALWT' in c]

# Sampling design variables
cols_to_drop += [c for c in df_dropped_targetNaN.columns if 'VESTR' in c or 'VEREP' in c or 'ESTRAT' in c]

# Interview month/day/time variables
cols_to_drop += [c for c in df_dropped_targetNaN.columns if c.startswith('INTV_')]

# Meta / administration variables
cols_to_drop += [c for c in df_dropped_targetNaN.columns if c.endswith('_A') or c.endswith('_E')]

# Year and file metadata
cols_to_drop += [c for c in ['YEAR', 'FILEDATE', 'QUARTER'] if c in df_dropped_targetNaN.columns]

# Remove duplicates in the list
cols_to_drop = list(set(cols_to_drop))

len(cols_to_drop), cols_to_drop[:20]

(7,
 ['WTPOUND2',
  'ANALWT2_C',
  'FILEDATE',
  'VESTR_C',
  'QUESTID2',
  'WTANSWER',
  'VEREP'])

drop_patterns = [
    r'^CASEID', r'^QUESTID', r'^PANEL', r'^VERSION', 
    r'^INTV', r'^FILE', r'^YEAR', 
    r'^WT', r'WGT', r'WEIGHT', r'ANALWT', 
    r'^VESTR', r'^VEREP',
    r'_A$', r'_E$', r'_ORIG$', r'_R$', r'_RC$'
]

import re

cols_to_drop = [
    c for c in df_dropped_targetNaN.columns
    if any(re.search(pat, c) for pat in drop_patterns)
]

len(cols_to_drop), cols_to_drop[:30]

(7,
 ['QUESTID2',
  'FILEDATE',
  'ANALWT2_C',
  'VESTR_C',
  'VEREP',
  'WTANSWER',
  'WTPOUND2'])

df_Removed_Columns = df_dropped_targetNaN.drop(columns=cols_to_drop, errors='ignore')

drop_patterns = [
    r'^CASEID', r'^QUESTID', r'^PANEL', r'^VERSION', 
    r'^INTV', r'^FILE', r'^YEAR', 
    r'^WT', r'WGT', r'WEIGHT', r'ANALWT', 
    r'^VESTR', r'^VEREP',
    r'_A$', r'_E$', r'_ORIG$', r'_R$', r'_RC$'
]

import re

cols_to_drop2 = [
    c for c in df_Removed_Columns.columns
    if any(re.search(pat, c) for pat in drop_patterns)
]

len(cols_to_drop2), cols_to_drop2[:30]

(0, [])

# Checking to see which variables have >99% missingess
na_percent = df_Removed_Columns.isna().mean() * 100

# Filter columns with > 99% missing
cols_over_99 = na_percent[na_percent > 99]

# Print count
print("Number of variables with >99% missingness:", len(cols_over_99))

# Print the variable names
print("\nVariables with >99% missingness:")
print(cols_over_99.index.tolist())

# Optional: show top 20 instead of the whole list
print("\nPreview (first 20):")
print(cols_over_99.head(20))

Number of variables with >99% missingness: 282

Variables with >99% missingness:
['EDUSCKEST', 'EDUSKPEST', 'HLCALLFG', 'HLCALL99', 'BKPOSTOB', 'BKOTHOF2', 'ALTOTFG', 'ALFQFLG', 'ALDYSFG', 'MJFQFLG', 'BLRECFL2', 'BLNT30C1', 'BLNT30C2', 'RSNOMRJ', 'RSNMRJMO', 'CCTOTFG', 'CCFQFLG', 'CRTOTFG', 'CRFQFLG', 'HRTOTFG', 'HRFQFLG', 'HALTOTFG', 'HALFQFLG', 'INHTOTFG', 'INHFQFLG', 'METOTFG', 'MEFQFLG', 'PNRNORXFG', 'TRQNORXFG', 'STMNORXFG', 'SEDNORXFG', 'SRCSEDNM2', 'SRCFRPNRNM', 'SRCFRTRQNM', 'SRCFRSEDNM', 'SRCCLFRPNR', 'SRCCLFRTRQ', 'SRCCLFRSED', 'OTCFLAG', 'OTDGNDLA', 'OTDGNDLB', 'OTDGNDLC', 'OTDGNDLD', 'OTDGNDLE', 'SUNTINSCV', 'SUNTENCV', 'SUNTNOOPN', 'YEATNDYR', 'YEHMSLYR', 'YESCHFLT', 'YESCHWRK', 'YESCHIMP', 'YESCHINT', 'YETCGJOB', 'YELSTGRD', 'YECIGFRNDOF2', 'YECIGNEXTYR2', 'YESTSCIG', 'YESTSMJ', 'YESTSALC', 'YESTSDNK', 'YEPCHKHW', 'YEPHLPHW', 'YEPCHORE', 'YEPLMTTV', 'YEPLMTSN', 'YEPGDJOB', 'YEPPROUD', 'YEYARGUP', 'YEYFGTSW', 'YEYFGTGP', 'YEYHGUN', 'YEYSELL', 'YEYSTOLE', 'YEYATTAK', 'YEPPKCIG', 'YEPMJEVR', 'YEPMJMO', 'YEPALDLY', 'YEGPKCIG', 'YEGMJEVR', 'YEGMJMO', 'YEGALDLY', 'YEFPKCIG', 'YEFMJEVR', 'YEFMJMO', 'YEFALDLY', 'YETLKNON', 'YETLKPAR', 'YETLKBGF', 'YETLKOTA', 'YETLKSOP', 'YEPRTDNG', 'YEPRBSLV', 'YEVIOPRV', 'YEDGPRGP', 'YESLFHLP', 'YEPRGSTD', 'YESCHACT', 'YECOMACT', 'YEFAIACT', 'YEOTHACT', 'YEDECLAS', 'YEDERGLR', 'YEDESPCL', 'YEPVNTYR', 'YERLGSVC', 'YERLGIMP', 'YERLDCSN', 'YERLFRND', 'YUSUITHK', 'YUCOSUITHK', 'YUSUIPLN', 'YUCOSUIPLN', 'SCHFELT', 'TCHGJOB', 'AVGGRADE', 'STNDSCIG', 'STNDSMJ', 'STNDALC', 'STNDDNK', 'PARCHKHW', 'PARHLPHW', 'PRCHORE2', 'PRLMTTV2', 'PARLMTSN', 'PRGDJOB2', 'PRPROUD2', 'ARGUPAR', 'YOFIGHT2', 'YOGRPFT2', 'YOHGUN2', 'YOSELL2', 'YOSTOLE2', 'YOATTAK2', 'PRPKCIG2', 'PRMJEVR2', 'PRMJMO', 'PRALDLY2', 'YFLPKCG2', 'YFLTMRJ2', 'YFLMJMO', 'YFLADLY2', 'FRDPCIG2', 'FRDMEVR2', 'FRDMJMON', 'FRDADLY2', 'TALKPROB', 'PRTALK3', 'PRBSOLV2', 'PREVIOL2', 'PRVDRGO2', 'GRPCNSL2', 'PREGPGM2', 'YTHACT2', 'DRPRVME3', 'ANYEDUC3', 'RLGATTD', 'RLGIMPT', 'RLGDCSN', 'RLGFRND', 'YUSUITHKYR', 'YUCOSUITHK2', 'YUSUIPLNYR', 'YUCOSUIPLN2', 'YODPREV', 'YODSCEV', 'YOLOSEV', 'YODPDISC', 'YODPLSIN', 'YODSLSIN', 'YOLSI2WK', 'YODPR2WK', 'YOWRHRS', 'YOWRDST', 'YOWRCHR', 'YOWRIMP', 'YODPPROB', 'YOWRPROB', 'YOWRAGE', 'YOWRDEPR', 'YOWRDISC', 'YOWRLSIN', 'YOWRPLSR', 'YOWRELES', 'YOWREMOR', 'YOWRGAIN', 'YOWRGROW', 'YOWRPREG', 'YOWRGNL2', 'YOWRLOSE', 'YOWRDIET', 'YOWRLSL2', 'YOWRSLEP', 'YOWRSMOR', 'YOWRENRG', 'YOWRSLOW', 'YOWRSLNO', 'YOWRJITT', 'YOWRJINO', 'YOWRTHOT', 'YOWRCONC', 'YOWRDCSN', 'YOWRNOGD', 'YOWRWRTH', 'YO_MDEA1', 'YO_MDEA2', 'YO_MDEA3', 'YO_MDEA4', 'YO_MDEA5', 'YO_MDEA6', 'YO_MDEA7', 'YO_MDEA8', 'YODSMMDE', 'YOPBINTF', 'YOPBDLYA', 'YOPBRMBR', 'YOPBAGE', 'YOPBNUM', 'YOPB2WK', 'YOPSHMGT', 'YOPSWORK', 'YOPSRELS', 'YOPSSOC', 'YOPSDAYS', 'YOSEEDOC', 'YOFAMDOC', 'YOOTHDOC', 'YOPSYCH', 'YOPSYMD', 'YOSOCWRK', 'YOCOUNS', 'YOOTHMHP', 'YONURSE', 'YORELIG', 'YOHERBAL', 'YOOTHHLP', 'YOTMTNOW', 'YORX12MO', 'YORXNOW', 'YORXHLP', 'YOTMTHLP', 'YMDELT', 'YMDEYR', 'YMDEAUD5YR', 'YMIUD5YANY', 'YMSUD5YANY', 'YMDERSUD5ANY', 'YMDESUD5ANYO', 'YTXMDEYR', 'YRXMDEYR', 'YMDETXRX', 'YDOCMDE', 'YOMDMDE', 'YPSY1MDE', 'YPSY2MDE', 'YSOCMDE', 'YCOUNMDE', 'YOMHMDE', 'YNURSMDE', 'YRELMDE', 'YHBCHMDE', 'YHLTMDE', 'YALTMDE', 'YMDEHPRX', 'YMDEHPO', 'YMDERXO2', 'YMDEHARX', 'YSDSHOME', 'YSDSWRK', 'YSDSREL', 'YSDSSOC', 'YSDSOVRL', 'MDEIMPY', 'YMDEIMAD5YR', 'YMIMI5YANY', 'YMIMR5YANY', 'YMIMS5YANY', 'CIRROSAGE', 'HEPBCAGE', 'HIVAIDSAG', 'GQTYPE2']

Preview (first 20):
EDUSCKEST     99.915804
EDUSKPEST     99.960118
HLCALLFG      99.953471
HLCALL99      99.953471
BKPOSTOB     100.000000
BKOTHOF2      99.519199
ALTOTFG       99.339729
ALFQFLG       99.295416
ALDYSFG       99.960118
MJFQFLG       99.550218
BLRECFL2      99.915804
BLNT30C1      99.884785
BLNT30C2      99.802805
RSNOMRJ       99.922451
RSNMRJMO      99.898079
CCTOTFG       99.995569
CCFQFLG       99.953471
CRTOTFG       99.995569
CRFQFLG       99.991137
HRTOTFG       99.986706
dtype: float64

# Look at how often these variables are non-missing among people with IRAMDEYR=1
df_99 = df_Removed_Columns.drop(columns=cols_over_99.index)
df_99.shape

(45133, 2347)

Split X and Y

'''Missing Values (NaNs):
-Many variables had missing values due to survey skip patterns or special response codes.
-For numeric variables, missing values were replaced with the median of the column. The median is robust to outliers and preserves the typical value without disproportionately influencing the results.
-For categorical variables, missing values were replaced with the mode (most frequent category), which preserves the dominant response and minimizes distortion.
-This ensures that all columns are usable for tree-based algorithms without introducing extreme bias, especially important for feature importance ranking.

Encoding Categorical Variables:
-Tree-based models like RandomForest cannot directly process string categories.
-We converted categorical variables to numeric codes using LabelEncoder, which assigns each category a unique integer. Importantly, RandomForest treats these codes as discrete categories, not as ordered values, so this transformation does not bias splits or feature importance.

Rationale:
-These steps allow us to use the full set of variables in a first-pass feature importance analysis with RandomForest.
-Alternative models like LightGBM or CatBoost can handle missing values natively, but median/mode imputation provides a simple, consistent approach for exploratory feature ranking.
-No other cleaning (e.g., dropping additional columns) is required for this stage, as we have already removed identifiers, survey weights, timestamps, and variables with extremely high missingness.'''


from sklearn.preprocessing import LabelEncoder

# Separate target
TARGET = "IRAMDEYR"
y = df_99[TARGET]
X = df_99.drop(columns=[TARGET])


# Identify numeric vs categorical
numeric_cols = X.select_dtypes(include=[np.number]).columns.tolist()
categorical_cols = X.select_dtypes(exclude=[np.number]).columns.tolist()

print(f"Numeric columns: {len(numeric_cols)}")
print(f"Categorical columns: {len(categorical_cols)}")

Numeric columns: 2346
Categorical columns: 0

# Impute missing values

# Numeric -> median
for col in numeric_cols:
    median_val = X[col].median()
    X[col] = X[col].fillna(median_val)

# Categorical -> mode
for col in categorical_cols:
    mode_val = X[col].mode()[0]
    X[col] = X[col].fillna(mode_val)

# Check if any NaNs remain
print("Any NaNs left in the dataset?", X.isna().any().any())

# Count of NaNs per column (should all be 0)
nan_counts = X.isna().sum()
print(nan_counts[nan_counts > 0])  # this will print nothing if there are no NaNs

total_nans = X.isna().sum().sum()
print("Total NaNs in the dataset:", total_nans)

Any NaNs left in the dataset? False
Series([], dtype: int64)
Total NaNs in the dataset: 0

# Encode categorical variables
le_dict = {}  # save encoders in case needed later
for col in categorical_cols:
    le = LabelEncoder()
    X[col] = le.fit_transform(X[col])
    le_dict[col] = le

# Combine X and y for later use
df_prepared = pd.concat([X, y], axis=1)

# Check if any NaNs remain
print("Any NaNs left in the dataset?", y.isna().any().any())

# Count of NaNs per column (should all be 0)
nan_counts = y.isna().sum()
print(nan_counts[nan_counts > 0])  # this will print nothing if there are no NaNs

total_nans = y.isna().sum().sum()
print("Total NaNs in the dataset:", total_nans)

Any NaNs left in the dataset? False
[]
Total NaNs in the dataset: 0

Feature Importance / Tiering

Correlation only works on linear continuous relationships. Our dataset is mostly categorical and ordinal, so Pearson correlation wouldn’t capture any useful signal. That’s why the papers and frameworks use model-based feature selection, like LASSO, ElasticNet, RFE, and Random Forest. Those methods work regardless of whether the relationship is linear or categorical.

# Removing target leakage features that was used to impute / recode the target IRAMDEYR
# IRAMDEYR → AMDEYR → AMDELT → ADSMMDEA → MDE symptom items + fallback items + timing

leakage_vars = [

    # Target + imputation
    "IRAMDEYR", "IIAMDEYR",

    # Pre-imputation past-year MDE
    "AMDEYR",

    # Lifetime MDE + imputation
    "AMDELT", "IRAMDELT", "IIAMDELT",

    # MDE with impairment (logically dependent on AMDEYR)
    "AMDEIMP", "IRAMDEIMP", "IIAMDEIMP",

    # DSM summary
    "ADSMMDEA",

    # Raw DSM symptom items (if present)
    "D_MDEA1", "D_MDEA2", "D_MDEA3", "D_MDEA4", "D_MDEA5",
    "D_MDEA6", "D_MDEA7", "D_MDEA8", "D_MDEA9",

    # Symptom recodes (one-digit suffixes, if present)
    "AD_MDEA1", "AD_MDEA2", "AD_MDEA3", "AD_MDEA4",
    "AD_MDEA5", "AD_MDEA6", "AD_MDEA7", "AD_MDEA8",

    # Symptom recodes (true DSM flags: two-digit suffix)
    "AD_MDEA11", "AD_MDEA21", "AD_MDEA31", "AD_MDEA41",
    "AD_MDEA51", "AD_MDEA61", "AD_MDEA71", "AD_MDEA81", "AD_MDEA91",

    # Fallback/gating items
    "ADPB2WK",
    "ADDPREV", "ADDSCEV", "ADLOSEV", "ADLSI2WK", "ADDPR2WK",
    "ADWRHRS", "ADWRDST", "ADWRCHR", "ADWRIMP", "ADDPPROB",

    # Soft-leakage variables
    "ARXMDEYR",   # received prescription meds
    "ATXMDEYR",   # received counseling/therapy
    "AHLTMDE",     # told by provider they have depression
]

print(leakage_vars)
len(leakage_vars)

['IRAMDEYR', 'IIAMDEYR', 'AMDEYR', 'AMDELT', 'IRAMDELT', 'IIAMDELT', 'AMDEIMP', 'IRAMDEIMP', 'IIAMDEIMP', 'ADSMMDEA', 'D_MDEA1', 'D_MDEA2', 'D_MDEA3', 'D_MDEA4', 'D_MDEA5', 'D_MDEA6', 'D_MDEA7', 'D_MDEA8', 'D_MDEA9', 'AD_MDEA1', 'AD_MDEA2', 'AD_MDEA3', 'AD_MDEA4', 'AD_MDEA5', 'AD_MDEA6', 'AD_MDEA7', 'AD_MDEA8', 'AD_MDEA11', 'AD_MDEA21', 'AD_MDEA31', 'AD_MDEA41', 'AD_MDEA51', 'AD_MDEA61', 'AD_MDEA71', 'AD_MDEA81', 'AD_MDEA91', 'ADPB2WK', 'ADDPREV', 'ADDSCEV', 'ADLOSEV', 'ADLSI2WK', 'ADDPR2WK', 'ADWRHRS', 'ADWRDST', 'ADWRCHR', 'ADWRIMP', 'ADDPPROB', 'ARXMDEYR', 'ATXMDEYR', 'AHLTMDE']





50

X_no_leakage = X.drop(columns=[col for col in leakage_vars if col != 'IRAMDEYR'], errors='ignore')
X_no_leakage.shape

(45133, 2315)

# Initialize the model
from sklearn.ensemble import RandomForestClassifier

rf_model_no_leakage = RandomForestClassifier(
    n_estimators=500,
    max_depth=None,
    random_state=42,
    n_jobs=-1
)

rf_model_no_leakage.fit(X_no_leakage, y)

RandomForestClassifier(n_estimators=500, n_jobs=-1, random_state=42)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

RandomForestClassifier

?Documentation for RandomForestClassifieriFitted

    <div class="estimator-table">
        <details>
            <summary>Parameters</summary>
            <table class="parameters-table">
              <tbody>

    <tr class="user-set">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('n_estimators',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">n_estimators&nbsp;</td>
        <td class="value">500</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('criterion',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">criterion&nbsp;</td>
        <td class="value">&#x27;gini&#x27;</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('max_depth',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">max_depth&nbsp;</td>
        <td class="value">None</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('min_samples_split',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">min_samples_split&nbsp;</td>
        <td class="value">2</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('min_samples_leaf',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">min_samples_leaf&nbsp;</td>
        <td class="value">1</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('min_weight_fraction_leaf',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">min_weight_fraction_leaf&nbsp;</td>
        <td class="value">0.0</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('max_features',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">max_features&nbsp;</td>
        <td class="value">&#x27;sqrt&#x27;</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('max_leaf_nodes',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">max_leaf_nodes&nbsp;</td>
        <td class="value">None</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('min_impurity_decrease',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">min_impurity_decrease&nbsp;</td>
        <td class="value">0.0</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('bootstrap',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">bootstrap&nbsp;</td>
        <td class="value">True</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('oob_score',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">oob_score&nbsp;</td>
        <td class="value">False</td>
    </tr>


    <tr class="user-set">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('n_jobs',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">n_jobs&nbsp;</td>
        <td class="value">-1</td>
    </tr>


    <tr class="user-set">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('random_state',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">random_state&nbsp;</td>
        <td class="value">42</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('verbose',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">verbose&nbsp;</td>
        <td class="value">0</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('warm_start',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">warm_start&nbsp;</td>
        <td class="value">False</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('class_weight',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">class_weight&nbsp;</td>
        <td class="value">None</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('ccp_alpha',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">ccp_alpha&nbsp;</td>
        <td class="value">0.0</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('max_samples',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">max_samples&nbsp;</td>
        <td class="value">None</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('monotonic_cst',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">monotonic_cst&nbsp;</td>
        <td class="value">None</td>
    </tr>

              </tbody>
            </table>
        </details>
    </div>
</div></div></div></div></div><script>function copyToClipboard(text, element) {
// Get the parameter prefix from the closest toggleable content
const toggleableContent = element.closest('.sk-toggleable__content');
const paramPrefix = toggleableContent ? toggleableContent.dataset.paramPrefix : '';
const fullParamName = paramPrefix ? `${paramPrefix}${text}` : text;

const originalStyle = element.style;
const computedStyle = window.getComputedStyle(element);
const originalWidth = computedStyle.width;
const originalHTML = element.innerHTML.replace('Copied!', '');

navigator.clipboard.writeText(fullParamName)
    .then(() => {
        element.style.width = originalWidth;
        element.style.color = 'green';
        element.innerHTML = "Copied!";

        setTimeout(() => {
            element.innerHTML = originalHTML;
            element.style = originalStyle;
        }, 2000);
    })
    .catch(err => {
        console.error('Failed to copy:', err);
        element.style.color = 'red';
        element.innerHTML = "Failed!";
        setTimeout(() => {
            element.innerHTML = originalHTML;
            element.style = originalStyle;
        }, 2000);
    });
return false;

}

document.querySelectorAll('.fa-regular.fa-copy').forEach(function(element) { const toggleableContent = element.closest('.sk-toggleable__content'); const paramPrefix = toggleableContent ? toggleableContent.dataset.paramPrefix : ''; const paramName = element.parentElement.nextElementSibling.textContent.trim(); const fullParamName = paramPrefix ? ${paramPrefix}${paramName} : paramName;

element.setAttribute('title', fullParamName);

});

EVERYTHING BELOW THIS IS WIP, AND SHOULD NOT BE EDITED

Tiering with Tier 0 + Tiers 1–4

'''Tier 0 serves a different role from Tiers 1-4.
Tier 0 is a master candidate list for Phase 2 (clinical review), so a simple “top 100 ranked features” is the right method.

Tiers 1-4 are performance-driven tiers that reflect different cumulative contributions of importance. They should remain based on percentage thresholds, because that structure tests how predictive performance scales with increasingly informative subsets.

So Tier 0 = fixed-K;
Tiers 1-4 = cumulative importance thresholds.'''

# Get feature importances from your Random Forest model (already leakage-clean)
importances = pd.Series(
    rf_model_no_leakage.feature_importances_,
    index=X_no_leakage.columns
).sort_values(ascending=False)

# Tier 0: Top-100 features
K = 100
tier0 = importances.head(K).index.tolist()

# Tiers 1–4 (cumulative thresholds)
cumulative_importance = importances.cumsum()

tier1 = importances[cumulative_importance <= 0.25].index.tolist()
tier2 = importances[cumulative_importance <= 0.45].index.tolist()
tier3 = importances[cumulative_importance <= 0.60].index.tolist()
tier4 = importances[cumulative_importance <= 0.70].index.tolist()

# Combine into dictionary
tiers = {
    'Tier 0 (Top 100)': tier0,
    'Tier 1': tier1,
    'Tier 2': tier2,
    'Tier 3': tier3,
    'Tier 4': tier4
}

# Print tier sizes
print("Feature counts per tier:")
for tier_name, features in tiers.items():
    print(f"{tier_name}: {len(features)} features")

# Plot feature importances for each tier
for tier_name, features in tiers.items():
    if len(features) == 0:
        continue
        
    plt.figure(figsize=(10, 12))
    sns.barplot(
        x=importances[features].values,
        y=importances[features].index,
        palette='viridis'
    )
    plt.title(f"{tier_name} Feature Importances ({len(features)} features)")
    plt.xlabel("Importance")
    plt.ylabel("Feature")
    plt.tight_layout()
    plt.show()

Feature counts per tier:
Tier 0 (Top 100): 100 features
Tier 1: 5 features
Tier 2: 12 features
Tier 3: 29 features
Tier 4: 52 features


C:\Users\agila\AppData\Local\Temp\ipykernel_16380\2494543076.py:47: FutureWarning: 

Passing `palette` without assigning `hue` is deprecated and will be removed in v0.14.0. Assign the `y` variable to `hue` and set `legend=False` for the same effect.

  sns.barplot(

png

C:\Users\agila\AppData\Local\Temp\ipykernel_16380\2494543076.py:47: FutureWarning: 

Passing `palette` without assigning `hue` is deprecated and will be removed in v0.14.0. Assign the `y` variable to `hue` and set `legend=False` for the same effect.

  sns.barplot(

png

C:\Users\agila\AppData\Local\Temp\ipykernel_16380\2494543076.py:47: FutureWarning: 

Passing `palette` without assigning `hue` is deprecated and will be removed in v0.14.0. Assign the `y` variable to `hue` and set `legend=False` for the same effect.

  sns.barplot(

png

C:\Users\agila\AppData\Local\Temp\ipykernel_16380\2494543076.py:47: FutureWarning: 

Passing `palette` without assigning `hue` is deprecated and will be removed in v0.14.0. Assign the `y` variable to `hue` and set `legend=False` for the same effect.

  sns.barplot(

png

C:\Users\agila\AppData\Local\Temp\ipykernel_16380\2494543076.py:47: FutureWarning: 

Passing `palette` without assigning `hue` is deprecated and will be removed in v0.14.0. Assign the `y` variable to `hue` and set `legend=False` for the same effect.

  sns.barplot(

png

Logistic Regression Baseline (With Weights + Metrics)

Re-Cleaning Dataset for Logistics Regression

df_dropped_targetNaN.shape

(45133, 2636)

drop_patterns_with_weight = [
    r'^CASEID', r'^QUESTID', r'^PANEL', r'^VERSION',
    r'^INTV', r'^FILE', r'^YEAR',
    r'^WT(?!ANSWER$)(?!POUND.*$)',  # Drop WT* except WTANSWER/WTPOUND (we drop separately)
    r'WGT', r'WEIGHT',
    r'^VESTR', r'^VEREP', r'^ESTRAT',
    r'_A$', r'_E$', r'_ORIG$', r'_R$', r'_RC$'
]

# Identify columns to drop without weight variable
cols_to_drop_with_weight = [
    c for c in df_dropped_targetNaN.columns
    if any(re.search(pat, c) for pat in drop_patterns_with_weight)
]

# Never drop ANALWT2_C
cols_to_drop_with_weight = [c for c in cols_to_drop if c != "ANALWT2_C"]

# Also manually drop these (we know they must go)
manual_dropping = [
    "WTANSWER", "WTPOUND2",  # useless weight components
    "QUESTID2", "FILEDATE",
    "VESTR_C", "VEREP"       # replicate/strata → DROP
]

cols_to_drop_with_weight += manual_dropping

# Remove duplicates
cols_to_drop_with_weight = list(set(cols_to_drop_with_weight))

len(cols_to_drop_with_weight), cols_to_drop_with_weight[:30]

(6, ['WTPOUND2', 'FILEDATE', 'VESTR_C', 'QUESTID2', 'WTANSWER', 'VEREP'])

df_removed_columns_LR = df_dropped_targetNaN.drop(columns=cols_to_drop_with_weight, errors='ignore')

# Checking to see which variables have >99% missingess
na_percent = df_removed_columns_LR.isna().mean() * 100

# Filter columns with > 99% missing
cols_over_99_LR = na_percent[na_percent > 99]

# Print count
print("Number of variables with >99% missingness:", len(cols_over_99_LR))

# Print the variable names
print("\nVariables with >99% missingness:")
print(cols_over_99_LR.index.tolist())

# Optional: show top 20 instead of the whole list
print("\nPreview (first 20):")
print(cols_over_99_LR.head(300))

Number of variables with >99% missingness: 282

Variables with >99% missingness:
['EDUSCKEST', 'EDUSKPEST', 'HLCALLFG', 'HLCALL99', 'BKPOSTOB', 'BKOTHOF2', 'ALTOTFG', 'ALFQFLG', 'ALDYSFG', 'MJFQFLG', 'BLRECFL2', 'BLNT30C1', 'BLNT30C2', 'RSNOMRJ', 'RSNMRJMO', 'CCTOTFG', 'CCFQFLG', 'CRTOTFG', 'CRFQFLG', 'HRTOTFG', 'HRFQFLG', 'HALTOTFG', 'HALFQFLG', 'INHTOTFG', 'INHFQFLG', 'METOTFG', 'MEFQFLG', 'PNRNORXFG', 'TRQNORXFG', 'STMNORXFG', 'SEDNORXFG', 'SRCSEDNM2', 'SRCFRPNRNM', 'SRCFRTRQNM', 'SRCFRSEDNM', 'SRCCLFRPNR', 'SRCCLFRTRQ', 'SRCCLFRSED', 'OTCFLAG', 'OTDGNDLA', 'OTDGNDLB', 'OTDGNDLC', 'OTDGNDLD', 'OTDGNDLE', 'SUNTINSCV', 'SUNTENCV', 'SUNTNOOPN', 'YEATNDYR', 'YEHMSLYR', 'YESCHFLT', 'YESCHWRK', 'YESCHIMP', 'YESCHINT', 'YETCGJOB', 'YELSTGRD', 'YECIGFRNDOF2', 'YECIGNEXTYR2', 'YESTSCIG', 'YESTSMJ', 'YESTSALC', 'YESTSDNK', 'YEPCHKHW', 'YEPHLPHW', 'YEPCHORE', 'YEPLMTTV', 'YEPLMTSN', 'YEPGDJOB', 'YEPPROUD', 'YEYARGUP', 'YEYFGTSW', 'YEYFGTGP', 'YEYHGUN', 'YEYSELL', 'YEYSTOLE', 'YEYATTAK', 'YEPPKCIG', 'YEPMJEVR', 'YEPMJMO', 'YEPALDLY', 'YEGPKCIG', 'YEGMJEVR', 'YEGMJMO', 'YEGALDLY', 'YEFPKCIG', 'YEFMJEVR', 'YEFMJMO', 'YEFALDLY', 'YETLKNON', 'YETLKPAR', 'YETLKBGF', 'YETLKOTA', 'YETLKSOP', 'YEPRTDNG', 'YEPRBSLV', 'YEVIOPRV', 'YEDGPRGP', 'YESLFHLP', 'YEPRGSTD', 'YESCHACT', 'YECOMACT', 'YEFAIACT', 'YEOTHACT', 'YEDECLAS', 'YEDERGLR', 'YEDESPCL', 'YEPVNTYR', 'YERLGSVC', 'YERLGIMP', 'YERLDCSN', 'YERLFRND', 'YUSUITHK', 'YUCOSUITHK', 'YUSUIPLN', 'YUCOSUIPLN', 'SCHFELT', 'TCHGJOB', 'AVGGRADE', 'STNDSCIG', 'STNDSMJ', 'STNDALC', 'STNDDNK', 'PARCHKHW', 'PARHLPHW', 'PRCHORE2', 'PRLMTTV2', 'PARLMTSN', 'PRGDJOB2', 'PRPROUD2', 'ARGUPAR', 'YOFIGHT2', 'YOGRPFT2', 'YOHGUN2', 'YOSELL2', 'YOSTOLE2', 'YOATTAK2', 'PRPKCIG2', 'PRMJEVR2', 'PRMJMO', 'PRALDLY2', 'YFLPKCG2', 'YFLTMRJ2', 'YFLMJMO', 'YFLADLY2', 'FRDPCIG2', 'FRDMEVR2', 'FRDMJMON', 'FRDADLY2', 'TALKPROB', 'PRTALK3', 'PRBSOLV2', 'PREVIOL2', 'PRVDRGO2', 'GRPCNSL2', 'PREGPGM2', 'YTHACT2', 'DRPRVME3', 'ANYEDUC3', 'RLGATTD', 'RLGIMPT', 'RLGDCSN', 'RLGFRND', 'YUSUITHKYR', 'YUCOSUITHK2', 'YUSUIPLNYR', 'YUCOSUIPLN2', 'YODPREV', 'YODSCEV', 'YOLOSEV', 'YODPDISC', 'YODPLSIN', 'YODSLSIN', 'YOLSI2WK', 'YODPR2WK', 'YOWRHRS', 'YOWRDST', 'YOWRCHR', 'YOWRIMP', 'YODPPROB', 'YOWRPROB', 'YOWRAGE', 'YOWRDEPR', 'YOWRDISC', 'YOWRLSIN', 'YOWRPLSR', 'YOWRELES', 'YOWREMOR', 'YOWRGAIN', 'YOWRGROW', 'YOWRPREG', 'YOWRGNL2', 'YOWRLOSE', 'YOWRDIET', 'YOWRLSL2', 'YOWRSLEP', 'YOWRSMOR', 'YOWRENRG', 'YOWRSLOW', 'YOWRSLNO', 'YOWRJITT', 'YOWRJINO', 'YOWRTHOT', 'YOWRCONC', 'YOWRDCSN', 'YOWRNOGD', 'YOWRWRTH', 'YO_MDEA1', 'YO_MDEA2', 'YO_MDEA3', 'YO_MDEA4', 'YO_MDEA5', 'YO_MDEA6', 'YO_MDEA7', 'YO_MDEA8', 'YODSMMDE', 'YOPBINTF', 'YOPBDLYA', 'YOPBRMBR', 'YOPBAGE', 'YOPBNUM', 'YOPB2WK', 'YOPSHMGT', 'YOPSWORK', 'YOPSRELS', 'YOPSSOC', 'YOPSDAYS', 'YOSEEDOC', 'YOFAMDOC', 'YOOTHDOC', 'YOPSYCH', 'YOPSYMD', 'YOSOCWRK', 'YOCOUNS', 'YOOTHMHP', 'YONURSE', 'YORELIG', 'YOHERBAL', 'YOOTHHLP', 'YOTMTNOW', 'YORX12MO', 'YORXNOW', 'YORXHLP', 'YOTMTHLP', 'YMDELT', 'YMDEYR', 'YMDEAUD5YR', 'YMIUD5YANY', 'YMSUD5YANY', 'YMDERSUD5ANY', 'YMDESUD5ANYO', 'YTXMDEYR', 'YRXMDEYR', 'YMDETXRX', 'YDOCMDE', 'YOMDMDE', 'YPSY1MDE', 'YPSY2MDE', 'YSOCMDE', 'YCOUNMDE', 'YOMHMDE', 'YNURSMDE', 'YRELMDE', 'YHBCHMDE', 'YHLTMDE', 'YALTMDE', 'YMDEHPRX', 'YMDEHPO', 'YMDERXO2', 'YMDEHARX', 'YSDSHOME', 'YSDSWRK', 'YSDSREL', 'YSDSSOC', 'YSDSOVRL', 'MDEIMPY', 'YMDEIMAD5YR', 'YMIMI5YANY', 'YMIMR5YANY', 'YMIMS5YANY', 'CIRROSAGE', 'HEPBCAGE', 'HIVAIDSAG', 'GQTYPE2']

Preview (first 20):
EDUSCKEST      99.915804
EDUSKPEST      99.960118
HLCALLFG       99.953471
HLCALL99       99.953471
BKPOSTOB      100.000000
                 ...    
YMIMS5YANY    100.000000
CIRROSAGE      99.725256
HEPBCAGE       99.058339
HIVAIDSAG      99.736335
GQTYPE2        99.811668
Length: 282, dtype: float64

df_99_LR = df_removed_columns_LR.drop(columns=cols_over_99_LR.index)
df_99_LR.shape

(45133, 2348)

# Remove target leakage features that was used to impute / recode the target IRAMDEYR

leakage_vars_LR = [

    # Target + imputation
    #"IRAMDEYR", "IIAMDEYR",
    "IIAMDEYR",

    # Pre-imputation past-year MDE
    "AMDEYR",

    # Lifetime MDE + imputation
    "AMDELT", "IRAMDELT", "IIAMDELT",

    # MDE with impairment (logically dependent on AMDEYR)
    "AMDEIMP", "IRAMDEIMP", "IIAMDEIMP",

    # DSM summary
    "ADSMMDEA",

    # Raw DSM symptom items (if present)
    "D_MDEA1", "D_MDEA2", "D_MDEA3", "D_MDEA4", "D_MDEA5",
    "D_MDEA6", "D_MDEA7", "D_MDEA8", "D_MDEA9",

    # Symptom recodes (one-digit suffixes, if present)
    "AD_MDEA1", "AD_MDEA2", "AD_MDEA3", "AD_MDEA4",
    "AD_MDEA5", "AD_MDEA6", "AD_MDEA7", "AD_MDEA8",

    # Symptom recodes (true DSM flags: two-digit suffix)
    "AD_MDEA11", "AD_MDEA21", "AD_MDEA31", "AD_MDEA41",
    "AD_MDEA51", "AD_MDEA61", "AD_MDEA71", "AD_MDEA81", "AD_MDEA91",

    # Fallback/gating items
    "ADPB2WK",
    "ADDPREV", "ADDSCEV", "ADLOSEV", "ADLSI2WK", "ADDPR2WK",
    "ADWRHRS", "ADWRDST", "ADWRCHR", "ADWRIMP", "ADDPPROB",

    # Soft-leakage variables
    "ARXMDEYR",   # received prescription meds
    "ATXMDEYR",   # received counseling/therapy
    "AHLTMDE",     # told by provider they have depression
]

# Drop only those leakage columns that actually exist in this PUF
cols_to_drop_leakaging = [c for c in leakage_vars_LR if c in df_99_LR.columns]
df_no_leakage_LR = df_99_LR.drop(columns=cols_to_drop_leakaging)
print("Dropped leakage columns:", cols_to_drop_leakaging)

Dropped leakage columns: ['IIAMDEYR', 'AMDEYR', 'AMDELT', 'IRAMDELT', 'IIAMDELT', 'AMDEIMP', 'IRAMDEIMP', 'IIAMDEIMP', 'ADSMMDEA', 'AD_MDEA1', 'AD_MDEA2', 'AD_MDEA3', 'AD_MDEA4', 'AD_MDEA5', 'AD_MDEA6', 'AD_MDEA7', 'AD_MDEA8', 'ADPB2WK', 'ADDPREV', 'ADDSCEV', 'ADLOSEV', 'ADLSI2WK', 'ADDPR2WK', 'ADWRHRS', 'ADWRDST', 'ADWRCHR', 'ADWRIMP', 'ADDPPROB', 'ARXMDEYR', 'ATXMDEYR', 'AHLTMDE']

df_no_leakage_LR.shape

(45133, 2317)

Define Target, Weight, and Base Feature Matrix

# Columns
TARGET_COL = "IRAMDEYR"
WEIGHT_COL = "ANALWT2_C"

# Sanity checks
assert TARGET_COL in df_no_leakage_LR.columns, "IRAMDEYR not found in df_no_leakage_LR"
assert WEIGHT_COL in df_no_leakage_LR.columns, "ANALWT2_C not found in df_no_leakage_LR"

# Raw target and weights
y_raw = df_no_leakage_LR[TARGET_COL]
w     = df_no_leakage_LR[WEIGHT_COL]

# Binary-encode target: 1 = past-year MDE, 0 = no MDE
y = (y_raw == 1).astype(int)

# Features = all columns except target + weight
X = df_no_leakage_LR.drop(columns=[TARGET_COL, WEIGHT_COL], errors='ignore')

print("Base X shape:", X.shape)
print("y positive rate:", y.mean().round(3))

Base X shape: (45133, 2315)
y positive rate: 0.115

Tier 1 and Subset Features

# Tier to evaluate
tier_name = "Tier 1"
tier_features = tiers[tier_name]

# Keep only those columns in X
X_tier = X[tier_features].copy()

print(f"{tier_name}: {len(tier_features)} features")
print("X_tier shape:", X_tier.shape)

Tier 1: 5 features
X_tier shape: (45133, 5)

Train / Test Split (with Weights Kept Aligned)

from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test, w_train, w_test = train_test_split(
    X_tier, y, w,
    test_size=0.2,
    random_state=42,
    stratify=y
)

print("Train shape:", X_train.shape, y_train.shape)
print("Test shape:",  X_test.shape,  y_test.shape)

Train shape: (36106, 5) (36106,)
Test shape: (9027, 5) (9027,)

Define Preprocessing (Median Impute + Scaling for Numeric)

'''Here we define numeric vs categorical and build a ColumnTransformer that:
-Imputes numeric with median, scales with StandardScaler
-Imputes categorical with mode, encodes with OneHotEncoder'''

'''Why median is better for NSDUH over mean imputation as NSDUH variables are mostly:
-Ordinal categories (0,1,2,3)
-Non-normal/long-tailed distributions
-Heavily skewed (e.g., hours, days, counts)
-Integer-coded survey recodes
-Mean imputation pushes values toward the wrong part of the distribution
-Median imputation better reflects the actual underlying "typical" respondent
-Median preserves the ordinal structure
-This prevents distortion of the predictors and reduces bias in logistic regression.'''

from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import StandardScaler
from sklearn.impute import SimpleImputer

# All numeric
numeric_cols = X_train.columns.tolist()
categorical_cols = []

numeric_transformer = Pipeline(steps=[
    ('imputer', SimpleImputer(strategy='median')),
    ('scaler',  StandardScaler())
])

preprocessor = ColumnTransformer(
    transformers=[
        ('num', numeric_transformer, numeric_cols)
    ],
    remainder='drop'
)

Building Logistic Regression Pipeline + CV Hyperparameter Tuning

from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import StratifiedKFold, GridSearchCV

log_reg_pipe = Pipeline(steps=[
    ('preprocess', preprocessor),
    ('model', LogisticRegression(
        max_iter=2000,
        solver='liblinear'
    ))
])

param_grid = {
    'model__C': [0.01, 0.1, 1.0, 10.0],
    'model__class_weight': [None, 'balanced']
}

cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

grid = GridSearchCV(
    estimator=log_reg_pipe,
    param_grid=param_grid,
    cv=cv,
    scoring='roc_auc',
    n_jobs=1,
    verbose=1
)

grid.fit(X_train, y_train)

print("Best params:", grid.best_params_)
print("Best CV ROC AUC:", grid.best_score_)

Fitting 5 folds for each of 8 candidates, totalling 40 fits
Best params: {'model__C': 1.0, 'model__class_weight': None}
Best CV ROC AUC: 0.8406564536348655

Refit best model on full training set with survey weights

# Refit best estimator WITH survey weights (correct way)
best_log_reg = grid.best_estimator_

best_log_reg.fit(
    X_train,
    y_train,
    model__sample_weight=w_train
)

Pipeline(steps=[('preprocess',
             ColumnTransformer(transformers=[(&#x27;num&#x27;,
                                              Pipeline(steps=[(&#x27;imputer&#x27;,
                                                               SimpleImputer(strategy=&#x27;median&#x27;)),
                                                              (&#x27;scaler&#x27;,
                                                               StandardScaler())]),
                                              [&#x27;ADPSDAYS&#x27;, &#x27;ASDSSOC2&#x27;,
                                               &#x27;ADPSHMGT&#x27;, &#x27;ADPSRELS&#x27;,
                                               &#x27;ADPSWORK&#x27;])])),
            (&#x27;model&#x27;,
             LogisticRegression(max_iter=2000, solver=&#x27;liblinear&#x27;))])</pre><b>In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. <br />On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.</b></div><div class="sk-container" hidden><div class="sk-item sk-dashed-wrapped"><div class="sk-label-container"><div class="sk-label fitted sk-toggleable"><input class="sk-toggleable__control sk-hidden--visually" id="sk-estimator-id-4" type="checkbox" ><label for="sk-estimator-id-4" class="sk-toggleable__label fitted sk-toggleable__label-arrow"><div><div>Pipeline</div></div><div><a class="sk-estimator-doc-link fitted" rel="noreferrer" target="_blank" href="https://scikit-learn.org/1.7/modules/generated/sklearn.pipeline.Pipeline.html">?<span>Documentation for Pipeline</span></a><span class="sk-estimator-doc-link fitted">i<span>Fitted</span></span></div></label><div class="sk-toggleable__content fitted" data-param-prefix="">
    <div class="estimator-table">
        <details>
            <summary>Parameters</summary>
            <table class="parameters-table">
              <tbody>

    <tr class="user-set">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('steps',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">steps&nbsp;</td>
        <td class="value">[(&#x27;preprocess&#x27;, ...), (&#x27;model&#x27;, ...)]</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('transform_input',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">transform_input&nbsp;</td>
        <td class="value">None</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('memory',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">memory&nbsp;</td>
        <td class="value">None</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('verbose',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">verbose&nbsp;</td>
        <td class="value">False</td>
    </tr>

              </tbody>
            </table>
        </details>
    </div>
</div></div></div><div class="sk-serial"><div class="sk-item sk-dashed-wrapped"><div class="sk-label-container"><div class="sk-label fitted sk-toggleable"><input class="sk-toggleable__control sk-hidden--visually" id="sk-estimator-id-5" type="checkbox" ><label for="sk-estimator-id-5" class="sk-toggleable__label fitted sk-toggleable__label-arrow"><div><div>preprocess: ColumnTransformer</div></div><div><a class="sk-estimator-doc-link fitted" rel="noreferrer" target="_blank" href="https://scikit-learn.org/1.7/modules/generated/sklearn.compose.ColumnTransformer.html">?<span>Documentation for preprocess: ColumnTransformer</span></a></div></label><div class="sk-toggleable__content fitted" data-param-prefix="preprocess__">
    <div class="estimator-table">
        <details>
            <summary>Parameters</summary>
            <table class="parameters-table">
              <tbody>

    <tr class="user-set">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('transformers',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">transformers&nbsp;</td>
        <td class="value">[(&#x27;num&#x27;, ...)]</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('remainder',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">remainder&nbsp;</td>
        <td class="value">&#x27;drop&#x27;</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('sparse_threshold',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">sparse_threshold&nbsp;</td>
        <td class="value">0.3</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('n_jobs',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">n_jobs&nbsp;</td>
        <td class="value">None</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('transformer_weights',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">transformer_weights&nbsp;</td>
        <td class="value">None</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('verbose',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">verbose&nbsp;</td>
        <td class="value">False</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('verbose_feature_names_out',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">verbose_feature_names_out&nbsp;</td>
        <td class="value">True</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('force_int_remainder_cols',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">force_int_remainder_cols&nbsp;</td>
        <td class="value">&#x27;deprecated&#x27;</td>
    </tr>

              </tbody>
            </table>
        </details>
    </div>
</div></div></div><div class="sk-parallel"><div class="sk-parallel-item"><div class="sk-item"><div class="sk-label-container"><div class="sk-label fitted sk-toggleable"><input class="sk-toggleable__control sk-hidden--visually" id="sk-estimator-id-6" type="checkbox" ><label for="sk-estimator-id-6" class="sk-toggleable__label fitted sk-toggleable__label-arrow"><div><div>num</div></div></label><div class="sk-toggleable__content fitted" data-param-prefix="preprocess__num__"><pre>[&#x27;ADPSDAYS&#x27;, &#x27;ASDSSOC2&#x27;, &#x27;ADPSHMGT&#x27;, &#x27;ADPSRELS&#x27;, &#x27;ADPSWORK&#x27;]</pre></div></div></div><div class="sk-serial"><div class="sk-item"><div class="sk-serial"><div class="sk-item"><div class="sk-estimator fitted sk-toggleable"><input class="sk-toggleable__control sk-hidden--visually" id="sk-estimator-id-7" type="checkbox" ><label for="sk-estimator-id-7" class="sk-toggleable__label fitted sk-toggleable__label-arrow"><div><div>SimpleImputer</div></div><div><a class="sk-estimator-doc-link fitted" rel="noreferrer" target="_blank" href="https://scikit-learn.org/1.7/modules/generated/sklearn.impute.SimpleImputer.html">?<span>Documentation for SimpleImputer</span></a></div></label><div class="sk-toggleable__content fitted" data-param-prefix="preprocess__num__imputer__">
    <div class="estimator-table">
        <details>
            <summary>Parameters</summary>
            <table class="parameters-table">
              <tbody>

    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('missing_values',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">missing_values&nbsp;</td>
        <td class="value">nan</td>
    </tr>


    <tr class="user-set">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('strategy',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">strategy&nbsp;</td>
        <td class="value">&#x27;median&#x27;</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('fill_value',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">fill_value&nbsp;</td>
        <td class="value">None</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('copy',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">copy&nbsp;</td>
        <td class="value">True</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('add_indicator',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">add_indicator&nbsp;</td>
        <td class="value">False</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('keep_empty_features',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">keep_empty_features&nbsp;</td>
        <td class="value">False</td>
    </tr>

              </tbody>
            </table>
        </details>
    </div>
</div></div></div><div class="sk-item"><div class="sk-estimator fitted sk-toggleable"><input class="sk-toggleable__control sk-hidden--visually" id="sk-estimator-id-8" type="checkbox" ><label for="sk-estimator-id-8" class="sk-toggleable__label fitted sk-toggleable__label-arrow"><div><div>StandardScaler</div></div><div><a class="sk-estimator-doc-link fitted" rel="noreferrer" target="_blank" href="https://scikit-learn.org/1.7/modules/generated/sklearn.preprocessing.StandardScaler.html">?<span>Documentation for StandardScaler</span></a></div></label><div class="sk-toggleable__content fitted" data-param-prefix="preprocess__num__scaler__">
    <div class="estimator-table">
        <details>
            <summary>Parameters</summary>
            <table class="parameters-table">
              <tbody>

    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('copy',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">copy&nbsp;</td>
        <td class="value">True</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('with_mean',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">with_mean&nbsp;</td>
        <td class="value">True</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('with_std',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">with_std&nbsp;</td>
        <td class="value">True</td>
    </tr>

              </tbody>
            </table>
        </details>
    </div>
</div></div></div></div></div></div></div></div></div></div><div class="sk-item"><div class="sk-estimator fitted sk-toggleable"><input class="sk-toggleable__control sk-hidden--visually" id="sk-estimator-id-9" type="checkbox" ><label for="sk-estimator-id-9" class="sk-toggleable__label fitted sk-toggleable__label-arrow"><div><div>LogisticRegression</div></div><div><a class="sk-estimator-doc-link fitted" rel="noreferrer" target="_blank" href="https://scikit-learn.org/1.7/modules/generated/sklearn.linear_model.LogisticRegression.html">?<span>Documentation for LogisticRegression</span></a></div></label><div class="sk-toggleable__content fitted" data-param-prefix="model__">
    <div class="estimator-table">
        <details>
            <summary>Parameters</summary>
            <table class="parameters-table">
              <tbody>

    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('penalty',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">penalty&nbsp;</td>
        <td class="value">&#x27;l2&#x27;</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('dual',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">dual&nbsp;</td>
        <td class="value">False</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('tol',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">tol&nbsp;</td>
        <td class="value">0.0001</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('C',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">C&nbsp;</td>
        <td class="value">1.0</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('fit_intercept',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">fit_intercept&nbsp;</td>
        <td class="value">True</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('intercept_scaling',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">intercept_scaling&nbsp;</td>
        <td class="value">1</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('class_weight',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">class_weight&nbsp;</td>
        <td class="value">None</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('random_state',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">random_state&nbsp;</td>
        <td class="value">None</td>
    </tr>


    <tr class="user-set">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('solver',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">solver&nbsp;</td>
        <td class="value">&#x27;liblinear&#x27;</td>
    </tr>


    <tr class="user-set">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('max_iter',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">max_iter&nbsp;</td>
        <td class="value">2000</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('multi_class',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">multi_class&nbsp;</td>
        <td class="value">&#x27;deprecated&#x27;</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('verbose',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">verbose&nbsp;</td>
        <td class="value">0</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('warm_start',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">warm_start&nbsp;</td>
        <td class="value">False</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('n_jobs',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">n_jobs&nbsp;</td>
        <td class="value">None</td>
    </tr>


    <tr class="default">
        <td><i class="copy-paste-icon"
             onclick="copyToClipboard('l1_ratio',
                      this.parentElement.nextElementSibling)"
        ></i></td>
        <td class="param">l1_ratio&nbsp;</td>
        <td class="value">None</td>
    </tr>

              </tbody>
            </table>
        </details>
    </div>
</div></div></div></div></div></div></div><script>function copyToClipboard(text, element) {
// Get the parameter prefix from the closest toggleable content
const toggleableContent = element.closest('.sk-toggleable__content');
const paramPrefix = toggleableContent ? toggleableContent.dataset.paramPrefix : '';
const fullParamName = paramPrefix ? `${paramPrefix}${text}` : text;

const originalStyle = element.style;
const computedStyle = window.getComputedStyle(element);
const originalWidth = computedStyle.width;
const originalHTML = element.innerHTML.replace('Copied!', '');

navigator.clipboard.writeText(fullParamName)
    .then(() => {
        element.style.width = originalWidth;
        element.style.color = 'green';
        element.innerHTML = "Copied!";

        setTimeout(() => {
            element.innerHTML = originalHTML;
            element.style = originalStyle;
        }, 2000);
    })
    .catch(err => {
        console.error('Failed to copy:', err);
        element.style.color = 'red';
        element.innerHTML = "Failed!";
        setTimeout(() => {
            element.innerHTML = originalHTML;
            element.style = originalStyle;
        }, 2000);
    });
return false;

}
document.querySelectorAll('.fa-regular.fa-copy').forEach(function(element) {
const toggleableContent = element.closest('.sk-toggleable__content');
const paramPrefix = toggleableContent ? toggleableContent.dataset.paramPrefix : '';
const paramName = element.parentElement.nextElementSibling.textContent.trim();
const fullParamName = paramPrefix ? ${paramPrefix}${paramName} : paramName;
element.setAttribute('title', fullParamName);

});

Predicting on test set (Probabilities + Class Labels)
# Probability for class 1 (MDE)
y_proba = best_log_reg.predict_proba(X_test)[:, 1]

# Hard labels using 0.5 threshold
y_pred = (y_proba >= 0.5).astype(int)
Compute Weighted Metrics: Accuracy, Sensitivity, Specificity, PPV, NPV, F1, F2, ROC AUC, PR AUC
from sklearn.metrics import (
    accuracy_score,
    recall_score,
    precision_score,
    f1_score,
    fbeta_score,
    roc_auc_score,
    precision_recall_curve,
    auc,
    confusion_matrix
)

# Confusion matrix with weights
tn, fp, fn, tp = confusion_matrix(
    y_test, y_pred,
    sample_weight=w_test
).ravel()

# Basic metrics
accuracy    = accuracy_score(y_test, y_pred, sample_weight=w_test)
sensitivity = recall_score(y_test, y_pred, sample_weight=w_test)          # TPR / recall
specificity = tn / (tn + fp) if (tn + fp) > 0 else np.nan

ppv = precision_score(y_test, y_pred, sample_weight=w_test)               # Positive Predictive Value
npv = tn / (tn + fn) if (tn + fn) > 0 else np.nan                         # Negative Predictive Value

f1 = f1_score(y_test, y_pred, sample_weight=w_test)
f2 = fbeta_score(y_test, y_pred, beta=2, sample_weight=w_test)

# ROC AUC (weighted)
roc_auc = roc_auc_score(y_test, y_proba, sample_weight=w_test)

# PR AUC (via precision-recall curve)
precisions, recalls, _ = precision_recall_curve(
    y_test, y_proba,
    sample_weight=w_test
)
pr_auc = auc(recalls, precisions)

print(f"=== {tier_name} — Weighted Test Metrics ===")
print(f"Accuracy:              {accuracy:.3f}")
print(f"Sensitivity (Recall):  {sensitivity:.3f}")
print(f"Specificity:           {specificity:.3f}")
print(f"PPV (Precision):       {ppv:.3f}")
print(f"NPV:                   {npv:.3f}")
print(f"F1:                    {f1:.3f}")
print(f"F2:                    {f2:.3f}")
print(f"ROC AUC:               {roc_auc:.3f}")
print(f"PR AUC:                {pr_auc:.3f}")
print("Confusion matrix (weighted):")
print(f"  TN: {tn:.1f},  FP: {fp:.1f}")
print(f"  FN: {fn:.1f},  TP: {tp:.1f}")
=== Tier 1 — Weighted Test Metrics ===
Accuracy:              0.972
Sensitivity (Recall):  0.635
Specificity:           1.000
PPV (Precision):       0.991
NPV:                   0.971
F1:                    0.774
F2:                    0.684
ROC AUC:               0.829
PR AUC:                0.836
Confusion matrix (weighted):
  TN: 47082423.8,  FP: 21888.1
  FN: 1426902.7,  TP: 2477163.9

Bootstrapped 95% CI for ROC AUC and F1
n_bootstrap = 500   # can bump to 1000 if laptop handles it
rng = np.random.RandomState(42)

roc_values = []
f1_values  = []

y_test_arr = y_test.to_numpy()
w_test_arr = w_test.to_numpy()

for i in range(n_bootstrap):
    # resample indices with replacement
    idx = rng.randint(0, len(y_test_arr), len(y_test_arr))
    
    y_b   = y_test_arr[idx]
    proba_b = y_proba[idx]
    pred_b  = y_pred[idx]
    w_b   = w_test_arr[idx]

    # ROC AUC
    try:
        roc_values.append(
            roc_auc_score(y_b, proba_b, sample_weight=w_b)
        )
    except ValueError:
        # can happen if bootstrap sample has only one class
        continue

    # F1
    f1_values.append(
        f1_score(y_b, pred_b, sample_weight=w_b)
    )

roc_ci = (np.percentile(roc_values, 2.5), np.percentile(roc_values, 97.5))
f1_ci  = (np.percentile(f1_values, 2.5),  np.percentile(f1_values, 97.5))

print(f"ROC AUC 95% CI: [{roc_ci[0]:.3f}, {roc_ci[1]:.3f}]")
print(f"F1      95% CI: [{f1_ci[0]:.3f},  {f1_ci[1]:.3f}]")
ROC AUC 95% CI: [0.795, 0.861]
F1      95% CI: [0.737,  0.807]

import numpy as np
import pandas as pd
from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import StandardScaler
from sklearn.impute import SimpleImputer
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import StratifiedKFold, GridSearchCV, train_test_split
from sklearn.metrics import (
    accuracy_score, recall_score, precision_score, f1_score, fbeta_score,
    roc_auc_score, precision_recall_curve, auc, confusion_matrix
)

def evaluate_tier(tier_name, tier_features, X, y, w):

    print(f"\n=== Running Logistic Regression for {tier_name} ({len(tier_features)} features) ===")

    # Subset X to this tier
    X_tier = X[tier_features].copy()

    # Train-test split
    X_train, X_test, y_train, y_test, w_train, w_test = train_test_split(
        X_tier, y, w,
        test_size=0.2,
        random_state=42,
        stratify=y
    )

    # Preprocessing (median impute + scaling)
    numeric_cols = X_train.columns.tolist()
    numeric_transformer = Pipeline(steps=[
        ('imputer', SimpleImputer(strategy='median')),
        ('scaler', StandardScaler())
    ])

    preprocessor = ColumnTransformer(
        transformers=[('num', numeric_transformer, numeric_cols)],
        remainder='drop'
    )

    # Logistic Regression Pipeline + CV
    log_reg_pipe = Pipeline(steps=[
        ('preprocess', preprocessor),
        ('model', LogisticRegression(
            max_iter=2000,
            solver='liblinear'
        ))
    ])

    param_grid = {
        'model__C': [0.01, 0.1, 1.0, 10.0],
        'model__class_weight': [None, 'balanced']
    }

    cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

    grid = GridSearchCV(
        estimator=log_reg_pipe,
        param_grid=param_grid,
        cv=cv,
        scoring='roc_auc',
        n_jobs=1,
        verbose=0
    )

    grid.fit(X_train, y_train)

    best_model = grid.best_estimator_

    # Refit using survey weights
    best_model.fit(X_train, y_train, model__sample_weight=w_train)

    # Predict
    y_proba = best_model.predict_proba(X_test)[:, 1]
    y_pred = (y_proba >= 0.5).astype(int)

    # Weighted metrics
    tn, fp, fn, tp = confusion_matrix(y_test, y_pred, sample_weight=w_test).ravel()

    metrics = {
        "Tier": tier_name,
        "Feature_Count": len(tier_features),
        "Accuracy": accuracy_score(y_test, y_pred, sample_weight=w_test),
        "Sensitivity": recall_score(y_test, y_pred, sample_weight=w_test),
        "Specificity": tn / (tn + fp),
        "PPV": precision_score(y_test, y_pred, sample_weight=w_test),
        "NPV": tn / (tn + fn),
        "F1": f1_score(y_test, y_pred, sample_weight=w_test),
        "F2": fbeta_score(y_test, y_pred, beta=2, sample_weight=w_test),
        "ROC_AUC": roc_auc_score(y_test, y_proba, sample_weight=w_test),
        "PR_AUC": auc(*precision_recall_curve(
            y_test, y_proba, sample_weight=w_test
        )[1::-1]),
        "TN": tn, "FP": fp, "FN": fn, "TP": tp
    }

    print(f"Finished {tier_name} ✓")

    return metrics
results = []

for tier_name, tier_features in tiers.items():
    metrics = evaluate_tier(tier_name, tier_features, X, y, w)
    results.append(metrics)

# Convert to DataFrame for comparison
results_df = pd.DataFrame(results)

results_df
=== Running Logistic Regression for Tier 0 (Top 100) (100 features) ===
Finished Tier 0 (Top 100) ✓

=== Running Logistic Regression for Tier 1 (5 features) ===
Finished Tier 1 ✓

=== Running Logistic Regression for Tier 2 (12 features) ===
Finished Tier 2 ✓

=== Running Logistic Regression for Tier 3 (29 features) ===
Finished Tier 3 ✓

=== Running Logistic Regression for Tier 4 (52 features) ===
Finished Tier 4 ✓




  
    
      
      Tier
      Feature_Count
      Accuracy
      Sensitivity
      Specificity
      PPV
      NPV
      F1
      F2
      ROC_AUC
      PR_AUC
      TN
      FP
      FN
      TP
    
  
  
    
      0
      Tier 0 (Top 100)
      100
      0.995128
      0.991536
      0.995425
      0.947268
      0.999296
      0.968897
      0.982354
      0.999676
      0.997569
      4.688882e+07
      215490.696164
      3.304401e+04
      3.871023e+06
    
    
      1
      Tier 1
      5
      0.971597
      0.634509
      0.999535
      0.991241
      0.970585
      0.773737
      0.683721
      0.829097
      0.835851
      4.708242e+07
      21888.125341
      1.426903e+06
      2.477164e+06
    
    
      2
      Tier 2
      12
      0.985758
      0.970803
      0.986997
      0.860880
      0.997554
      0.912543
      0.946628
      0.997614
      0.982910
      4.649183e+07
      612482.848786
      1.139883e+05
      3.790078e+06
    
    
      3
      Tier 3
      29
      0.985917
      0.971832
      0.987085
      0.861812
      0.997640
      0.913522
      0.947637
      0.998514
      0.987752
      4.649595e+07
      608366.414519
      1.099692e+05
      3.794097e+06
    
    
      4
      Tier 4
      52
      0.986358
      0.969399
      0.987763
      0.867826
      0.997439
      0.915805
      0.947226
      0.998456
      0.987261
      4.652790e+07
      576409.973102
      1.194673e+05
      3.784599e+06
    
  


Results Summary
What you found:
Tier 1 (5 features) works but misses a lot of true MDE cases:

Sensitivity = 0.63 → It catches 63% of real depression cases.
Specificity = 1.00 → It never falsely marks healthy cases as depression.
It’s “ultra-conservative”—great for avoiding false alarms, but not good for screening.

Tier 2–4 (12–52 features) perform extremely well:

Sensitivity ≈ 0.97 → They catch 97% of all real MDE cases.
Specificity ≈ 0.987 → They correctly identify healthy people 98.7% of the time.
AUC ≈ 0.998 → Near-perfect discrimination.
These models perform at the same level as the best clinical screeners used in psychiatry research, and in some aspects even better.

Tier 0 (100 features) is the maximum-possible-performance tier:

It achieves nearly perfect separation: AUC 0.9997
But it uses too many predictors to be practical.

A model using about 29–52 variables (Tier 3–4) can detect past-year major depression with extremely high accuracy: around 97% sensitivity and 99% specificity. These tiers offer nearly the same predictive value as using 100 features, but with far fewer inputs—meaning they can be adapted into clinical workflows.
Final verdict: Using only 12–30 variables, the model can detect depression cases with about 97% sensitivity and 99% specificity. This is comparable to or better than many psychiatric screening tools. The smallest tier (5 variables) is much more conservative and misses cases, so it isn’t suitable for screening. The larger tiers show that adding more contextual variables dramatically improves ability to identify patients with depression.Using only 12–30 variables, the model can detect depression cases with about 97% sensitivity and 99% specificity. This is comparable to or better than many psychiatric screening tools. The smallest tier (5 variables) is much more conservative and misses cases, so it isn’t suitable for screening. The larger tiers show that adding more contextual variables dramatically improves ability to identify patients with depression.
Detailed Technical Analysis
Tier 0 (Top 100 features)
Best possible / upper bound:

Accuracy: 0.995
Sensitivity: 0.992
Specificity: 0.995
ROC AUC: 0.9997
PR AUC: 0.9975

Interpretation:

The model almost perfectly separates depressed vs. non-depressed.
This shows that NSDUH data contains very strong signal for detecting IRAMDEYR, even after removing leakage.
But this is not deployable—too large, not clinically feasible.

Tier 1 (5 features)

Accuracy: 0.972
Sensitivity: 0.635
Specificity: 0.9995
ROC AUC: 0.829

Interpretation:

Very high specificity → almost no false positives
Low sensitivity → misses many real cases
It’s not good enough for a screening use case but fine for confirmatory use.
Why?: With only 5 predictors, the model doesn’t capture enough variation in clinical presentations.

Tier 2 (12 features)

Sensitivity: 0.971
Specificity: 0.987
ROC AUC: 0.9976

Interpretation:

Massive jump in sensitivity (+33% over Tier 1).
This shows that around 10–15 features are enough to fully capture the underlying construct.
Tier 2 is already extremely strong.

Tier 3 (29 features)

Sensitivity: 0.972
Specificity: 0.987
ROC AUC: 0.9985

Interpretation:

Performance increases slightly over Tier 2.
This is likely the clinically optimal tradeoff between performance and complexity.

Tier 4 (52 features)

Sensitivity: 0.969
Specificity: 0.988
ROC AUC: 0.9984

Interpretation:

Very similar to Tier 3.
You gain almost no performance beyond ~30 features.

	Tier	Feature_Count	Accuracy	Sensitivity	Specificity	PPV	NPV	F1	F2	ROC_AUC	PR_AUC	TN	FP	FN	TP
0	Tier 0 (Top 100)	100	0.995128	0.991536	0.995425	0.947268	0.999296	0.968897	0.982354	0.999676	0.997569	4.688882e+07	215490.696164	3.304401e+04	3.871023e+06
1	Tier 1	5	0.971597	0.634509	0.999535	0.991241	0.970585	0.773737	0.683721	0.829097	0.835851	4.708242e+07	21888.125341	1.426903e+06	2.477164e+06
2	Tier 2	12	0.985758	0.970803	0.986997	0.860880	0.997554	0.912543	0.946628	0.997614	0.982910	4.649183e+07	612482.848786	1.139883e+05	3.790078e+06
3	Tier 3	29	0.985917	0.971832	0.987085	0.861812	0.997640	0.913522	0.947637	0.998514	0.987752	4.649595e+07	608366.414519	1.099692e+05	3.794097e+06
4	Tier 4	52	0.986358	0.969399	0.987763	0.867826	0.997439	0.915805	0.947226	0.998456	0.987261	4.652790e+07	576409.973102	1.194673e+05	3.784599e+06