Hello,
I am working on a gaussian process model in which I am adding structure from a correlation matrix (vcv in code below) to a model. The model has been tested and validated using test data, and while sampling is quick for small datasets, it is impossibly slow when I try to scale up (increase Nspp, see code below). I am looking for any ideas how to make the model more computationally efficient so I can use it for larger datasets.
My colleagues and I think that the model may be struggling with iterating through a larger matrix, but that’s just our current guess, so alternative hypotheses are most welcome.
Any suggestions for how to improve the model efficiency would be greatly appreciated!
Below is code for the Stan model:
functions {
matrix lambda_vcv(matrix vcv, real lambda, real sigma){
matrix[rows(vcv),cols(vcv)] local_vcv;
matrix[rows(vcv),cols(vcv)] sigma_mat;
local_vcv = vcv * lambda;
for(i in 1:rows(local_vcv))
local_vcv[i,i] = vcv[i,i];
sigma_mat = diag_matrix(rep_vector(sigma, rows(vcv)));
return(sigma_mat * local_vcv * sigma_mat);
}
}
data {
int<lower=1> N;
int<lower=1> Nspp;
int<lower=1, upper= Nspp > species[N];
vector[N] ypred; // response
vector[N] year; // predictor
matrix[Nspp, Nspp] Vphy; //
// Priors
real a_z_prior_mu;
real a_z_prior_sigma;
real lam_interceptsa_prior_alpha;
real lam_interceptsa_prior_beta;
real sigma_interceptsa_prior_mu;
real sigma_interceptsa_prior_sigma;
real b_z_prior_mu;
real b_z_prior_sigma;
real lam_interceptsb_prior_alpha;
real lam_interceptsb_prior_beta;
real sigma_interceptsb_prior_mu;
real sigma_interceptsb_prior_sigma;
real sigma_y_prior_mu;
real sigma_y_prior_sigma;
}
parameters {
real<lower=0> sigma_y;
real<lower=0, upper=1> lam_interceptsa;
real<lower=0> sigma_interceptsa;
real<lower=0, upper=1> lam_interceptsb;
real<lower=0> sigma_interceptsb;
vector[Nspp] b; // slope
real b_z;
vector[Nspp] a; // intercept
real a_z;
}
model {
real yhat[N];
for(i in 1:N){
yhat[i] =
a[species[i]] + b[species[i]] * year[i]; }
a ~ multi_normal(rep_vector(a_z,Nspp), lambda_vcv(Vphy, lam_interceptsa, sigma_interceptsa));
b ~ multi_normal(rep_vector(b_z,Nspp), lambda_vcv(Vphy, lam_interceptsb, sigma_interceptsb));
ypred ~ normal(yhat, sigma_y);
// Prior
a_z ~ normal(a_z_prior_mu, a_z_prior_sigma);
lam_interceptsa ~ beta(lam_interceptsa_prior_alpha, lam_interceptsa_prior_beta);
sigma_interceptsa ~ normal(sigma_interceptsa_prior_mu, sigma_interceptsa_prior_sigma);
b_z ~ normal(b_z_prior_mu, b_z_prior_sigma);
lam_interceptsb ~ beta(lam_interceptsb_prior_alpha, lam_interceptsb_prior_beta);
sigma_interceptsb ~ normal(sigma_interceptsb_prior_mu, sigma_interceptsb_prior_sigma);
sigma_y ~ normal(sigma_y_prior_mu, sigma_y_prior_sigma);
}
The R code to generate our test data is the following:
Nyear <- 10
Nstudy <- 30
Nspp <- 1000
N <- Nspp * Nstudy * Nyear
m <- 0.6
lam <- 0.7 # lambda
sigma_y <- 0.01 # sigma_y
sigma_interceptsb <- 0.1
#generate a tree
spetree <- pbtree(n = Nspp, nsim = 1, b = 1, complete = FALSE, scale = 1)
spetree$tip.label <- paste("s", 1:Nspp, sep = "")
scaledtree <- rescale(spetree, model = "lambda", lam)
slopez <- fastBM(scaledtree, a = m, mu = 0, sig2 = sigma_interceptsb ^ 2)
phylosig(x = slopez, tree = spetree, method = "lambda")
vcvtree <- vcv(spetree, corr = TRUE)
dat <- data.frame(matrix(NA, N, 1))
names(dat) <- c("year")
dat$year <- c(1:Nyear)
dat$study <- rep(c(1:Nstudy), each = Nspp)
dat$species <- rep(1:Nspp, Nstudy)
dat$slopes <-rnorm(Nspp, 5, 3)
dat$smpl.int <- rtruncnorm(Nspp, 0, 365, 188, 20)
dat$doy <- as.integer(dat$year*dat$slopes + dat$smpl.int)
phypriors <- list(a_z_prior_mu = 150,
a_z_prior_sigma = 50,
astudy_prior_mu = 0,
astudy_prior_sigma = 50,
lam_interceptsa_prior_alpha = 1,
lam_interceptsa_prior_beta = 1,
sigma_interceptsa_prior_mu = 40,
sigma_interceptsa_prior_sigma = 20,
b_z_prior_mu = 0,
b_z_prior_sigma = 10,
lam_interceptsb_prior_alpha = 1,
lam_interceptsb_prior_beta = 1,
sigma_interceptsb_prior_mu = 5,
sigma_interceptsb_prior_sigma = 5,
sigma_y_prior_mu = 20,
sigma_y_prior_sigma = 10
)
simu_inits <- function(chain_id) {
a_z.temp <- rnorm(n = Nspp, mean = phypriors[["a_z_prior_mu"]], sd = phypriors[["a_z_prior_sigma"]])
b_z.temp <- rnorm(n = Nspp, mean = phypriors[["b_z_prior_mu"]], sd = phypriors[["b_z_prior_sigma"]])
return(append(list(a = a_z.temp,
b_year = b_z.temp),
phypriors))
}
datalist <- append(list(N = nrow(dat),
Nspp = length(unique(dat$species)),
ypred = dat$doy,
species = dat$species,
year = dat$year,
Vphy = vcvtree), phypriors)
Thanks!
