Stats Assignment

BDS Toolbox

Source: Abe

1: Bias-Variance Trade-off

Q1

In your own words: what is irreducible error? In your own words: explain the bias-variance trade-off

Irreducible error is a part of the total error that cannot be eliminated by any model, even models that capture true relationships perfectly. This erros is due to random variation, measurement errors, unobserved factors, etc..

The bias-variance trade-off is the balance between a model being too simple and too flexible. A simple model has high bias and can miss existing patterns (underfitting), while an overly flexible model has high variance and can fit data too closely, including noise (overfitting).

Q2

The following code makes a test and train variable in the df dataset, and fits a very complex and very basic model.

df <- data.frame(time = 1:10)
df <- df %>% mutate(train_score = 5 + 3 * scale(time) + rnorm(length(time)), 
                    test_score = 5 + 3 * scale(time) + rnorm(length(time)))

mod_basic <- lm(train_score ~ 1, data = df)
mod_complex <- lm(train_score ~ 1 + poly(time, 9), data = df)

Use your R skills to make the plots on slide 14, 15 and 21 from the lecture slides.

df_long <- df %>%
  select(time, train_score, test_score) %>%
  pivot_longer(
    cols = c(train_score, test_score),
    names_to = "set",
    values_to = "score"
  )

# Plot similar to slide 15
ggplot(df_long, aes(x = time, y = score, colour = set)) +
  geom_point() +
  theme_bw()

# Predictions of the basic model
df$pred_basic <- predict(mod_basic)

# Plot of the basic model on the test observations
ggplot(df, aes(x = time, y = test_score)) +
  geom_point(color = "red") +
  geom_line(aes(y = pred_basic)) +
  labs(title = "Basic model") +
  theme_bw()

# Predictions of the complex model
df$pred_complex <- predict(mod_complex)

# Plot similar to slide 21
ggplot(df, aes(x = time, y = test_score)) +
  geom_point(color = "red") +
  geom_line(aes(y = pred_complex)) +
  labs(title = "Complex model") +
  theme_bw()

mse_basic <- mean((df$test_score - df$pred_basic)^2)
mse_complex <- mean((df$test_score - df$pred_complex)^2)

data.frame(
  model = c("Basic", "Complex"),
  test_MSE = c(mse_basic, mse_complex)
)

2: Linear Model

# oefenwebDatabase::close_connections()
con <- oefenwebDatabase::connect()

get data:

user_domain_ratings <- DBI::dbGetQuery(con, "SELECT * FROM user_domain_ratings")

make test and train data:

n <- nrow(user_domain_ratings)
test_ids <- sample(1:n, round(n/5))
prop.table(table((1:n) %in% test_ids))
#> 
#>    FALSE     TRUE 
#> 0.800001 0.199999

test_data <- user_domain_ratings[test_ids, ]
dim(user_domain_ratings)
#> [1] 392132      7

user_domain_ratings <- user_domain_ratings[!(1:n) %in% test_ids, ]
dim(user_domain_ratings)
#> [1] 313706      7

do a simple linear regression

mod <- lm(rating ~ num_played, data = user_domain_ratings)

summary(mod)
#> 
#> Call:
#> lm(formula = rating ~ num_played, data = user_domain_ratings)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -92.418  -8.087   1.103   8.081  40.955 
#> 
#> Coefficients:
#>               Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) -8.041e+00  2.296e-02 -350.27   <2e-16 ***
#> num_played   3.387e-03  3.684e-05   91.94   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 11.6 on 313704 degrees of freedom
#> Multiple R-squared:  0.02624,    Adjusted R-squared:  0.02624 
#> F-statistic:  8454 on 1 and 313704 DF,  p-value: < 2.2e-16

Q3

Significant, yes. But how would you interpret this result? Is this a good model?

user_domain_ratings$y_pred <- predict(mod)

user_domain_ratings %>% 
  slice_sample(prop = .2) %>% 
  ggplot(aes(x = 100 * round(num_played / 100), y = round(rating))) + 
  geom_count(alpha = .2, colour = "blue") +
  geom_line(aes(y = y_pred)) + 
  xlab("Number of items made") + 
  ylab("Player rating") +
  theme_bw()

It seems the relationship is not described by a straight line very well. So the this linear model is probably too simple for this data.

Q4

By just using the num_played variable can you make a better model? Make sure to evaluate on the test data and make sure that you don’t overfit. [See 3.6.5 Non-linear Transformations of the Predictors]

# Keep observations where log(num_played) is defined
user_domain_ratings_log <- user_domain_ratings %>%
  filter(num_played > 0) %>%
  mutate(num_played_log = log(num_played))

