  In this paper we recast the problem of missing values in the covariates of a
regression model as a latent Gaussian Markov random field (GMRF) model in a
fully Bayesian framework. Our proposed approach is based on the definition of
the covariate imputation sub-model as a latent effect with a GMRF structure. We
show how this formulation works for continuous covariates and provide some
insight on how this could be extended to categorical covariates.
  The resulting Bayesian hierarchical model naturally fits within the
integrated nested Laplace approximation (INLA) framework, which we use for
model fitting. Hence, our work fills an important gap in the INLA methodology
as it allows to treat models with missing values in the covariates.
  As in any other fully Bayesian framework, by relying on INLA for model
fitting it is possible to formulate a joint model for the data, the imputed
covariates and their missingness mechanism. In this way, we are able to tackle
the more general problem of assessing the missingness mechanism by conducting a
sensitivity analysis on the different alternatives to model the non-observed
covariates.
  Finally, we illustrate the proposed approach with two examples on modeling
health risk factors and disease mapping. Here, we rely on two different
imputation mechanisms based on a typical multiple linear regression and a
spatial model, respectively. Given the speed of model fitting with INLA we are
able to fit joint models in a short time, and to easily conduct sensitivity
analyses.
