A Bayesian Approach to Spherical Factor Analysis for Binary Data

The formulation of classical factor models for binary data in terms of random utility functions described in Section 2 lends itself naturally to extensions to more general embedding spaces. In particular, we can generalize the model by embedding the positions $\bm{\beta}_{i}$, $\bm{\psi}_{j}$ and $\bm{\zeta}_{j}$ into a connected Riemannian manifold $\mathcal{D}$ equipped with a metric $\rho:\mathcal{D}\times\mathcal{D}\to M\subseteq\mathbb{R}^{+}$, and then rewriting the utility functions in (2) as

$\displaystyle U_{+}(\bm{\psi}_{j},\bm{\beta}_{i})$

$\displaystyle=-\left\{\rho\left(\bm{\psi}_{j},\bm{\beta}_{i}\right)\right\}^{2}+\epsilon_{i,j},$

$\displaystyle U_{-}(\bm{\zeta}_{j},\bm{\beta}_{i})$

$\displaystyle=-\left\{\rho\left(\bm{\zeta}_{j},\bm{\beta}_{i}\right)\right\}^{2}+\nu_{i,j}.$

(3)

Assuming that the errors $\epsilon_{i,j}$ and $\nu_{i,j}$ are such that their differences $\upsilon_{i,j}=\nu_{i,j}-\epsilon_{i,j}$ are independent with cumulative distribution function $G_{j}$, this formulation leads to a likelihood function of the form

$$p\left(\mbox{\boldmath$Y$}\mid\{\bm{\psi}_{j}\}_{j=1}^{J},\{\bm{\zeta}_{j}\}_{j=1}^{J},\{\bm{\beta}_{i}\}_{i=1}^{I}\right)=\prod_{i=1}^{I}\prod_{j=1}^{J}\theta_{i,j}^{y_{i,j}}\left(1-\theta_{i,j}\right)^{1-y_{i,j}},$$

where $Y$ is the $I\times J$ matrix of observations whose rows correspond to $\mbox{\boldmath$y$}_{i}^{T}$, and

$\displaystyle\theta_{i,j}$

$\displaystyle=G_{j}\left(e\left(\bm{\psi}_{j},\bm{\zeta}_{j},\bm{\beta}_{i}\right)\right),$

$\displaystyle e\left(\bm{\psi}_{j},\bm{\zeta}_{j},\bm{\beta}_{i}\right)$

$\displaystyle=\left\{\rho\left(\bm{\zeta}_{j},\bm{\beta}_{i}\right)\right\}^{2}-\left\{\rho\left(\bm{\psi}_{j},\bm{\beta}_{i}\right)\right\}^{2}.$

(4)

In this paper, we focus on the case in which $\mathcal{D}$ corresponds to the unit hypersphere in $K+1$ dimensions, $\mathcal{S}^{K}$, which is equipped with its geodesic distance $\rho_{K}$, the shortest angle separating two points measured over the great circle connecting them. There are two alternative ways to describe such a space. The first one uses a hyperspherical coordinate system and relies on a vector of $K$ angles, $\mbox{\boldmath$\phi$}=(\phi_{1},\ldots,\phi_{K})\in[-\pi,\pi]\times[-\pi/2,\pi/2]^{K-1}$. The second one uses the fact that $\mathcal{S}^{K}$ is a submanifold of $\mathbb{R}^{K+1}$, and relies on a vector of $K+1$ Cartesian coordinates $\mbox{\boldmath$x$}=(x_{1},\ldots,x_{K+1})$ subject to the constraint $\left\|\mbox{\boldmath$x$}\right\|=1$. The two representations are connected through the transformation

	$\displaystyle x_{1}=$	$\displaystyle\cos\phi_{1}\cos\phi_{2}\cos\phi_{3}\cdots\cos\phi_{K-1},$		(5)
	$\displaystyle x_{2}=$	$\displaystyle\sin\phi_{1}\cos\phi_{2}\cos\phi_{3}\cdots\cos\phi_{K-1},$
	$\displaystyle x_{3}=$	$\displaystyle\sin\phi_{2}\cos\phi_{3}\dots\cos\phi_{K-1}$
	$\displaystyle\vdots$
	$\displaystyle x_{K}=$	$\displaystyle\sin\phi_{K-1}\cos\phi_{K},$
	$\displaystyle x_{K+1}=$	$\displaystyle\sin\phi_{K}.$

While the hyperspherical coordinate representation tends to be slightly more interpretable and directly reflects the true dimensionality of the embedding space, Cartesian coordinates often result in more compact expressions. For example the distance metric $\rho_{K}$ can be simply written as $\rho_{K}(\mbox{\boldmath$x$},\mbox{\boldmath$z$})=\arccos\left(\mbox{\boldmath$x$}^{T}\mbox{\boldmath$z$}\right)$. Furthermore, the Cartesian coordinate representation will be key to the development of our computational approaches. Notation-wise, in the sequel, we adopt the convention of using Greek letters when representing points in $\mathcal{S}^{K}$ through their angular representation, and using Roman letters when representing them through embedded Cartesian coordinates.

Lastly, we discuss the selection of the link function. $G_{j}$ must account for the fact that, when $\mathcal{S}^{K}$ is used as embedding space and $\rho_{K}$ as the distance metric, the function $e\left(\bm{\psi}_{j},\bm{\zeta}_{j},\bm{\beta}_{i}\right)$ introduced in (4) has as its range the interval $[-\pi^{2},\pi^{2}]$. While various choices are possible, we opt to work with the cumulative distribution of a shifted and scaled beta distribution,

