ProTuS Analysis

1. Packages and data import

# Install packages if needed:
# install.packages("readxl")
# install.packages("dplyr")
# install.packages("tidyr")
# install.packages("ggplot2")
# install.packages("lme4")

library(readxl)
library(dplyr)
library(tidyr)
library(ggplot2)
library(lme4)

# Change file path if necessary
file <- "/Users/mette/Desktop/NTNU/Learning analytics/DataSet-ProTuS.xlsx"

learners <- read_excel(file, sheet = "Learners") %>%
  select(where(~ !all(is.na(.))))

actions_raw <- read_excel(file, sheet = "Users actions") %>%
  select(where(~ !all(is.na(.))))

action_types_raw <- read_excel(file, sheet = "Action type") %>%
  select(where(~ !all(is.na(.))))

2. Initial data inspection

dim(learners)

[1] 116   5

dim(actions_raw)

[1] 13247    13

names(learners)

[1] "id"             "year"           "masteryID"      "ELOrating"      "iteraction_num"

names(actions_raw)

 [1] "idUsers_action"        "learner"               "lesson"                "type"                  "menu_item"             "learning_content_name" "category"              "Elo_rating_student"    "Elo_rating_content"   
[10] "session"               "changed_item"          "learning_content_type" "time"                 

nrow(learners)

[1] 116

n_distinct(actions_raw$learner)

[1] 81

sum(duplicated(actions_raw))

[1] 0

colSums(is.na(actions_raw))

       idUsers_action               learner                lesson                  type             menu_item learning_content_name              category    Elo_rating_student    Elo_rating_content               session 
                    0                     0                     0                     1                     1                     1                     1                     1                     1                     1 
         changed_item learning_content_type                  time 
                    1                     1                     1 

The Learners sheet contains 116 learners. Not all learners are represented in the provided action log.

learner_action_coverage <- learners %>%
  mutate(has_action_records = id %in% actions_raw$learner) %>%
  count(has_action_records)

print(learner_action_coverage)

# A tibble: 2 × 2
  has_action_records     n
  <lgl>              <int>
1 FALSE                 35
2 TRUE                  81

3. Data cleaning

Create a clean lookup table for action types.

action_lookup <- action_types_raw %>%
  select(idaction_type, action_type) %>%
  filter(!is.na(idaction_type)) %>%
  distinct(idaction_type, .keep_all = TRUE)

# Verify that each action type ID occurs only once
action_lookup %>%
  count(idaction_type) %>%
  filter(n > 1)

Remove the single malformed action record with a missing action type and add the action labels.

actions_clean <- actions_raw %>%
  filter(!is.na(type)) %>%
  left_join(
    action_lookup,
    by = c("type" = "idaction_type"),
    relationship = "many-to-one"
  ) %>%
  arrange(learner, time, idUsers_action)

Verify the cleaned dataset.

data_summary <- tibble(
  measure = c(
    "Learners in Learners sheet",
    "Learners represented in action log",
    "Raw action records",
    "Clean action records"
  ),
  n = c(
    nrow(learners),
    n_distinct(actions_raw$learner),
    nrow(actions_raw),
    nrow(actions_clean)
  )
)

print(data_summary)

# A tibble: 4 × 2
  measure                                n
  <chr>                              <int>
1 Learners in Learners sheet           116
2 Learners represented in action log    81
3 Raw action records                 13247
4 Clean action records               13246

# Verify that action IDs are unique after cleaning
stopifnot(sum(duplicated(actions_clean$idUsers_action)) == 0)

4. Elo events and exercise data

Check which action types contain non-zero Elo values.

elo_events <- actions_clean %>%
  filter(Elo_rating_student != 0)

elo_events %>%
  count(type, action_type, sort = TRUE) %>%
  print(n = Inf)

# A tibble: 1 × 3
   type action_type          n
  <dbl> <chr>            <int>
1    24 exercise started  2269

nrow(elo_events)

[1] 2269

n_distinct(elo_events$learner)

[1] 62

summary(elo_events$Elo_rating_student)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1232    1332    1368    1394    1455    1672 

All observed non-zero Elo values occur on coding exercise events (type 24).

