What smoothing method to visualise changing Bernoulli success probability?

Question

Suppose I am running an experiment that goes as follows:

1. I am going to show participants in my experiment a picture, and I tel lthem to look at it for between 0 and 50 seconds.
2. I ask each participant the same question concerning the picture, and record whether they are able to answer the question correctly or not.

My hypothesis is that if the time is too short they won't be able to sufficiently remember the picture, and if it is too long they will get bored and forget the picture, thus, I think that there is an optimal time (between 0-50 seconds) for which recall probability peaks.

Given that my collected data is essentially pairs of points x = time looking at picture and y = 0 or 1, how would I using curve smoothing methods to visualise the relationship between time and recall success? It would be nice if the smoother had the property that is locally well calibrated.

Example data in R simulating my hypothesis:

set.seed(0)
x <- seq(0, 50, 0.1)
prob <- - x*(x - 50)/1000
tmp_df <- data.frame(x, prob)
tmp_df <- tmp_df[rep(1:nrow(tmp_df), times = 10),]
tmp_df$success <- rbinom(n = nrow(tmp_df), size = 1, prob = tmp_df$prob)
tmp_df$keep <- rbinom(n = nrow(tmp_df), size = 1, prob = tmp_df$prob)

For example, the default method = 'gam' smoothing in ggplot2 seems to reproduce the inverse quadratic relationship pretty well, and seems to be quite well calibrated:

ggplot(tmp_df[tmp_df$keep == 1, ], aes(x = x, y = success)) +
    geom_point(position = position_jitter(width = 0.2, height = 0.2), size = 0.1) +
    geom_smooth()

It seems to work quite well even if due to bad experimental design I did not sample the time space uniformly:

ggplot(tmp_df, aes(x = x, y = success)) +
    geom_point(position = position_jitter(width = 0.2, height = 0.2), size = 0.1) +
    geom_smooth()

But, one issue already is that the confidence interval near $x = 0$ and $x = 50$ goes into the negative y-axis.

Answer 1

You should fit a binomial GAM:

library(mgcv)
fit <- gam(success ~ s(x), data = tmp_df, family = binomial)
tmp_df <- cbind(tmp_df, predict(fit, se.fit = TRUE))
tmp_df$lwr <- tmp_df$fit - 2 * tmp_df$se.fit
tmp_df$upr <- tmp_df$fit + 2 * tmp_df$se.fit

library(boot)
tmp_df$fit <- inv.logit(tmp_df$fit)
tmp_df$lwr <- inv.logit(tmp_df$lwr)
tmp_df$upr <- inv.logit(tmp_df$upr)



ggplot(tmp_df[tmp_df$keep == 1, ], aes(x = x, y = success)) +
  geom_point(position = position_jitter(width = 0.2, height = 0.2), size = 0.1) +
  geom_line(aes(y = fit), col = "blue") +
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.2)