# Fit improved model
mod_log <- lm(rating ~ num_played_log,
              data = user_domain_ratings_log)

summary(mod_log)
#> 
#> Call:
#> lm(formula = rating ~ num_played_log, data = user_domain_ratings_log)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -42.570  -7.156   0.152   7.704  39.319 
#> 
#> Coefficients:
#>                 Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)    -13.70240    0.06358  -215.5   <2e-16 ***
#> num_played_log   1.68937    0.01279   132.1   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 10.63 on 278294 degrees of freedom
#> Multiple R-squared:  0.05901,    Adjusted R-squared:  0.059 
#> F-statistic: 1.745e+04 on 1 and 278294 DF,  p-value: < 2.2e-16

# Add predictions and residuals
user_domain_ratings_log$y_pred <- predict(mod_log)
user_domain_ratings_log$resid <- residuals(mod_log)

# Plot model
user_domain_ratings_log %>%
  slice_sample(prop = .2) %>%
  ggplot(aes(x = num_played_log, y = rating)) +
  geom_count(alpha = .2, colour = "blue") +
  geom_line(aes(y = y_pred)) +
  xlab("Log number of items made") +
  ylab("Player rating") +
  theme_bw()

# Plot residuals
user_domain_ratings_log %>%
  slice_sample(prop = .2) %>%
  ggplot(aes(x = y_pred, y = resid)) +
  geom_point(alpha = .2, colour = "blue") +
  geom_hline(yintercept = 0, lty = 2) +
  xlab("Fitted Values") +
  ylab("Residuals") +
  theme_bw()

# Test data
test_data_log <- test_data %>%
  filter(num_played > 0) %>%
  mutate(num_played_log = log(num_played))

# Original model on the same test observations
pred_original <- predict(mod, newdata = test_data_log)

# Log model
pred_log <- predict(mod_log, newdata = test_data_log)

# Compare test MSE
mse_original <- mean((test_data_log$rating - pred_original)^2, na.rm = TRUE)

mse_log <- mean((test_data_log$rating - pred_log)^2, na.rm = TRUE)

data.frame(
  model = c("Original", "Log"),
  test_MSE = c(mse_original, mse_log)
)

We improved the model by taking log. We fitted the linear model and inspected the residuals.

3. Multiple Linear Regression

Get more data:

users <- DBI::dbGetQuery(con, "SELECT id AS new_user_id, gender, grade, date_of_birth FROM users")
user_domain_ratings <- left_join(user_domain_ratings, users)

## Joining with `by = join_by(new_user_id)`

And extend the model:

mod <- lm(rating ~ num_played + 
            as.factor(domain_id) + 
            gender + 
            grade,
          data = user_domain_ratings)
summary(mod)
#> 
#> Call:
#> lm(formula = rating ~ num_played + as.factor(domain_id) + gender + 
#>     grade, data = user_domain_ratings)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -62.482  -3.297  -0.087   3.349  37.222 
#> 
#> Coefficients:
#>                          Estimate Std. Error  t value Pr(>|t|)    
#> (Intercept)            -2.461e+01  6.858e-02 -358.852  < 2e-16 ***
#> num_played              2.117e-03  2.113e-05  100.173  < 2e-16 ***
#> as.factor(domain_id)2   3.961e-01  5.081e-02    7.796 6.43e-15 ***
#> as.factor(domain_id)3  -3.896e+00  5.102e-02  -76.367  < 2e-16 ***
#> as.factor(domain_id)4  -1.352e+01  5.098e-02 -265.242  < 2e-16 ***
#> as.factor(domain_id)5   4.439e+00  5.394e-02   82.288  < 2e-16 ***
#> as.factor(domain_id)15  1.576e+01  5.140e-02  306.595  < 2e-16 ***
#> as.factor(domain_id)16  1.252e+01  5.371e-02  233.103  < 2e-16 ***
#> as.factor(domain_id)60  3.177e+00  5.227e-02   60.780  < 2e-16 ***
#> genderm                 1.084e+00  2.622e-02   41.348  < 2e-16 ***
#> grade                   2.234e+00  8.276e-03  269.915  < 2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 5.976 on 208087 degrees of freedom
#>   (105608 observations deleted due to missingness)
#> Multiple R-squared:  0.7212, Adjusted R-squared:  0.7212 
#> F-statistic: 5.382e+04 on 10 and 208087 DF,  p-value: < 2.2e-16

Q5

Based on the time differene between last_played and the date_of_birth calculate the age of each player. Now improve the model by using multiple variables. What is the best model that you can build. This is how do you get the higher explained variance? Again, make sure to evaluate on the test data and make sure that you don’t overfit.

