Lecture 03 - Demo Notebook¶

We recommend using Noto for this lecture tutorial, where we've already installed the dependencies of the pymer4 package and statsmodels.

We extended the data with extra features. The feature description is found here.

In [1]:
%matplotlib inline
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

from sklearn.model_selection import train_test_split
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import root_mean_squared_error

# Import the linear regression model class
from pymer4.models import Lm

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

# Import Gaussian modeling
import statsmodels.formula.api as smf

# Data directory
DATA_DIR = "./../../data"
In [2]:
# Parse the data
df = pd.read_csv('{}/aggregated_extended_fc.csv'.format(DATA_DIR))
df = df.fillna('NaN')
list(df.columns)
display(df)
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
0 0 1.9 2334.4 2951.8 0.850000 16.8 38.1 1392.858333 790.762032 55068.387500 ... 2334.4 2951.8 55068.387500 0.006061 0.968519 0.137436 3.0 2.0 16.8 38.1
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
3 3 2.8 2737.1 8185.5 3.691250 13.5 46.4 20203.590260 656.052901 27596.864484 ... 2737.1 8185.5 27596.864484 0.011879 0.500000 0.151223 4.6 3.2 13.5 46.4
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
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
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
284 294 2.2 2452.4 4623.1 18.147762 36.4 71.3 3614.055952 853.195566 16834.900000 ... 2452.4 4623.1 16834.900000 0.086765 0.800216 0.191907 3.3 1.4 36.4 71.3
285 296 0.9 1643.2 1932.4 0.000000 0.4 31.2 926.916667 616.918475 -12860.522222 ... 1643.2 1932.4 -12860.522222 0.121485 0.359663 0.028497 1.2 0.9 0.4 31.2
286 297 1.4 2718.6 360.3 0.180000 2.0 15.3 346.437500 122.017326 0.000000 ... 2718.6 360.3 0.000000 0.000000 0.100000 0.029412 0.7 0.6 2.0 15.3
287 298 0.9 0.1 1954.9 0.597368 15.9 3.3 350.758333 266.095738 0.000000 ... 0.1 1954.9 0.000000 0.000000 0.311397 0.087702 0.0 0.0 15.9 3.3

288 rows × 49 columns

Linear Regression¶

We will fit a simple linear regression with the pymer4 library with one variable.

In [3]:
# Initialize linear regression model using 1 predictor (time_in_prob) and sample data
model = Lm("grade ~ mu_speed_playback_mean", data=df)

# Fit it
print(model.fit())

Formula: grade~mu_speed_playback_mean

Family: gaussian	 Estimator: OLS

Std-errors: non-robust	CIs: standard 95%	Inference: parametric 

Number of observations: 288	 R^2: 0.019	 R^2_adj: 0.016

Log-likelihood: -463.538 	 AIC: 931.076	 BIC: 938.402

Fixed effects:

                        Estimate  2.5_ci  97.5_ci     SE   DF  T-stat  P-val  \
Intercept                  3.524   3.068    3.981  0.232  286  15.201  0.000   
mu_speed_playback_mean     0.555   0.095    1.015  0.234  286   2.377  0.018   

                        Sig  
Intercept               ***  
mu_speed_playback_mean    *

Let's visualize the fit of our linear regression model.

In [4]:
# fit a regression line
x_pred = df.mu_speed_playback_mean
y_pred = model.fits

plt.figure()
plt.scatter(df.mu_speed_playback_mean, df.grade)
plt.plot(x_pred,y_pred, color='red')
plt.xlabel("Mean playback speed")
plt.ylabel("Grade")
plt.show()
No description has been provided for this image

Two-Variable Regression¶

Next, we will do a regression with two variables. Which of these variables has the larger influence on grade?

In [5]:
# Linear regression with two variables

# Initialize model using 2 predictor (time_in_problem) and sample data
regression_two_variables = Lm("grade ~ ch_time_in_prob_sum + wa_num_subs_perc_correct", data=df)

