Multivariate Probit

Definition

The multivariate probit model jointly analyzes $p$ binary outcome variables for the same subject by positing an underlying $p$ -dimensional continuous latent response vector $\mathbf{y}_i$ distributed multivariate-normally. Subject $i$ is observed to succeed on outcome $j$ ( $u_{ij}=1$ ) whenever the $j$ -th latent variable exceeds zero. Correlation among binary outcomes is captured via shared latent factors rather than a free $p\times p$ correlation matrix, making the model tractable for $p>2$ .

Key Ideas

Latent response model: $\mathbf{y}_i = \mathbf{B}\mathbf{x}_i + \mathbf{\Lambda}\boldsymbol{\theta} + \boldsymbol{\zeta}$ , where $\mathbf{x}_i$ is a $q$ -vector of covariates, $\mathbf{B}$ is $p\times q$ , $\boldsymbol{\theta} \sim \mathcal{N}(0, I_m)$ is an $m$ -dimensional latent factor ( $m < p$ ), $\mathbf{\Lambda}$ is the $p\times m$ factor loading matrix, and $\boldsymbol{\zeta} \sim \mathcal{N}(0, d^2 I_p)$ with $d_j^2 = 1 - \sum_k \lambda_{jk}^2$ .
Response probability: Conditional on $\boldsymbol{\theta}$ , outcomes are independent: $P(u_{ij}=1|\boldsymbol{\theta}) = \Phi(z_{ij})$ where $z_{ij} = (\beta_j \mathbf{x}_i + \lambda_j \boldsymbol{\theta})/d_j$ .
Marginal probability: Integrating out $\boldsymbol{\theta}$ gives the probability of pattern $\mathbf{u}_i$ : $P(\mathbf{u}_i) = \int L_i(\boldsymbol{\theta}) \prod_k \phi(\theta_k)\,d\boldsymbol{\theta}$ , approximated by $m$ -fold Gauss-Hermite quadrature.
All $2^p$ pattern probabilities are estimable from the model for any covariate vector; bootstrapped standard errors obtainable.
Single joint test: Covariate effects are tested simultaneously for all $p$ outcomes, accounting for their inter-correlation — more powerful than $p$ separate univariate tests.

How It Works

Estimation (Bock-Gibbons 1996; Gibbons-Wilcox-Gok 1998):

Expectation-Maximization (EM) initialization: Run EM cycles starting from $\lambda_{jk}=\mathbf{1}_{k=j}$ and $d_j=1$ to obtain suitable $\Lambda$ values.
Fisher scoring: Fix $m(m+1)/2$ elements of $\Lambda$ in lower-triangular form to resolve rotational indeterminacy; score the remaining parameters with standard errors from the information matrix (Aitchison-Silvey constrained approach or direct fixing).
Quadrature: Gauss-Hermite with $Q \geq 5$ points per dimension approximates the $m$ -fold integral over $\boldsymbol{\theta}$ .
Empirical prior (Bock-Aitkin 1981): Replace Gaussian quadrature weights with nonparametric posterior-density estimates $g(\mathbf{T}_k) \propto \sum_i P(\mathbf{u}_i|\mathbf{T}_k)A(\mathbf{T}_k)$ ; iterate to joint convergence.
Sampling weights: $\log L = \sum_i r_i \log P(\mathbf{u}_i)$ where $r_i$ are observation weights; derivatives scale accordingly.

Model selection: Change in log-likelihood (deviance) as factors are added; chi-squared approximation satisfactory for deviance differences even with ungrouped data. Absolute fit testable only when covariates are excluded (response patterns can be grouped into $2^p$ cells).

Identification: $pq + pm - m(m+1)/2$ estimable parameters — equals the saturated model mean parameters plus pairwise correlations. Lower-triangular $\Lambda$ (set $\lambda_{jk}=0$ for $k>j$ ) removes rotational indeterminacy.

Why It Matters

