<jupyter_start><jupyter_text>The most common traits of substance abuseIn 2017, nearly 20 million American adults (almost 7% of the population) battled a substance abuse problem, costing the US more than $740 billion annually in lost workplace productivity, healthcare expenses, and crime-related costs, according to the [American Addiction Centers](https://americanaddictioncenters.org/rehab-guide/addiction-statistics).Substance abuse has a cure, but many people are too ashamed to reach out for help, which makes the problem even harder to address.What if we could tell whether someone has a potential problem with drugs without asking them explicitly? Maybe we could look out for common traits among people who suffer with substance abuse problems and provide aid in advance. Organizations could target ad campaigns at the right people and right moment and increase their effectiveness.But what would these common traits be? 1. Problem understandingThis study will try to answer the following questions:1. Which age group shows the largest share of people with a possible substance problem? And which gender?2. Which kind of professional is more vulnerable to substance abuse?3. Are people with kids or partners less likely to have a substance problem?4. Are religious people less likely to have a substance problem?5. Does making healthy choices (e.g., not smoking or eating plant-based diets) correlate negatively with substance abuse?And last but not least, how accurately can we predict substance abuse with just this bunch of information?This notebook focus on the technical aspects of the study. Please visit [my blog post](https://brunoarine.com/blog/the-most-common-traits-of-substance-abuse/) where I discuss my findings in an accessible language. 2. Data understandingTo tackle the previous questions, I'm going to use the [OkCupid Profile Dataset](https://github.com/rudeboybert/JSE_OkCupid), a file containing answers from 60k users covering a number of topics; from marital status, job, and kids, to alcohol and drugs consumption.<jupyter_code>import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from scipy.stats import chi2_contingency
df = pd.read_csv("../data/profiles.csv")<jupyter_output><empty_output><jupyter_text>Selecting the useful variablesThese are the variables that will help me address the proposed questions:- diet- drinks- drugs- job- offspring- pets- religion- sex- smokes- statusLet's ditch all the rest:- body_type- education- ethnicity- height- income- last_online- location- orientation- sign- speaks- essay from 0 to 9<jupyter_code>df.drop(["body_type", "education", "ethnicity", "height", "income", "last_online", "location", "orientation",
        "sign", "speaks", "essay0", "essay1", "essay2", "essay3", "essay4", "essay5", "essay6",
        "essay7", "essay8", "essay9"], axis=1, inplace=True)<jupyter_output><empty_output><jupyter_text>Variable descriptionsThe following details are also provided in the dataset codebook.Missing data is blank. Some variables have two factors, which are denoted in this codebook by a semicolon (e.g. "graduated from; two-year college) though they are found without punctuation in the dataset. Details for specific variables are found in parenthesis.**diet** - mostly/strictly; anything, vegetarian, vegan, kosher, halal, other**drinks** - very often, often, socially, rarely, desperately, not at all**drugs** - never, sometimes, often**job** - student, art/music/writing, banking/finance, administration, technology, construction, education, entertainment/media, management, hospitality, law, medicine, military, politics/government, sales/marketing, science/engineering, transportation, unemployed, other, rather not say, retire**offspring** - has a kid, has kids, doesnt have a kid, doesn't want kids; ,and/,but might want them, wants them, doesnt want any, doesnt want more**orientation** - straight, gay, bisexual**pets** - has dogs, likes dogs, dislikes dogs; and has cats, likes cats, dislikes cats**religion** - agnosticism, atheism, Christianity, Judaism, Catholicism, Islam, Hinduism, Buddhism, Other; and very serious about it, and somewhat serious about it, but not too serious about it, and laughing about it**sex** - m, f**smokes** - yes, sometimes, when drinking, trying to quit, no**status** - single, seeing someone, married, in an open relationship Let's check if the actual data type is the same as their description in the codebook.<jupyter_code>df.info()<jupyter_output><class 'pandas.core.frame.DataFrame'>
RangeIndex: 59946 entries, 0 to 59945
Data columns (total 11 columns):
 #   Column     Non-Null Count  Dtype 
---  ------     --------------  ----- 
 0   age        59946 non-null  int64 
 1   diet       35551 non-null  object
 2   drinks     56961 non-null  object
 3   drugs      45866 non-null  object
 4   job        51748 non-null  object
 5   offspring  24385 non-null  object
 6   pets       40025 non-null  object
 7   religion   39720 non-null  object
 8   sex        59946 non-null  object
 9   smokes     54434 non-null  object
 10  status     59946 non-null  object
dtypes: int64(1), object(10)
memory usage: 5.0+ MB<jupyter_text>Missing values<jupyter_code>df.isnull().mean()<jupyter_output><empty_output><jupyter_text>The database has 60k rows, but there's a huge number of missing values, especially in the `offspring` column. People have their own personal reasons not to answer a survey, like taboo or fear of exposition. That makes it a MNAR (Missing Not At Random) scenario, which is quite difficult to handle. I'm not a huge fan of imputing the missing values with estimates in this case (even with multivariate methods like KNN) because doing so will likely create some sort of bias.<jupyter_code>df.dropna(axis=0).shape<jupyter_output><empty_output><jupyter_text>That will would leave us only ~15% of original data, which is still beter than coming up with estimates (and risking affect the $\chi^2$ tests later).Nonetheless, there's no need to drop all the missing rows right now. While analyzing the variables individually, I can drop *just the missing values that belong to the variables being analyzed.* That way, I still got a bit of statistical power for chi-squared tests. After the individual analyzes are done, then all missing values will be dropped before training the machine learning model.In sum: the NaNs-containing rows are going to be dropped shortly before I evaluate the final model. Until then, only the missing values in the target variable will be dropped, along with the missing values in each of the column being analyzed along Section 4 - Evaluation.  3. Preparing the data  Defining and creating the target variableFirst, we need to clearly define what substance abuse is. According to [Wikipedia](https://en.wikipedia.org/wiki/Substance_abuse):> Substance abuse, also known as drug abuse, is use of a drug in amounts or by methods which are harmful to the individual or others. (...) Drugs most often associated with this term include: alcohol, amphetamines, barbiturates, benzodiazepines, cannabis, cocaine, hallucinogens, methaqualone, and opioids.The OkCupid dataset isn't precise in terms of which drugs its users consume, but it gives us a few clues about quantity: - The `drinks` variable has 6 possible values: very often, often, socially, rarely, desperately, not at all- The `drugs` variable has 3 possible values: never, sometimes, oftenI'm going to assume that users who drink "very often", or use drugs "often" are likely candidates for the group of users with a substance problem. I dismissed "desperately" entries because this sounds more like an answer to signal you're cool.I'm going to name it as `abuse`, and convert the possible values of `drinks` and `drugs` the variables `alcohol_problem`, `drugs_problem`, and `both_problems`, which will either be 0 or 1. But before that, I'm going to drop the missing values in both columns because I want the target column to be devoid of NaNs.<jupyter_code>df.dropna(axis=0, subset=["drinks", "drugs"], inplace=True)
df.shape<jupyter_output><empty_output><jupyter_text>Good, dropping the null values for the target variable left us with 75% of the original data, which is acceptable for me.<jupyter_code># Creates alcohol_problem column with values:
# 1 if drinks = "very often", 0 if anything else.
# Not considering the "desperately" entries improved the model accuracy at the end of this notebook,
# which kinda proves my point.
df["alcohol_problem"] = df['drinks'].str.contains("very often")

# Creates drugs_problem column with values:
# 1 if drinks = "often"
# 0 if anything else
df["drugs_problem"] = df['drugs'].str.contains("often")

# Users have an abuse problem (abuse = 1) if they have either an alcohol or drug abuse problem;
df['abuse'] = (df['alcohol_problem'] | df['drugs_problem'])<jupyter_output><empty_output><jupyter_text>A user may have a drug abuse problem, an alcohol abuse problem, or both. I want to discriminate each so that I can make a stacked bar plot in the analysis section.<jupyter_code>df['both_problems'] = (df["drugs_problem"] & df["alcohol_problem"])
df['alcohol_problem_only'] = (df['alcohol_problem'] & ~df['drugs_problem'])
df['drugs_problem_only'] = (~df['alcohol_problem'] & df['drugs_problem'])<jupyter_output><empty_output><jupyter_text>Age  Several populational traits have a potentially strong correlation with age (e.g. having kids). Therefore I want to create a variable called `age_group` so that I can assess whether my hypotheses hold for all ages when analyzing the plots.<jupyter_code>age_group_labels = ["18 to 25", "25 to 29", "30 to 39", "40+"]
df["age_group"] = pd.cut(df['age'], bins=[18,25,30,40,100], labels=age_group_labels)<jupyter_output><empty_output><jupyter_text>The decision upon age intervals was such that each bin contained roughly the same number of users.<jupyter_code>df["age_group"].value_counts()<jupyter_output><empty_output><jupyter_text>Gender<jupyter_code># Creates a female variable, so that gender can be either 0 or 1 to maintain the same format as
# the other variables
df["female"] = df["sex"].replace({"m": 0, "f": 1})<jupyter_output><empty_output><jupyter_text>Jobs<jupyter_code># Create an "unemployed" variable to reduce the number of features for the machine learning model
df["unemployed"] = df["job"].copy()
df.loc[df["job"].notna(), "unemployed"] = (df["job"] == "unemployed") | (df["job"] == "rather not say")<jupyter_output><empty_output><jupyter_text>Status The `okcupid_codebook.txt` file tells that all missing values are presented as an empty string in the CSV file. That's not the case with `status`. (I thought it was strange that such a column didn't have one single missing value in a dataset with 60k answers!) The missing values here are stated as "unknown". I'll replace those with NaNs.<jupyter_code>df['status'].replace("unknown", np.nan, inplace=True)<jupyter_output><empty_output><jupyter_text>According to the `okcupid_codebook.txt` file, the possible values of `status` are single, seeing someone, married, and in an open relationship. I'll replace the value "available" with the latter since it's more explainative.<jupyter_code>df["status"] = df["status"].replace({"available": "in an open relationship"})

