Bayesian spatial models are widely used to analyse data that arise in scientific disciplines such as health, ecology, and the environment. Traditionally, Markov chain Monte Carlo (MCMC) methods have been used to fit these type of models. However, these are highly computationally intensive methods that present a wide range of issues in terms of convergence and can become infeasible in big data problems. The integrated nested Laplace approximation (INLA) method is a computational less-intensive alternative to MCMC that allows us to perform approximate Bayesian inference in latent Gaussian models such as generalised linear mixed models and spatial and spatio-temporal models. This approach can be used in combination with the stochastic partial differential equation (SPDE) approach to analyse geostatistical data that have been collected at particular sites to predict the spatial process underlying the data as well as to assess the effect of covariates and model other sources of variability. Here we demonstrate how to fit a Bayesian spatial model using the INLA and SPDE approaches applied to freely available data of malaria prevalence and risk factors in Mozambique. We show how to fit and interpret the model to predict malaria risk and assess the effect of covariates using the R-INLA package, and provide the R code necessary to reproduce the results or to use it in other spatial applications.