Construct the exercise-level data. Elo change is calculated between consecutive observed exercise events within each learner.

exercise_data <- actions_clean %>%
  filter(type == 24) %>%
  arrange(learner, time, idUsers_action) %>%
  group_by(learner) %>%
  mutate(
    exercise_number = row_number(),
    previous_Elo = lag(Elo_rating_student),
    Elo_gain = Elo_rating_student - previous_Elo
  ) %>%
  ungroup()

The first observed exercise event for each learner has no preceding observed Elo value and therefore has no calculable Elo change.

exercise_data %>%
  summarise(
    n_exercises = n(),
    n_learners = n_distinct(learner),
    n_Elo_transitions = sum(!is.na(Elo_gain)),
    mean_gain = mean(Elo_gain, na.rm = TRUE),
    median_gain = median(Elo_gain, na.rm = TRUE),
    sd_gain = sd(Elo_gain, na.rm = TRUE),
    min_gain = min(Elo_gain, na.rm = TRUE),
    max_gain = max(Elo_gain, na.rm = TRUE),
    positive = sum(Elo_gain > 0, na.rm = TRUE),
    zero = sum(Elo_gain == 0, na.rm = TRUE),
    negative = sum(Elo_gain < 0, na.rm = TRUE)
  )

5. Construct preparatory engagement variables

The relevant action types are:

• 22 = coding example
• 23 = coding challenge
• 24 = coding exercise

Preparatory engagement is defined as example and challenge engagement occurring after the preceding observed exercise and before the current exercise.

actions_sequence <- actions_clean %>%
  arrange(learner, time, idUsers_action) %>%
  group_by(learner) %>%
  mutate(
    exercises_so_far = cumsum(type == 24),
    upcoming_exercise = if_else(
      type == 24,
      exercises_so_far,
      exercises_so_far + 1
    )
  ) %>%
  ungroup()

Count example and challenge engagement events before each exercise.

Repeated engagement with the same content is retained as additional observed engagement.

preparation <- actions_sequence %>%
  filter(type %in% c(22, 23)) %>%
  group_by(learner, upcoming_exercise) %>%
  summarise(
    examples_before = sum(type == 22),
    challenges_before = sum(type == 23),
    .groups = "drop"
  ) %>%
  mutate(
    preparatory_before = examples_before + challenges_before
  )

Join the preparatory engagement variables to the exercise-level data. Exercises without observed example or challenge engagement in the preceding interval receive a count of zero.

exercise_data <- exercise_data %>%
  left_join(
    preparation,
    by = c(
      "learner",
      "exercise_number" = "upcoming_exercise"
    )
  ) %>%
  mutate(
    examples_before = replace_na(examples_before, 0),
    challenges_before = replace_na(challenges_before, 0),
    preparatory_before = replace_na(preparatory_before, 0)
  )

6. Final analysis dataset

Exclude the first observed exercise for each learner because it has no preceding observed exercise from which an Elo change can be calculated.

Create binary variables used in the robustness analyses.

analysis_data <- exercise_data %>%
  filter(!is.na(Elo_gain)) %>%
  mutate(
    any_examples = examples_before > 0,
    any_preparation = preparatory_before > 0,
    Elo_increased = Elo_gain > 0
  )

Check the final analytical sample.

analysis_sample <- analysis_data %>%
  summarise(
    exercise_transitions = n(),
    learners = n_distinct(learner)
  )

print(analysis_sample)

# A tibble: 1 × 2
  exercise_transitions learners
                 <int>    <int>
1                 2207       57

Five of the 62 learners with observed exercise events have only one observed exercise and therefore cannot contribute to the analysis of between-exercise Elo changes.

exercise_counts <- exercise_data %>%
  group_by(learner) %>%
  summarise(
    n_exercises = n(),
    n_Elo_transitions = sum(!is.na(Elo_gain)),
    .groups = "drop"
  ) %>%
  arrange(n_exercises)

exercise_counts %>%
  filter(n_exercises == 1)

7. Descriptive statistics

