%matplotlib inline
# Principal package imports
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import numpy as np
import scipy as sc

# Scikit-learn package imports
from sklearn import feature_extraction, model_selection
from sklearn.metrics import root_mean_squared_error, roc_auc_score

# PyBKT package imports
import random
old_randint = random.randint
random.randint = lambda a, b: old_randint(a, int(b))
from pyBKT.models import Model
random.randint = old_randint

# Import the lmm model class
from pymer4.models import Lmer

DATA_DIR = "./../../data/"

assistments = pd.read_csv(DATA_DIR + 'assistments.csv', low_memory=False).dropna()
assistments.head()

print("Number of unique students in the dataset:", len(set(assistments['user_id'])))
print("Number of unique skills in the dataset:", len(set(assistments['skill_name'])))

Number of unique students in the dataset: 4151
Number of unique skills in the dataset: 110

def create_iterator(data):
    '''
    Create an iterator to split interactions in data into train and test, with the same student not appearing in two diverse folds.
    :param data:        Dataframe with student's interactions.
    :return:            An iterator.
    '''    
    # Both passing a matrix with the raw data or just an array of indexes works
    X = np.arange(len(data.index)) 
    # Groups of interactions are identified by the user id (we do not want the same user appearing in two folds)
    groups = data['user_id'].values 
    return model_selection.GroupShuffleSplit(n_splits=1, train_size=.8, test_size=0.2, random_state=0).split(X, groups=groups)

skills_subset = ['Circle Graph', 'Venn Diagram', 'Mode', 'Division Fractions', 'Finding Percents', 'Area Rectangle']
data = assistments[assistments['skill_name'].isin(skills_subset)]

print("Skill set:", set(data['skill_name']))
print("Number of unique students in the subset:", len(set(data['user_id'])))
print("Number of unique skills in the subset:", len(set(data['skill_name'])))

Skill set: {'Finding Percents', 'Mode', 'Circle Graph', 'Area Rectangle', 'Venn Diagram', 'Division Fractions'}
Number of unique students in the subset: 1527
Number of unique skills in the subset: 6

# Initialize the model
model = Model(seed=0)

# Fit the model on the entire dataset
%time model.fit(data=data)

predictions = model.predict(data=data)[['user_id', 'skill_name', 'correct', 'correct_predictions']]

# Rename the dataframe columns as per instructions
predictions.columns = ['user_id', 'skill_name', 'y_true', 'y_pred_bkt']

CPU times: user 60.8 ms, sys: 6.57 ms, total: 67.4 ms
Wall time: 77.3 ms

predictions.head()

def avg_y_by_x(x, y):
    '''
    Compute average learning curve and number of students over the number of opportunities. 
    x is the number of opportunities.
    y the success rates of the users (can be predicted success rate or true success rate).
    '''
    # Transform lists into arrays
    x = np.array(x)
    y = np.array(y)

    # Sort the integer id representing the number of opportunities in increasing order
    xs = sorted(list(set(x)))

    # Supporting lists to store the:
    # - xv: integer identifier of the number of opportunities
    # - yv: average value across students at that number of opportunities
    # - lcb and ucb: lower and upper confidence bound
    # - n_obs: number of observartions present at that number of opportunities (on per-skill plots, it is the #students)
    xv, yv, lcb, ucb, n_obs = [], [], [], [], []

    # For each integer identifier of the number of opportunities 0, ...
    for v in xs:
        ys = [y[i] for i, e in enumerate(x) if e == v] # We retrieve the values for that integer identifier
        if len(ys) > 0: 
            xv.append(v) # Append the integer identifier of the number of opportunities
            yv.append(sum(ys) / len(ys)) # Append the average value across students at that number of opportunities
            n_obs.append(len(ys)) # Append the number of observartions present at that number of opportunities

            
            # Prepare data for confidence interval computation
            unique, counts = np.unique(ys, return_counts=True)
            counts = dict(zip(unique, counts))

            if 0 not in counts:
                counts[0] = 0
            if 1 not in counts:
                counts[1] = 0

            # Calculate the 95% confidence intervals
            ci = sc.stats.beta.interval(0.95, 0.5 + counts[0], 0.5 + counts[1])
            lcb.append(ci[0])
            ucb.append(ci[1])

    return xv, yv, lcb, ucb, n_obs