# Create a "married" variable to reduce the number of features for the machine learning model
df["married"] = df["status"].str.contains("married")<jupyter_output><empty_output><jupyter_text>Offspring  Whoever answered "doesn't have kids" (and might or might not want them), or simply "wants kids" will have the `offspring` variable set to zero, while everything else will be replaced with ones.<jupyter_code>df['offspring'].replace(['doesn&rsquo;t have kids',
                         'doesn&rsquo;t have kids, and doesn&rsquo;t want any',
                         'doesn&rsquo;t have kids, but might want them',
                         'doesn&rsquo;t have kids, but wants them',
                         'might want kids',
                         'doesn&rsquo;t want kids',
                         'wants kids'], 0, inplace=True)

df['offspring'].replace(['has a kid',
                         'has a kid, and might want more',
                         'has a kid, and wants more',
                         'has a kid, but doesn&rsquo;t want more',
                         'has kids',
                         'has kids, and might want more',
                         'has kids, and wants more',
                         'has kids, but doesn&rsquo;t want more',], 1, inplace=True)<jupyter_output><empty_output><jupyter_text>Pets<jupyter_code># 1 if person has dogs (or cats), 0 otherwise
df['has_dogs'] = df['pets'].str.contains('has dogs')
df['has_cats'] = df['pets'].str.contains('has cats')