$\displaystyle G_{j}(z)=G_{\kappa_{j}}(z)$

$\displaystyle=\int_{-\pi^{2}}^{z}\frac{1}{2\pi^{2}}\frac{\Gamma(2\kappa_{j})}{\Gamma(\kappa_{j})\Gamma(\kappa_{j})}\Bigg{(}\frac{\pi^{2}+z}{2\pi^{2}}\Bigg{)}^{\kappa_{j}-1}\Bigg{(}\frac{\pi^{2}-z}{2\pi^{2}}\Bigg{)}^{\kappa_{j}-1},$

$\displaystyle z$

$\displaystyle\in[-\pi^{2},\pi^{2}].$

(6)

Note that $1/\sqrt{\kappa_{j}}$, which controls the dispersion of the density $g_{j}$ associated with $G_{j}$, plays an analogous role to the one that $\sigma_{j}$ played in Section 2. In fact, note that

$\displaystyle\lim_{\kappa_{j}\to\infty}\frac{g_{\kappa_{j}}(z)}{\sqrt{\frac{2\kappa_{j}+1}{2\pi^{5}}}\exp\left\{-\frac{2\kappa_{j}+1}{\pi^{4}}z^{2}\right\}}=1.$

(7)

This observation, together with the fact that $\rho_{K}(\mbox{\boldmath$x$},\mbox{\boldmath$z$})=\arccos\left(\mbox{\boldmath$x$}^{T}\mbox{\boldmath$z$}\right)\approx\left\|\mbox{\boldmath$x$}-\mbox{\boldmath$z$}\right\|$ for small values of $\left\|\mbox{\boldmath$x$}-\mbox{\boldmath$z$}\right\|$, make it clear that the projection of the likelihood of a spherical model in $\mathcal{S}^{K}$ on its tangent bundle is very close to a version of the Euclidean factor model in $\mathbb{R}^{K}$ described in (1) in which the link function is the probit link. We expand on this connection in Section 3.3

Finally, a short note on the connection between the likelihood functions associated with two models defined in $\mathcal{S}^{K}$ and $\mathcal{S}^{K-1}$, respectively. Let $\bm{\beta},\bm{\psi}\in\mathcal{S}^{K}$ and $\bm{\beta}^{*},\bm{\psi}^{*}\in\mathcal{S}^{K-1}$ be (the angular representations of) two pairs of points in spherical spaces of dimension $K$ and $K-1$, respectively. It is easy to verify that, when $\beta_{K}=\psi_{K}=0$ as well as $\beta_{k}=\beta_{k}^{*}$ and $\psi_{k}=\psi_{k}^{*}$ for all $k=1,\ldots,K-1$, then $\rho_{K}(\bm{\beta},\bm{\psi})=\rho_{K-1}(\bm{\beta}^{*},\bm{\psi}^{*})$. Hence, the likelihood function of a spherical factor model in which the latent positions are embedded in $\mathcal{S}^{K}$ and for which $\beta_{i,K}=\psi_{j,K}=\zeta_{j,K}=0$ for all $i$ and $j$, is identical to the likelihood function that would be obtained by directly fitting a model in which the latent positions are embedded in $\mathcal{S}^{K-1}$. This observation will play an important role in the next Section, which is focused on the defining prior distributions for the latent positions.

As we discussed in the introduction, a key challenge involved in the implementation of the spherical factor models we just described is selecting priors for the latent positions $\bm{\psi}_{1},\ldots,\bm{\psi}_{J}$ and $\bm{\zeta}_{1},\ldots,\bm{\zeta}_{J}$ and $\bm{\beta}_{1},\ldots,\bm{\beta}_{J}$. This challenge is specially daunting in settings where estimating the dimension $K$ of the latent space is of interest.

To illustrate the issues involved, consider the use of the von Mises–Fisher family as priors for the latent positions. The Hausdorff density with respect to the uniform measure on the sphere for the von Mises–Fisher distribution on $\mathbb{R}^{d}$ is given by

$\displaystyle p_{\mathcal{H}}(\mbox{\boldmath$x$}\mid\bm{\eta},\omega)$

$\displaystyle=\frac{\omega^{d/2-1}}{(2\pi)^{d/2}I_{d/2-1}(\omega)}\exp\left\{\omega\bm{\eta}^{T}\mbox{\boldmath$x$}\right\},$

$\displaystyle\left\|\mbox{\boldmath$x$}\right\|$

$\displaystyle=1,$

where $I_{\nu}(\cdot)$ is the modified Bessel function of order $\nu$, $\bm{\eta}$ satisfies $\left\|\bm{\eta}\right\|=1$ and represents the mean direction, and $\omega$ is a scalar precision parameter. The von Mises–Fisher distribution is a natural prior in this setting because it is the spherical analogue to the multivariate Gaussian distribution with diagonal, homoscedastic covariance matrix that is sometimes used as a prior for the latent positions in traditional factor models. Indeed, note that