def plot_learning_curve(skill_name):
    '''
    Plot learning curve using BKT model for skill `skill_name`. 
    '''
    preds = predictions[predictions['skill_name'] == skill_name] # Retrieve predictions for the current skill

    xp = []
    yp = {}
    for col in preds.columns: # For y_true and and y_pred_bkt columns, initialize an empty list for curve values
        if 'y_' in col:
            yp[col] = []

    for user_id in preds['user_id'].unique(): # For each user
        user_preds = preds[preds['user_id'] == user_id] # Retrieve the predictions on the current skill for this user 
        xp += list(np.arange(len(user_preds))) # The x-axis values go from 0 to |n_opportunities|-1
        for col in preds.columns: 
            if 'y_' in col: # For y_true and and y_pred_bkt columns
                yp[col] += user_preds[col].tolist() # The y-axis value is the success rate for this user at that opportunity

    fig, axs = plt.subplots(2, 1, gridspec_kw={'height_ratios': [3, 2]}) # Initialize the plotting figure

    lines = []
    for col in preds.columns:
        if 'y_' in col: # For y_true and and y_pred_bkt columns   
            x, y, lcb, ucb, n_obs = avg_y_by_x(xp, yp[col]) # Calculate mean and 95% confidence intervals for success rate 
            y = [1-v for v in y] # Transform success rate in error rate
            if col == 'y_true': # In case of ground-truth data, we also show the confidence intervals
                axs[0].fill_between(x, lcb, ucb, alpha=.1)
            model_line, = axs[0].plot(x, y, label=col) # Plot the curve
            lines.append(model_line) # Store the line to then set the legend    

    # Make decorations for the learning curve plot
    axs[0].set_title(skill_name)
    axs[0].legend(handles=lines)
    axs[0].set_ylabel('Error')
    axs[0].set_ylim(0, 1)
    axs[0].set_xlim(0, None)

    # Plot the number of observations per number of opportunities bars and make decorations
    axs[1].set_xlabel('#Opportunities')
    axs[1].bar([i for i in range(len(n_obs))], n_obs)
    axs[1].set_ylabel('#Observations')
    axs[1].set_ylim(0, 750)
    axs[1].set_xlim(0, None)

    # Plot the learning curve and the bar plot 
    plt.show()

for skill_name in skills_subset[:3]: # For the first three skills under consideration
    plot_learning_curve(skill_name)

for skill_name in skills_subset[3:]: # For the remaining three skills under consideration    
    plot_learning_curve(skill_name)

# Data processing
# Number of attempts before current
def preprocess_data(data):
    data['aux'] = 1
    data['prev_attempts'] = data.sort_values('order_id').groupby(['user_id', 'skill_name'])['aux'].cumsum() -1

    # Number of correct and incorrect attempts before current attempt
    data['correct_aux'] = data.sort_values('order_id').groupby(['user_id', 'skill_name'])['correct'].cumsum()
    data['before_correct_num'] = data.sort_values('order_id').groupby(['user_id', 'skill_name'])['correct_aux'].shift(periods=1, fill_value=0)
    data['before_wrong_num'] = data['prev_attempts'] - data['before_correct_num']
    return data

data = preprocess_data(data)
data.head()

# Obtain indexes
train_index, test_index = next(create_iterator(data))
# Split the data
X_train, X_test = data.iloc[train_index], data.iloc[test_index]

# Initialize and fit the model
model = Lmer("correct ~ (1|user_id) + (1|skill_name) + (0 + prev_attempts|skill_name)", data=X_train, family='binomial')
%time model.fit() 
# Compute predictions
X_test['afm_predictions'] = model.predict(data=X_test, verify_predictions=False)
X_test.head()

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: correct~(1|user_id)+(1|skill_name)+(0+prev_attempts|skill_name)

Family: binomial	 Inference: parametric

Number of observations: 40258	 Groups: {'user_id': 1221.0, 'skill_name': 6.0}

Log-likelihood: -16797.782 	 AIC: 33603.565

Random effects:

                       Name     Var     Std
user_id         (Intercept) 2.56000 1.60000
skill_name      (Intercept) 0.68300 0.82700
skill_name.1  prev_attempts 0.00500 0.06900

No random effect correlations specified

Fixed effects:

CPU times: user 17 s, sys: 4.89 s, total: 21.9 s
Wall time: 30.5 s