descriptive_statistics <- analysis_data %>%
  summarise(
    n_observations = n(),
    n_learners = n_distinct(learner),

    mean_Elo_gain = mean(Elo_gain),
    sd_Elo_gain = sd(Elo_gain),
    median_Elo_gain = median(Elo_gain),
    min_Elo_gain = min(Elo_gain),
    max_Elo_gain = max(Elo_gain),

    mean_examples = mean(examples_before),
    median_examples = median(examples_before),
    max_examples = max(examples_before),

    mean_challenges = mean(challenges_before),
    median_challenges = median(challenges_before),
    max_challenges = max(challenges_before),

    mean_preparatory = mean(preparatory_before),
    median_preparatory = median(preparatory_before),
    max_preparatory = max(preparatory_before),

    pct_examples_before = mean(examples_before > 0) * 100,
    pct_challenges_before = mean(challenges_before > 0) * 100,
    pct_any_preparation = mean(preparatory_before > 0) * 100
  )

print(descriptive_statistics)

# A tibble: 1 × 19
  n_observations n_learners mean_Elo_gain sd_Elo_gain median_Elo_gain min_Elo_gain max_Elo_gain mean_examples median_examples max_examples mean_challenges median_challenges max_challenges mean_preparatory median_preparatory
           <int>      <int>         <dbl>       <dbl>           <dbl>        <dbl>        <dbl>         <dbl>           <int>        <int>           <dbl>             <int>          <int>            <dbl>              <int>
1           2207         57          2.11        9.65               0          -62          156         0.225               0           25          0.0435                 0             19            0.269                  0
# ℹ 4 more variables: max_preparatory <int>, pct_examples_before <dbl>, pct_challenges_before <dbl>, pct_any_preparation <dbl>

Descriptive statistics by type of preparatory engagement.

preparation_summary <- analysis_data %>%
  mutate(
    preparation_group = case_when(
      examples_before == 0 & challenges_before == 0 ~ "None",
      examples_before > 0 & challenges_before == 0 ~ "Examples only",
      examples_before == 0 & challenges_before > 0 ~ "Challenges only",
      examples_before > 0 & challenges_before > 0 ~ "Examples + challenges"
    ),
    preparation_group = factor(
      preparation_group,
      levels = c(
        "None",
        "Examples only",
        "Challenges only",
        "Examples + challenges"
      )
    )
  ) %>%
  group_by(preparation_group) %>%
  summarise(
    n = n(),
    mean_Elo_gain = mean(Elo_gain),
    median_Elo_gain = median(Elo_gain),
    sd_Elo_gain = sd(Elo_gain),
    .groups = "drop"
  )

print(preparation_summary)

# A tibble: 4 × 5
  preparation_group         n mean_Elo_gain median_Elo_gain sd_Elo_gain
  <fct>                 <int>         <dbl>           <dbl>       <dbl>
1 None                   2043          1.71               0        8.76
2 Examples only           127          5.77               0        9.95
3 Challenges only          11         17.5                0       46.7 
4 Examples + challenges    26          9                  0       18.5 

8. Primary inferential analysis

Because each learner contributes multiple exercise transitions, linear mixed-effects models are used with a random intercept for learner.

Exercise number is included as a control for practice volume.

H1a: Example engagement

model_H1a <- lmer(
  Elo_gain ~ examples_before + exercise_number + (1 | learner),
  data = analysis_data,
  REML = FALSE
)

summary(model_H1a)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: Elo_gain ~ examples_before + exercise_number + (1 | learner)
   Data: analysis_data

      AIC       BIC    logLik -2*log(L)  df.resid 
  16215.0   16243.5   -8102.5   16205.0      2202 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-6.6634 -0.2307 -0.1696 -0.1127 15.9841 

Random effects:
 Groups   Name        Variance Std.Dev.
 learner  (Intercept)  3.225   1.796   
 Residual             88.840   9.426   
Number of obs: 2207, groups:  learner, 57

Fixed effects:
                 Estimate Std. Error t value
(Intercept)      2.507127   0.409561   6.121
examples_before  1.153487   0.173649   6.643
exercise_number -0.006196   0.005340  -1.160

Correlation of Fixed Effects:
            (Intr) exmpl_