$$\lim_{\omega\to\infty}\frac{\frac{\omega^{d/2-1}}{(2\pi)^{d/2}I_{d/2-1}(\omega)}\exp\left\{\omega\bm{\eta}^{T}\mbox{\boldmath$x$}\right\}}{\left(\frac{\omega}{2\pi}\right)^{d/2}\exp\left\{-\frac{\omega}{2}(\mbox{\boldmath$x$}-\bm{\eta})^{T}(\mbox{\boldmath$x$}-\bm{\eta})\right\}}=1.$$

We are interested in the behavior of the prior on

$$\theta_{i,j}=G_{\kappa_{j}}\left(\left\{\rho_{K}\left(\bm{\zeta}_{j},\bm{\beta}_{i}\right)\right\}^{2}-\left\{\rho_{K}\left(\bm{\psi}_{j},\bm{\beta}_{i}\right)\right\}^{2}\right)$$

(recall Equation (4)) induced by a von Mises–Fisher prior with mean direction $\mathbf{0}$ and precision $\omega$ for $\bm{\beta}_{i}$, and independent von Mises–Fisher priors with mean direction $\mathbf{0}$ and precision $\tau$ for both $\bm{\psi}_{j}$ and $\bm{\zeta}_{j}$. Because, $\bm{\psi}_{j}$ and $\bm{\zeta}_{j}$ are assigned the same prior distributions, it is clear from a simple symmetry argument that $\mathrm{E}\left(\theta_{i,j}\right)=1/2$ for all values of $\omega$, $\tau$ and $\kappa_{j}$. On the other hand, Figure 1 shows the value $\mathrm{Var}\left(\theta_{i,j}\right)$ as a function of the embedding dimension $K$ for various combinations of $\omega$, $\tau$ and $\kappa_{j}$. Note that, in every case, it is apparent that $\lim_{K\to\infty}\mathrm{Var}\left(\theta_{i,j}\right)\to 0$, which implies that $\theta_{i,j}$ converges in distribution to a point mass at $1/2$ as the number of latent dimensions grows.

This behavior, which is clearly unappealing, is driven by two features. First, the von Mises–Fisher prior has a single scale parameter for all its dimensions. Secondly, the surface area of the $K$ dimensional sphere, $\frac{2\pi^{K/2}}{\Gamma(K/2)}$, tends to zero as $K\to\infty$. In fact, not only the von Mises-Fisher distribution, but also other widely used distributions on the sphere such as the Fisher-Bingham and the Watson distributions suffer from the same drawback. Dryden (2005) proposes to overcome this phenomenon by scaling the (common) concentration parameter of the von Mises-Fisher according to the number of dimensions. However, this approach is unappealing in our setting for two reasons. First, it does not appear to overcome the degeneracy issues just discussed. Secondly, under this approach changes in $K$ would fundamentally change the nature and structure of the prior distribution. In the sequel, we discuss a novel class of flexible priors distributions on $\mathcal{S}^{K}$ that allow for different scales across various dimensions and can be used to construct priors for $\theta_{i,j}$ that do not concentrate on a point mass as the number of dimensions increase.

We build our new class of prior distributions in terms of the vector of angles $\bm{\phi}=(\phi_{1},\ldots,\phi_{K})\in[-\pi,\pi]\times[-\pi/2,\pi/2]^{K-1}$ associated with the hyperspherical coordinate system. In particular, we assign each component of $\phi$ a (scaled) von Mises distribution centered at 0, so that

$\displaystyle p(\bm{\phi}\mid\bm{\omega})=\left(\frac{1}{2\pi}\right)^{K}2^{K-1}\frac{1}{I_{0}(\omega_{1})}\exp\left\{\omega_{1}\cos\phi_{1}\right\}\prod_{k=2}^{K}\frac{1}{I_{0}(\omega_{k})}\exp\left\{\omega_{k}\cos 2\phi_{k}\right\},$

(8)

where $\bm{\omega}=(\omega_{1},\ldots,\omega_{K})$ is a vector of dimension-specific precisions. We call this prior the spherical von Mises distribution, and denote it by $\mbox{\boldmath$\phi$}\sim\mbox{SvM}(\mbox{\boldmath$\omega$})$.

One key feature of this prior is that the parameters $\omega_{1},\ldots,\omega_{K}$ control the concentration of each marginal distribution around their mean. In particular, when $\omega_{K}\to\infty$, the marginal distribution of $\phi_{K}$ becomes a point mass at zero, and $p(\bm{\phi}\mid\bm{\omega})$ concentrates all its mass in a greater nested hypersphere of dimension $K-1$ (see Figure 2 for an illustration).

A second key feature of this class of priors is that

$\displaystyle\lim_{\mbox{\boldmath$\omega$}\to\mathbf{\infty}}\frac{\left(\frac{1}{2\pi}\right)^{K}2^{K-1}\frac{1}{I_{0}(\omega_{1})}\exp\left\{\omega_{1}\cos\phi_{1}\right\}\prod_{k=2}^{K}\frac{1}{I_{0}(\omega_{k})}\exp\left\{\omega_{k}\cos 2\phi_{k}\right\}}{\left(\frac{1}{2\pi}\right)^{K/2}2^{K-1}\left\{\prod_{k=1}^{K}\omega_{k}\right\}^{1/2}\exp\left\{-\frac{1}{2}\left[\omega_{1}\phi_{1}^{2}+4\sum_{k=2}^{K}\omega_{k}\phi_{k}^{2}\right]\right\}}=1.$

(9)