# Make a copy of pets so that all NaNs are at the right position, then
# has_pets = 1 when user either has a dog or a cat
df['has_pets'] = df['pets']
df.loc[df['pets'].notna(), 'has_pets'] = (df['has_dogs'] == 1) | (df['has_cats'] == 1)<jupyter_output><empty_output><jupyter_text>Religion The `religion` variable has several possible values plus a few modifiers ("laughing about it, "somewhat serious about it", etc). I thought of creating a new variable out of these modifiers called `religious`, which is going to be 1 if:- `religion` is not atheism or agnosticism (for obvious reasons), and- the `religion` modifier is "very serious about it".<jupyter_code># religious is 1 if "very serious", 0 otherwise
df['religious'] = df['religion'].str.contains('very serious')

# Conditions: religious remains 1 IF original value is 1 AND not atheist or agnostic
cond_religious = df['religious'] == 1
cond_notna = df['religious'].notna()
cond_not_atheist = ~(df['religion'].str.contains("agnosticism") | df['religion'].str.contains("atheism"))

# Replace only rows where values are not nulls. If I don't do that, condition testing will turn
# missing values into boolean, because NaN doesn't contain strings "agnosticism" or "atheism"
df.loc[cond_notna, "religious"] = cond_religious & cond_not_atheist<jupyter_output><empty_output><jupyter_text>Healthy variables Smokers will be smokers, independently of frequency, so setting the variable to 0 who answer "no" and 1 for everyone else.<jupyter_code>df["smokes"] = df["smokes"].replace({"yes": 1,
                                    "sometimes": 1,
                                    "when drinking": 1,
                                    "trying to quit": 1,
                                    "no": 0})<jupyter_output><empty_output><jupyter_text>The "diet" variable has 6 possible values (anything, vegetarian, vegan, kosher, halal, and other) and 2 modifiers (mostly/strictly), which might create unnecessary dummy variables that won't add value to the model. So let's boil it down to diet type without the nuance.<jupyter_code>df.loc[df['diet'].notna(), 'vegetarian'] = df['diet'].str.contains("vegetarian") | df['diet'].str.contains("vegan")<jupyter_output><empty_output><jupyter_text>4. Evaluation Auxiliary functions:<jupyter_code>def chi_test(df, column):
    """Calculates the p-value of chi-squared tests for a chosen column. The function also
    performs the test for all age group in the DataFrame.
    """
    df_dropna = df.dropna(axis=0, subset=[column, 'abuse'])
    print("Chi-squared test")
    print("----------------")
    for group in age_group_labels:
        df_group = df_dropna[df_dropna["age_group"] == group]
        cont_table = pd.crosstab(df_group[column], df_group['abuse'])
        c, p, dof, expected = chi2_contingency(cont_table)
        print(f"Age group: {group}, p = {p:.2g}, chi2 = {c:.2g}, n = {len(df_group[column])}")
  
def vert_barplot(df, bars, label, groupby, title=None, color=None, figsize=(9,5)):
    """Groups the DataFrame by a specified label (like "age_group), and plots a pretty bar chart
    with custom parameters.
    
    Parameters
    ----------
    
    df : DataFrame
    bars : list of tuples
        Contains pairs of columns labels and values to be plotted.
        Example: [("offspring", 0), ("offspring, 1")]
    groupby : str
        Column label by which the DataFrame is going to be grouped
    title : str
        x-axis label
    color : list
        List of colors to be passed to the plotting method.
    figsize : tuple
        Figure size
    """
    
    df_dropna = df.dropna(axis=0, subset=[x[0] for x in bars].append("abuse"))
    table_mean = []
    for i in range(len(bars)):
        col_name = bars[i][0]
        col_value = bars[i][1]
        table_mean.append(df_dropna[df_dropna[col_name] == col_value]
                                .groupby(groupby)["abuse"]
                                .mean()
                                .to_frame())
    table_p = pd.concat(table_mean, axis=1)
    
    table_percent = table_p * 100  # percentages should be multiplied by 100
    
    
    table_percent.columns = label
    table_percent.plot.bar(figsize=figsize,
                           color=color,
                           rot=0,
                           capsize=0)
    plt.ylabel("% of users with substance problem")
    plt.xlabel(title)
    sns.despine()
    plt.tight_layout()  # prevents excess white space in the borders
    plt.savefig(f"../results/{bars[0][0]}.png")
    plt.show()
  
def horz_barplot(df, column, figsize=(10,10), barwidth=0.75, legend_height=1.05):
    """Groups the DataFrame by a specified label (like "age_group), and plots a pretty horizontal
    bar chart with custom parameters.
    
    Parameters
    ----------
    
    df : DataFrame
    column : str
        Column label to be stacked
    figsize : tuple
        Figure size
    barwidth : float
        Bar width passed to Matplotlib
    legend_height : float
        Legend box position above the chart (relative pos)
    """
    df_dropna = df.dropna(axis=0, subset=[column, "abuse"])
    column_stacked = df_dropna.groupby(column).mean()[['alcohol_problem_only',
                                                   'drugs_problem_only',
                                                   'both_problems',
                                                   'abuse']].sort_values('abuse')*100
    sns.set_palette("rocket_r")
    column_stacked.drop('abuse', axis=1, inplace=True)
    column_stacked.plot(kind='barh', stacked=True, figsize=figsize, width=barwidth, linewidth=0)
    plt.ylabel("")
    plt.xlabel("% of users with possible substance problems")
    plt.gca().invert_yaxis()
    plt.legend(loc='upper center', 
               bbox_to_anchor=(0.5, legend_height),
              ncol=3, fancybox=True, shadow=True)
    sns.despine()
    plt.tight_layout()
    plt.savefig(f"../results/{column}.png")
    plt.show()<jupyter_output><empty_output><jupyter_text>Answering the questions<jupyter_code>df["abuse"].value_counts()<jupyter_output><empty_output><jupyter_text>Less than 3% of the users in this dataset have a probable substance problem. This is a specially imbalanced dataset, so any attempts to predict an outcome will likely be impacted. Techniques to deal with it can be found at the end of the document.  Question 1 - Which age group shows the largest percentage of people with a possible substance addiction problem? And does substance abuse disproportionately affect any gender?<jupyter_code>gender_colors = sns.color_palette("colorblind").as_hex()
color=[gender_colors[0], gender_colors[4]]
vert_barplot(df=df,
             bars=[("female", 0), ("female", 1)],
             label=["male", "female"],
             groupby="age_group",
             title="age group",
             color=color)
chi_test(df, "female")<jupyter_output>Chi-squared test
----------------
Age group: 18 to 25, p = 0.58, chi2 = 0.31, n = 10472
Age group: 25 to 29, p = 0.66, chi2 = 0.19, n = 12611
Age group: 30 to 39, p = 0.0011, chi2 = 11, n = 12734
Age group: 40+, p = 0.016, chi2 = 5.8, n = 8202<jupyter_text>Question 2 - Which kind profession is more vulnerable to substance abuse?<jupyter_code>horz_barplot(df, "job")<jupyter_output><empty_output><jupyter_text>Question 3 - Are people with kids or partners less likely to have a substance abuse problem?<jupyter_code>horz_barplot(df, "status", figsize=(10,4), barwidth=0.6, legend_height=1.3)<jupyter_output><empty_output><jupyter_text>Young adults are much less likely to have kids, and that will impact the outcome of the correlation study with substance abuse. It would be more reasonable to observe differences by age groups.<jupyter_code>table_offspring = pd.concat([df["age_group"].value_counts(), (1-df.groupby("age_group")["offspring"].mean())*100], axis=1)
table_offspring.columns = ["no. of users", "percentage without kids"]
table_offspring.index.rename("age group", inplace=True)
table_offspring
vert_barplot(df=df,
             bars=[("offspring", 1), ("offspring",0)], 
             label=["with kids", "without kids"],
             groupby="age_group",
             title="age group",
             color=sns.color_palette("Paired"))
chi_test(df, "offspring")
vert_barplot(df=df,
             bars=[("has_cats", 1), ("has_dogs",0), ("has_pets", 0)], 
             label=["cat owners", "dog owners", "no pets"],
             groupby="age_group",
             title="age group",
             color=sns.color_palette("Paired"))
chi_test(df, "has_cats")
chi_test(df, "has_dogs")<jupyter_output>Chi-squared test
----------------
Age group: 18 to 25, p = 0.89, chi2 = 0.02, n = 7112
Age group: 25 to 29, p = 0.13, chi2 = 2.3, n = 8496
Age group: 30 to 39, p = 0.59, chi2 = 0.29, n = 8604
Age group: 40+, p = 0.9, chi2 = 0.016, n = 5880<jupyter_text>Question 4 - Are religious people less likely to have a substance abuse problem?<jupyter_code>vert_barplot(df=df,
             bars=[("religious", 1), ("religious",0)], 
             label=["takes religion seriously", "atheist or cool about it"],
             groupby="age_group",
             title="age group",
             color=sns.color_palette("Paired"))
chi_test(df, "religious")<jupyter_output>Chi-squared test
----------------
Age group: 18 to 25, p = 7.5e-06, chi2 = 20, n = 7099
Age group: 25 to 29, p = 0.00058, chi2 = 12, n = 8496
Age group: 30 to 39, p = 0.0058, chi2 = 7.6, n = 8744
Age group: 40+, p = 0.046, chi2 = 4, n = 6391<jupyter_text>Question 5 - Does making healthy choices (e.g., not smoking or eating plant-based diets) correlate negatively with substance abuse?<jupyter_code>palette = sns.color_palette("colorblind").as_hex()
color = [palette[2], palette[1]]
vert_barplot(df=df,
             bars=[("vegetarian", 1), ("vegetarian",0)], 
             label=["vegetarian", "non-vegetarian"],
             groupby="age_group",
             title="age group",
             color=color)
chi_test(df, "vegetarian")
vert_barplot(df=df,
             bars=[("smokes", 1), ("smokes",0)], 
             label=["smoker", "non-smoker"],
             groupby="age_group",
             title="age group",
             color=sns.color_palette("rocket_r", n_colors=2))
chi_test(df, "smokes")<jupyter_output>Chi-squared test
----------------
Age group: 18 to 25, p = 4.3e-68, chi2 = 3e+02, n = 9966
Age group: 25 to 29, p = 1.1e-55, chi2 = 2.5e+02, n = 12004
Age group: 30 to 39, p = 9.1e-36, chi2 = 1.6e+02, n = 12284
Age group: 40+, p = 7.3e-10, chi2 = 38, n = 8012<jupyter_text>Model evaluation Let's clear up all unused columns and drop all rows with missing values.<jupyter_code>df_useful = df.drop(["diet","pets","religion","age_group", "alcohol_problem",
                    "alcohol_problem_only", "drugs_problem", "drugs_problem_only",
                    "both_problems", "has_pets", "job", "status", "sex","drinks", "drugs"], axis=1)

