[n.b. I’ve posted a similar question with the brms tag inquiring about a brms-specific solution, but I’m still curious if my more-raw-stan-centered proposal below makes any sense]
I have space-time data (x,y,time) from the neuroimaging field where the data is collected hierarchically such that there are multiple “trials” of data from each of multiple individuals, and I’d like to achieve inference on the 2D+1D smoothed hypersurface averaged-across-trials-and-individuals while permitting systematic deviations from this surface in each individual (and trials just being random noise deviates from the individual’s surface), ideally with shrinkage applied to the magnitude of the individual deviations.
I have experience with implementing Hierarchical 1D GPs in this manner, but I don’t know how to add the spatial dimensions; would it be as simple as something akin to the the automatic relevance determination approach such that I’m still computing point-wise distances but using separate parameters to encode the influence of time temporal distance vs the influence of spatial distance?
Did you have a look at Rob Tranguccis repo with multi-output GPs? The code is probably a bit outdated and could perhaps be more efficient. But maybe it’s still helpful to you?
Good call. I see that the gp_3dim_ncp.stan model specifically seems to implement the isotropy within space and anisotropy between space and time that I’d be looking for (if I’m even using those terms correctly). Very different from how I was thinking of the ARD analogy though, so I’ll have to work on understanding how this works…
Space time models should use a spatiotemporal covariances function. Or a kronecker product structure. Resist the urge to think of it as a 3D space model.
@rtrangucci Does your like imply agreement that the code at gp_3dim_ncp.stan implements the Kronecker product approach that Dan suggests? Just wanted to confirm before I dive into adding hierarchy.
I am really interested in spatio-temporal GPS and specifically I am using INLA to infer model parametres about rainfall events. I would really like to use STAN (both HMC or maybe future Laplace approximation when it will be ready, I ping @charlesm93 if there are news) to extend my INLA model(s), thanks to STAN flexibility and clarity in model building. Is there a way to know your advances in the 2D+1D GPs’ topic?
Rob’s code implements the non-hierarchical scenario as Dan suggests it should be done. Presumably you’re looking either for a write up or performance improvements (i.e. you’re happy with the structure of the model and just want it to sample faster)? I think I remember seeing Rob already published on the Kronecker product stuff, so you could search for his pubs. For speed, you could look at the GPU acceleration stuff (cholesky decompose, etc) and maybe also use the new reduce_sum for the final likelihood part. I’m completely naive about approximate Bayes, but generally wary; have you looked at the splines stuff as an alternative venue for approximation? (i.e. instead of approximating at the level of the Bayesian computation, approximating in the structure of the model and keeping full Bayes)
Yes, I have a non-hierarchical model scenario and I would like to extend my model with censoring (which is something I cannot do with INLA AFAIK) and have more control over generate quantities in space and time.
I never tried the “spline stuff”: have you a link with a simple model to start playing with? Firstly I would generate an artificial rain event and try to recover the parameters.
So speed is important, but firstly I would like to have a model running and to build piece-wise to what I have in mind. I started with INLA on single and real rainfall events, because I found the method in literature. Now I feel that is the moment to take a leap toward generalization (non-linear trend and censoring low values).
By the way i completely missed Kroenker functions in Stan user manual…
It’s time to RT(F)M, i know :-)
thanks @anon75146577 and @mike-lawrence . I will write down some examples starting again from here, solve with INLA and with Stan (maybe starting from a simple model before complicate it too much) and put everything on discourse for rewiev to see if there is a future for my model.
@fabio There is a usable prototype for the embedded Laplace approximation. See this arxiv paper (which needs to be updated).
You do have constrains on which likelihood you can use: for now bernoulli with a logit link, and Poisson with a log link; but the code is expandable. Where you get the most flexibility is in specifying the covariance matrix. And you get to use HMC on the hyperparameters.
One option is use Template Model Builder (TMB) and \texttt{tmbstan}. TMB allows incorporate a GRF by the SPDE method. The problem here is that you need build a template in TMB (.cpp code).
Another problem is to find a correct convergence of the posterior distribution and good mixed of chains for the hyperparameters of the spatial field (\tau and \kappa) with NUTS. I don’t know if the parameterization of the SPDE method can affect this (\tau-\kappa or \rho-\sigma).
The work of @charlesm93 should be a great help for all us if we want use \texttt{Stan} for spatial models.