exampls_bfr -0.157       
exercs_nmbr -0.441  0.077

H1b: Preparatory engagement

model_H1b <- lmer(
  Elo_gain ~ preparatory_before + exercise_number + (1 | learner),
  data = analysis_data,
  REML = FALSE
)

summary(model_H1b)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: Elo_gain ~ preparatory_before + exercise_number + (1 | learner)
   Data: analysis_data

      AIC       BIC    logLik -2*log(L)  df.resid 
  16194.1   16222.6   -8092.1   16184.1      2202 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-6.6916 -0.2305 -0.1698 -0.1128 15.8619 

Random effects:
 Groups   Name        Variance Std.Dev.
 learner  (Intercept)  2.572   1.604   
 Residual             88.207   9.392   
Number of obs: 2207, groups:  learner, 57

Fixed effects:
                    Estimate Std. Error t value
(Intercept)         2.415102   0.390273   6.188
preparatory_before  1.101194   0.135372   8.135
exercise_number    -0.005961   0.005248  -1.136

Correlation of Fixed Effects:
            (Intr) prprt_
prprtry_bfr -0.156       
exercs_nmbr -0.471  0.079

9. Model diagnostics

Residual diagnostics are examined for the primary mixed-effects model.

par(mfrow = c(1, 2))

plot(
  fitted(model_H1b),
  resid(model_H1b),
  xlab = "Fitted values",
  ylab = "Residuals",
  main = "Residuals vs fitted"
)

abline(h = 0, lty = 2)

qqnorm(
  resid(model_H1b),
  main = "Normal Q-Q plot"
)

qqline(resid(model_H1b))

par(mfrow = c(1, 1))

The residual diagnostics indicate substantial deviation from normality. The Elo-change distribution also contains a large proportion of zero changes. Logistic mixed-effects models are therefore used below as a complementary robustness analysis.

10. Sensitivity analysis

The largest observed preparatory-engagement count is 44 events. A sensitivity analysis excludes this single extreme observation to assess its influence on the continuous mixed-model estimates.

analysis_sensitivity <- analysis_data %>%
  filter(preparatory_before < 44)

model_H1a_sensitivity <- lmer(
  Elo_gain ~ examples_before + exercise_number + (1 | learner),
  data = analysis_sensitivity,
  REML = FALSE
)

model_H1b_sensitivity <- lmer(
  Elo_gain ~ preparatory_before + exercise_number + (1 | learner),
  data = analysis_sensitivity,
  REML = FALSE
)

summary(model_H1a_sensitivity)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: Elo_gain ~ examples_before + exercise_number + (1 | learner)
   Data: analysis_sensitivity

      AIC       BIC    logLik -2*log(L)  df.resid 
  16177.2   16205.7   -8083.6   16167.2      2201 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-6.7118 -0.2449 -0.1847 -0.1219 16.1169 

Random effects:
 Groups   Name        Variance Std.Dev.
 learner  (Intercept)  2.508   1.584   
 Residual             87.845   9.373   
Number of obs: 2206, groups:  learner, 57

Fixed effects:
                 Estimate Std. Error t value
(Intercept)      2.598682   0.388928   6.682
examples_before  0.674130   0.192743   3.498
exercise_number -0.007559   0.005233  -1.445

Correlation of Fixed Effects:
            (Intr) exmpl_
exampls_bfr -0.170       
exercs_nmbr -0.474  0.085

summary(model_H1b_sensitivity)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: Elo_gain ~ preparatory_before + exercise_number + (1 | learner)
   Data: analysis_sensitivity

      AIC       BIC    logLik -2*log(L)  df.resid 
  16173.1   16201.5   -8081.5   16163.1      2201 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-6.7167 -0.2427 -0.1822 -0.1216 15.9958 

Random effects:
 Groups   Name        Variance Std.Dev.
 learner  (Intercept)  2.375   1.541   
 Residual             87.724   9.366   
Number of obs: 2206, groups:  learner, 57

Fixed effects:
                    Estimate Std. Error t value
(Intercept)         2.543141   0.385531   6.596
preparatory_before  0.696274   0.171450   4.061
exercise_number    -0.007171   0.005216  -1.375