# Get age
user_domain_ratings <- user_domain_ratings %>%
  mutate(
    age = as.numeric(
      difftime(as.Date(last_played),
               as.Date(date_of_birth),
               units = "days")
    ) / 365.25
  )

# Log nr played
user_domain_ratings <- user_domain_ratings %>%
  filter(num_played > 0) %>%
  mutate(num_played_log = log(num_played))

# Fit MR
mod_multiple <- lm(
  rating ~ num_played_log +
    as.factor(domain_id) +
    gender +
    grade +
    age,
  data = user_domain_ratings
)

summary(mod_multiple)
#> 
#> Call:
#> lm(formula = rating ~ num_played_log + as.factor(domain_id) + 
#>     gender + grade + age, data = user_domain_ratings)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -40.512  -2.965   0.112   3.077  35.028 
#> 
#> Coefficients:
#>                          Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)            -36.398326   0.083117 -437.92   <2e-16 ***
#> num_played_log           1.447754   0.009403  153.97   <2e-16 ***
#> as.factor(domain_id)2    0.705331   0.045797   15.40   <2e-16 ***
#> as.factor(domain_id)3   -3.196597   0.047457  -67.36   <2e-16 ***
#> as.factor(domain_id)4  -12.431808   0.047802 -260.07   <2e-16 ***
#> as.factor(domain_id)5    5.186042   0.048634  106.64   <2e-16 ***
#> as.factor(domain_id)15  17.055612   0.047450  359.44   <2e-16 ***
#> as.factor(domain_id)16  14.108016   0.049447  285.32   <2e-16 ***
#> as.factor(domain_id)60   3.298280   0.046802   70.47   <2e-16 ***
#> genderm                  1.139305   0.023971   47.53   <2e-16 ***
#> grade                    0.658617   0.011106   59.30   <2e-16 ***
#> age                      1.575779   0.010459  150.66   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 5.291 on 195371 degrees of freedom
#>   (82913 observations deleted due to missingness)
#> Multiple R-squared:  0.7593, Adjusted R-squared:  0.7593 
#> F-statistic: 5.602e+04 on 11 and 195371 DF,  p-value: < 2.2e-16

# Fit interaction between num_played_log and domain_id
mod_interaction <- lm(
  rating ~ num_played_log * as.factor(domain_id) +
    gender +
    grade +
    age,
  data = user_domain_ratings
)

summary(mod_interaction)
#> 
#> Call:
#> lm(formula = rating ~ num_played_log * as.factor(domain_id) + 
#>     gender + grade + age, data = user_domain_ratings)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -42.253  -2.752   0.139   2.944  32.842 
#> 
#> Coefficients:
#>                                        Estimate Std. Error  t value Pr(>|t|)
#> (Intercept)                           -30.12081    0.15678 -192.117   <2e-16
#> num_played_log                          0.36648    0.02596   14.117   <2e-16
#> as.factor(domain_id)2                  -3.33629    0.20280  -16.451   <2e-16
#> as.factor(domain_id)3                 -13.47280    0.19152  -70.345   <2e-16
#> as.factor(domain_id)4                 -25.29986    0.18422 -137.333   <2e-16
#> as.factor(domain_id)5                   6.09094    0.18512   32.903   <2e-16
#> as.factor(domain_id)15                 11.43503    0.18856   60.645   <2e-16
#> as.factor(domain_id)16                 10.12749    0.18190   55.675   <2e-16
#> as.factor(domain_id)60                 -7.86759    0.19150  -41.084   <2e-16
#> genderm                                 1.08894    0.02319   46.961   <2e-16
#> grade                                   0.64798    0.01074   60.341   <2e-16
#> age                                     1.56701    0.01012  154.910   <2e-16
#> num_played_log:as.factor(domain_id)2    0.69716    0.03610   19.312   <2e-16
#> num_played_log:as.factor(domain_id)3    1.92065    0.03479   55.205   <2e-16
#> num_played_log:as.factor(domain_id)4    2.44692    0.03322   73.661   <2e-16
#> num_played_log:as.factor(domain_id)5   -0.31824    0.03327   -9.566   <2e-16
#> num_played_log:as.factor(domain_id)15   0.98347    0.03481   28.252   <2e-16
#> num_played_log:as.factor(domain_id)16   0.57237    0.03528   16.222   <2e-16
#> num_played_log:as.factor(domain_id)60   2.03351    0.03376   60.238   <2e-16
#>                                          
#> (Intercept)                           ***
#> num_played_log                        ***
#> as.factor(domain_id)2                 ***
#> as.factor(domain_id)3                 ***
#> as.factor(domain_id)4                 ***
#> as.factor(domain_id)5                 ***
#> as.factor(domain_id)15                ***
#> as.factor(domain_id)16                ***
#> as.factor(domain_id)60                ***
#> genderm                               ***
#> grade                                 ***
#> age                                   ***
#> num_played_log:as.factor(domain_id)2  ***
#> num_played_log:as.factor(domain_id)3  ***
#> num_played_log:as.factor(domain_id)4  ***
#> num_played_log:as.factor(domain_id)5  ***
#> num_played_log:as.factor(domain_id)15 ***
#> num_played_log:as.factor(domain_id)16 ***
#> num_played_log:as.factor(domain_id)60 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 5.11 on 195364 degrees of freedom
#>   (82913 observations deleted due to missingness)
#> Multiple R-squared:  0.7755, Adjusted R-squared:  0.7754 
#> F-statistic: 3.748e+04 on 18 and 195364 DF,  p-value: < 2.2e-16

