  It has become increasingly common to collect high-dimensional binary response
data; for example, with the emergence of new sampling techniques in ecology. In
smaller dimensions, multivariate probit (MVP) models are routinely used for
inferences. However, algorithms for fitting such models face issues in scaling
up to high dimensions due to the intractability of the likelihood, involving an
integral over a multivariate normal distribution having no analytic form.
Although a variety of algorithms have been proposed to approximate this
intractable integral, these approaches are difficult to implement and/or
inaccurate in high dimensions. Our main focus is in accommodating
high-dimensional binary response data with a small to moderate number of
covariates. We propose a two-stage approach for inference on model parameters
while taking care of uncertainty propagation between the stages. We use the
special structure of latent Gaussian models to reduce the highly expensive
computation involved in joint parameter estimation to focus inference on
marginal distributions of model parameters. This essentially makes the method
embarrassingly parallel for both stages. We illustrate performance in
simulations and applications to joint species distribution modeling in ecology.
