Lecture 4 - Demo Notebook

Prelimilaries: Imports and stuff

%matplotlib inline
import numpy as np 
import pandas as pd
import matplotlib.pyplot as plt

from sklearn.model_selection import train_test_split
from sklearn import tree
from sklearn.ensemble import RandomForestClassifier
from sklearn.neighbors import KNeighborsClassifier
from sklearn.metrics import roc_auc_score, balanced_accuracy_score
from sklearn.preprocessing import MinMaxScaler, normalize
from sklearn.linear_model import LogisticRegression

import statsmodels.api as sm
import statsmodels.formula.api as smf

from scipy.spatial.distance import pdist, cdist, squareform

# Data directory
DATA_DIR = "./../../data"

Section 1: Pre-processing

Read the data

# Parse the aggregated data frame
df_lq = pd.read_csv('{}/aggregated_extended_fc.csv'.format(DATA_DIR))
ts = pd.read_csv('{}/time_series_extended_fc.csv'.format(DATA_DIR))

Clean the data

We remove inactive students that did not click during weekdays and weekend for the first 5 weeks of the semester.

def remove_inactive_students(df, ts):
    df = df.fillna('NaN')
    
    #find all users weeks with 0 clicks on weekends and 0 clicks on weekdays during the first weeks of the semester
    df_first = ts[ts.week < 5]
    rows = np.where(np.logical_and(df_first.ch_total_clicks_weekend==0, df_first.ch_total_clicks_weekday == 0).to_numpy())[0]
    df_zero = df_first.iloc[rows,:]
    dropusers = np.unique(df_zero.user)

    ts = ts[ts.user.isin(dropusers)==False]
    df = df[df.user.isin(dropusers)==False]
    return df, ts

df_lq, ts = remove_inactive_students(df_lq, ts)
# print(df_lq.columns)

display(df_lq)

	user	ch_num_sessions	ch_time_in_prob_sum	ch_time_in_video_sum	ch_ratio_clicks_weekend_day	ch_total_clicks_weekend	ch_total_clicks_weekday	ch_time_sessions_mean	ch_time_sessions_std	bo_delay_lecture	...	time_in_problem	time_in_video	lecture_delay	content_anticipation	mean_playback_speed	relative_video_pause	submissions	submissions_correct	clicks_weekend	clicks_weekday
1	1	3.4	1698.4	9227.8	0.567500	4.0	179.4	3068.720238	1257.504407	-2883.367738	...	1698.4	9227.8	-2883.367738	0.009091	1.122014	0.361389	9.0	5.9	4.0	179.4
2	2	5.3	2340.6	10801.3	26.562274	94.6	129.2	1750.289268	1024.134043	10027.216667	...	2340.6	10801.3	10027.216667	0.237488	0.807090	0.272210	6.1	3.0	94.6	129.2
4	4	2.5	3787.3	7040.0	1.543889	58.4	64.9	3373.908333	1363.320365	-914.633333	...	3787.3	7040.0	-914.633333	0.290421	0.846794	0.196403	0.3	0.1	58.4	64.9
5	5	2.5	2568.9	3718.9	0.009677	31.2	50.8	1753.647500	1190.793589	12406.195238	...	2568.9	3718.9	12406.195238	0.036241	1.090196	0.308608	3.8	2.4	31.2	50.8
6	6	4.2	5475.2	9711.6	0.476263	105.1	46.7	20410.677619	1561.548415	-13723.616667	...	5475.2	9711.6	-13723.616667	0.278603	0.932590	0.303720	1.9	1.2	105.1	46.7
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
272	281	2.2	2600.3	3071.8	0.213467	20.6	122.6	2429.026667	336.303590	-69623.272700	...	2600.3	3071.8	-69623.272700	0.026691	1.717685	0.086820	10.9	6.7	20.6	122.6
276	285	1.8	2111.4	2130.4	0.236626	7.0	21.2	1378.954762	544.967013	23125.010000	...	2111.4	2130.4	23125.010000	0.001667	0.839917	0.061321	5.1	2.6	7.0	21.2
277	286	1.9	2224.8	2408.1	0.326239	20.5	46.7	1744.663333	664.335107	4812.464286	...	2224.8	2408.1	4812.464286	0.023393	1.137915	0.171036	6.5	4.1	20.5	46.7
281	291	4.0	6898.3	5657.2	0.080266	131.5	15.7	2323.840952	1623.425778	-5626.100000	...	6898.3	5657.2	-5626.100000	0.157089	1.048739	0.192380	3.2	1.7	131.5	15.7
283	293	3.5	8127.5	113.4	0.632304	28.9	20.6	7963.627500	1001.514794	0.000000	...	8127.5	113.4	0.000000	0.010000	0.246698	0.065734	5.0	3.4	28.9	20.6