# Test data
test_data <- left_join(test_data, users)

test_data <- test_data %>%
  mutate(
    age = as.numeric(
      difftime(as.Date(last_played),
               as.Date(date_of_birth),
               units = "days")
    ) / 365.25
  ) %>%
  filter(num_played > 0) %>%
  mutate(num_played_log = log(num_played))

# Predictions
pred_multiple <- predict(mod_multiple, newdata = test_data)
pred_interaction <- predict(mod_interaction, newdata = test_data)

mse_multiple <- mean(
  (test_data$rating - pred_multiple)^2,
  na.rm = TRUE
)

mse_interaction <- mean(
  (test_data$rating - pred_interaction)^2,
  na.rm = TRUE
)

data.frame(
  model = c("Multiple", "Interaction"),
  test_MSE = c(mse_multiple, mse_interaction)
)

Adding multiple predictors considerably improved the model. It explained about 75.9% of the rating variance. We then also added an interaction term, whcih performed better on the test data. Its test MSE was 26.03, compared with 27.79 for the model without the interaction, so it seems to be the better model. Since the improvement is present on the test data, there is no indication of overfitting.

Q6

What happened with the effect of grade after including the effect of gender? [Tip; check page 101]

# Without gender
mod_grade <- lm(
  rating ~ num_played_log +
    as.factor(domain_id) +
    grade,
  data = user_domain_ratings
)

summary(mod_grade)
#> 
#> Call:
#> lm(formula = rating ~ num_played_log + as.factor(domain_id) + 
#>     grade, data = user_domain_ratings)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -49.138  -3.076   0.007   3.166  35.512 
#> 
#> Coefficients:
#>                          Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)            -29.528816   0.056244 -525.01   <2e-16 ***
#> num_played_log           1.745162   0.007504  232.57   <2e-16 ***
#> as.factor(domain_id)2    0.403149   0.039521   10.20   <2e-16 ***
#> as.factor(domain_id)3   -3.477518   0.041624  -83.55   <2e-16 ***
#> as.factor(domain_id)4  -12.714497   0.042203 -301.27   <2e-16 ***
#> as.factor(domain_id)5    5.363509   0.043087  124.48   <2e-16 ***
#> as.factor(domain_id)15  16.738204   0.041364  404.65   <2e-16 ***
#> as.factor(domain_id)16  14.131678   0.043038  328.36   <2e-16 ***
#> as.factor(domain_id)60   2.573178   0.041038   62.70   <2e-16 ***
#> grade                    1.849150   0.006869  269.19   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 5.575 on 274660 degrees of freedom
#>   (3626 observations deleted due to missingness)
#> Multiple R-squared:  0.7409, Adjusted R-squared:  0.7409 
#> F-statistic: 8.725e+04 on 9 and 274660 DF,  p-value: < 2.2e-16

# With gender
mod_grade_gender <- lm(
  rating ~ num_played_log +
    as.factor(domain_id) +
    gender +
    grade,
  data = user_domain_ratings
)

