Gaussian process (GP) models are computationally demanding for large datasets. Much work has been done to avoid expensive matrix operations that arise in parameter estimation with larger datasets via sparse and/or reduced rank covariance matrices (Datta et al. 2016 provide a nice review). What follows is an implementation of a spatial Gaussian predictive process Poisson GLM in Stan, following...