import warnings
# Ignore warnings related to future Pandas changes (eventual syntax recommendations)
warnings.simplefilter(action='ignore', category=FutureWarning)
# from google.colab import files
import pandas as pd
import numpy as np

df = pd.read_csv("healthcare-dataset-stroke-data.csv")

# I will prepare and clean the original dataset by converting to dataframe,
# handling missing data, converting data types, and renaming column
# names for consistency in capitalization.

# After each step, I will display the head of the dataframe.

# Replace missing BMIs with the mean BMI
mean_bmi = df["bmi"].mean()
df["bmi"].fillna(mean_bmi, inplace = True)
df.head()

# Check how many "Unknown" rows there are to determine how to handle
unknown_count = df["smoking_status"].value_counts().get("Unknown", 0)
print(unknown_count/len(df)) # as a proportion
# Because unknown is in nearly 1/3 of the dataframe's rows, we want to replace
# it (case by case) with data that makes sense

0.3021526418786693

# Group based on gender and heart disease, find mode for each category
# Similarly to in HW 3, replace each Unknown in a given category with the corresponding mode
df["smoking_status"] = df["smoking_status"].replace("Unknown", np.nan)

male_disease = (df["gender"] == "Male") & (df["heart_disease"] == 1)
male_disease_mode = df[male_disease]["smoking_status"].mode()[0]
df.loc[male_disease, "smoking_status"] = df[male_disease]["smoking_status"].fillna(male_disease_mode)

female_disease = (df["gender"] == "Female") & (df["heart_disease"] == 1)
female_disease_mode = df[female_disease]["smoking_status"].mode()[0]
df.loc[female_disease, "smoking_status"] = df[female_disease]["smoking_status"].fillna(female_disease_mode)

male_no_disease = (df["gender"] == "Male") & (df["heart_disease"] == 0)
male_no_disease_mode = df[male_no_disease]["smoking_status"].mode()[0]
df.loc[male_no_disease, "smoking_status"] = df[male_no_disease]["smoking_status"].fillna(male_no_disease_mode)

female_no_disease = (df["gender"] == "Female") & (df["heart_disease"] == 0)
female_no_disease_mode = df[female_no_disease]["smoking_status"].mode()[0]
df.loc[female_no_disease, "smoking_status"] = df[female_no_disease]["smoking_status"].fillna(female_no_disease_mode)

df.head()

# The "ever_married" column consists of Yes/No, so I will replace these with 1/0 respectively
df.loc[df["ever_married"] == "Yes", "ever_married"] = 1
df.loc[df["ever_married"] == "No", "ever_married"] = 0

# This will make it easier for us to analyze this column alongside other
# columns with 0/1 values, such as "hypertension" and "heart_disease."

df.head()

# Enforce that all column names should be strictly lowercase for consistency
# Since "Residence_type" is the only name that doesn't follow this
df = df.rename(columns = {"Residence_type": "residence_type"})
df.head()

df.to_csv("cleaned_stroke_dataset.csv", index = False) # saving new df for convenient access
# files.download("cleaned_stroke_dataset.csv")

# we will load our cleaned dataset into a new dataframe
new_df = pd.read_csv("cleaned_stroke_dataset.csv")

# we want to learn more about our dataset that we have first, before conducting hypothesis testing
# let's first check if our data has properly been cleaned by seeing if there is a deviation in the missing values
# additionally, we want to check if the type of the value is correct
new_df.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 5110 entries, 0 to 5109
Data columns (total 12 columns):
 #   Column             Non-Null Count  Dtype  
---  ------             --------------  -----  
 0   id                 5110 non-null   int64  
 1   gender             5110 non-null   object 
 2   age                5110 non-null   float64
 3   hypertension       5110 non-null   int64  
 4   heart_disease      5110 non-null   int64  
 5   ever_married       5110 non-null   int64  
 6   work_type          5110 non-null   object 
 7   residence_type     5110 non-null   object 
 8   avg_glucose_level  5110 non-null   float64
 9   bmi                5110 non-null   float64
 10  smoking_status     5110 non-null   object 
 11  stroke             5110 non-null   int64  