Somewhat informally, this means that for large values of the precision parameters, the prior behaves like a zero-mean multivariate normal distribution with covariance matrix $\mbox{diag}\{1/\omega_{1},1/(4\omega_{2}),\ldots,1/(4\omega_{K-1})\}$.

Finally, we note that the density spherical von Mises distribution can be written in terms of the Cartesian coordinates $\mbox{\boldmath$x$}=(x_{1},\ldots,x_{K+1})^{T}$ associated with $\bm{\phi}$ (recall Equation (5)). More specifically, the density of the Hausdorff measure with respect to the uniform distribution on the sphere associated with (8) is given by

$$p(\mbox{\boldmath$x$}\mid\bm{\omega})=\frac{1}{\prod_{k=1}^{K}\sqrt{\sum_{t=1}^{k+1}x_{t}^{2}}}\frac{1}{2\pi I_{o}(\omega_{1})}\exp\left\{\omega_{1}\frac{x_{1}}{\sqrt{x_{1}^{2}+x_{2}^{2}}}\right\}\\ \left\{\prod_{k=2}^{K}\frac{1}{\pi I_{o}(\omega_{k})}\right\}\exp\left\{-\sum_{k=2}^{K}\omega_{k}\left(2\frac{x_{k+1}^{2}}{\sum_{t=1}^{k+1}x_{t}^{2}}-1\right)\right\},\quad\quad\mbox{\boldmath$x$}^{T}\mbox{\boldmath$x$}=1.$$

(10)

See the online supplementary materials for the derivation.

Coming back to our spherical factor model, we use the spherical von Mises distribution as priors for the $\bm{\psi}_{j}$s, $\bm{\zeta}_{j}$s and $\bm{\beta}_{i}$s. In particular, we set

	$\displaystyle\bm{\psi}_{i}$	$\displaystyle\sim\mbox{SvM}(\tau,2^{2}\tau,3^{2}\tau,\ldots,K^{2}\tau)$		(11)
	$\displaystyle\bm{\zeta}_{j}$	$\displaystyle\sim\mbox{SvM}(\tau,2^{2}\tau,3^{2}\tau,\ldots,K^{2}\tau)$		(12)
	$\displaystyle\bm{\beta}_{j}$	$\displaystyle\sim\mbox{SvM}(\omega,2^{2}\omega,3^{2}\omega,\ldots,K^{2}\omega)$		(13)

independently across all $i=1,\ldots,I$ and $j=1,\ldots,J$. By setting up the priors on the latent positions to have increasing marginal precisions, we can avoid the degeneracy issues that arose in the case of the von Mises-Fisher priors. This is demonstrated in Figure 3, which shows the value of $\mathrm{Var}\left(\theta_{i,j}\right)$ as function of the latent dimension $K$ for various combinations of $\omega$, $\tau$ and $\kappa_{j}$ under (11-13). Note that, in every case, the variance of the prior decreases slightly with the addition of the first few dimensions, but then seems to quickly stabilize around a constant that is determined by the value of the hyperparameters. Figure 4 explores in more detail the shape of the implied prior distribution on $\theta_{i,j}$, as well as the effect of various hyperparameters. Note that the implied prior can either be unimodal (which typically happens for relatively large values $\omega$ and $\tau$), or trimodal (which happens for low values of $\omega$ and $\tau$). The intuition behind the formulation in Equations (11-13) comes from the fact that these polynomial rates ensures that $\mathrm{Var}\left(\theta_{i,j}\right)$ converges to a constant with the number of dimensions under the Gaussian approximation in Equation (9).

In addition to avoiding degeneracy issues, using an increasing sequence for the marginal precision of all three sets of $\bm{\psi}_{j}$s, $\bm{\zeta}_{j}$s and $\bm{\beta}_{i}$s means that the prior encourages the model to concentrate most of the mass on a embedded sub-sphere of $\mathcal{S}^{K}$ (recall the discussion at the end of Section 3.1). Hence, we can think of $K$ as representing the highest intrinsic dimension allowed by our prior, and as the combination of $\omega$ and $\tau$ as controlling the “effective” prior dimension of the space, $K^{*}\leq K$. Another implication of this prior structure is that we can think of our prior as encouraging a decomposition of the (spherical) variance that is reminiscent of the principal nested sphere analysis introduced in Jung et al. (2012).

In Section 3.1 we argued that the spherical model of dimension $K$ contains all other spherical models of lower dimension as special cases. Similarly, the Euclidean space of dimension $K+1$ includes all spherical spaces of dimension $K$ or lower. It is also true that the $K$-dimensional Euclidean model discussed in Section 2 (with a probit link) can be seen as a limit case of our spherical model on $\mathcal{S}^{K}$.

The argument relies on projecting the spherical model onto the tangent bundle of the manifold, and making $\omega\to\infty$, $\tau\to\infty$ and $\kappa_{j}\to\infty$ while keeping $\omega/\tau$ and $\omega/\kappa_{j}$ constant for all $j$. As mentioned in previous sections, under these circumstances, (i) the link function $G_{\kappa_{j}}$ defined in Equation (6) projects onto the probit link function, (ii) the geodesic distance between the original latent positions converges to the Euclidean distance between their projections onto the tangent space, and (iii) the spherical von-Mises priors defined in Equations (11-13) project onto Gaussian priors.