Correlation of Fixed Effects:
            (Intr) prprt_
prprtry_bfr -0.180       
exercs_nmbr -0.483  0.096

Compare the primary and sensitivity estimates.

fixef(model_H1a)

    (Intercept) examples_before exercise_number 
    2.507127432     1.153486703    -0.006195978 

fixef(model_H1a_sensitivity)

    (Intercept) examples_before exercise_number 
    2.598682080     0.674129935    -0.007559097 

fixef(model_H1b)

       (Intercept) preparatory_before    exercise_number 
       2.415101609        1.101193953       -0.005960955 

fixef(model_H1b_sensitivity)

       (Intercept) preparatory_before    exercise_number 
       2.543141007        0.696273513       -0.007170756 

11. Logistic robustness analysis

The logistic robustness analysis examines whether engagement is associated with the probability of a positive Elo change.

Here, Elo_increased = TRUE when Elo_gain > 0. Zero and negative changes are therefore classified as no positive increase.

Event-count models

model_H1a_logistic <- glmer(
  Elo_increased ~ examples_before + exercise_number + (1 | learner),
  data = analysis_data,
  family = binomial
)

model_H1b_logistic <- glmer(
  Elo_increased ~ preparatory_before + exercise_number + (1 | learner),
  data = analysis_data,
  family = binomial
)

summary(model_H1a_logistic)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: Elo_increased ~ examples_before + exercise_number + (1 | learner)
   Data: analysis_data

      AIC       BIC    logLik -2*log(L)  df.resid 
   1582.3    1605.1    -787.2    1574.3      2203 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.5193 -0.3867 -0.3129 -0.2578  5.0246 

Random effects:
 Groups  Name        Variance Std.Dev.
 learner (Intercept) 0.4086   0.6392  
Number of obs: 2207, groups:  learner, 57

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -2.004839   0.144031 -13.919  < 2e-16 ***
examples_before  0.209197   0.047729   4.383 1.17e-05 ***
exercise_number -0.002765   0.001929  -1.434    0.152    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) exmpl_
exampls_bfr -0.200       
exercs_nmbr -0.430  0.099

summary(model_H1b_logistic)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: Elo_increased ~ preparatory_before + exercise_number + (1 | learner)
   Data: analysis_data

      AIC       BIC    logLik -2*log(L)  df.resid 
   1579.9    1602.7    -786.0    1571.9      2203 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.2700 -0.3859 -0.3130 -0.2587  5.0790 

Random effects:
 Groups  Name        Variance Std.Dev.
 learner (Intercept) 0.3986   0.6313  
Number of obs: 2207, groups:  learner, 57

Fixed effects:
                    Estimate Std. Error z value Pr(>|z|)    
(Intercept)        -2.016464   0.143615 -14.041  < 2e-16 ***
preparatory_before  0.196596   0.042392   4.638 3.52e-06 ***
exercise_number    -0.002578   0.001928  -1.337    0.181    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) prprt_
prprtry_bfr -0.212       
exercs_nmbr -0.437  0.122

Any engagement versus no engagement

model_H1a_logistic_any <- glmer(
  Elo_increased ~ any_examples + exercise_number + (1 | learner),
  data = analysis_data,
  family = binomial
)

model_H1b_logistic_any <- glmer(
  Elo_increased ~ any_preparation + exercise_number + (1 | learner),
  data = analysis_data,
  family = binomial
)

summary(model_H1a_logistic_any)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: Elo_increased ~ any_examples + exercise_number + (1 | learner)
   Data: analysis_data

      AIC       BIC    logLik -2*log(L)  df.resid 
   1564.3    1587.1    -778.2    1556.3      2203 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.2968 -0.3799 -0.2985 -0.2489  4.9066 

Random effects:
 Groups  Name        Variance Std.Dev.
 learner (Intercept) 0.4212   0.649   
Number of obs: 2207, groups:  learner, 57

