  We develop Bayesian predictive stacking for geostatistical models, where the
primary inferential objective is to provide inference on the latent spatial
random field and conduct spatial predictions at arbitrary locations. We exploit
analytically tractable posterior distributions for regression coefficients of
predictors and the realizations of the spatial process conditional upon process
parameters. We subsequently combine such inference by stacking these models
across the range of values of the hyper-parameters. We devise stacking of means
and posterior densities in a manner that is computationally efficient without
resorting to iterative algorithms such as Markov chain Monte Carlo (MCMC) and
can exploit the benefits of parallel computations. We offer novel theoretical
insights into the resulting inference within an infill asymptotic paradigm and
through empirical results showing that stacked inference is comparable to full
sampling-based Bayesian inference at a significantly lower computational cost.
