GAM predictions using empirical bayes - function template

gam_bayes_predict.R

# Template for doing empirical Bayes prediction from a GAM with mgcv
# CJ Brown 2021-08-23
#
# See Wood 2017 and ?predict.gam for more info

#First of all its often easiest to create your prediction dataframe with
# visreg, e.g. this makes a df across all values of NEOLI_Total
# and "has_fee": 
xfit <-  visreg(tour_mod, xvar = "NEOLI_Total", partial = FALSE, 
                plot = FALSE, by = "has_fee")
newdat1 <- xfit$fit %>% filter((NEOLI_Total %% 1) == 0)

#Function template for the prediction.
# you will need to modify it, see below 
gam_bayes_predict <- function(model, newdata, nsim = 1000){
  rmvn <- function(n, mu, sig) { ## MVN random deviates
    L <- mroot(sig)
    m <- ncol(L)
    t(mu + L %*%  matrix(rnorm(m * n), m, n))
  }
  #1000 draws of the parameters from covariance matrix 
  br <- rmvn(1000, coef(model), model$Vp)
  lp1 <- predict(model, newdata = newdata, type = "lpmatrix")
  
  #Simulate posterior
  res <- matrix(NA, nrow = nrow(newdata),
                                   ncol = nsim)
  for (i in 1:nsim){
    val1 <- lp1 %*% br[i,]  
    #MODIFY HERE for the particular non-linear function
    # you want to use
    res[, i] <- 10^val1
  }
  
  #MODIFY here for calculating intervals
  #e.g. for probability >X
  #X <- 5
  # apply(res>X, 1, sum)/nsim
  
  #Or for just returning probabilty quantiles on res: 
  return(
  data.frame(t(apply(res, 1, quantile, probs = c(0.025, 0.25, 0.5, 0.75, 0.975)))) %>%
    cbind(newdata)
}

@maxcampb will find this useful