Hi everyone,
I have been using Gaussian processes to model the north-south migration of GPS-tagged animals over time. Due to the number of observations (from several hundred to several thousand obs. per individual, meaning around 40000 obs. in total) I have been using the Hilbert space basis function approximation.
Currently, my most elaborate model is a heteroskedastic GP with a slow trend and yearly cyclic trend (adapting the code from here and here) - this model samples without issue.
As a next step, I would like to extend the model to incorporate the individual ID of each animal - although individuals generally follow a similar migration pattern, there is variation in the exact timing and overall distance that they migrate. I am thinking that a way to model this would be some sort of hierarchical GP, a bit like the one in this tutorial. - each individual would have its own GP (planning to use Matérn 1/2 kernel for these) , but would contribute to an overall average migration trajectory (which would presumably be modelled by the slow trend + cyclic GP as in my current model).
However, I am struggling to code a Hilbert basis function version of this type of model. In particular, it’s unclear to me exactly how the indexing would work - presumably I’d need to have separate basis functions and spectral densities for each individual, and combine these using append_row in normal_id_glm? Am I thinking along the right lines here, or am I barking up the wrong tree?
To further complicate things, the GPS data for each individual have different start and end points (not all individuals were tagged at once, tags were collected at different times) and thus quite different lengths. Also, data from different individuals may be spaced differently (e.g. some GPS tags return daily fixes, while others only return weekly fixes).
Any help with this problem would be greatly appreciated!
Heteroskedastic model code (no individual ID):
functions{
vector diagSPD_EQ(real alpha, real rho, real L, int M) {
return sqrt((alpha) * sqrt(2*pi()) * rho * exp(-0.5*(rho*pi()/2/L)^2 * linspaced_vector(M, 1, M)^2));
}
vector diagSPD_periodic(real alpha, real rho, int M) {
real a = 1/rho^2;
array[M] int one_to_M;
for (m in 1:M) one_to_M[m] = m;
vector[M] q = sqrt(alpha * 2 / exp(a) * to_vector(modified_bessel_first_kind(one_to_M, a)));
return append_row(q,q);
}
matrix PHI(int N, int M, real L, vector x) {
return sin(diag_post_multiply(rep_matrix(pi()/(2*L) * (x+L), M), linspaced_vector(M, 1, M)))/sqrt(L);
}
matrix PHI_periodic(int N, int M, real w0, vector x) {
matrix[N,M] mw0x = diag_post_multiply(rep_matrix(w0*x, M), linspaced_vector(M, 1, M));
return append_col(cos(mw0x), sin(mw0x));
}
}
data{
int N; // Number of observations
array[N] real y; // Northing
array[N] real x; // Time (day)
real<lower=0> c_f1; // factor c to determine the boundary value L
int<lower=1> M_f1; // number of basis functions for smooth function
int<lower=1> J_f2; // number of cos and sin functions for periodic function
real<lower=0> c_g1; // factor c to determine the boundary value L
int<lower=1> M_g1; // number of basis functions for smooth function
}
transformed data {
// Standardize data
real xmean = mean(x);
real ymean = mean(y);
real xsd = sd(x);
real ysd = sd(y);
vector[N] xn = (to_vector(x) - xmean)/xsd;
vector[N] yn = (to_vector(y) - ymean)/ysd;
// Basis functions for f1
real L_f1 = c_f1*max(xn);
matrix[N,M_f1] PHI_f1 = PHI(N, M_f1, L_f1, xn);
// Basis functions for f2
real period_year = 365.25/xsd;
matrix[N,2*J_f2] PHI_f2 = PHI_periodic(N, J_f2, 2*pi()/period_year, xn);
// Basis functions for g1
real L_g1 = c_g1*max(xn);
matrix[N,M_g1] PHI_g1 = PHI(N, M_g1, L_g1, xn);
// Concatenated basis functions
matrix[N,M_f1+2*J_f2] PHI_f = append_col(PHI_f1, PHI_f2);
}
parameters {
vector[M_f1] beta_f1; // the basis functions coefficients for f1
vector[2*J_f2] beta_f2; // the basis functions coefficients for f2
vector[M_g1] beta_g1; // the basis functions coefficients for g1
real<lower=0> lengthscale_f1; // lengthscale for f1
real<lower=0> lengthscale_f2; // lengthscale for f2
real<lower=0> lengthscale_g1; // lengthscale for g1
real<lower=0> sigma_f1; // scale of f1
real<lower=0> sigma_f2; // scale of f2
real<lower=0> sigma_g1; // scale of g1
real intercept_g;
}
model {
// spectral densities for f1, f2, and g1
vector[M_f1] diagSPD_f1 = diagSPD_EQ(sigma_f1, lengthscale_f1, L_f1, M_f1);
vector[2*J_f2] diagSPD_f2 = diagSPD_periodic(sigma_f2, lengthscale_f2, J_f2);
vector[M_f1] diagSPD_g1 = diagSPD_EQ(sigma_g1, lengthscale_g1, L_g1, M_g1);
// priors
beta_f1 ~ normal(0, 1);
beta_f2 ~ normal(0, 1);
beta_g1 ~ normal(0, 1);
lengthscale_f1 ~ exponential(1);
lengthscale_f2 ~ normal(0, 0.1);
lengthscale_g1 ~ exponential(1);
sigma_f1 ~ exponential(1);
sigma_f2 ~ exponential(1);
sigma_g1 ~ exponential(1);
intercept_g ~ normal(0, 0.1);
// Model
yn ~ normal_id_glm(PHI_f,
0.0,
append_row(diagSPD_f1 .* beta_f1, diagSPD_f2 .* beta_f2),
exp(intercept_g + PHI_g1 * (diagSPD_g1 .* beta_g1)));
}
generated quantities {
vector[N] f1;
vector[N] f2;
vector[N] f;
vector[N] sigma;
vector[N] log_lik;
{
// spectral densities for f1, f2, and g1
vector[M_f1] diagSPD_f1 = diagSPD_EQ(sigma_f1, lengthscale_f1, L_f1, M_f1);
vector[2*J_f2] diagSPD_f2 = diagSPD_periodic(sigma_f2, lengthscale_f2, J_f2);
vector[M_f1] diagSPD_g1 = diagSPD_EQ(sigma_g1, lengthscale_g1, L_g1, M_g1);
// functions scaled back to original scale
f1 = (0.0 + PHI_f1 * (diagSPD_f1 .* beta_f1))*ysd;
f2 = (PHI_f2 * (diagSPD_f2 .* beta_f2))*ysd;
f = f1 + f2 + ymean;
sigma = exp(intercept_g + PHI_g1 * (diagSPD_g1 .* beta_g1))*ysd;
// log_liks for loo
for (n in 1:N) log_lik[n] = normal_lpdf(y[n] | f[n], sigma[n]);
}
}