  Motivated by problems from neuroimaging in which existing approaches make use
of "mass univariate" analysis which neglects spatial structure entirely, but
the full joint modelling of all quantities of interest is computationally
infeasible, a novel method for incorporating spatial dependence within a
(potentially large) family of model-selection problems is presented. Spatial
dependence is encoded via a Markov random field model for which a variant of
the pseudo-marginal Markov chain Monte Carlo algorithm is developed and
extended by a further augmentation of the underlying state space. This approach
allows the exploitation of existing unbiased marginal likelihood estimators
used in settings in which spatial independence is normally assumed thereby
facilitating the incorporation of spatial dependence using non-spatial
estimates with minimal additional development effort.
  The proposed algorithm can be realistically used for analysis of %smaller
subsets of large image moderately sized data sets such as $2$D slices of whole
$3$D dynamic PET brain images or other regions of interest. Principled
approximations of the proposed method, together with simple extensions based on
the augmented spaces, are investigated and shown to provide similar results to
the full pseudo-marginal method. Such approximations and extensions allow the
improved performance obtained by incorporating spatial dependence to be
obtained at negligible additional cost. An application to measured PET image
data shows notable improvements in revealing underlying spatial structure when
compared to current methods that assume spatial independence.