df_useful.dropna(axis=0, inplace=True)<jupyter_output><empty_output><jupyter_text>There are a few categorical variables remaining. Let's transform them into dummy variables.<jupyter_code>df_useful.info()
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import precision_score, roc_auc_score, f1_score, recall_score, plot_roc_curve
from sklearn.utils import resample
from sklearn.pipeline import make_pipeline
from sklearn.model_selection import StratifiedKFold, cross_val_score<jupyter_output><empty_output><jupyter_text>To ensure convergence, I could create a pipeline and pre-treat the data with StandardScaler. But I want the data to be as pure as possible, and the LogisticRegression without penalty. All of that so to ensure that I get **interpretability** from the LR coefficients.<jupyter_code>model = LogisticRegression(penalty="none",  max_iter=200)

y = df_useful["abuse"]
X = df_useful.drop(["abuse"], axis=1)

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)<jupyter_output><empty_output><jupyter_text>Our dataset is severely imbalanced and all attempts to train the model will be impacted unless we deal with it first. I decided to upsample the training set, which means randomly replicating positive cases until we have a 50/50 class split in the data set.<jupyter_code># concatenate our training data back together
df_train = pd.concat([X_train, y_train], axis=1)

# separate minority and majority classes
not_abuse = df_train[df_train["abuse"] == 0]
abuse =  df_train[df_train["abuse"] == 1]