dtypes: float64(3), int64(5), object(4)
memory usage: 479.2+ KB

# now, let's find the summary statistics (count, mean, std, min, max, quartiles) of our dataset
new_df.describe()

# lets count the number of features and entries
new_df.shape

(5110, 12)

# we want to check which one of the cateogircal features are over-represented
# are there more of one gender?
display(new_df['gender'].value_counts())

# are there more of one work type?
display(new_df['work_type'].value_counts())

# are there more of one residence type?
display(new_df['residence_type'].value_counts())

# are there more of one smoking status?
display(new_df['smoking_status'].value_counts())

gender
Female    2994
Male      2115
Other        1
Name: count, dtype: int64

work_type
Private          2925
Self-employed     819
children          687
Govt_job          657
Never_worked       22
Name: count, dtype: int64

residence_type
Urban    2596
Rural    2514
Name: count, dtype: int64

smoking_status
never smoked       3406
formerly smoked     915
smokes              789
Name: count, dtype: int64

# now, we want to see if there are any outliers in our dataset in terms of age, avg_glucose_level, and bmi
# this is because these are the only numerical variables that are not binary
# best way to do this will be to see it in a boxplot (the outliers are dots)

import matplotlib.pyplot as plt

# for age
# reference: lec 8, slide 23
fig, ax = plt.subplots()
# makes the box plot
ax.boxplot(new_df['age'], tick_labels = ['Age'])
# labels
ax.set_xlabel('Category')
ax.set_ylabel('Years')
ax.set_title('Observation of Outliers for Age')
plt.show()

# for avg_glucose_level
# reference: lec 8, slide 23
fig1, ax1 = plt.subplots()
# makes the box plot
ax1.boxplot(new_df['avg_glucose_level'], tick_labels = ['Avg. Glucose Level'])
# labels
ax1.set_xlabel('Category')
ax1.set_ylabel('mg/dL')
ax1.set_title('Observation of Outliers for Avg. Glucose Level')
plt.show()

# for bmi
# reference: lec 8, slide 23
fig2, ax2 = plt.subplots()
# makes the box plot
ax2.boxplot(new_df['bmi'], tick_labels = ['BMI'])
# labels
ax2.set_xlabel('Category')
ax2.set_ylabel('kg/m^2')
ax2.set_title('Observation of Outliers for BMI')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

glucose_stroke = new_df[new_df['stroke'] == 1]['avg_glucose_level']
glucose_no_stroke = new_df[new_df['stroke'] == 0]['avg_glucose_level']

# Plotting the histogram for the mean glucose level of those who had a stroke

plt.hist(glucose_stroke)
plt.title("Distribution of Average Glucose Level For Those Who Had a Stroke")
plt.xlabel("Average Glucose Level")
plt.ylabel("Frequency")

# Show plot
plt.show()

# Plotting the histogram for the mean glucose level of those who did not have a stroke

plt.hist(glucose_no_stroke)
plt.title("Distribution of Average Glucose Level For Those Who Did Not Have a Stroke")
plt.xlabel("Average Glucose Level")
plt.ylabel("Frequency")

# Show plot
plt.show()

from scipy.stats import mannwhitneyu

res = mannwhitneyu(glucose_stroke, glucose_no_stroke, alternative="greater")

res.pvalue

np.float64(1.8201836355446622e-09)

# Plot a boxplot
box = plt.boxplot(
    [glucose_stroke, glucose_no_stroke],tick_labels=["Stroke", "No Stroke"],
    patch_artist=True
)

# Color the boxes
box['boxes'][0].set_facecolor('lightgreen')
box['boxes'][1].set_facecolor('skyblue')

plt.title("Average Glucose Levels For Those Who Had a Stroke vs Did Not Have a Stroke")
plt.xlabel("Groups")
plt.ylabel("Average Glucose Levels")

# Simple legend
plt.legend(["Stroke", "No Stroke"])