def plot_metrics(df_preds, model):

    # Compute RMSE and AUC per skill
    skill_rmse = df_preds.groupby(['skill_name']).apply(lambda df_preds: root_mean_squared_error(df_preds.correct, df_preds[model]), include_groups=False)
    print('RMSE', np.mean(skill_rmse), np.std(skill_rmse))

    skill_auc = df_preds.groupby(['skill_name']).apply(lambda df_preds: roc_auc_score(df_preds.correct, df_preds[model]), include_groups=False)
    print('AUC', np.mean(skill_auc), np.std(skill_auc))

    # Create overall RMSE and RMSE per skill data frames
    skills_all = ['Skills']*len(skills_subset)
    df_overall_rmse = pd.DataFrame(list(zip(skills_all, skill_rmse)),
                   columns =['x', 'RMSE'])
    df_skill_rmse = pd.DataFrame(list(zip(skills_subset, skill_rmse)),
                   columns =['x', 'RMSE'])

    # Create overall AUC and AUC per skill data frames
    df_overall_auc = pd.DataFrame(list(zip(skills_all, skill_auc)),
                   columns =['x', 'AUC'])
    df_skill_auc = pd.DataFrame(list(zip(skills_subset, skill_auc)),
                   columns =['x', 'AUC'])

    # Two bar plots for RMSE: first one with std (mean RMSE over all skills), then one bar plot with a bar for each specific skill
    fig, axes = plt.subplots(1, 2, figsize=(10,5))
    sns.barplot(ax=axes[0],x='x', y='RMSE', data=df_overall_rmse, estimator=np.mean, errorbar='sd')
    axes[0].set_title('Overall RMSE')
    axes[0].set_xlabel('')

    sns.barplot(ax=axes[1],x='x', y='RMSE', data=df_skill_rmse, estimator=np.mean, errorbar='sd')
    plt.xticks(rotation=90)
    axes[1].set_title('RMSE across skills')
    axes[1].set_xlabel('')

    plt.show()

    # Two bar plots for AUC: first one with std (mean AUC over all skills), then one bar plot with a bar for each specific skill
    fig, axes = plt.subplots(1, 2, figsize=(10,5))
    sns.barplot(ax=axes[0],x='x', y='AUC', data=df_overall_auc, estimator=np.mean, errorbar='sd')
    axes[0].set_title('Overall AUC')
    axes[0].set_xlabel('')

    sns.barplot(ax=axes[1],x='x', y='AUC', data=df_skill_auc, estimator=np.mean, errorbar='sd')
    plt.xticks(rotation=90)
    axes[1].set_title('AUC across skills')
    axes[1].set_xlabel('')

    plt.show()

plot_metrics(X_test, 'afm_predictions')

RMSE 0.4187859248639723 0.06896549818413097
AUC 0.6051593380658353 0.06187187963009105

# Initialize and fit the model
model = Lmer("correct ~ (1|user_id) + (1|skill_name) + (0 + before_correct_num|skill_name) + (0 + before_wrong_num|skill_name)", data=X_train, family='binomial')
%time model.fit() 
# Compute predictions
X_test['pfa_predictions'] = model.predict(data=X_test, verify_predictions=False)
X_test.head()

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: correct~(1|user_id)+(1|skill_name)+(0+before_correct_num|skill_name)+(0+before_wrong_num|skill_name)

Family: binomial	 Inference: parametric

Number of observations: 40258	 Groups: {'user_id': 1221.0, 'skill_name': 6.0}

Log-likelihood: -16385.969 	 AIC: 32781.939

Random effects:

                            Name     Var     Std
user_id              (Intercept) 1.74800 1.32200
skill_name           (Intercept) 0.69900 0.83600
skill_name.1  before_correct_num 0.02600 0.16200
skill_name.2    before_wrong_num 0.00000 0.01000

No random effect correlations specified

Fixed effects:

CPU times: user 31.3 s, sys: 4.24 s, total: 35.5 s
Wall time: 49.4 s

plot_metrics(X_test, 'pfa_predictions')

RMSE 0.40175375765684646 0.06679636660411878
AUC 0.7263985371288011 0.09235565673977136

df_preds = pd.DataFrame()

# Train a BKT model for each skill
for skill in skills_subset:
    print("--{}--".format(skill))
    X_train_skill = X_train[X_train['skill_name'] == skill]
    X_test_skill = X_test[X_test['skill_name'] == skill]
    # Initialize and fit the model
    model = Model(seed=0)
    %time model.fit(data=X_train_skill) 
    preds = model.predict(data=X_test_skill) [['user_id', 'order_id', 'skill_name', 'correct', 'prev_attempts',
       'before_correct_num', 'before_wrong_num', 'afm_predictions', 'pfa_predictions', 'correct_predictions']]
    df_preds = pd.concat([df_preds, preds], ignore_index=True)

X_test = df_preds
X_test.columns = ['user_id', 'order_id', 'skill_name', 'correct', 'prev_attempts',
       'before_correct_num', 'before_wrong_num', 'afm_predictions', 'pfa_predictions', 'bkt_predictions']
X_test = X_test.dropna(subset=['afm_predictions', 'pfa_predictions', 'bkt_predictions'])
X_test.head()

--Circle Graph--
CPU times: user 8.02 ms, sys: 2.58 ms, total: 10.6 ms
Wall time: 16.7 ms
--Venn Diagram--
CPU times: user 15.8 ms, sys: 2.95 ms, total: 18.8 ms
Wall time: 33.4 ms
--Mode--
CPU times: user 3.89 ms, sys: 68 μs, total: 3.96 ms
Wall time: 4.08 ms
--Division Fractions--
CPU times: user 5.35 ms, sys: 320 μs, total: 5.67 ms
Wall time: 6.64 ms
--Finding Percents--
CPU times: user 8.2 ms, sys: 596 μs, total: 8.8 ms
Wall time: 9.91 ms
--Area Rectangle--
CPU times: user 9.98 ms, sys: 587 μs, total: 10.6 ms
Wall time: 11.6 ms