# upsample minority
abuse_upsampled = resample(abuse,
                          replace=True, # sample with replacement
                          n_samples=len(not_abuse), # match number in majority class
                          random_state=42) # reproducible results

# combine majority and upsampled minority
upsampled = pd.concat([not_abuse, abuse_upsampled])

y_train = upsampled.abuse
X_train = upsampled.drop('abuse', axis=1)
model.fit(X_train, y_train)<jupyter_output><empty_output><jupyter_text>Taking the exponent of model's coefficients will return the odds ratio for each feature. For instance, if you're female you're half as likely to have a substance abuse problem, whereas being unenmployed makes it 4 to 5 times more likely to have a substance abuse problem.<jupyter_code>coef_table = pd.DataFrame(model.coef_[0], index=X_train.columns, columns=["coef"])
np.exp(coef_table)
palette = sns.color_palette("colorblind").as_hex()
colors = [palette[3] if x < 0 else palette[2] for x in coef_table.values]
pos = np.arange(len(coef_table))
plt.figure(figsize=(10,5))
plt.barh(pos, coef_table["coef"], color=colors)
plt.yticks(pos, coef_table.index)
plt.axvline(x=0, color='.2')
sns.despine()
plt.tight_layout()
plt.savefig("../results/coef.png")
plt.show()
model.predict_proba(X_test)[:,1]
print("Accuracy:",model.score(X_test, y_test))
print("Recall:",recall_score(y_test, model.predict(X_test)))
print("ROC AUC:",roc_auc_score(y_test, model.predict_proba(X_test)[:,1]))<jupyter_output>Accuracy: 0.7796810395747195
Recall: 0.825
ROC AUC: 0.8442604355716878<jupyter_text>According to the metrics, the model performed better than I initially expected, even dropping all missing values in the table and working with 15% of the original data.<jupyter_code>plot_roc_curve(model, X_test, y_test)
1 - (y_test[y_test == 0] == model.predict(X_test[y_test == 0])).mean()<jupyter_output><empty_output><jupyter_text>The ROC curve shows that the model is decent at detecting positive cases, but it still needs a bit of work concerning false positives. At its optimal decision threshold (near to the curve knee), the model has 82.5% probability of correctly detecting positive cases and 22.2% probability of false alarm, which would be slightly annoying in a real world application, to say the least.  How other models will perform?We could try other kinds of models and see how they perform in comparison the Logistic Regression (like the Random Forest).<jupyter_code>from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import GridSearchCV
# The more estimators the better, until a certain point. Let's start with 500.
forest = RandomForestClassifier(n_estimators=500, max_depth=3)