234 rows × 49 columns

display(ts)

	user	week	ch_num_sessions	ch_time_in_prob_sum	ch_time_in_video_sum	ch_total_clicks_weekend	ch_total_clicks_weekday	bo_delay_lecture	ma_competency_anti	ma_content_anti	mu_speed_playback_mean	mu_frequency_action_relative_video_pause	wa_num_subs	wa_num_subs_correct	wa_num_subs_perc_correct
1	1	0	7.0	326.0	15525.0	40.0	227.0	4492.833333	0.090909	0.090909	1.319288	0.345528	7.0	4.0	0.571429
2	2	0	4.0	1224.0	12209.0	1.0	258.0	-8998.000000	0.060606	0.060606	1.000000	0.230415	13.0	8.0	0.615385
4	4	0	4.0	1294.0	12037.0	140.0	46.0	-9146.333333	0.181818	0.181818	1.184140	0.267606	3.0	1.0	0.333333
5	5	0	2.0	1324.0	4440.0	119.0	0.0	0.000000	0.060606	0.060606	1.894958	0.407080	0.0	0.0	0.000000
6	6	0	3.0	1773.0	14462.0	102.0	144.0	-90548.166667	0.121212	0.121212	1.000000	0.413793	5.0	3.0	0.600000
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
2864	281	9	0.0	0.0	0.0	0.0	0.0	0.000000	0.000000	0.000000	0.000000	0.000000	0.0	0.0	0.000000
2868	285	9	0.0	0.0	0.0	0.0	0.0	0.000000	0.000000	0.000000	0.000000	0.000000	0.0	0.0	0.000000
2869	286	9	0.0	0.0	0.0	0.0	0.0	0.000000	0.000000	0.000000	0.000000	0.000000	0.0	0.0	0.000000
2873	291	9	0.0	0.0	0.0	0.0	0.0	0.000000	0.000000	0.000000	0.000000	0.000000	0.0	0.0	0.000000
2875	293	9	1.0	0.0	0.0	14.0	0.0	0.000000	0.000000	0.000000	0.000000	0.000000	0.0	0.0	0.000000

2340 rows × 15 columns

Prepare data for classification

Add a pass/fail label

# We first add a column to the dataframe containing the outcome variable
# compute pass/fail label
df_lq['passed'] = df_lq.grade >= 4
df_lq['passed'] = df_lq['passed'].astype(int)

Remove "bad" features and Split Data

# We then split the data in a train-test split (stratified by the outcome variable)
X = df_lq.drop(['user','grade', 'gender', 'category', 'year', 'passed'], axis=1)
y = df_lq['passed']
X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, random_state=0, stratify=y) # split train and validation data set

Print pass/fail proportions

# The class proportions in train and validation sets are the same, thanks to the stratification on y
print(y_train.value_counts(normalize=True))
print(y_val.value_counts(normalize=True))

passed
1    0.604278
0    0.395722
Name: proportion, dtype: float64
passed
1    0.595745
0    0.404255
Name: proportion, dtype: float64

Define Evaluation Metrics (will see later in the slides)

def compute_scores(clf, X_train, y_train, X_test, y_test, roundnum = 3):
    clf.fit(X_train, y_train)
    y_pred = clf.predict(X_test)
    accuracy =  balanced_accuracy_score(y_test, y_pred)
    
    y_pred_proba = clf.predict_proba(X_test)[:,1]
    auc = roc_auc_score(y_test, y_pred_proba)
    
    return round(accuracy,roundnum), round(auc,roundnum)

Section 2: Decision Trees

Compute a decision tree of max depth 2 over all the features