This observation is important for at least three reasons. Firstly, it can guide prior elicitation, specially for the precision parameters $\omega$ and $\tau$. Secondly, spherical models can be seen as interpolating (in terms of complexity) between Euclidean models of adjacent dimensions. In particular, note that the likelihood for Euclidean models involves $\mathcal{O}\left(\{I+2J\}K\right)$ parameters (since the $\sigma_{j}$s are not identifiable and typically fixed) while the likelihood of the spherical model relies on $\mathcal{O}\left(\{I+2J\}K+J\right)$ parameters, and finally, $\mathcal{O}\left(\{I+2J\}K\right)<\mathcal{O}\left(\{I+2J\}K+J\right)<\mathcal{O}\left(\{I+2J\}\{K+1\}\right)$. Thirdly, it enhances the interpretability of the model. In particular, the reciprocal of the precision, $\frac{1}{\omega}$ can be interpreted as a measure of sphericity of the latent space that can be directly contrasted across datasets.

The likelihood function for the spherical factor model discussed in Section 3.1 is invariant to simultaneous rotations of all latent positions. However, the structure of the priors in (11-13) induces weak identifiability for the $\mbox{\boldmath$\psi$}_{j}$s, $\mbox{\boldmath$\zeta$}_{j}$s, and $\bm{\beta}_{i}$s in the posterior distribution. Reflections, which are not addressed by our prior formulation, can be accounted for by fixing the octant of the hypersphere to which a particular $\bm{\beta}_{i}$ belongs. In our implementation, this constraint is enforced a posteriori by post-processing the samples generated by the Markov chain Monte Carlo algorithm described in Section 4. On the other hand, while the scale parameters $\sigma_{1},\ldots,\sigma_{J}$ are not identifiable in the Euclidean model, the analogous precisions $\kappa_{1},\ldots,\kappa_{J}$ are identifiable in the spherical model. This is because $\mathcal{S}^{K}$ has finite measure, and therefore the likelihood is not invariant to rescalings of the latent positions. More generally, partial Procrustes analysis (which accounts for invariance to translations and rotations, but retains the scale) can be used to postprocess the samples to ensure identifiability.

Completing the specification of the model requires that we assign hyperpriors to $\omega$, $\tau$ and $\kappa_{1},\ldots,\kappa_{J}$. The precisions $\omega$ and $\tau$ are assigned independent Gamma distributions, $\omega\sim\mbox{Gam}(a_{\omega},b_{\omega})$ and $\tau\sim\mbox{Gam}(a_{\tau},b_{\tau})$. Similarly, the concentration parameters associated with the link function are assumed to be conditionally independent and given a common prior, $\kappa_{j}\sim\mbox{Gam}(c,\lambda)$, where $\lambda$ is in turn given a conditionally conjugate Gamma hyperprior, $\lambda\sim\mbox{Gam}(a_{\lambda},b_{\lambda})$. The parameters $a_{\omega}$, $b_{\omega}$, $a_{\tau}$, $b_{\tau}$, $a_{\lambda}$, $b_{\lambda}$ and $c$ for these hyperpriors are assigned to strongly favor configurations in which $\beta_{i,1}\in[-\pi/2,\pi/2]$ (which is consistent with the assumption that a Euclidean model is approximately correct), and so that the induced prior distribution on $\theta_{i,j}$ is, for large $K$, close to a $\mbox{Beta}(1/2,1/2)$ distribution (the proper, reference prior for the probability of a Bernoulli distribution). Various combinations of hyperparameters satisfy these requirements, most of which lead to priors on $\theta_{i,j}$ that are trimodal. We recommend picking one in which the hyperpriors are not very concentrated (e.g., see Figure 5) and studying the sensitivity of the analysis to the prior choice. In our experience, inferences tend to be robust to reasonable changes in the hyperparameters (see Section 5.1.1).

We conducted a simulation study involving four distinct scenarios to evaluate our spherical model. For each scenario, we generate one data set consisting of $J=700$ items and $I=100$ subjects (similar in size to the roll call data from the U.S. Senate presented in Section 5.2). In our first three scenarios, the data is simulated from spherical factor models on $\mathcal{S}^{2}$, $\mathcal{S}^{3}$ and $\mathcal{S}^{5}$, respectively. In all cases, the subject-specific latent positions, as well as the item-specific latent positions, are sampled from spherical von-Mises distributions where all component-wise precisions are equal to 2. Note that this data generation mechanism is slightly different from the model we fit to the data (in which the precision of the components increase with the index of the dimension). For the fourth scenario, data was simulated from a Euclidean probit model (recall Section 2) in which the intercepts $\mu_{1},\ldots,\mu_{J}$, the discrimination parameters $\bm{\alpha}_{1},\ldots,\bm{\alpha}_{J}$, and latent traits $\bm{\beta}_{1},\ldots,\bm{\beta}_{I}$ are sampled from standard Gaussian distributions.

In each scenario, both spherical and Euclidean probit factor models of varying dimension are fitted to the simulated datasets. The left column of Figure 6 shows the value of the Deviance Information Criteria (DIC) (Spiegelhalter et al., 2002, 2014; Gelman et al., 2014) as a function of the dimension $K$ of the fitted model’s latent space for each of the four scenarios in our simulation. In our models, the DIC is defined as:

$$DIC=\ell\left(\mathrm{E}\{\bm{\Theta}\mid\mbox{data}\}\right)-2\mathrm{Var}\left(\ell\left(\bm{\Theta}\right)\mid\mbox{data}\right),$$

where $\bm{\Theta}=[\theta_{i,j}]$ is the matrix of probabilities given by $\theta_{i,j}=G_{\kappa_{j}}\left(\left\{\rho(\bm{\zeta}_{j},\bm{\beta}_{i})\right\}^{2}-\left\{\rho(\bm{\psi}_{j},\bm{\beta}_{i})\right\}^{2}\right)$ for the spherical model and $\theta_{i,j}=\Phi\left(\mu_{j}+\alpha_{j}\beta_{i}\right)$ for the Euclidean probit model, $\ell\left(\bm{\Theta}\right)=\sum_{i=1}^{I}\sum_{j=1}^{J}\left\{y_{i,j}\log\theta_{i,j}+(1-y_{i,j})\log(1-\theta_{i,j})\right\}$, and the expectations and variances are computed with respect to the posterior distribution of the parameters (which are in turn approximated from the samples generated by our Markov chain Monte Carlo algorithm). The first term in the DIC can be understood as a goodness-of-fit measure, while the second one can be interpreted as a measure of model complexity.

From Figure 6, we can see that DIC is capable of identifying the presence of sphericity in the underlying latent space, as well as its correct dimension. In Scenarios 1 to 3, the optimal model under DIC is correctly identified as a spherical model with the right dimension. Furthermore, as we would expect, the optimal Euclidean model in each case has one additional dimension (i.e., the dimension of the space corresponds to the lowest-dimensional Euclidean space in which the true spherical latent space can be embedded). For Scenario 4, DIC correctly selects a Euclidean model in three dimensions, followed very closely by a spherical space of the same dimension. Again, these results match our previous discussion about the ability of a spherical model to approximate a Euclidean model.

The second column of Figure 6 shows the in-sample predictive accuracy of the models in each scenario. This accuracy is closely linked to the goodness-of-fit component of the DIC. Note that, from the point of view of this metric, both the spherical and Euclidean models have about the same performance. Furthermore, in both cases, as the dimension increases, the predictive accuracy increases sharply at first, but quickly plateaus once we reach the right model dimension, generating a clear elbow in the graph. Similar elbows can be seen in the third column of Figure 6, which shows the amount of variance associated with each component of a principal nested spheres decomposition (Jung et al., 2012) of the subject’s latent positions for the $\mathcal{S}^{10}$. This decomposition can be interpreted as a version of principal component analysis for data on an spherical manifold. Note, however, that the elbow in Scenario 4 is much less sharp than the elbows we observed in Scenarios 1, 2 and 3.

To get a better sense of the relationship between the spherical and Euclidean models, Figure 7 presents histograms of the posterior samples for the hyperparameters $\omega$, $\tau$ and $1/\lambda$ of our spherical model for Scenario 2 (left column) and Scenario 4 (right column). We show these two scenarios because the true dimension of the latent space is 3 in both cases. Note that the posterior distributions of these hyperparameters concentrate on very different values depending on whether the truth corresponds to a Euclidean or a spherical model. Furthermore, in all cases the posterior distributions are clearly different from the priors.

Finally, Figure 8 presents posterior estimates for the subject’s latent positions for Scenario 1. To facilitate comparisons between the truth and the estimates, we plot each dimension separately. While the first dimension seems to be reconstructed quite accurately (the coverage rate for the 95% credible intervals shown in the left panel is 98%, with an approximate 95% confidence interval of $(0.953,1.000)$), the reconstruction of the second component is less so (the coverage in this case is 68%, with an approximate 95% confidence interval of $(0.578,0.782)$).

Factor models are widely used for the analysis of voting records from deliberative bodies (e.g., see Jackman, 2001, Clinton et al., 2004 and Bafumi et al., 2005). In this context, the embedding space provides an abstract representation for the different policy positions, and legislator-specific latent traits are interpreted as their revealed preferences over one or more policy dimensions (e.g., economic vs. social issues in a two dimensional setting, see Moser et al., 2020) or as their ideological preferences in a liberal-conservative scale (in one-dimensional settings, see for example Jessee, 2012).

While Euclidean policy spaces often provide an appropriate description of legislator’s decision-making process, they might not be appropriate in all settings. The idea that spherical policy spaces might be appealing alternatives to Euclidean ones dates back at least to Weisberg (1974), who provides a number of examples and notes that “circular shapes may be expected for alliance structures and for vote coalitions where extremists of the left and right coalesce for particular purposes”. Spherical policy spaces also provide a mechanism to operationalize the so-called “horseshoe theory” (Pierre, 2002; Taylor, 2006, pg. 118), which asserts that the far-left and the far-right are closer to each other than they are to the political center, in an analogous way to how the opposite ends of a horseshoe are close to each other.

