Add gps for other covariates and use Hilbert equations for predictions

[nicholasjclark]

Can add gp() function as in brms by including an attribute for the linear predictor matrix (non_mgcv_terms). Will then only need to update the mvgam_predict function and make sure all uses of the Xp matrix have the attribute included before calling mvgam_predict.

For gps, can use existing stan functions and improve efficiency of predictions using basis function weights:

# Evaluate Laplacian eigenfunction for a given GP basis function
phi = function(boundary, m, centred_covariate) { 
  1 / sqrt(boundary) * sin((m * pi)/(2 * boundary) * (centred_covariate + boundary)) 
}

# Evaluate eigenvalues for a given GP basis function
lambda = function(boundary, m) {
  ((m * pi)/(2 * boundary))^2 
}

# Spectral density squared exponential Gaussian Process kernel
spd = function(alpha_gp, rho_gp, eigenvalues) { 
  (alpha_gp^2) * sqrt(2 * pi) * rho_gp * exp(-0.5 * (rho_gp^2) * (eigenvalues^2)) 
}

library(mvgam)
simdat = sim_mvgam(n_series = 1, trend_model = 'GP')
mod = mvgam(y ~ 1,
            trend_model = 'GP',
            data = simdat$data_train,
            run_model = T)
code(mod)

object <- mod
new_covariate <- 1:100
covariate <- 'time'

# Need original covariate from observed data (centred)
covariate_cent <- unique(object$obs_data[[covariate]]) - 
  mean(unique(object$obs_data[[covariate]]))

# Need to know number of original GP basis functions
# (best to store this with the model as could be confusing when multiple
# GPs are used in one model)
num_basis_line <- object$model_file[grep('num_gp_basis = ', object$model_file)]
num_gp_basis <- as.numeric(unlist(regmatches(num_basis_line, 
                                        gregexpr("[[:digit:]]+", num_basis_line))))

# Get draws of GP parameters, including the basis function weights
gp_pars <- as.data.frame(object, variable = 'trend_params')

# Get vector of eigenvalues of covariance matrix
eigenvalues <- vector()
for(m in 1:num_gp_basis){
  eigenvalues[m] <- lambda(boundary = (5.0/4) * 
                             (max(covariate_cent) - min(covariate_cent)),
         m = m)
}

# Get vector of eigenfunctions
eigenfunctions <- matrix(NA, nrow = length(new_covariate),
                         ncol = num_gp_basis)
for(m in 1:num_gp_basis){
  eigenfunctions[, m] <- phi(boundary = (5.0/4) * (max(covariate_cent) - min(covariate_cent)),
                           m = m,
                           centred_covariate = new_covariate - 
                             mean(unique(object$obs_data[[covariate]])))
}

# Compute diagonal of covariance matrix
gp_preds <- matrix(NA, nrow = NROW(gp_pars), ncol = NROW(eigenfunctions))
for(i in 1:NROW(gp_pars)){
  diag_SPD <- sqrt(spd(alpha_gp = gp_pars$`alpha_gp[1]`[i], 
                       rho_gp = gp_pars$`rho_gp[1]`[i], 
                       sqrt(eigenvalues)))
  b_gp <- unlist(gp_pars[i, grep('b_gp', colnames(gp_pars))])
  gp_preds[i, ] <- (diag_SPD * b_gp) %*% t(eigenfunctions)
}

# Compute GP draw
plot(gp_preds[1,], type = 'l',
     ylim = range(gp_preds))
for(i in 2:50){
  lines(gp_preds[i,])
}

# Draws should be identical
trends <- mvgam:::mcmc_chains(object$model_output, 'trend')
plot(trends[1,], type = 'l')
lines(gp_preds[1,], col = 'red')

[nicholasjclark]

Implemented now but requires further testing. Other aspects still needed are:

1. Allow updating of priors
2. Work out the scaling function from brms
3. Implement tests for continuous by variables
4. Add gp params to variables
5. Add to the dynamic wrapper and set as default
6. Ensure documentation is clear, with k=10 and c=5/4 as defaults
7. Write a vignette to support the functions

[nicholasjclark]

All done