clf = tree.DecisionTreeClassifier(max_depth=2, random_state=0, criterion='entropy')
accuracy, auc = compute_scores(clf, X_train, y_train, X_val, y_val)
print("Decision tree. Balanced Accuracy = {}, AUC = {}".format(accuracy, auc))

Decision tree. Balanced Accuracy = 0.577, AUC = 0.602

Visualize the decision tree

plt.figure(figsize=(20, 10))
tree.plot_tree(clf, feature_names=X_train.columns);

Does depth improves perfromance ?

# We can change the max depth
accuracy_list = []
auc_list = []
for depth in range(1,len(X_train.columns)):
    clf = tree.DecisionTreeClassifier(max_depth=depth, random_state=0, criterion='entropy')
    accuracy, auc = compute_scores(clf, X_train, y_train, X_val, y_val)
    accuracy_list.append(accuracy)
    auc_list.append(auc)
    # print("Decision tree. Depth = {}, Balanced Accuracy = {}, AUC = {}".format(depth, accuracy, auc))
x = list(range(1,len(X_train.columns)))
plt.plot(x, accuracy_list, 'r', label = 'accuraزy')
plt.plot(x, auc_list, 'b', label = 'auc')
plt.xlabel("Decision tree Depth")
plt.ylabel("EvaluationMetrics")
plt.ylim([0,1])
plt.legend()
plt.show()

Section 3: Random Forests

Next, we will use a random forest classifier instead of a decision tree.

rf = RandomForestClassifier(n_estimators=100, random_state=0, criterion='entropy') # create a Random Forest
accuracy, auc = compute_scores(rf, X_train, y_train, X_val, y_val)
print("Random Forest. Balanced Accuracy = {}, AUC = {}".format(accuracy, auc))

Random Forest. Balanced Accuracy = 0.586, AUC = 0.585

For a single tree, in fact, keeping a low depth is necessary to avoid overfitting and to reduce the variance. Random forests, instead, can have a higher depth, and consequently a lower bias, since the variance is reduced in the aggregation step.

In this case, decision trees seem to perform better than random forests. A reason for this behavior could be that the single tree is already very "stable", i.e. it will change a little in response to little changes in the data. If this was the case, the submodels in the ensemble forest would be all very similar to the single tree, if they were allowed to choose among all the features at every split. Since, though, only a random subset of features is considered at each split, some subtrees would choose bad splits and have overall bad performances.

Section 4: K-Nearest Neighbors

We only use the euclidean distance since all our features are numerical

feature = 'ch_time_in_prob_sum'

# Compute the pairwise distance matrix for all the elements of the training set
X_train_dist = squareform(pdist(X_train[feature].to_numpy().reshape(-1,1), metric='euclidean'))

# Compute the distance between all elements of the training set and of the validation set
X_val_dist = cdist(X_val[feature].to_numpy().reshape(-1,1), X_train[feature].to_numpy().reshape(-1,1), metric='euclidean')

X_train_dist

array([[   0. ,  181.9,  932.5, ..., 1821.3, 1884.6, 2220.8],
       [ 181.9,    0. ,  750.6, ..., 1639.4, 2066.5, 2402.7],
       [ 932.5,  750.6,    0. , ...,  888.8, 2817.1, 3153.3],
       ...,
       [1821.3, 1639.4,  888.8, ...,    0. , 3705.9, 4042.1],
       [1884.6, 2066.5, 2817.1, ..., 3705.9,    0. ,  336.2],
       [2220.8, 2402.7, 3153.3, ..., 4042.1,  336.2,    0. ]])

print('Training set size:', X_train.shape)
print('Validation set size:', X_val.shape)
print('Training pairwise distances size:', X_train_dist.shape)
print('Validation distances size:', X_val_dist.shape)

Training set size: (187, 44)
Validation set size: (47, 44)
Training pairwise distances size: (187, 187)
Validation distances size: (47, 187)

knn = KNeighborsClassifier(n_neighbors=5, metric='precomputed')

accuracy, auc = compute_scores(knn, X_train_dist, y_train, X_val_dist, y_val)
print("k-nearest neighbors. Balanced Accuracy = {}, AUC = {}".format(accuracy, auc))

k-nearest neighbors. Balanced Accuracy = 0.533, AUC = 0.548

Section 5: Logistic regression

We normalize the data data using the MinMaxScaler such that all the features are on the same scale.