In this Section we use our spherical factor models to investigate the geometry of policy spaces in two different U.S. Senates, the 102^nd (which met between January 3, 1991 and January 3, 1993, during the last two years of George H. W. Bush’s presidency) and the 115^th (which met between January 3, 2017 and January 3, 2019, during the first two years of Donald Trump’s presidency). The data we analyze consists of the record of roll call votes, and it is publicly available from https://voteview.com/. Before presenting our results, we provide some additional background into the data and the application. Roll call votes are those in which each legislators votes ”yea” or ”nay” as their names are called by the clerk, so that the names of senators voting on each side are recorded. It is important to note that not all votes fall into this category. Voice votes (in which the position of individual legislator are not recorded) are fairly common. This means that roll call votes represent a biased sample from which to infer legislator’s preferences. Nonetheless, roll call data is widely used in political science. The U.S. Senate is composed of two representatives from each of the 50 States in the Union. The total number of Senators included in each raw dataset might be slightly higher because of vacancies, which are typically filled as they arise. While turnover in membership is typically quite low, these changes lead to voting records with (ignorable) missing values on votes that came up during the period in which a Senator was not part of the chamber. Additionally, missing values can also occur because of temporary absences from the chamber, or because of explicit or implicit abstentions. As is common practice, we exclude from our analysis any senators that missed more than 40% of the votes cast during a given session, which substantially reduces the number missing values. The remaining ones, which are a small percentage of the total number of votes, were treated in our analysis as if missing completely at random. While this assumption is not fully supported by empirical evidence (e.g., see Rodríguez & Moser, 2015), it is extremely common in applications. Furthermore, the relatively low frequency of missing values in these datasets suggests that deviations from it will likely have a limited impact in our results. See table 1 for a summary of the features of the two post processed datasets.

As in Section 5.1, we fit both Euclidean and spherical factor models of varying dimensions to each of these two datasets. The left column of Figure 13 shows the DIC values associated with these models. For the 102^nd Senate, we can see that a Euclidean model of dimension 3 seems to provide the best fit overall while, among the spherical models, a model of dimension 3 also seems to outperform. On the other hand, for the 115^th Senate, a two-dimensional spherical model seems to provide the best fit to the data, closely followed by a circular (one-dimensional spherical) model. The right column of Figure 13 presents the principal nested sphere decomposition associated with the posterior mean configuration estimated by an eight-dimensional spherical model on each of the two Senates. The differences are substantial. In the case of the 115^th Senate (for which DIC suggests that the geometry of the latent space is indeed spherical), the first component of the decomposition explains more than 95% of the variability in the latent space, and the first three components explain close to 99%. On the other hand, for the 102^nd Senate, the first three spherical dimensions explain less than 70% of the total variability.

At a high level, these results are consistent with the understanding that scholars of American politics have of contemporaneous congressional voting patterns. While congressional voting in the U.S. has tended to be at least two-dimensional for most of its history (with one dimension roughly aligning with economic issues, and the other(s) corresponding to a combination of various social and cultural issues), there is clear evidence that it has become more and more unidimensional over the last 40 years (e.g., see Poole & Rosenthal, 2000). Increasing polarization has also been well documented in the literature (e.g., see Jessee, 2012 and Hare & Poole, 2014). The rise of extreme factions within both the Republican and Democratic parties willing to vote against their more mainstream colleagues, however, is a relatively new phenomenon that is just starting to be documented (see Ragusa & Gaspar, 2016, Lewis, 2019a, Lewis, 2019b) and can explain the dominance of the spherical model for the 115^th Senate voting data.

Figure 14 presents histograms of the posterior samples for the hyperparameters $\omega$, $\tau$ and $1/\lambda$ for each of the two Senates. Note that, as with the simulations, the posterior distributions differ substantially from the priors. Furthermore, the model prefers relatively smaller values of $\omega$ and relatively larges values of $\tau$ and $1/\lambda$ for the 115^th Senate when compared with the 102^nd.

To conclude this Section, we focus our attention on the one-dimensional versions of the Euclidean and spherical models. While one-dimensional models are not preferred by the DIC criteria in our datasets, researchers often use them in practice because the resulting ranking of legislators often reflects their ordering in a liberal-conservative (left-right) scale. Figure 15 compares the median rank order of legislators estimated under the one-dimensional Euclidean and one-dimensional spherical (circular) models for the 102^nd (left panel) and the 115^th (right panel) Senates. To generate the ranks under the circular model, we “unwrap” the circle and compute the ranks on the basis of the associated angles (which live in the interval $[-\pi,\pi]$). For the 102^nd Senate (where DIC would prefer the one-dimensional Euclidean model over the circular model), the median ranks of legislators under both models are very similar. This result provides support to our argument that the spherical model can provide a very good approximation to the Euclidean model when there is no evidence of sphericity in the data. On the other hand, for the 115^th (where the circular model would be preferred over the one-dimensional Euclidean model), the ranking of Republican legislators differs substantially between both models. Digging a little bit deeper, we can see that the five legislators for which the ranks differ the most are Rand Paul (KY), Mike Lee (UT), Jeff Flake (AZ), Bob Corker (TN) and John Neely Kennedy (LA), all known for being hard line conservatives, but also for bucking their party and voting with (left wing) Democrats on some specific issues.