Fixed effects:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -2.113753   0.148887 -14.197  < 2e-16 ***
any_examplesTRUE  1.315045   0.203150   6.473 9.59e-11 ***
exercise_number  -0.002168   0.001934  -1.121    0.262    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) a_TRUE
any_xmpTRUE -0.288       
exercs_nmbr -0.430  0.126

summary(model_H1b_logistic_any)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: Elo_increased ~ any_preparation + exercise_number + (1 | learner)
   Data: analysis_data

      AIC       BIC    logLik -2*log(L)  df.resid 
   1559.8    1582.6    -775.9    1551.8      2203 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.3139 -0.3796 -0.2960 -0.2475  4.9103 

Random effects:
 Groups  Name        Variance Std.Dev.
 learner (Intercept) 0.4215   0.6492  
Number of obs: 2207, groups:  learner, 57

Fixed effects:
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)         -2.140276   0.149844 -14.283  < 2e-16 ***
any_preparationTRUE  1.354819   0.197975   6.843 7.74e-12 ***
exercise_number     -0.001952   0.001933  -1.010    0.313    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) a_TRUE
any_prpTRUE -0.305       
exercs_nmbr -0.433  0.139

12. Tables for the paper

Table 1: Analytical sample and Elo changes

table1 <- analysis_data %>%
  summarise(
    `Exercise transitions` = n(),
    `Learners` = n_distinct(learner),
    `Mean Elo change` = round(mean(Elo_gain), 2),
    `SD Elo change` = round(sd(Elo_gain), 2),
    `Median Elo change` = median(Elo_gain),
    `Minimum` = min(Elo_gain),
    `Maximum` = max(Elo_gain),
    `Positive, n (%)` = paste0(
      sum(Elo_gain > 0),
      " (",
      round(mean(Elo_gain > 0) * 100, 1),
      "%)"
    ),
    `Zero, n (%)` = paste0(
      sum(Elo_gain == 0),
      " (",
      round(mean(Elo_gain == 0) * 100, 1),
      "%)"
    ),
    `Negative, n (%)` = paste0(
      sum(Elo_gain < 0),
      " (",
      round(mean(Elo_gain < 0) * 100, 1),
      "%)"
    )
  )

print(table1)

# A tibble: 1 × 10
  `Exercise transitions` Learners `Mean Elo change` `SD Elo change` `Median Elo change` Minimum Maximum `Positive, n (%)` `Zero, n (%)` `Negative, n (%)`
                   <int>    <int>             <dbl>           <dbl>               <dbl>   <dbl>   <dbl> <chr>             <chr>         <chr>            
1                   2207       57              2.11            9.65                   0     -62     156 272 (12.3%)       1852 (83.9%)  83 (3.8%)        

Table 2: Elo change by preparatory engagement

table2 <- analysis_data %>%
  mutate(
    preparation_group = case_when(
      examples_before == 0 & challenges_before == 0 ~ "None",
      examples_before > 0 & challenges_before == 0 ~ "Examples only",
      examples_before == 0 & challenges_before > 0 ~ "Challenges only",
      examples_before > 0 & challenges_before > 0 ~ "Examples + challenges"
    ),
    preparation_group = factor(
      preparation_group,
      levels = c(
        "None",
        "Examples only",
        "Challenges only",
        "Examples + challenges"
      )
    )
  ) %>%
  group_by(preparation_group) %>%
  summarise(
    n = n(),
    mean_Elo_change = round(mean(Elo_gain), 2),
    SD = round(sd(Elo_gain), 2),
    median_Elo_change = median(Elo_gain),
    .groups = "drop"
  )

print(table2)

# A tibble: 4 × 5
  preparation_group         n mean_Elo_change    SD median_Elo_change
  <fct>                 <int>           <dbl> <dbl>             <dbl>
1 None                   2043            1.71  8.76                 0
2 Examples only           127            5.77  9.95                 0
3 Challenges only          11           17.4  46.7                  0
4 Examples + challenges    26            9    18.5                  0

Table 3: Primary linear mixed-effects models

Extract coefficients and 95% Wald confidence intervals directly from the fitted models.

coef_H1a <- as.data.frame(coef(summary(model_H1a)))
coef_H1b <- as.data.frame(coef(summary(model_H1b)))