forest.fit(X_train, y_train)
print("Random Forest auROC:", roc_auc_score(y_test, forest.predict_proba(X_test)[:,1]))<jupyter_output>Random Forest auROC: 0.8414322444041138<jupyter_text>Why did it perform so badly? Maybe doing a cross-validated GridSearch to see if things can be improved.<jupyter_code>parameters = {"n_estimators": [200,400,500,600,800,900,1000],
              "criterion": ["gini", "entropy"],
              "max_depth": [None, 5, 10, 15, 30, 60],
           }
clf = GridSearchCV(forest, parameters, scoring='roc_auc', verbose=3, n_jobs=5)
clf.fit(X_train, y_train)
for params, mean, std in zip(clf.cv_results_['params'],
                             clf.cv_results_['mean_test_score'],
                             clf.cv_results_['std_test_score']):
    print("{} {:.3f} +/- {:.3f}".format(params,mean,std))<jupyter_output>{'criterion': 'gini', 'max_depth': None, 'n_estimators': 200} 0.953 +/- 0.004
{'criterion': 'gini', 'max_depth': None, 'n_estimators': 400} 0.953 +/- 0.003
{'criterion': 'gini', 'max_depth': None, 'n_estimators': 500} 0.953 +/- 0.003
{'criterion': 'gini', 'max_depth': None, 'n_estimators': 600} 0.953 +/- 0.003
{'criterion': 'gini', 'max_depth': None, 'n_estimators': 800} 0.953 +/- 0.003
{'criterion': 'gini', 'max_depth': None, 'n_estimators': 900} 0.953 +/- 0.003
{'criterion': 'gini', 'max_depth': None, 'n_estimators': 1000} 0.953 +/- 0.003
{'criterion': 'gini', 'max_depth': 5, 'n_estimators': 200} 0.836 +/- 0.006
{'criterion': 'gini', 'max_depth': 5, 'n_estimators': 400} 0.836 +/- 0.006
{'criterion': 'gini', 'max_depth': 5, 'n_estimators': 500} 0.836 +/- 0.005
{'criterion': 'gini', 'max_depth': 5, 'n_estimators': 600} 0.836 +/- 0.006
{'criterion': 'gini', 'max_depth': 5, 'n_estimators': 800} 0.836 +/- 0.006
{'criterion': 'gini', 'max_depth': 5, 'n_estimators': 900} 0.835 +/- 0.005
{'c[...]<jupyter_text>Weird. The standard deviation for all models are interestingly small. And the model performs great whatever the parameters. Even the number of trees doesn't significantly influence the model's performance, so that we can get away with a pretty decent score with just 200 estimators. But as we saw before, the model performs badly when tested against the validation set. What could be going wrong?<jupyter_code>from sklearn.metrics import plot_confusion_matrix
# Plot non-normalized confusion matrix
titles_options = [("Confusion matrix, without normalization", None),
                  ("Normalized confusion matrix", 'true')]
for title, normalize in titles_options:
    disp = plot_confusion_matrix(forest, X_test, y_test,
                                # display_labels=class_names,
                                 cmap=plt.cm.Blues,
                                 normalize=normalize)
    disp.ax_.set_title(title)

    print(title)
    print(disp.confusion_matrix)

plt.show()<jupyter_output>Confusion matrix, without normalization
[[1418  235]
 [  23   17]]
Normalized confusion matrix
[[0.85783424 0.14216576]
 [0.575      0.425     ]]<jupyter_text>The problem is that Random Forests are based on Decision Trees, which are highly sensible to unbalanced classes. The models are probably _overfitting_, and that would explain why the models score so high in the cross-validation GridSearch. In that case, I'm gonna use the most realistic value (with a limited max_depth), and see how it performs against the Logistic Regression.<jupyter_code>forest = RandomForestClassifier(n_estimators=500, max_depth=3)
forest.fit(X_train, y_train)

fig, ax = plt.subplots()
plot_roc_curve(model, X_test, y_test, ax=ax, color="orange")
plot_roc_curve(forest, X_test, y_test, ax=ax, color="blue")
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1122976122976123, 'fraction_numerical': 0.026964026964026964, 'mean_word_length': 3.9948080553807426, 'pattern_counts': {'":': 11, '<': 142, '<?xml version=': 0, '>': 142, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '37', 'dir': 'jupyter-structured-clean-dedup', 'id': '385591', 'license': 'CC0-1.0', 'n_tokens_mistral': 10374, 'n_tokens_neox': 9401, 'n_words': 3335, 'path': 'notebooks/project.ipynb', 'repo_name': 'brunoarine/vegetarian-market-study', 'size': '31746'}