# Fit regression model
regression_two_variables.fit()
print(regression_two_variables.coefs)

Formula: grade~ch_time_in_prob_sum+wa_num_subs_perc_correct

Family: gaussian	 Estimator: OLS

Std-errors: non-robust	CIs: standard 95%	Inference: parametric 

Number of observations: 288	 R^2: 0.110	 R^2_adj: 0.104

Log-likelihood: -449.516 	 AIC: 905.031	 BIC: 916.020

Fixed effects:

                          Estimate    2.5_ci   97.5_ci        SE   DF  \
Intercept                 3.410119  3.148091  3.672148  0.133123  285   
ch_time_in_prob_sum       0.000157  0.000094  0.000220  0.000032  285   
wa_num_subs_perc_correct  0.716132  0.035683  1.396581  0.345700  285   

                             T-stat     P-val  Sig  
Intercept                 25.616335  0.000000  ***  
ch_time_in_prob_sum        4.921856  0.000001  ***  
wa_num_subs_perc_correct   2.071542  0.039208    *

The two variables have very different scales: one is time in seconds {0, inf} and one is percentage {0, 1}. Therefore, we cannot directly compare them.

To make the coefficients comparable, we will standardize them by computing the z-score.

In [6]:
# compute z-score for time in problem and percentage correct
df['time_in_prob_z'] = (df.ch_time_in_prob_sum - df.ch_time_in_prob_sum.mean())/df.ch_time_in_prob_sum.std()
df['percentage_correct_z'] = (df.wa_num_subs_perc_correct - df.wa_num_subs_perc_correct.mean())/df.wa_num_subs_perc_correct.std()

# re-compute the regression
# Initialize model using 2 predictor (time_in_problem) and sample data
regression_comparable = Lm("grade ~ time_in_prob_z + percentage_correct_z", data=df)

# fit the regression model
regression_comparable.fit()
print(regression_comparable.coefs)

Formula: grade~time_in_prob_z+percentage_correct_z

Family: gaussian	 Estimator: OLS

Std-errors: non-robust	CIs: standard 95%	Inference: parametric 

Number of observations: 288	 R^2: 0.110	 R^2_adj: 0.104

Log-likelihood: -449.516 	 AIC: 905.031	 BIC: 916.020

Fixed effects:

                      Estimate    2.5_ci   97.5_ci        SE   DF     T-stat  \
Intercept             4.048611  3.914245  4.182978  0.068265  285  59.307699   
time_in_prob_z        0.346488  0.207923  0.485054  0.070398  285   4.921856   
percentage_correct_z  0.145832  0.007266  0.284398  0.070398  285   2.071542   

                         P-val  Sig  
Intercept             0.000000  ***  
time_in_prob_z        0.000001  ***  
percentage_correct_z  0.039208    *

We observe that the time in problem attribute has a larger impact on grade than percentage correct.

Another option is to use a MinMax Scaling, i.e. we are not standardizing the features but normalizing them between 0 and 1. We don't need to apply the scaler to the percentage of correct solution as this features is naturally between 0 and 1.