summary(mod_grade_gender)
#> 
#> Call:
#> lm(formula = rating ~ num_played_log + as.factor(domain_id) + 
#>     gender + grade, data = user_domain_ratings)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -48.799  -3.151   0.028   3.279  34.652 
#> 
#> Coefficients:
#>                          Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)            -30.687850   0.077953 -393.67   <2e-16 ***
#> num_played_log           1.905494   0.009384  203.07   <2e-16 ***
#> as.factor(domain_id)2    0.764419   0.048320   15.82   <2e-16 ***
#> as.factor(domain_id)3   -2.706608   0.049965  -54.17   <2e-16 ***
#> as.factor(domain_id)4  -11.827097   0.050272 -235.26   <2e-16 ***
#> as.factor(domain_id)5    5.628703   0.051224  109.88   <2e-16 ***
#> as.factor(domain_id)15  17.228766   0.050062  344.15   <2e-16 ***
#> as.factor(domain_id)16  14.573354   0.052083  279.81   <2e-16 ***
#> as.factor(domain_id)60   3.663696   0.049322   74.28   <2e-16 ***
#> genderm                  1.035967   0.025292   40.96   <2e-16 ***
#> grade                    1.742593   0.008934  195.04   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 5.592 on 195831 degrees of freedom
#>   (82454 observations deleted due to missingness)
#> Multiple R-squared:  0.7312, Adjusted R-squared:  0.7312 
#> F-statistic: 5.328e+04 on 10 and 195831 DF,  p-value: < 2.2e-16

data.frame(
  model = c("Without gender", "With gender"),
  grade_coefficient = c(
    coef(mod_grade)["grade"],
    coef(mod_grade_gender)["grade"]
  )
)

After adding gender to the model, the effect of grade becomes smaller because some of the effect of grades is now accounted for by gender, but the change is quite small.

Q7

Now who would you explain this model to a collegue without a data science background? Would you want to change the model to more simplified model?

I would explain that the model predicts player rating using practice, domain, gender, grade and age. The effect of practice can also differ between domains. I would keep the more complex model for prediction because it performs better on the test data, but use the simpler model if interpretability is more important.

4. Generalized LM

Bonus: Run the following code, inspect the models, see if you can understand / improve / visualise some results. You can hand in your code, but this is not required.

GLM is more of the same, but you need different tools! Lets see if can we predict the use of the question-mark button. Get some data

logs <- DBI::dbGetQuery(con, "SELECT * FROM log_records_1")

logs <- logs %>% mutate(
  q = 1 * (answer == "¿")
)

mod_basic <- glm(q ~ show_coins, data = logs, family = binomial())

summary(mod_basic)
#> 
#> Call:
#> glm(formula = q ~ show_coins, family = binomial(), data = logs)
#> 
#> Coefficients:
#>             Estimate Std. Error  z value Pr(>|z|)    
#> (Intercept) -3.45921    0.01160 -298.239   <2e-16 ***
#> show_coins   0.11167    0.01237    9.031   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 560696  on 1910170  degrees of freedom
#> Residual deviance: 560612  on 1910169  degrees of freedom
#> AIC: 560616
#> 
#> Number of Fisher Scoring iterations: 6

mod0 <- glm(q ~ 1, data = logs, family = binomial())
table(predict(mod0, type = "response"))
#> 
#> 0.0335179415908896 
#>            1910171

anova(mod_basic, mod0, test="Chisq")

## Analysis of Deviance Table
## 
## Model 1: q ~ show_coins
## Model 2: q ~ 1
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1   1910169     560612                          
## 2   1910170     560696 -1  -83.751 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

logs$thbe <- logs$new_user_domain_rating - logs$item_rating
mod_improved <- glm(q ~ show_coins + thbe + as.factor(difficulty), data = logs, family = binomial())
summary(mod_improved)
#> 
#> Call:
#> glm(formula = q ~ show_coins + thbe + as.factor(difficulty), 
#>     family = binomial(), data = logs)
#> 
#> Coefficients:
#>                         Estimate Std. Error  z value Pr(>|z|)    
#> (Intercept)            -2.896742   0.019868 -145.802  < 2e-16 ***
#> show_coins              0.126375   0.012596   10.033  < 2e-16 ***
#> thbe                   -0.646660   0.006803  -95.061  < 2e-16 ***
#> as.factor(difficulty)1 -0.081523   0.015740   -5.179 2.23e-07 ***
#> as.factor(difficulty)2  0.043604   0.016447    2.651  0.00802 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 560696  on 1910170  degrees of freedom
#> Residual deviance: 540656  on 1910166  degrees of freedom
#> AIC: 540666
#> 
#> Number of Fisher Scoring iterations: 7

data.frame(
  model = c("Basic", "Improved"),
  AIC = c(AIC(mod_basic), AIC(mod_improved))
)

logs$pred <- predict(mod_improved, type = "response")

logs %>%
  group_by(difficulty) %>%
  summarise(
    predicted_probability = mean(pred, na.rm = TRUE)
  ) %>%
  ggplot(aes(
    x = as.factor(difficulty),
    y = predicted_probability
  )) +
  geom_col() +
  xlab("Difficulty") +
  ylab("Predicted probability of question-mark use") +
  theme_bw()

oefenwebDatabase::close_connections()
#> [[1]]
#> [1] TRUE