ci_H1a <- as.data.frame(
  confint(
    model_H1a,
    parm = "beta_",
    method = "Wald"
  )
)

ci_H1b <- as.data.frame(
  confint(
    model_H1b,
    parm = "beta_",
    method = "Wald"
  )
)

Create Table 3.

The intercepts are omitted from the paper table so that the table focuses on the predictors relevant to H1a and H1b.

table3_H1a <- data.frame(
  Model = "H1a",
  Predictor = rownames(coef_H1a),
  Beta = coef_H1a$Estimate,
  SE = coef_H1a$`Std. Error`,
  t = coef_H1a$`t value`,
  CI_low = ci_H1a$`2.5 %`,
  CI_high = ci_H1a$`97.5 %`
)

table3_H1b <- data.frame(
  Model = "H1b",
  Predictor = rownames(coef_H1b),
  Beta = coef_H1b$Estimate,
  SE = coef_H1b$`Std. Error`,
  t = coef_H1b$`t value`,
  CI_low = ci_H1b$`2.5 %`,
  CI_high = ci_H1b$`97.5 %`
)

table3 <- bind_rows(
  table3_H1a,
  table3_H1b
) %>%
  filter(Predictor != "(Intercept)") %>%
  mutate(
    Predictor = recode(
      Predictor,
      examples_before = "Example engagement events",
      preparatory_before = "Preparatory engagement events",
      exercise_number = "Exercise number"
    ),
    across(
      c(Beta, SE, t, CI_low, CI_high),
      ~ round(.x, 3)
    )
  )

print(table3)

  Model                     Predictor   Beta    SE      t CI_low CI_high
1   H1a     Example engagement events  1.153 0.174  6.643  0.813   1.494
2   H1a               Exercise number -0.006 0.005 -1.160 -0.017   0.004
3   H1b Preparatory engagement events  1.101 0.135  8.135  0.836   1.367
4   H1b               Exercise number -0.006 0.005 -1.136 -0.016   0.004

Table 4: Logistic robustness analysis

Function for extracting odds ratios, Wald confidence intervals, and p-values.

extract_logistic_result <- function(model, predictor, label) {

  model_summary <- coef(summary(model))

  beta <- model_summary[predictor, "Estimate"]
  se <- model_summary[predictor, "Std. Error"]
  p_value <- model_summary[predictor, "Pr(>|z|)"]

  data.frame(
    Predictor = label,
    OR = exp(beta),
    CI_low = exp(beta - 1.96 * se),
    CI_high = exp(beta + 1.96 * se),
    p = p_value
  )
}

Create Table 4.

table4 <- bind_rows(
  extract_logistic_result(
    model_H1a_logistic,
    "examples_before",
    "Example-event count"
  ),
  extract_logistic_result(
    model_H1b_logistic,
    "preparatory_before",
    "Preparatory-event count"
  ),
  extract_logistic_result(
    model_H1a_logistic_any,
    "any_examplesTRUE",
    "Any example engagement"
  ),
  extract_logistic_result(
    model_H1b_logistic_any,
    "any_preparationTRUE",
    "Any preparatory engagement"
  )
) %>%
  mutate(
    across(
      c(OR, CI_low, CI_high),
      ~ round(.x, 3)
    ),
    p = if_else(
      p < .001,
      "< .001",
      format(round(p, 3), nsmall = 3)
    )
  )

print(table4)

                   Predictor    OR CI_low CI_high      p
1        Example-event count 1.233  1.123   1.354 < .001
2    Preparatory-event count 1.217  1.120   1.323 < .001
3     Any example engagement 3.725  2.501   5.547 < .001
4 Any preparatory engagement 3.876  2.629   5.714 < .001

13. Figures for the paper

Figure 1: Distribution of Elo rating changes

ggplot(analysis_data, aes(x = Elo_gain)) +
  geom_histogram(
    binwidth = 5,
    boundary = 0,
    color = "black"
  ) +
  labs(
    title = "Distribution of Elo rating changes",
    subtitle = paste0(
      "Changes between consecutive observed exercise events (N = ",
      nrow(analysis_data),
      ")"
    ),
    x = "Elo rating change",
    y = "Number of exercise transitions"
  ) +
  theme_minimal()