scaler = MinMaxScaler()
scaler.fit(X_train)

X_train_scaled = scaler.transform(X_train)
X_val_scaled  = scaler.transform(X_val)

clf = LogisticRegression(random_state=0)
accuracy, auc = compute_scores(clf, X_train_scaled, y_train, X_val_scaled, y_val)
print("Logistic Regression. Balanced Accuracy = {}, AUC = {}".format(accuracy, auc))

Logistic Regression. Balanced Accuracy = 0.559, AUC = 0.626

Section 6: Time Series - Your Turn

Build a classifier that can predict whether students pass the course after half of the course (5 weeks). You will need to use the data frame ts for this task. You can use kNN, RF, or decision tree. Train your model on the training data and predict on the test data.

We first drop all the weeks before week 5 and then perform the train/test split.

# Consider only data up to the 5th week
ts = ts[ts.week <= 5]

# Train-test split done on the users, so that all the rows corresponding to one user go into the same set.
users = ts.user.unique()
y = df_lq.passed
users_train, users_val, y_train, y_val = train_test_split(users, y, test_size=0.2, random_state=0, stratify=y)

X_train = ts[ts.user.isin(users_train)]
X_val = ts[ts.user.isin(users_val)]

# Sort indexes to make label arrays consistent with the data
y_train = y_train.sort_index()
y_val = y_val.sort_index()

Decision tree/Random forest

For decision tree and RF, we need to aggregate the features to be able to feed them directly in the mode. We use the mean as an aggregation function.

# Aggregate features

groups_train = X_train.drop('week', axis = 1).groupby('user', as_index = False)
groups_val = X_val.drop('week', axis = 1).groupby('user', as_index = False)

standard_train = groups_train.std()
averages_train = groups_train.mean()
standard_val = groups_val.std()
averages_val = groups_val.mean()

X_train_aggregate = pd.concat([standard_train, averages_train], axis=1)
X_val_aggregate = pd.concat([standard_val, averages_val], axis=1)
X_train_aggregate.head()

	user	ch_num_sessions	ch_time_in_prob_sum	ch_time_in_video_sum	ch_total_clicks_weekend	ch_total_clicks_weekday	bo_delay_lecture	ma_competency_anti	ma_content_anti	mu_speed_playback_mean	...	ch_total_clicks_weekend	ch_total_clicks_weekday	bo_delay_lecture	ma_competency_anti	ma_content_anti	mu_speed_playback_mean	mu_frequency_action_relative_video_pause	wa_num_subs	wa_num_subs_correct	wa_num_subs_perc_correct
0	2	3.391165	2109.754962	6234.491952	82.771976	103.319730	100718.818602	0.099941	0.099941	0.000000	...	120.000000	176.833333	24461.694444	0.140258	0.140258	1.000000	0.332633	10.000000	4.833333	0.441964
1	5	1.095445	2254.745433	2119.963773	57.595139	38.301001	42682.471658	0.074596	0.074596	0.452140	...	52.000000	25.833333	-114.222222	0.060401	0.060401	1.395326	0.379019	2.500000	1.333333	0.166667
2	6	1.505545	2872.580959	4712.799423	70.164806	57.094366	38052.014201	0.093191	0.093191	0.000000	...	126.500000	56.166667	-22872.694444	0.197672	0.197672	1.000000	0.380767	1.333333	0.833333	0.211111
3	7	1.760682	2224.075868	3908.428094	73.841497	78.863807	70834.903484	0.045196	0.045196	0.000000	...	133.833333	80.500000	-77474.511111	0.149894	0.149894	1.000000	0.380495	3.666667	1.500000	0.338384
4	8	1.366260	3821.751924	2660.061704	110.199214	182.871540	33286.469085	0.110765	0.110765	0.027152	...	234.333333	143.000000	3480.708333	0.210803	0.210803	1.011085	0.444130	2.666667	1.500000	0.291667

5 rows × 28 columns

clf = tree.DecisionTreeClassifier(max_depth=3, min_samples_leaf=5, min_samples_split=2, 
                                  random_state=0, criterion='entropy')

accuracy, auc = compute_scores(clf, X_train_aggregate, y_train, X_val_aggregate, y_val)
print("Decision tree. Balanced Accuracy = {}, AUC = {}".format(accuracy, auc))