Improves on univariate analysis: A single joint test of covariate effects across all $p$ outcomes is statistically more powerful than $p$ separate regressions and correctly handles the multiple-testing problem via the model's correlation structure.
Distinct from generalized estimating equations (GEE): Full likelihood enables marginal likelihood model comparison and exact pattern probability inference; GEE is semiparametric and provides no pattern probabilities.
Distinct from multinomial probit: Multinomial probit models choosing one alternative from $p$ mutually exclusive options (Independence of Irrelevant Alternatives (IIA) avoidance); multivariate probit models $p$ simultaneous binary outcomes that need not be mutually exclusive.
Pattern probability interpretation: The $2^p$ estimated probabilities quantify joint service-use profiles conditional on covariates — directly relevant for cost projection and policy simulation.
Empirical prior robustness: Bock-Aitkin nonparametric prior is a built-in sensitivity check against the multivariate-normality assumption with minimal additional computation.

Sparse Factor Analytic Probit Extension (Hahn-Carvalho-Scott 2012)

For high-dimensional binary outcomes ( $p$ large relative to $n$ ), Hahn, Carvalho, and Scott (2012) replace the unstructured $\Sigma$ with a low-rank sparse decomposition: $\Sigma = BVB^\top + I$ , where $B$ is a $p \times k$ lower-triangular loading matrix with spike-and-slab priors on each entry, and $V$ is a free $k \times k$ positive-definite factor covariance.

Key extensions relative to Chib-Greenberg (1998):

Identification: Lower-triangular $B$ with positive diagonal (Geweke-Zhou 1996) instead of correlation-constrained $\Sigma$
Sparsity: Spike-and-slab prior $b_{js} \mid \nu_s, q_s \sim q_s \cdot \mathcal{N}(0, \nu_s) + (1-q_s) \cdot \delta_0(b_{js})$ with $\text{Beta}(1,1)$ on $q_s$ for automatic multiplicity correction (Scott-Berger 2006)
Correlated factors: Free $V$ (rather than $V = I_k$ ) absorbs inter-factor correlation, fostering a sparser $B$
Partisanship factor: Adds a scalar factor $\gamma_i$ with truncated-normal identification ( $R^+/D^-$ ) for a known interpretable dimension

See Sparse Factor Probit and Hahn, Carvalho, Scott (2012).

MCMC for High-Dimensional Pick Any/J Data (Edwards-Allenby 2003)

"Pick any/ $J$ " surveys — where respondents check all that apply from $J$ items — yield multivariate binomial responses. The identification challenge is that the $J$ -dimensional latent normal must have unit variances (i.e., a correlation matrix $R$ , not a covariance matrix $\Sigma$ ), but directly sampling $R$ in Markov Chain Monte Carlo (MCMC) is computationally hard for $J > 4$ .

Edwards and Allenby (2003) resolve this by ignoring the restriction during sampling:

Fit an unrestricted covariance matrix $\Omega$ under an inverted-Wishart prior; draw $(\mathbf{y}_k, \Omega, \gamma)$ from the standard three-block Gibbs sampler.
Postprocess each draw: let $C = \mathrm{diag}(\Omega)^{-1/2}$ , then $\hat\alpha = C\gamma$ and $\hat R = C\Omega C'$ .

Key results: (i) $\hat R$ is guaranteed positive definite; (ii) mean squared error (MSE) is lower than tetrachoric correlations at all tested sample sizes (Tables 1–3); (iii) tetrachoric correlations fail entirely when any cell count is zero, while the Bayesian estimator is always identified through the prior. The method handles $J \geq 50$ easily and unlocks standard multivariate analysis (principal component analysis (PCA), canonical correlation, conditional normals) on the estimated latent correlation matrix.

The general principle: For any non-identified likelihood with a proper prior, the posterior concentrates on the identified subspace. MCMC chains on individual non-identified parameters wander indefinitely, but postprocessed functions of those draws that are identified converge normally. Identification shifts from a binary (likelihood identified or not) to a continuous concern about the degree of posterior concentration.

See Edwards and Allenby (2003).

Open Questions

Scalability: Gauss-Hermite quadrature is exponential in $m$ ; for $m>3$ or $p\gg 5$ , sparse grid or Monte Carlo integration alternatives are needed.
Choice of $m$ : deviance testing is sensitive to sample size; information criteria may be preferable for large $n$ .
Sensitivity of sparse factor results to the selection rule for roll-call votes (e.g., conditioning on close votes) is an open applied question.
Standard errors of derived quantities (eigenvalues, eigenvectors, canonical correlations) from the Edwards-Allenby postprocessed draws are not directly available; delta-method or bootstrap extension needed.