rmse_bkt = root_mean_squared_error(X_test['bkt_predictions'],X_test['correct'])
rmse_afm = root_mean_squared_error(X_test['afm_predictions'],X_test['correct'])
rmse_pfa = root_mean_squared_error(X_test['pfa_predictions'],X_test['correct'])

rmse = [rmse_bkt, rmse_afm, rmse_pfa]
models = ['BKT', 'AFM', 'PFA']

plt.bar(models, rmse)
plt.ylabel('RMSE')

Text(0, 0.5, 'RMSE')

auc_bkt = roc_auc_score(X_test['correct'], X_test['bkt_predictions'])
auc_afm = roc_auc_score(X_test['correct'], X_test['afm_predictions'])
auc_pfa = roc_auc_score(X_test['correct'], X_test['pfa_predictions'])

auc = [auc_bkt, auc_afm, auc_pfa]
models = ['BKT', 'AFM', 'PFA']

plt.bar(models, auc)
plt.ylabel('AUC')

Text(0, 0.5, 'AUC')

	user_id	order_id	skill_name	correct
0	64525	33022537	Box and Whisker	1
1	64525	33022709	Box and Whisker	1
2	70363	35450204	Box and Whisker	0
3	70363	35450295	Box and Whisker	1
4	70363	35450311	Box and Whisker	0

	user_id	order_id	skill_name	correct	correct_predictions	state_predictions	aux	prev_attempts	correct_aux	before_correct_num	before_wrong_num	afm_predictions
53382	53167	26451283	Mode	1	0.50000	0.00105	1	0	1	0	0	0.85507
53383	53167	26451302	Mode	1	0.00000	0.00000	1	1	2	1	0	0.86928
157253	53167	32517498	Division Fractions	0	0.50000	0.00905	1	0	0	0	0	0.62272
157254	53167	32517627	Division Fractions	1	0.00000	0.00000	1	1	1	0	1	0.64312
157255	53167	32517648	Division Fractions	1	0.00000	0.00000	1	2	2	1	1	0.66302

	user_id	order_id	skill_name	correct	correct_predictions	state_predictions	aux	prev_attempts	correct_aux	before_correct_num	before_wrong_num	afm_predictions	pfa_predictions
53382	53167	26451283	Mode	1	0.50000	0.00105	1	0	1	0	0	0.85507	0.83721
53383	53167	26451302	Mode	1	0.00000	0.00000	1	1	2	1	0	0.86928	0.87290
157253	53167	32517498	Division Fractions	0	0.50000	0.00905	1	0	0	0	0	0.62272	0.61971
157254	53167	32517627	Division Fractions	1	0.00000	0.00000	1	1	1	0	1	0.64312	0.62047
157255	53167	32517648	Division Fractions	1	0.00000	0.00000	1	2	2	1	1	0.66302	0.65994

	user_id	order_id	skill_name	correct	prev_attempts	before_correct_num	before_wrong_num	afm_predictions	pfa_predictions	bkt_predictions
0	64525	28186893	Circle Graph	1	0	0	0	0.48266	0.46224	0.50000
1	64525	28187093	Circle Graph	1	1	1	0	0.49251	0.48999	0.00000
2	64525	32413158	Circle Graph	1	2	2	0	0.50236	0.51780	0.00000
3	64525	33022751	Circle Graph	0	3	3	0	0.51221	0.54551	0.00000
4	64525	33023039	Circle Graph	1	4	3	1	0.52205	0.54598	0.00000

Lecture 7 - Demo Notebook¶

BKT - Learning Curves¶

Additive Factors Model (AFM) and Performance Factors Analysis (PFA)¶

Model Comparison on Subset¶

	user_id	skill_name	y_true	y_pred_bkt
3957	14	Circle Graph	0	0.50000
3958	14	Circle Graph	1	0.00000
3959	14	Circle Graph	0	0.00000
3960	14	Circle Graph	0	0.00000
3961	14	Circle Graph	0	0.00000

	user_id	order_id	skill_name	correct	correct_predictions	state_predictions	aux	prev_attempts	correct_aux	before_correct_num	before_wrong_num
3957	14	21617623	Circle Graph	0	0.50000	0.01498	1	0	0	0	0
3958	14	21617632	Circle Graph	1	0.00000	0.00000	1	1	1	0	1
3959	14	21617641	Circle Graph	0	0.00000	0.00000	1	2	1	1	1
3960	14	21617650	Circle Graph	0	0.00000	0.00000	1	3	1	1	2
3961	14	21617659	Circle Graph	0	0.00000	0.00000	1	4	1	1	3