plt.show()

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
# load our cleaned dataset into a new dataframe
new_df = pd.read_csv("cleaned_stroke_dataset.csv")

#isolate the dataframe to only individuals who have suffered a stroke
stroke_df = new_df[new_df['stroke'] == 1]
#bmis of the group of individuals who formerly smoked
bmi_former = new_df[(new_df['stroke'] == 1) &
 (new_df['smoking_status'] == 'formerly smoked')]['bmi']
#bmis of the group of individuals who never smoked
bmi_never = new_df[(new_df['stroke'] == 1) &
 (new_df['smoking_status'] == 'never smoked')]['bmi']
#bmis of the group of individuals who current smoke
bmi_smokes = new_df[(new_df['stroke'] == 1) &
 (new_df['smoking_status'] == 'smokes')]['bmi']

print("Testing for normality of the BMIs across the three groups")
from scipy import stats

print(stats.shapiro(bmi_former))
print(stats.shapiro(bmi_never))
print(stats.shapiro(bmi_smokes))

Testing for normality of the BMIs across the three groups
ShapiroResult(statistic=np.float64(0.9210027177418553), pvalue=np.float64(0.00017487815375835016))
ShapiroResult(statistic=np.float64(0.9080854687207409), pvalue=np.float64(1.78520786822663e-07))
ShapiroResult(statistic=np.float64(0.967286938347141), pvalue=np.float64(0.2669849457546159))

# Plotting the histogram for the BMIs of former smokers who suffered a stroke

plt.hist(bmi_former)
plt.title("Distribution of BMIs for Former Smokers Who Suffered a Stroke")
plt.xlabel("BMI (kg/m\u00b2)")
plt.ylabel("Frequency")

# Show plot
plt.show()

# Plotting the histogram for the BMIs of never smokers who suffered a stroke

plt.hist(bmi_never)
plt.title("Distribution of BMIs for Never Smokers Who Suffered a Stroke")
plt.xlabel("BMI (kg/m\u00b2)")
plt.ylabel("Frequency")

# Show plot
plt.show()

# Plotting the histogram for the BMIs of current smokers who suffered a stroke

plt.hist(bmi_smokes)
plt.title("Distribution of BMIs for Current Smokers Who Suffered a Stroke")
plt.xlabel("BMI (kg/m\u00b2)")
plt.ylabel("Frequency")

# Show plot
plt.show()

#Perform a Kruskal-Wallis test
statistic, p_value = stats.kruskal(bmi_former, bmi_never, bmi_smokes)

print(f"The p-value is: {p_value:.4e}")

The p-value is: 3.9884e-01

colors = ['skyblue', 'salmon', 'lightgreen', 'gold']

ax = new_df.groupby('smoking_status')['bmi'].median().plot(
    kind='bar',
    color=colors
)

plt.title('Median BMIs across Smoking Groups')
plt.xlabel('Smoking Groups')
plt.ylabel('Median BMI (kg/m\u00b2)')
ax.set_xticklabels([label.get_text().capitalize() for label in ax.get_xticklabels()],
                   rotation=45)
ax.legend(title='Smoking Status')

plt.show()

import seaborn as sns



df = stroke_df[["smoking_status", "bmi"]].dropna()

order = ["never smoked", "formerly smoked", "smokes"]
df = df[df["smoking_status"].isin(order)]

groups = [
    df[df["smoking_status"] == cat]["bmi"]
    for cat in order
]

stat, p = stats.kruskal(*groups)

sns.set(style="whitegrid", context="talk")
plt.figure(figsize=(10, 6))

sns.violinplot(
    data=df,
    x="smoking_status",
    y="bmi",
    order=order,
    palette="Set2",
    inner=None
)

sns.boxplot(
    data=df,
    x="smoking_status",
    y="bmi",
    order=order,
    width=0.25,
    showcaps=True,
    boxprops={'facecolor':'white', 'zorder':3},
    showfliers=False,
    whiskerprops={'linewidth':2}
)

sns.stripplot(
    data=df,
    x="smoking_status",
    y="bmi",
    order=order,
    color="black",
    alpha=0.35,
    jitter=True
)

