Skip to content

Instantly share code, notes, and snippets.

Show Gist options
  • Save gavinsimpson/e3587067817536a467dec091d479840b to your computer and use it in GitHub Desktop.
Save gavinsimpson/e3587067817536a467dec091d479840b to your computer and use it in GitHub Desktop.
Predicting from a model fitted with `bam()` when `discrete = TRUE` was used to both fit and predict from the model with `terms` in the `predict` call doesn't include the model intercept
# Weird behaviour of predict.bam()
# packages
library("mgcv")
# simulate some data
set.seed(42)
df <- gamSim(1, n = 10000)
# fit a bam with and without discrete = TRUE and a gam
m_bam <- bam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = df, method = "fREML")
m_dam <- bam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = df, method = "fREML",
discrete = TRUE)
m_gam <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = df, method = "REML")
# data to predict at
new_df <- as.data.frame(lapply(m_gam$var.summary, `[`, 2))
sms <- c("s(x0)", "s(x1)", "s(x2)")
p_bam_d <- predict(m_bam, newdata = new_df, terms = sms)
p_bam <- predict(m_bam, newdata = new_df, terms = sms, discrete = FALSE)
p_dam_d <- predict(m_dam, newdata = new_df, terms = sms)
p_dam <- predict(m_dam, newdata = new_df, terms = sms, discrete = FALSE)
p_gam <- predict(m_gam, newdata = new_df, terms = sms)
output <- data.frame(fun = rep(c("bam", "gam"), times = c(4,1)),
discrete_fit = c(FALSE, FALSE, TRUE, TRUE, FALSE),
discrete_predict = c(TRUE, FALSE, TRUE, FALSE, FALSE),
with_terms = c(p_bam_d, p_bam, p_dam_d, p_dam, p_gam))
# now with exclude instead of terms, to see if there is a difference there
p_bam_d <- predict(m_bam, newdata = new_df, exclude = "s(x3)")
p_bam <- predict(m_bam, newdata = new_df, exclude = "s(x3)", discrete = FALSE)
p_dam_d <- predict(m_dam, newdata = new_df, exclude = "s(x3)")
p_dam <- predict(m_dam, newdata = new_df, exclude = "s(x3)", discrete = FALSE)
p_gam <- predict(m_gam, newdata = new_df, exclude = "s(x3)")
output <- transform(output,
with_exclude = c(p_bam_d, p_bam, p_dam_d, p_dam, p_gam))
output
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment