Expand to include all gp types supported in brms

[nicholasjclark]

Very recent updates to {brms} on Github now allow for a broader range of kernels, which work for both exact and approximate GPs (including two versions of Matern). At present {mvgam} only supports squared exponential, AND it only allows for one-dimensional GP effects. But it is likely possible to let {brms} do all the work to generate the eigenfunction design matrices for any GPs that it supports and to insert those directly when preparing predictions for {mvgam} models:

library(brms)
library(mvgam)

# simulate some nonlinear data
dat <- mgcv::gamSim(1, n = 30, scale = 2)

# initiate a brms GP model using the 'mock' backend so it doesn't actually fit
brms_mock <- brm(y ~ gp(x1, k = 10, scale = FALSE) - 1, data = dat,
                 mock_fit = 1, backend = "mock", rename = FALSE)

# Reduce the size of the brmsfit object
str(brms_mock)
brms_mock$opencl <- NULL
brms_mock$data.name <- NULL
brms_mock$algorithm <- NULL
brms_mock$backend <- NULL
brms_mock$stan_args <- NULL
brms_mock$model <- NULL
brms_mock$stan_funs <- NULL
brms_mock$threads <- NULL
brms_mock$prior <- NULL
brms_mock$family <- NULL

# Pull out the GP design matrix (eigenfunctions)
brms_design <- standata(brms_mock)$Xgp_1

# Now get eigenfunctions from an equivalent mvgam model to show they are 
# identical
mvgam_pre <- mvgam(y ~ gp(x1, k = 10, scale = FALSE) - 1,
                   data = dat,
                   family = gaussian(),
                   run_model = FALSE)
mvgam_design <- standata(mvgam_pre)$X
summary(brms_design - mvgam_design)

# Now for prediction, we need to evaluate newdata
newdat <- data.frame(y = 0, x1 = runif(50, 0, 1))
brms_newdesign <- standata(brms_mock, newdata = newdat,
                           internal = TRUE)$Xgp_1
mvgam_newdesign <- mvgam:::obs_Xp_matrix(newdata = newdat, 
                                         mgcv_model = mvgam_pre$mgcv_model)

summary(brms_newdesign - mvgam_newdesign)

# if names can be tracked appropriately, should be able to pull out the 
# brms design matrix for any approximate gp terms that brms supports; it is just a matter
# of ensuring the appropriate parameters are used in the Stan code
brms_mock2 <- brm(y ~ gp(x1, x2, k = 10, scale = FALSE) - 1, data = dat,
                 mock_fit = 1, backend = "mock", rename = FALSE)
dim(standata(brms_mock2)$Xgp_1)