Latent variable models for ordinal data represent a useful tool in different fields of research in which the constructs of interest are not directly observable. In such models, problems related to the integration of the likelihood function can arise since analytical solutions do not exist. Numerical approximations, like the widely used Gauss Hermite (GH) quadrature, are generally applied to solve these problems. However, GH becomes unfeasible as the number of latent variables increases. Thus, alternative solutions have to be found. In this paper, we propose an extended version of the Laplace method for approximating the integrals, known as fully exponential Laplace approximation. It is computational feasible also in presence of many latent variables, and it is more accurate than the classical Laplace method.