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!