  Quantifying spatial and/or temporal associations in multivariate geolocated
data of different types is achievable via spatial random effects in a Bayesian
hierarchical model, but severe computational bottlenecks arise when spatial
dependence is encoded as a latent Gaussian process (GP) in the increasingly
common large scale data settings on which we focus. The scenario worsens in
non-Gaussian models because the reduced analytical tractability leads to
additional hurdles to computational efficiency. In this article, we introduce
Bayesian models of spatially referenced data in which the likelihood or the
latent process (or both) are not Gaussian. First, we exploit the advantages of
spatial processes built via directed acyclic graphs, in which case the spatial
nodes enter the Bayesian hierarchy and lead to posterior sampling via routine
Markov chain Monte Carlo (MCMC) methods. Second, motivated by the possible
inefficiencies of popular gradient-based sampling approaches in the
multivariate contexts on which we focus, we introduce the simplified manifold
preconditioner adaptation (SiMPA) algorithm which uses second order information
about the target but avoids expensive matrix operations. We demostrate the
performance and efficiency improvements of our methods relative to alternatives
in extensive synthetic and real world remote sensing and community ecology
applications with large scale data at up to hundreds of thousands of spatial
locations and up to tens of outcomes. Software for the proposed methods is part
of R package 'meshed', available on CRAN.
