  Efficient Bayesian inference remains a computational challenge in
hierarchical models. Simulation-based approaches such as Markov Chain Monte
Carlo methods are still popular but have a large computational cost. When
dealing with the large class of Latent Gaussian Models, the INLA methodology
embedded in the R-INLA software provides accurate Bayesian inference by
computing deterministic mixture representation to approximate the joint
posterior, from which marginals are computed. The INLA approach has from the
beginning been targeting to approximate univariate posteriors. In this paper we
lay out the development foundation of the tools for also providing joint
approximations for subsets of the latent field. These approximations inherit
Gaussian copula structure and additionally provide corrections for skewness.
The same idea is carried forward also to sampling from the mixture
representation, which we now can adjust for skewness.

Computation