medians = df.groupby("smoking_status")["bmi"].median()

for i, cat in enumerate(order):
    plt.scatter(i, medians[cat], color="red", s=120, zorder=5, label="Median" if i == 0 else "")

plt.title("BMI Distribution by Smoking Status (Kruskal–Wallis Test)", weight="bold")
plt.xlabel("Smoking Status")
plt.ylabel("BMI")

plt.text(
    1,
    max(df["bmi"]) * 0.98,
    f"Kruskal-Wallis p = {p:.3f}",
    ha='center',
    fontsize=12,
    bbox=dict(facecolor='white', alpha=0.6, edgecolor='gray')
)

plt.legend()
plt.tight_layout()
plt.show()

from scipy.stats import chi2_contingency

# Contingency table for hypertension and stroke
contingency_table = pd.crosstab(new_df["hypertension"], new_df["stroke"])
contingency_table

# Bar chart of raw counts
contingency_table.plot(kind="bar")

plt.title("Stroke Counts by Hypertension Status")
plt.xlabel("Hypertension")
plt.ylabel("Count")
plt.legend(title="Stroke", labels=["No Stroke (0)", "Stroke (1)"])
plt.xticks(rotation=0)

plt.show()

# Perform chi-squared test
chi2, p, dof, expected = chi2_contingency(contingency_table)

print("Chi-squared statistic:", chi2)
print("Degrees of freedom:", dof)
print("P-value:", p)
print("Expected frequencies:")
print(expected)

Chi-squared statistic: 81.6053682482931
Degrees of freedom: 1
P-value: 1.661621901511823e-19
Expected frequencies:
[[4387.2665362  224.7334638]
 [ 473.7334638   24.2665362]]

# Bar chart showing stroke proportions by hypertension status
proportions = contingency_table.div(contingency_table.sum(axis=1), axis=0)

proportions.plot(kind="bar")

plt.title("Hypertension vs Stroke")
plt.xlabel("Hypertension")
plt.ylabel("Proportion")
plt.legend(title="Stroke", labels=["No Stroke (0)", "Stroke (1)"])
plt.xticks(rotation=0)

plt.show()

# Stroke rate comparison
stroke_rate_no_htn = contingency_table.loc[0,1] / contingency_table.loc[0].sum()
stroke_rate_htn = contingency_table.loc[1,1] / contingency_table.loc[1].sum()

print("Stroke rate (no hypertension):", stroke_rate_no_htn)
print("Stroke rate (hypertension):", stroke_rate_htn)

Stroke rate (no hypertension): 0.03967909800520382
Stroke rate (hypertension): 0.13253012048192772

from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score, classification_report
import warnings
warnings.filterwarnings('ignore')

# initializing X and Y
# id and stroke are not features to train the model on
# ref: HW4 2.2
X = new_df.drop(columns = ['id', 'stroke'])
Y = new_df['stroke'].astype(int)

# ref: lec 13, slide 26-27
# use one-hot encoding to create numeric features
X = pd.get_dummies(X, columns=["gender", "residence_type", "work_type", "smoking_status"])

# ref: HW4 2.2
seed = 42
test_size = 0.2

# split the data 0.8 and 0.2 for train and test
X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size = test_size, random_state = seed)

# normalizing the train and test sets
s = StandardScaler()
# ref: lec 18, slide 31
X_train_scaled = s.fit_transform(X_train)
X_test_scaled = s.transform(X_test)

model = RandomForestClassifier(random_state = seed)

# ref: HW4 4.1
np.random.seed(seed)
# ref: lec 16, slide 46
model.fit(X_train_scaled, y_train)

# prints the accuracy of the Random Forest model
# also displays the classification report of the model
# ref: lec 12, slide 50
y_pred = model.predict(X_test_scaled)
accuracy = accuracy_score(y_test, y_pred)
print(f"Accuracy of Random Forest: {accuracy:.3f}")
print(classification_report(y_test, y_pred, target_names=['Bad Rating (False)', 'Good Rating (True)']))

