Summary The integrated nested Laplace approximation (INLA) is a convenient way to obtain approximations to the posterior marginals for parameters in Bayesian hierarchical models when the latent effects can be expressed as a Gaussian Markov random field. In addition, its implementation in the R-INLA package for the R statistical software provides an easy way to fit models using the INLA in practice in a fraction of the time that other computer-intensive methods (e.g. Markov chain Monte Carlo methods) take to fit the same model. Although the INLA provides a fast approximation to the marginals of the model parameters, it is difficult to use it with models that are not implemented in R-INLA. It is also difficult to make multivariate posterior inference on the parameters of the model as the INLA focuses on the posterior marginals and not the joint posterior distribution. We describe how to use the INLA within the Metropolis–Hastings algorithm to fit complex spatial models and to estimate the joint posterior distribution of a small number of parameters. We illustrate the benefits of this new method with two examples. In the first, a spatial econometrics model with two auto-correlation parameters (for the response and the error term) is considered. This model is not currently available in R-INLA, and multivariate inference is often required to assess dependence between the two spatial auto-correlation parameters in the model. Furthermore, the estimation of spillover effects is based on the joint posterior distribution of a spatial auto-correlation parameter and a covariate coefficient. In the second example, a multivariate spatial model for several diseases is proposed for disease mapping. This model includes a shared specific spatial effect as well as disease-specific spatial effects. Dependence on the shared spatial effect is modulated via disease-specific weights. By inspecting the joint posterior distribution of these weights it is possible to assess which diseases have a similar spatial pattern.