In [7]:
scaler = MinMaxScaler()
scaler.fit((df['ch_time_in_prob_sum']).to_numpy().reshape(-1,1))
df['time_in_problem_norm'] = scaler.transform((df['ch_time_in_prob_sum']).to_numpy().reshape(-1,1))
df.head()
Out[7]:
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 ... content_anticipation mean_playback_speed relative_video_pause submissions submissions_correct clicks_weekend clicks_weekday time_in_prob_z percentage_correct_z time_in_problem_norm
0 0 1.9 2334.4 2951.8 0.850000 16.8 38.1 1392.858333 790.762032 55068.387500 ... 0.006061 0.968519 0.137436 3.0 2.0 16.8 38.1 -0.212453 -0.354933 0.209856
1 1 3.4 1698.4 9227.8 0.567500 4.0 179.4 3068.720238 1257.504407 -2883.367738 ... 0.009091 1.122014 0.361389 9.0 5.9 4.0 179.4 -0.501291 1.835571 0.152682
2 2 5.3 2340.6 10801.3 26.562274 94.6 129.2 1750.289268 1024.134043 10027.216667 ... 0.237488 0.807090 0.272210 6.1 3.0 94.6 129.2 -0.209637 0.438664 0.210414
3 3 2.8 2737.1 8185.5 3.691250 13.5 46.4 20203.590260 656.052901 27596.864484 ... 0.011879 0.500000 0.151223 4.6 3.2 13.5 46.4 -0.029567 0.077082 0.246058
4 4 2.5 3787.3 7040.0 1.543889 58.4 64.9 3373.908333 1363.320365 -914.633333 ... 0.290421 0.846794 0.196403 0.3 0.1 58.4 64.9 0.447380 -1.190914 0.340468

5 rows × 52 columns

In [8]:
# re-compute the regression
# Initialize model using 2 predictor (time_in_problem) and sample data
regression_scaled = Lm("grade ~ time_in_problem_norm + wa_num_subs_perc_correct", data=df)

# fit the regression model
regression_scaled.fit()
print(regression_scaled.coefs)

Formula: grade~time_in_problem_norm+wa_num_subs_perc_correct

Family: gaussian	 Estimator: OLS

Std-errors: non-robust	CIs: standard 95%	Inference: parametric 

Number of observations: 288	 R^2: 0.110	 R^2_adj: 0.104

Log-likelihood: -449.516 	 AIC: 905.031	 BIC: 916.020

Fixed effects:

                          Estimate    2.5_ci   97.5_ci        SE   DF  \
Intercept                 3.410119  3.148091  3.672148  0.133123  285   
time_in_problem_norm      1.750411  1.050396  2.450426  0.355640  285   
wa_num_subs_perc_correct  0.716132  0.035683  1.396581  0.345700  285   

                             T-stat     P-val  Sig  
Intercept                 25.616335  0.000000  ***  
time_in_problem_norm       4.921856  0.000001  ***  
wa_num_subs_perc_correct   2.071542  0.039208    *

Generalized Linear Models¶

Now that we have successfully experimented with regression, we are interested in predicting whether a student will pass/fail a course.

In [9]:
import statsmodels.api as sm
import statsmodels.formula.api as smf

# compute pass/fail label
df['passed'] = df.grade >= 4
df['passed'] = df['passed'].astype(int)

# logistic regression
mod1 = smf.glm(formula='passed ~ wa_num_subs_perc_correct', data=df, family=sm.families.Binomial()).fit()
print(mod1.summary())

                 Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:                 passed   No. Observations:                  288
Model:                            GLM   Df Residuals:                      286
Model Family:                Binomial   Df Model:                            1
Link Function:                  Logit   Scale:                          1.0000
Method:                          IRLS   Log-Likelihood:                -193.28
Date:                Tue, 28 Jul 2026   Deviance:                       386.57
Time:                        16:21:31   Pearson chi2:                     288.
No. Iterations:                     4   Pseudo R-squ. (CS):           0.003213
Covariance Type:            nonrobust                                         
============================================================================================
                               coef    std err          z      P>|z|      [0.025      0.975]
--------------------------------------------------------------------------------------------
Intercept                    0.2520      0.201      1.251      0.211      -0.143       0.647
wa_num_subs_perc_correct     0.5718      0.596      0.960      0.337      -0.595       1.739
============================================================================================

Mixed Effect Models¶

Sometimes, we might deal with correlated samples. For example, our data set contains students from different origin. We might hypothesize that students having a similar background (category) behave more similar as they come from the same school/education system.

We will therefore now run a model with a random intercept, i.e. we will assume that the intercept depends on the category the students are in.