# displays all the feature importances
# ref: HW4 5.1
importances = model.feature_importances_
feature_importance_df = pd.DataFrame({
    'Feature': X.columns,
    'Importance': importances
})

feature_importance_df = feature_importance_df.sort_values(by = 'Importance', ascending = False)

# this shows the top 5 important features
print("Feature Importances:")
print(feature_importance_df[:5])

Accuracy of Random Forest: 0.939
                    precision    recall  f1-score   support

Bad Rating (False)       0.94      1.00      0.97       960
Good Rating (True)       0.00      0.00      0.00        62

          accuracy                           0.94      1022
         macro avg       0.47      0.50      0.48      1022
      weighted avg       0.88      0.94      0.91      1022

Feature Importances:
             Feature  Importance
4  avg_glucose_level    0.273913
0                age    0.241790
5                bmi    0.228814
1       hypertension    0.028308
2      heart_disease    0.024642

# now, we will use the age, glucose level, and bmi as determining features
# ref: HW4 2.2
X = new_df[['age', 'avg_glucose_level', 'bmi']]
Y = new_df['stroke'].astype(int)

# split the data 0.8 and 0.2 for train and test
X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size = test_size, random_state = seed)

# ref: HW4 3.1
# got warnings without the class_weight arg
# this is because there were not many people that had a stroke compared to those that did not
# ref: https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html
model = LogisticRegression(random_state=seed, max_iter = 1000, class_weight = 'balanced')

# ref: HW4 4.1
np.random.seed(seed)
model.fit(X_train, y_train)

# prints the accuracy of the Logistic Regression model
# also displays the classification report of the model
# ref: lec 12, slide 50
y_pred = model.predict(X_test)
accuracy = accuracy_score(y_test, y_pred)
print(f"Accuracy of Logistic Regression: {accuracy:.3f}")
print(classification_report(y_test, y_pred, target_names=['Bad Rating (False)', 'Good Rating (True)']))

Accuracy of Logistic Regression: 0.745
                    precision    recall  f1-score   support

Bad Rating (False)       0.98      0.74      0.84       960
Good Rating (True)       0.17      0.81      0.28        62

          accuracy                           0.74      1022
         macro avg       0.58      0.77      0.56      1022
      weighted avg       0.93      0.74      0.81      1022

# now we can use our trained model to predict stroke risk
def predict_risk(age, avg_glucose_level, bmi):
    # this is the user input
    input_df = pd.DataFrame({
        "age": [age],
        "avg_glucose_level": [avg_glucose_level],
        "bmi": [bmi]
    })
    # this is just predict
    # ref: HW4 4.1
    pred = model.predict(input_df)

    # this just reprints what the user put in
    print("User Input:")
    print("Age:", age)
    print("Average Glucose Level:", avg_glucose_level)
    print("BMI:", bmi)

    # this is just saying if there is a high or low stroke risk
    if pred[0] == 1:
        print("\nResult: High Stroke Risk\n")
    else:
        print("\nResult: Low Stroke Risk\n")

# demo
predict_risk(age = 67, avg_glucose_level = 228.69, bmi = 36.6)

predict_risk(age = 40, avg_glucose_level = 288.9, bmi = 25)

User Input:
Age: 67
Average Glucose Level: 228.69
BMI: 36.6

Result: High Stroke Risk

User Input:
Age: 40
Average Glucose Level: 288.9
BMI: 25

Result: Low Stroke Risk

import plotly.graph_objects as go

