If I understand the problem well, you have a nice linear model and will just need to juggle some indices. I made a small gist with fake data here, giving an example of how you might do this. Basically you want varying effects in the response to the environmental changes, even though those changes are shared by all the individuals in the sample. The key line is here:
python
species_effects = pm.Normal('species_effect', dims=['isotope', 'species'])
time_effects = pm.Normal('time_effect', dims=['isotope', 'time'])
env_beta = pm.Normal('env_betas', dims=['isotope', 'species', 'env_factors'])
mu_long = (species_effects[i_idx, j_idx] +
time_effects[i_idx, t_idx] +
(env_beta[i_idx, j_idx, :] * X_data[t_idx, :]).sum(axis=-1)) #einsum by hand
- Give all your parameters the dimensions you wrote down in your mathematical notation (i deviated somewhat from yours, allowing the response to vary by both isotope and species)
- Make “long” indices using e.g.
pd.factorizeto map the compact (I, J) and (I, T) parameter arrays to “long”(I * J * T, )arrays (forces everything to line up)
a. Note the use of “fancy-indexing” by passing two arrays – PyMC indexing matches Numpy - For matrix multiplication, stretch out all indices except one, then multiply and sum away that last index.
In my example, I allowed the response to environmental change to vary by isotope-species, but you could make that anything you like (e.g., isotope-species-individual) by adding indexes.
To be honest I’m a bit confused on how to implement GPs when there are multiple dimensions of variation (hierarchical GP?), so you would need to tag in an expert on those if that’s the way you want to go. For example, I’ve read this very old thread that uses block covariance matrices to stitch together many individuals in a dataset with (individual, time), but it seems like there has to be a better way. @bwengals ?