In [10]:
# Initialize model instance using 1 predictor with random intercepts and slopes
model = Lmer("passed ~ (1|category) + wa_num_subs_perc_correct", data=df, family='binomial')

# Fit it
print(model.fit())

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: passed~(1|category)+wa_num_subs_perc_correct

Family: binomial	 Inference: parametric

Number of observations: 288	 Groups: {'category': 5.0}

Log-likelihood: -191.991 	 AIC: 389.982

Random effects:

                 Name   Var    Std
category  (Intercept)  0.08  0.283

No random effect correlations specified

Fixed effects:

                          Estimate  2.5_ci  97.5_ci     SE     OR  OR_2.5_ci  \
(Intercept)                  0.184  -0.300    0.667  0.247  1.201      0.741   
wa_num_subs_perc_correct     0.661  -0.526    1.848  0.606  1.937      0.591   

                          OR_97.5_ci   Prob  Prob_2.5_ci  Prob_97.5_ci  \
(Intercept)                    1.949  0.546        0.426         0.661   
wa_num_subs_perc_correct       6.348  0.660        0.372         0.864   

                          Z-stat  P-val Sig  
(Intercept)                0.744  0.457      
wa_num_subs_perc_correct   1.092  0.275

Regression for Time-Series Data¶

Next, we analyze the data over time as we are dealing with time series interactions. First, we create a dataframe containing information about the user.

In [11]:
# parse the necessary data frames
df_ui = (df.loc[:,['user','grade','gender','category','year']]).copy()

# compute pass/fail label
df_ui['passed'] = (df_ui['grade'] >= 4).astype(int)
display(df_ui)
user grade gender category year passed
0 0 4.50 NaN NaN Y2-2018-19 1
1 1 4.50 M Suisse.Autres Y2-2018-19 1
2 2 5.25 M Suisse.PAM Y2-2018-19 1
3 3 4.50 F Suisse.Autres Y2-2018-19 1
4 4 4.75 F France Y2-2018-19 1
... ... ... ... ... ... ...
283 293 5.25 M France Y3-2019-20 1
284 294 5.25 F France Y3-2019-20 1
285 296 6.00 F France Y3-2019-20 1
286 297 5.00 M Suisse.PAM Y3-2019-20 1
287 298 4.00 M France Y3-2019-20 1

288 rows × 6 columns

Next, we parse the the data.

In [12]:
df_byweek = pd.read_csv('{}/fc_long_extended.csv'.format(DATA_DIR))
display(df_byweek)
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
0 0 0 4.0 5682.0 6417.0 12.0 102.0 -24339.200000 0.015152 0.015152 1.539474 0.315217 8.0 4.0 0.500000
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
3 3 0 11.0 3517.0 26500.0 10.0 141.0 -33102.111111 0.045455 0.045455 1.000000 0.301887 17.0 10.0 0.588235
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
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
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
2876 294 9 1.0 2263.0 0.0 0.0 8.0 0.000000 0.000000 0.000000 0.000000 0.000000 0.0 0.0 0.000000
2877 296 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
2878 297 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
2879 298 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

2880 rows × 15 columns

In [13]:
df_byuser = df_byweek.sort_values(by=['user', 'week']).reset_index(drop=True)

We can now run a model over the time series data of the students. In our first task, we are interested in predicting course grade early on during the semester. This type of information can be useful for an instructor in order to be able to provide intervention to struggling students. We will use again the category as a random effect. We will need to train a separate model for each week (i.e. predicting after 1 week of the course, after 2 weeks of the course, after 3 weeks, etc.). However, we will use the same equation for all models.

Step 1: We will write a function that aggregates the features for all weeks.

In [14]:
def aggregate_features(df_ui, df_byuser, week_nr):

    df_weeknr = df_byuser[df_byuser['week'] < week_nr]
    df_return = df_weeknr.groupby(['user']).mean()
    df_return['user'] = df_return.index
    
    # Return df with aggregated features
    df_return = df_return.set_index('user').join(df_ui.set_index('user'))
    df_return.reset_index()
    
    return df_return