# ref: HW4 6.2.3
def visualize_3d_plotly_interactive(X, Y, model):

    # ref: HW5 4.1
    glucose_values = np.linspace(X["avg_glucose_level"].min(), X["avg_glucose_level"].max(), 100)
    bmi_values = np.linspace(X["bmi"].min(), X["bmi"].max(), 100)

    # ref: https://numpy.org/doc/2.3/reference/generated/numpy.meshgrid.html
    # makes the values into tuples
    glucosev, bmiv = np.meshgrid(glucose_values, bmi_values)

    # solving for the decision boundary
    # since it is 3d, the model coeficients and interecept are stored in
    # the .coef_ and .intercept_ variables
    # ref: https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html
    agev = -((model.coef_[0][1] * glucosev + model.coef_[0][2] * bmiv + model.intercept_[0]) / model.coef_[0][0])

    # ref: HW4 6.2.3
    fig = go.Figure()

    # plotting the points without stroke
    fig.add_trace(go.Scatter3d(
        x = X[Y == 0]['avg_glucose_level'],
        y = X[Y == 0]['bmi'],
        z = X[Y == 0]['age'],
        mode = 'markers',
        marker = dict(size = 3, color = 'green', opacity = 0.3),
        name = 'No Stroke'
    ))

    # plotting the points with stroke
    fig.add_trace(go.Scatter3d(
        x = X[Y == 1]['avg_glucose_level'],
        y = X[Y == 1]['bmi'],
        z = X[Y == 1]['age'],
        mode = 'markers',
        marker = dict(size = 3, color = 'red', opacity = 0.7),
        name = 'Stroke'
    ))

    # plotting the decision boundary
    # ref: https://plotly.com/python/3d-surface-plots/
    fig.add_trace(go.Surface(
        x = glucosev,
        y = bmiv,
        z = agev,
        name = 'Decision Boundary',
        opacity = 0.5,
        # was getting a bar on the side that was not necessary
        showscale = False
    ))

    # updating the layout
    fig.update_layout(
        title = 'Logistic Regression Decision Boundary',
        scene = dict(
            xaxis_title = 'Avg. Glucose Level (mg/dL)',
            yaxis_title = 'BMI (kg/m^2)',
            zaxis_title = 'Age (Years)',
            camera = dict(
                eye = dict(x = 2.0, y = 2.0, z = 1.5)
            )
        ),
        width = 900,
        height = 700,
        showlegend = True
    )

    fig.show()

visualize_3d_plotly_interactive(X, Y, model)

	id	age	hypertension	heart_disease	ever_married	avg_glucose_level	bmi	stroke
count	5110.000000	5110.000000	5110.000000	5110.000000	5110.000000	5110.000000	5110.000000	5110.000000
mean	36517.829354	43.226614	0.097456	0.054012	0.656164	106.147677	28.893237	0.048728
std	21161.721625	22.612647	0.296607	0.226063	0.475034	45.283560	7.698018	0.215320
min	67.000000	0.080000	0.000000	0.000000	0.000000	55.120000	10.300000	0.000000
25%	17741.250000	25.000000	0.000000	0.000000	0.000000	77.245000	23.800000	0.000000
50%	36932.000000	45.000000	0.000000	0.000000	1.000000	91.885000	28.400000	0.000000
75%	54682.000000	61.000000	0.000000	0.000000	1.000000	114.090000	32.800000	0.000000
max	72940.000000	82.000000	1.000000	1.000000	1.000000	271.740000	97.600000	1.000000

Predicting Stroke Risk based on Health and Demographic Factors¶

Spring 2026 Data Science Project¶

Contributions¶

Introduction¶

Data Curation¶

Exploratory Data Analysis¶

Primary Analysis¶

Visualization¶

Insights and Conclusions¶

	id	gender	age	hypertension	heart_disease	ever_married	work_type	Residence_type	avg_glucose_level	bmi	smoking_status	stroke
0	9046	Male	67.0	0	1	Yes	Private	Urban	228.69	36.600000	formerly smoked	1
1	51676	Female	61.0	0	0	Yes	Self-employed	Rural	202.21	28.893237	never smoked	1
2	31112	Male	80.0	0	1	Yes	Private	Rural	105.92	32.500000	never smoked	1
3	60182	Female	49.0	0	0	Yes	Private	Urban	171.23	34.400000	smokes	1
4	1665	Female	79.0	1	0	Yes	Self-employed	Rural	174.12	24.000000	never smoked	1

stroke	0	1
hypertension
0	4429	183
1	432	66