Abnormally large confidence interval with binomial gam when p->1 at max of predictor range

Question

I am running a gam (mgcv in R) to model a non-linear effect of time on a binomial reponse (positive or negative sample). This is a minimal example of such a model:

gam<- gam(pcr ~ s(rel_time, k=7), family = binomial(link="logit"), data=df)

For one of the species there are exclusively positive samples the last quarter of the predictor range (rel_time), and therefore the predicted real values of p is extremely close to one, as expected. However, then I also get abnormally large predicted confidence intervals (see figure of predicted values, with jittered observations as well, below).

This is how i have calculated predicted values from the model:

new_data <- data.frame(rel_time= seq(from=-33, to=710, by=1))

gam_pred <- bind_cols(new_data, as.data.frame(predict(gam, newdata = new_data, se.fit = TRUE)))

test3_pred$fit_r<- plogis(test3_pred$fit)
test3_pred$up<- plogis(test3_pred$fit+2*test3_pred$se.fit)
test3_pred$down<- plogis(test3_pred$fit-2*test3_pred$se.fit)

It seems that Wilson confidence intervals should be the way to go (Confidence interval around binomial estimate of 0 or 1, but I have no clue how to calculate this when predicting from a gam.

I have experimented with different link functions, but with no major effect. If I run a simple glm, the problem disappears, but clearly this is no option. By the way: If I do not specify k to be smaller than 8, the estimates of p also becomes a strange, and drops down to around 0.2 dispite only postive values:

As should be visible from the figure, the sampling interval in longer towards the second half of the sampling period. I am not convinced that this change, when k is different, is due to the same issue.

I can share the data if needed.

I would appreciate if someone can provide advice on how to deal with this.

Answer 1

The confidence interval is being formed assuming a Gaussian approximation to the posterior distribution and I think you've hit one of the settings where this approximation can be not very good at all.

Luckily, mgcv provides an alternative Metropolis Hastings sampler which alternates proposals from a Student's t distribution, with random walk proposals from a shrunken version of the Bayesian parameter covariance matrix.

The details of the sampler are available at ?gam.mh and you can follow the example there for an illustration.

I have written a translation of the example that uses functionality in my gratia package that simplifies doing the posterior sampling required and formats it in a tidy way: https://gavinsimpson.github.io/gratia/articles/posterior-simulation.html#metropolis-hastings-sampler

The salient figure is

which shows much better credible intervals using the Metropolis Hastings sampler compared with the Gaussian approximation.

Before increasing k you should address the issue of smoothing autocorrelated data, as smoothers will fit the autocorrelation which may or may not be what you want. The time looks irregular here so that raises some problems currently for mgcv but brms should be able to handle it with a 1d spatial correlation structure if you want to stay within the penalized spline world. Also, with brms you won't need to do extra posterior sampling because it will already have done that for you using Stan as brms fits models using Hamiltonian Monte Carlo.