Step 2: We will split the data into a training and test set (20% users in the test set, stratified by pass/fail label). In our case, data stratification refers to choosing a sample with the same ratio of pass/fail as the initial dataset, so our training set and our test set are both representative of our original population. If you are interested, you can read more about stratifying test sets here.

In [15]:
# perform train/test split
df_week5 = aggregate_features(df_ui, df_byuser, 5)
df_train5, df_test5 = train_test_split(df_week5, test_size=0.2, random_state=0,  stratify=df_week5['passed'])

df_week10 = aggregate_features(df_ui, df_byuser, 10)
df_train10, df_test10 = train_test_split(df_week10, test_size=0.2, random_state=0,  stratify=df_week10['passed'])

Step 3: We will now train our model on the training data for 5 and 10 weeks.

In [16]:
# Train a multi-regression model for weeks 5 and 10
# Initialize model instance using 1 predictor with random intercepts and slopes
model5 = Lmer("passed ~ (1|category) +  wa_num_subs_perc_correct", data=df_train5, family='binomial')
model10 = Lmer("passed ~ (1|category) +  wa_num_subs_perc_correct", data=df_train10, family='binomial')

# Fit the models
print(model5.fit())
print(model10.fit())

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: passed~(1|category)+wa_num_subs_perc_correct

Family: binomial	 Inference: parametric

Number of observations: 230	 Groups: {'category': 5.0}

Log-likelihood: -152.134 	 AIC: 310.269

Random effects:

                 Name    Var    Std
category  (Intercept)  0.157  0.396

No random effect correlations specified

Fixed effects:

                          Estimate  2.5_ci  97.5_ci     SE     OR  OR_2.5_ci  \
(Intercept)                  0.172  -0.428    0.772  0.306  1.188      0.652   
wa_num_subs_perc_correct     0.463  -0.694    1.620  0.590  1.589      0.500   

                          OR_97.5_ci   Prob  Prob_2.5_ci  Prob_97.5_ci  \
(Intercept)                    2.164  0.543        0.395         0.684   
wa_num_subs_perc_correct       5.052  0.614        0.333         0.835   

                          Z-stat  P-val Sig  
(Intercept)                0.561  0.575      
wa_num_subs_perc_correct   0.784  0.433      
Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: passed~(1|category)+wa_num_subs_perc_correct

Family: binomial	 Inference: parametric

Number of observations: 230	 Groups: {'category': 5.0}

Log-likelihood: -151.683 	 AIC: 309.366

Random effects:

                 Name    Var    Std
category  (Intercept)  0.162  0.402

No random effect correlations specified

Fixed effects:

                          Estimate  2.5_ci  97.5_ci     SE     OR  OR_2.5_ci  \
(Intercept)                  0.091  -0.506    0.688  0.305  1.095      0.603   
wa_num_subs_perc_correct     0.880  -0.527    2.286  0.718  2.410      0.591   

                          OR_97.5_ci   Prob  Prob_2.5_ci  Prob_97.5_ci  \
(Intercept)                    1.989  0.523        0.376         0.665   
wa_num_subs_perc_correct       9.835  0.707        0.371         0.908   

                          Z-stat  P-val Sig  
(Intercept)                0.298  0.766      
wa_num_subs_perc_correct   1.226  0.220

Step 4: We predict on the test data and check the Accuracy

In [17]:
# predict on the test data for weeks 5, 10
predictions5 = model5.predict(df_test5, verify_predictions=False)
rmse5 = root_mean_squared_error(df_test5['passed'], predictions5)

predictions10 = model10.predict(df_test10, verify_predictions=False)
rmse10 = root_mean_squared_error(df_test10['passed'], predictions10)

print(rmse5)
print(rmse10)

0.5025608612692153
0.5024946963345416