predict.bam() produces different results when data is subset

Apologies in advance since I cannot provide a reproducible example due to the immense size of my model. I'll do my best to describe my situation fully, hopefully this will be sufficient.

My model looks like this:

model <- binary_outcome ~ s(height, side, bottom, top, by=level) + s(side_angle, height_angle, by=level) + fixed_eff + level


I'm trying to use splines to model non-linear relationships in the data, splitting them by level. I am fairly confident that this model structure makes sense. Due to data set size issues, I estimate a different model for each year.

I then fit the model:

fit <- bam(model, data=d, family=binomial(link="probit"), gc.level=1, nthreads=14, control=gam.control(trace=TRUE))


The first thing I've noticed is that the size of fit does not vary in proportion by d. The largest data sets sometimes create the smallest model sizes. Inspecting this further it seems to be that fit$smooths has an object X0 that drives the size of fit. The largest objects have nrow(X0) == nrow(d) whereas the smaller objects have nrow(X0) << nrow(d). This leads me to the biggest issue which is that the smaller fit objects tend to produce much worse predictions (measuring by RMSE) than the larger objects when I do: d$pred <- predict.bam(fit, newdata=d, type="response")


Now I've found a solution which seems to produce predictions in line with what I would expect, but it seems strange that it should work, and I would love to know why it's working. It involves subsetting d by level, applying d$$level <- factor(d$$level) and then running predict.bam() again on the subsetted data:

d_l <- d[d$$level == level,] d_l$$level <- factor(d$$level) d_l$$pred <- predict.bam(fit, newdata=d_l, type="response")


This seems odd since the data is essentially the same, the only difference really is that the level factor now only has one level since the data was subset in line 1.

Happy to elaborate or test out any suggestions/ideas. Thanks!

1 Answer

It appears when I do

d\$pred <- predict.gam(fit, newdata=d, type="response")


The predictions are in line with doing

d_l <- d[d$$level == level,] d_l$$level <- factor(d$$level) d_l$$pred <- predict.bam(fit, newdata=d_l, type="response")


Still unclear why this would be different.