Skip to content

Instantly share code, notes, and snippets.

@cbrown5
Created August 22, 2021 23:24
Show Gist options
  • Save cbrown5/9e3a24c323dfc8fc5781a2a058c183a9 to your computer and use it in GitHub Desktop.
Save cbrown5/9e3a24c323dfc8fc5781a2a058c183a9 to your computer and use it in GitHub Desktop.
GAM predictions using empirical bayes - function template
# 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)
}
@cbrown5
Copy link
Author

cbrown5 commented Aug 22, 2021

@maxcampb will find this useful

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment