The Integrated Nested Laplace Approximations (INLA) method has become a commonly used tool for researchers and practitioners to perform approximate Bayesian inference for various fields of applications. It has become essential to incorporate more complex models and expand the method’s capabilities with more features. In this dissertation, we contribute to the INLA method in different aspects. First, we present a new framework, INLA$^+$, based on dense matrices to perform approximate Bayesian inference. An application of the new approach is fitting disease-mapping models for count data with complex interactions. When the precision matrix is dense, the new approach scales better than the existing INLA method and utilizes the power of multiprocessors on shared and distributed memory architectures in today’s computational resources. Second, we propose an adaptive technique to improve gradient estimation for the convex gradient-based optimization framework in INLA. We propose a simple limited-memory technique for improving the accuracy of the numerical gradient of the marginal posterior of the hyperparameter by exploiting a coordinate transformation of the gradient and the history of previously taken descent directions. Third, we extend the commonly utilized Bayesian spatial model in disease mapping, known as the Besag model, into a non-stationary spatial model. This new model considers variations in spatial dependency among a predetermined number of sub-regions. The model incorporates multiple precision parameters, which enable different intensities of spatial dependence in each sub-region. To avoid overfitting and enhance generalization, we derive a joint penalized complexity prior for these parameters. These contributions expand the capabilities of the INLA method, improving its scalability, accuracy, and flexibility for a wider range of applications.