Figure 2: Mean Elo change by preparatory engagement

plot_data <- analysis_data %>%
  mutate(
    preparation_group = case_when(
      examples_before == 0 & challenges_before == 0 ~ "None",
      examples_before > 0 & challenges_before == 0 ~ "Examples only",
      examples_before == 0 & challenges_before > 0 ~ "Challenges only",
      examples_before > 0 & challenges_before > 0 ~ "Examples + challenges"
    ),
    preparation_group = factor(
      preparation_group,
      levels = c(
        "None",
        "Examples only",
        "Challenges only",
        "Examples + challenges"
      )
    )
  ) %>%
  group_by(preparation_group) %>%
  summarise(
    n = n(),
    mean = mean(Elo_gain),
    SE = sd(Elo_gain) / sqrt(n),
    lower = mean - 1.96 * SE,
    upper = mean + 1.96 * SE,
    .groups = "drop"
  ) %>%
  mutate(
    group_label = paste0(
      preparation_group,
      "\n(n = ",
      n,
      ")"
    )
  )

ggplot(
  plot_data,
  aes(x = group_label, y = mean)
) +
  geom_point(size = 3) +
  geom_errorbar(
    aes(
      ymin = lower,
      ymax = upper
    ),
    width = 0.15
  ) +
  geom_hline(
    yintercept = 0,
    linetype = "dashed"
  ) +
  labs(
    title = "Mean Elo change by preparatory engagement",
    subtitle = "Points show group means; error bars show 95% confidence intervals",
    x = "Engagement between consecutive exercise events",
    y = "Mean Elo rating change"
  ) +
  theme_minimal()

Figure 3: Predicted probability of a positive Elo change

Predictions are shown for 0–8 preparatory events because observations become sparse at higher engagement counts.

prediction_data <- data.frame(
  preparatory_before = 0:8,
  exercise_number = median(analysis_data$exercise_number)
)

prediction_data$predicted_probability <- predict(
  model_H1b_logistic,
  newdata = prediction_data,
  type = "response",
  re.form = NA
)

ggplot(
  prediction_data,
  aes(
    x = preparatory_before,
    y = predicted_probability
  )
) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  scale_y_continuous(
    labels = scales::percent_format(accuracy = 1)
  ) +
  labs(
    title = "Predicted probability of a positive Elo change",
    subtitle = paste0(
      "Exercise position held at its median; N = ",
      nrow(analysis_data),
      " transitions"
    ),
    x = "Preparatory engagement events before exercise",
    y = "Predicted probability"
  ) +
  theme_minimal()

Figure 4: Positive Elo changes with and without preparation

any_prep_plot <- analysis_data %>%
  mutate(
    preparation = if_else(
      preparatory_before > 0,
      "Any preparation",
      "No preparation"
    ),
    preparation = factor(
      preparation,
      levels = c(
        "No preparation",
        "Any preparation"
      )
    )
  ) %>%
  group_by(preparation) %>%
  summarise(
    n = n(),
    positive_n = sum(Elo_gain > 0),
    positive_pct = mean(Elo_gain > 0) * 100,
    .groups = "drop"
  ) %>%
  mutate(
    label = paste0(
      preparation,
      "\n(n = ",
      n,
      ")"
    )
  )

print(any_prep_plot)

# A tibble: 2 × 5
  preparation         n positive_n positive_pct label                       
  <fct>           <int>      <int>        <dbl> <chr>                       
1 No preparation   2043        219         10.7 "No preparation\n(n = 2043)"
2 Any preparation   164         53         32.3 "Any preparation\n(n = 164)"

ggplot(
  any_prep_plot,
  aes(x = label, y = positive_pct)
) +
  geom_col(width = 0.6) +
  geom_text(
    aes(
      label = paste0(
        round(positive_pct, 1),
        "%"
      )
    ),
    vjust = -0.5
  ) +
  labs(
    title = "Positive Elo changes by preparatory engagement",
    x = NULL,
    y = "Exercise transitions with positive Elo change (%)"
  ) +
  theme_minimal()