Let $\mu_{j}\sim N(0,1/2)$, $\alpha_{j,k}\sim N(0,1/2)$, $\beta_{j,k}\sim N(0,6/[\pi k]^{2})$, and $z_{i,j}(K)=\mu_{j}+\sum_{k=1}^{K}\alpha_{j,k}\beta_{i,k}$. Because $z_{i,j}(K)$ is the sum of $K+1$ random variables with finite second moments, a simple application of the central limit theorem indicates that $\{z_{i,j}(K)-\mathrm{E}(z_{i,j}(K))\}/\mathrm{Var}(z_{i,j}(K))$ converges in distribution to standard normal distribution as $K\to\infty$. Now, note that $(z_{i,j}(K))=0$ by construction, and that,

	$\displaystyle\mathrm{Var}\left\{z_{i,j}(K)\right\}$	$\displaystyle=\mathrm{Var}(\mu_{j})+\sum_{k=1}^{K}\mathrm{Var}(\alpha_{j,k})\mathrm{Var}(\beta_{i,k})$
		$\displaystyle=\frac{1}{2}+\frac{1}{2}\sum_{k=1}^{K}\frac{6}{\pi^{2}}\frac{1}{k^{2}}=\frac{1}{2}+\frac{1}{2}\frac{6}{\pi^{2}}\frac{\pi^{2}}{6}=1$

As a consequence, $\theta_{i,j}=\Phi\left(z_{i,j}(K)\right)$ converges in distribution to a uniform distribution in $[0,1]$.

The gradient with respect to the log of Hausdorff measure may appear forbidding to derive, especially the prior component of the full conditional because the expression of the gradient changes with the predetermined maximum dimension. Hence, we develop an automatic method to facilitate the computation of the gradient. Interestingly, the GHMC algorithm projects all the gradient to the tangent space through the momentum update, which allows the gradient computed with or without the constraint to be the same. Nevertheless, it is generally recommended to derive the gradients under the unit norm constraint since they are invariant to the reparameterization with respect to the constraint.

$\displaystyle\nabla\log p_{\mathcal{H}}(\mbox{\boldmath$x$}_{\beta_{i}}\mid\dots)=\nabla\log p_{\mathcal{H}_{l}}+\nabla\log p_{\mathcal{H}_{\pi}}+\nabla\log p_{J},$

(14)

where $\log p_{\mathcal{H}_{l}}$, $\log p_{\mathcal{H}_{\pi}}$ and $\nabla\log p_{J}$ is the gradient of log-likelihood, log of prior and log of Jacobian respectively with respect to the Hausdorff measure. For $\emph{K}=1$, the model becomes the circular factor model and we focus on the derivation for $K>1$ in this paper. Let $x_{\beta_{i,1}},\dots,x_{\beta_{i,K}}$ be the independent random variable and $\mbox{\boldmath$x$}_{\beta_{i,K+1}}$ be the dependent random variable. As a result, the gradient associated with the last dimension $x_{\beta_{i,K+1}}$ is 0. The first $K$ components of the gradient from the prior and Jacobian are shown as follows,

$$\nabla\log p_{\mathcal{H}_{\pi}}+\nabla\log p_{\mathcal{H}_{J}}=4\omega_{K}\begin{bmatrix}x_{1}\\ x_{2}\\ \vdots\\ x_{K}\end{bmatrix}+\omega_{1}x_{2}(x_{1}^{2}+x_{2}^{2})^{-\frac{3}{2}}\begin{bmatrix}x_{2}\\ -x_{1}\\ 0\\ \dots\end{bmatrix}\\ +4\begin{bmatrix}x_{1}\\ x_{2}\\ \vdots\\ x_{K}\end{bmatrix}\circ\Sigma_{\text{-}2}\begin{bmatrix}\omega_{2}x_{3}^{2}/\left(\sum\limits_{t=1}^{3}{x_{t}^{2}}\right)^{2}\\ \vdots\\ \omega_{K-2}x_{K-1}^{2}/\left(\sum\limits_{t=1}^{K-1}{x_{t}^{2}}\right)^{2}\\ \omega_{K-1}x_{K}^{2}/\left(\sum\limits_{t=1}^{K}{x_{t}^{2}}\right)^{2}\end{bmatrix}-4\begin{bmatrix}0\\ 0\\ \frac{\omega_{2}x_{3}}{\sum\limits_{t=1}^{3}{x_{t}^{2}}}\\ \vdots\\ \frac{\omega_{K-1}x_{K}}{\sum\limits_{t=1}^{K}{x_{t}^{2}}}\end{bmatrix}-\begin{bmatrix}x_{1}\\ x_{2}\\ \vdots\\ x_{K}\end{bmatrix}\circ\Sigma_{\text{-}1}\begin{bmatrix}\frac{1}{\sum\limits_{t=1}^{2}{x_{t}^{2}}}\\ \frac{1}{\sum\limits_{t=1}^{3}{x_{t}^{2}}}\\ \vdots\\ \frac{1}{\sum\limits_{t=1}^{K}{x_{t}^{2}}}\\ \end{bmatrix},$$

(15)

where $\Sigma_{-t}$ is an upper triangular matrix of size $\emph{K}-t$ without the first $t$ columns and all non-zero entries are 1. As a special case when $K=2$, $\Sigma_{-2}$ becomes a 0 matrix and $\Sigma_{-1}$ becomes scalar 1. Once the maximum dimension is specified, $\Sigma_{-1}$, $\Sigma_{-2}$ can be generated and all the gradients can then be vectored easily during the implementation without explicitly deriving the expression for different maximum dimension. Recall that, $\mbox{\boldmath$x$}_{\beta_{i}},\mbox{\boldmath$x$}_{\zeta_{j}}$ and $\mbox{\boldmath$x$}_{\psi_{j}}$ share the same prior and hence the corresponding gradient expression will also be the same.