Decision tree. Balanced Accuracy = 0.506, AUC = 0.536

rf = RandomForestClassifier(n_estimators=100, max_depth=5, random_state=0, criterion='entropy') # create a Random Forest

accuracy, auc = compute_scores(rf,X_train_aggregate, y_train, X_val_aggregate, y_val)
print("Random Forest. Balanced Accuracy = {}, AUC = {}".format(accuracy, auc))

Random Forest. Balanced Accuracy = 0.533, AUC = 0.609

K-Nearest Neighbors

For kNN, we don't need to aggregate the features. We directly compute the pairwise Euclidean distance between users, separately for each feature. We then normalize the pairwise distance matrices such that we can sum them up and feed them directly into the model.

feature = ['ch_time_in_prob_sum', 'ma_competency_anti','bo_delay_lecture']

X_train_vectors = X_train.groupby('user')[feature].agg(list)
X_val_vectors = X_val.groupby('user')[feature].agg(list)

X_train_vectors.head()

	ch_time_in_prob_sum	ma_competency_anti	bo_delay_lecture
user
2	[1224.0, 1310.0, 6185.0, 3195.0, 5545.0, 2515.0]	[0.0606060606060606, 0.0535714285714285, 0.12,...	[-8998.0, -22764.666666666668, 228700.0, -3375...
5	[1324.0, 266.0, 4700.0, 4014.0, 4881.0, 0.0]	[0.0606060606060606, 0.0, 0.0, 0.1351351351351...	[0.0, -67829.0, 67143.66666666667, 0.0, 0.0, 0.0]
6	[1773.0, 5152.0, 6475.0, 8562.0, 6753.0, 10076.0]	[0.1212121212121212, 0.1071428571428571, 0.22,...	[-90548.16666666669, -46688.0, 0.0, 0.0, 0.0, ...
7	[6410.0, 902.0, 5776.0, 4517.0, 1540.0, 3827.0]	[0.1212121212121212, 0.1071428571428571, 0.2, ...	[-96044.5, -100976.5, 0.0, -187836.66666666663...
8	[1527.0, 4861.0, 5211.0, 8722.0, 4512.0, 12462.0]	[0.1515151515151515, 0.1071428571428571, 0.12,...	[-26506.333333333332, 67210.33333333333, 0.0, ...

def normalize(distance_matrix):
    range_matrix = np.max(distance_matrix) - np.min(distance_matrix)
    normalized = (distance_matrix - np.min(distance_matrix)) / range_matrix
    return normalized

# Compute the pairwise distance matrix for all the elements of the training set, by computing
# the distances between the vectors, for each of the features selected, and summing up
X_train_dist = sum(map(lambda x: normalize(squareform(pdist(x[1].tolist(), metric='euclidean'))),
                       X_train_vectors.items()))

# Same thing but between all elements of the training set and of the validation set
X_val_dist = sum(map(lambda x: normalize(cdist(x[0][1].tolist(), x[1][1].tolist(), metric='euclidean')), 
                     zip(X_val_vectors.items(), X_train_vectors.items())))

X_train_dist

array([[0.        , 0.63479635, 0.79863403, ..., 0.92701674, 0.79609908,
        0.86466875],
       [0.63479635, 0.        , 0.75841214, ..., 0.57145625, 0.54598916,
        0.74351221],
       [0.79863403, 0.75841214, 0.        , ..., 0.97851606, 0.88172496,
        0.43477215],
       ...,
       [0.92701674, 0.57145625, 0.97851606, ..., 0.        , 0.34186977,
        0.81487542],
       [0.79609908, 0.54598916, 0.88172496, ..., 0.34186977, 0.        ,
        0.72684563],
       [0.86466875, 0.74351221, 0.43477215, ..., 0.81487542, 0.72684563,
        0.        ]])

knn = KNeighborsClassifier(n_neighbors=5, metric='precomputed')

accuracy, auc = compute_scores(knn, X_train_dist, y_train, X_val_dist, y_val)
print("k-nearest neighbors. Balanced Accuracy = {}, AUC = {}".format(accuracy, auc))

k-nearest neighbors. Balanced Accuracy = 0.523, AUC = 0.544