1 Introduction
Understanding how the brain represents behaviorally relevant information such as speech is a fundamental question in auditory neuroscience. Speech signals are inherently multidimensional, with subtle variations in frequency, intensity, and timing information distinguishing different categories (Caclin et al., Reference Caclin, McAdams, Smith and Winsberg2005, Reference Caclin, Brattico, Tervaniemi, Näätänen, Morlet, Giard and McAdams2006). These acoustic subtleties need to be robustly encoded by the auditory system and mapped onto existing speech representations (e.g., in native speakers) or emergent representations (e.g., in non-native learners). A major goal for auditory neuroscientists is to decipher how these attributes are differentially represented in the brain and how these neural representations are shaped by different auditory experiences across the lifespan. One critical step to elucidate these questions is to quantify the degree of similarity (or dissimilarity) between neural representations of sounds within the same sensory-perceptual representational space. This step is critical because the degree of dissimilarity between neural representations of speech sounds must reflect linguistically relevant differences across languages. For instance, native speakers of Mandarin Chinese rely on subtle syllable-level pitch changes to decode word meanings, exhibiting a more nuanced neural encoding of pitch contours (or lexical tones) than native speakers of languages that are not tonal (e.g., Bidelman et al., Reference Bidelman, Krishnan and Gandour2011; Feng, Gan, et al., Reference Feng, Gan, Llanos, Meng, Wang, Wong and Chandrasekaran2021; Krishnan et al., Reference Krishnan, Gandour and Bidelman2010; Raizada et al., Reference Raizada, Tsao, Liu, Holloway, Ansari and Kuhl2010).
Multidimensional scaling (MDS; Borg & Groenen, Reference Borg and Groenen2005; Carroll & Arabie, Reference Carroll and Arabie1980; Torgerson, Reference Torgerson1952, Reference Torgerson1958) and its extension to individual difference scaling (INDSCAL; Carroll & Chang, Reference Carroll and Chang1970) are two dimensionality reduction techniques that have been extensively used to quantify the degree of similarity between neural representations of sensory stimuli captured with modern neuroimaging technologies. Given data on dissimilarity matrices between pairs of stimuli (e.g., in our work, speech sounds), MDS and INDSCAL can extract lower-dimensional latent features from which a denoised version of the observed dissimilarities can be reconstructed. In the classical literature, these quantities are usually described in geometric terms: the point coordinates of the stimuli in the low-dimensional space form the configuration, and the coordinate axes are called dimensions. In our Bayesian formulation, we treat these coordinates as unobserved random variables endowed with a prior and a posterior distribution, and we refer to them as latent features to emphasize this latent-variable interpretation. Throughout this article, latent features therefore denote the coordinate vectors representing the stimuli in the low-dimensional space, endowed with a probability model. More specifically, the latent features represent the original stimuli as points in a low-dimensional space such that the pairwise (usually Euclidean) distances computed in the reduced space (i.e., from the latent features) match the observed pairwise dissimilarities as well as possible.
Standard MDS methods, when applied to a single dissimilarity matrix, extract common latent features that allow the reconstruction of that matrix. Since they are restricted to handling two-way pairwise dissimilarity data, when multiple subject-specific dissimilarity matrices are available, they are often reduced to a single group-level matrix via complete pooling (e.g., by averaging the distances element-wise) and analyzed with standard MDS, or analyzed by running separate MDS for each subject, followed by post-hoc alignment of the resulting configurations across subjects. In contrast, INDSCAL can accommodate three-way data comprising pairwise dissimilarities for different individuals, allowing for the extraction of individual latent features and denoised distances by re-weighting the common latent features for each subject (Carroll & Chang, Reference Carroll and Chang1970). INDSCAL can therefore be used to simultaneously achieve dimensionality reduction, data visualization, and, when appropriate, meaningful interpretation of the inferred latent features. MDS-type models have also been extended to alternating least squares scaling (ALSCAL; Takane et al., Reference Takane, Young and De Leeuw1977; Young et al., Reference Young, Takane and Lewyckyj1978) and proximity scaling (PROXSCAL; De Leeuw & Heiser, Reference De Leeuw and Heiser1980), which minimize the standardized residual sum of squares (STRESS; Kruskal, Reference Kruskal1964a, Reference Kruskal1964b; Shepard, Reference Shepard1962), a measure of scaled difference between observed distances and their model reconstructed values. Originally introduced for continuous data, such extensions made MDS more broadly applicable to include non-metric data as well as missing values. MDS has gradually transformed from a descriptive geometric optimization method into a rigorous statistical inference framework, shifting the objective function from minimizing an empirical fit measure, like Kruskal’s STRESS, to maximizing a formal likelihood function, which combines the underlying true distance in the latent MDS space and a random error term (Ramsay, Reference Ramsay1977; Takane & Carroll, Reference Takane and Carroll1981). Asymmetric Mahalanobis distance-based formulations have also been considered (Bradlow & Schmittlein, Reference Bradlow and Schmittlein2000). Furthermore, MDS-type models based on latent utility functions have been extensively developed for binary choice and ordered preference data for psychometric analysis in marketing applications. Here, the utilities are constructed using choice and individual-specific location vectors in the MDS latent space. Two main strategies exist for this construction: the “vector” approach uses the inner product, and the “unfolding” approach uses the squared Euclidean distances (DeSarbo & Cho, Reference DeSarbo and Cho1989; DeSarbo et al., Reference DeSarbo, Kim, Wedel and Fong1998, Reference DeSarbo, Park and Scott2008; Fong et al., Reference Fong, DeSarbo, Park and Scott2010; Park et al., Reference Park, DeSarbo and Liechty2008; Zinnes & Wolff, Reference Zinnes and Wolff1977, etc.).
For the rest of the article, we focus mainly on MDS methods for metric-valued continuous data, aligning with the requirements of our motivating application. More recently, such methods have seen widespread use in the neuroscience literature, becoming a staple technique over the last few decades. The notable appeal of MDS in this field may lie in its ability to process distance or dissimilarity matrices directly, thereby obviating the often-involved step of complex modeling of the raw, high-dimensional data sets. Di Liberto et al. (Reference Di Liberto, O’Sullivan and Lalor2015), Feng et al. (Reference Feng, Yi and Chandrasekaran2019), Feng, Gan, et al. (Reference Feng, Gan, Llanos, Meng, Wang, Wong and Chandrasekaran2021), Khalighinejad et al. (Reference Khalighinejad, da Silva and Mesgarani2017), Llanos et al. (Reference Llanos, German, Gnanateja and Chandrasekaran2021), Mesgarani et al. (Reference Mesgarani, Cheung, Johnson and Chang2014), and Zinszer et al. (Reference Zinszer, Anderson, Kang, Wheatley and Raizada2016) have used MDS to align acoustic and neural representations of sounds in a common space for comparison. In this body of literature, a higher degree of alignment between acoustic and neural representations is interpreted as a signature of a more robust neural encoding of stimulus signals. Here, by “alignment” we mean the correspondence between the space in which acoustic stimuli can be represented and the latent space in which the brain represents them. These works suggest that the degree of alignment between neural and/or acoustic representations of sounds can be investigated using a reduced set of latent features derived from the primary acoustic cues used to perceive sounds. For example, the perception of Mandarin tones seems to be based on two major acoustic dimensions: pitch height (high vs. low pitch) and pitch direction (rising vs. falling pitch) (Chandrasekaran, Gandour, et al., Reference Chandrasekaran, Gandour and Krishnan2007; Gandour, Reference Gandour1978). Some work (Gandour & Harshman, Reference Gandour and Harshman1978; Gandour, Reference Gandour1983) has also suggested a possibly relevant third dimension, namely, the magnitude of the slope in the pitch contour. In the last decades, MDS and INDSCAL have also become popular tools to assess the degree of structural alignment between acoustic and neural representations of speech sounds (Chandrasekaran, Gandour, et al., Reference Chandrasekaran, Gandour and Krishnan2007; Chandrasekaran, Krishnan, et al., Reference Chandrasekaran, Krishnan and Gandour2007; Raizada et al., Reference Raizada, Tsao, Liu, Holloway, Ansari and Kuhl2010). Much of the attractiveness of these tools lies in their ability to represent global dissimilarity patterns between complex neural signals in a way that is easy to visualize and interpret. While MDS is not exclusively conceived as a technique for dimensionality reduction, using a reduced number of dimensions, prior neuroscience work (Feng, Li, et al., Reference Feng, Li, Hsu, Wong, Chou and Chandrasekaran2021; Llanos et al., Reference Llanos, German, Gnanateja and Chandrasekaran2021; Mesgarani et al., Reference Mesgarani, Cheung, Johnson and Chang2014) has been able to capture fine-grained linguistically-relevant differences in sound processing, which suggests that theoretical principles of dimensionality reduction inform the representation of speech sounds in the brain.
Current metric MDS approaches in neuroscience are, however, often used as mere visualization tools. Indeed, joint modeling and uncertainty quantification across groups (e.g., native vs. non-native listeners) are not addressed by the existing methods. Furthermore, most methods do not natively support inference on the number of latent features.
The Bayesian paradigm possesses the machinery to resolve these issues, though they have yet to all be comprehensively addressed in the current literature. Following the foundational work by Oh and Raftery (Reference Oh and Raftery2001), who established a framework for performing inference and uncertainty quantification on latent coordinates via Markov chain Monte Carlo (MCMC), the Bayesian MDS literature has expanded to include scalable extensions for hyperbolic geometries (Liu et al., Reference Liu, Lubold, Raftery and McCormick2024; Praturu & Sharpee, Reference Praturu and Sharpee2022) and massive-scale phylogenetic applications (Holbrook et al., Reference Holbrook, Lemey, Baele, Dellicour, Brockmann, Rambaut and Suchard2021). However, these models generally assume population homogeneity, rendering them insufficient for characterizing the complex neuroplasticity and varied speech-sound encoding found across diverse language groups. While domain-specific adaptations have emerged, such as the bioinformatics model by Nguyen and Holmes (Reference Nguyen and Holmes2017), they are often restricted to a single latent dimension and fail to account for individual or group-specific variation in the underlying features. Even hierarchical approaches designed to capture heterogeneity, such as the model for orchestral interpretations of Beethoven symphonies by Yanchenko and Hoff (Reference Yanchenko and Hoff2020), situate variability at the level of denoised distances rather than within the unobserved latent features that induce those distances. Consequently, existing methods lack the multi-dimensional flexibility and feature-level hierarchical structure necessary to quantify the nuanced, group-specific latent representations required for our motivating neuroscience applications.
Indeed, inference about the latent features is challenging due to identifiability issues; while the reconstructed distances are identifiable, the latent features are unique up to transformations, that is, rotations, translations, and reflections. While these indeterminacies may be overlooked for simple visualization, anchoring the model becomes essential for meaningful interpretation and comparison. Specifically, in Bayesian frameworks, the posterior remains invariant under these transformations, necessitating explicit strategies for valid posterior summarization and uncertainty quantification. Common approaches to resolve this symmetry include imposition of strong priors to favor a specific orientation (e.g., DeSarbo et al., Reference DeSarbo, Kim, Wedel and Fong1998), or utilizing post-hoc alignment to map posterior draws to a common reference configuration (e.g., Oh & Raftery, Reference Oh and Raftery2001); however, such strategies may not fully resolve more complex symmetries, such as signed permutations of the latent dimensions. Alternatively, employing hard coordinate constraints to lock the configuration during estimation can provide a more concrete mathematical resolution (Bradlow & Schmittlein, Reference Bradlow and Schmittlein2000; Park et al., Reference Park, DeSarbo and Liechty2008); however, this rigidity may lead to sampler inefficiency, leading to poor MCMC mixing and slow convergence. Lin and Fong (Reference Lin and Fong2019) addressed identifiability by minimizing a loss function post-MCMC, allowing the sampler to move freely, avoiding the bottlenecks of hard constraints, while ensuring that the final latent features are unique and interpretable.
Building on these prior works, while addressing some of their aforementioned limitations, we propose a flexible Bayesian mixed MDS model for studying dissimilarity across subjects and groups. Our proposal allows performing biologically interpretable inference and straightforward uncertainty quantification for both individual and group-specific distances across stimuli, taking into account the individual idiosyncrasies as well as the heterogeneity between native and non-native Mandarin listeners at the level of the unobserved features. Our proposal thus takes into account two different sources of heterogeneity and does this at the latent feature level and not at the latent distance level. Moreover, by weighting the latent features, our proposal allows the latent axes to be identified uniquely and not be affected by rotation invariance, in contrast to standard MDS and in line with the classical INDSCAL. Importantly, this allows inferring biologically interpretable latent dimensions when such interpretation is meaningful.
Our proposed approach can also infer the effective number of relevant features in a semi-automated data-adaptive way. This is achieved by using cumulative shrinkage priors (Bhattacharya & Dunson, Reference Bhattacharya and Dunson2011) that impose increasing shrinkage on higher-dimensional latent features. By starting with a generous number of features, and then eliminating redundant higher-dimensional features which are shrunk toward near-zero values on-the-fly using an adaptive MCMC algorithm (following Roberts & Rosenthal, Reference Roberts and Rosenthal2009), this approach obviates the need for computationally intensive methods, such as refitting models for different latent dimensions or executing trans-dimensional MCMC steps (e.g., Richardson & Green, Reference Richardson and Green1997). Finally, we employ a post-processing procedure, adapting to recent developments to address similar issues in Bayesian latent factor models (LFMs; Papastamoulis & Ntzoufras, Reference Papastamoulis and Ntzoufras2022), which maintains excellent sampling efficiency while fully resolving the identifiability issues, including translations and signed permutations of the values of the latent dimensions across MCMC samples, enabling meaningful inference on the latent features coherently across groups and subjects.
We apply the proposed model to assess group-level differences in the neural representation of Mandarin Chinese tones between native Mandarin speakers and monolingual native English speakers. Neural representations of Mandarin tones were extracted from a dataset of frequency-following responses (FFRs) to Mandarin tones (Llanos et al., Reference Llanos, Xie and Chandrasekaran2017). The FFR is a scalp-recorded electrophysiological brain component that reflects phase-locked activity from ensembles of neurons along the central auditory nervous system. When the brain is stimulated with a periodic sound, neurons in the auditory system synchronize their oscillatory activity by firing at the same phase of each cycle in the stimulus waveform. This synchronized phase-locked activity is aggregated by the scalp-recorded FFR, which thus mimics the temporal structure of the sound with a high degree of fidelity. Prior cross-linguistic work (Krishnan et al., Reference Krishnan, Xu, Gandour and Cariani2005; Llanos et al., Reference Llanos, Xie and Chandrasekaran2017; Reetzke et al., Reference Reetzke, Xie, Llanos and Chandrasekaran2018) has shown that the FFR reflects language experience-dependent plasticity. Specifically, Mandarin lexical tones are more faithfully represented in the FFR of native speakers of Mandarin Chinese, relative to native speakers of English (e.g., Krishnan et al., Reference Krishnan, Xu, Gandour and Cariani2005). Therefore, Mandarin listeners are expected to convey bigger differences between FFRs to Mandarin tones than non-native speakers of Mandarin Chinese (Llanos et al., Reference Llanos, Xie and Chandrasekaran2017; Reetzke et al., Reference Reetzke, Xie, Llanos and Chandrasekaran2018). Motivated by these neuroscientific experiments and related prior literature, our proposed mixed MDS methodology allows scientists to map the observed pairwise FFR distances to a lower-dimensional common feature space, simultaneously enabling the evaluation of heterogeneity among language groups and subjects, presenting a principled statistical approach for comparing groups and individuals, while also allowing the visualization of the geometry. Notably, in our formulation, the latent axes of the reduced space are all uniquely identifiable, which makes inference and interpretation of the features biologically very meaningful.
An alternative approach to obtaining a low-dimensional representation of the FFR time series might be through traditional LFMs. See, e.g., Aguilar and West (Reference Aguilar and West2000), albeit in a financial time-series context. LFMs assume that observed time series can be represented as linear combinations of latent sources plus noise. However, these assumptions may not hold when relationships are nonlinear, nonstationary, or fundamentally distance-based, as is often considered to be the case with the geometrically complex patterns of brain activity. Additionally, while LFMs can recover latent components, they usually lack meaningful spatial arrangement, making intuitive visualization challenging. In contrast, MDS does not rely on any such restrictive assumptions on a specific generative model for the observed time series. By operating directly on pairwise dissimilarities, MDS can be applied to a wide range of distance and dissimilarity matrices and provides low-dimensional embeddings that summarize complex similarity structures. It provides interpretable, low-dimensional visual embeddings that preserve large-scale geometric patterns, making it especially well-suited for exploratory and visualization-driven neurobiological analyses, which likely contributed to its popularity in such contexts. These features are particularly important for our application as well, where the focus is on characterizing and visualizing differences in underlying auditory processing mechanisms from pairwise dissimilarities between the corresponding FFRs. In particular, our proposal is tailored to dissimilarity data with multiple sources of heterogeneity (subjects and groups) and, motivated by the auditory neuroscientific application, targets low-dimensional Euclidean representations for dissimilarities.
The rest of this article is organized as follows. Section 2 provides additional background on our motivating scientific experiments. Section 3 details our novel multi-group mixed MDS models. Section 4 outlines statistical and computational challenges and solution strategies. In particular, Section 4.1 reports the proposed MCMC algorithm, Section 4.2 describes a strategy to select the number of features adaptively, and Section 4.3 outlines the post-processing algorithm to solve the identifiability issues and allow inference on the latent features. Section 5 discusses the results of the proposed method applied to some synthetic numerical experiments. Section 6 presents the results of the proposed method applied to the aforementioned two language groups’ neural dissimilarities data. Section 7 contains concluding remarks.
2 Neural dissimilarities data
The neural dissimilarity matrices used in our analysis come from the previously published auditory neuroscience study by Llanos et al. (Reference Llanos, Xie and Chandrasekaran2017). The brain activities of
$n=28$
subjects,
$n_{1}=14$
Mandarin speakers and
$n_{2} = 14$
English (non-Mandarin) speakers, were recorded under exposure to different Mandarin tones. Individual distance matrices between stimuli were computed as Euclidean distances between FFRs. As a tonal language, Mandarin Chinese has four syllabic pitch contours or tones that are used to convey different lexical meanings. For instance, the syllable “ma” can be interpreted as “mother,” “hemp,” “horse,” or “scold” depending on whether it is pronounced with high level (T1), low-rising (T2), low-dipping (T3), or high-falling (T4) tones, respectively. These four tones, pronounced by native Mandarin speakers, represented the stimuli. Figure 1 provides context for the associated neural signals underlying our distance matrices. As FFRs reflect fast, phase-locked oscillations (often at hundreds of Hz), they are traditionally presented in the time domain to demonstrate neural tracking of fine temporal structure. In such representations, meaningful information is encoded in the rate of oscillation (neural pitch tracking) rather than peak amplitude. While qualitative cues exist, such as T3’s slower overall rate and T4’s terminal deceleration, these differences are naturally subtle.
Tone neural data: Mean (intra-tone cross-measurements) FFR in the time domain from one Mandarin-speaking and one non-Mandarin-speaking listener. T1, T2, T3, and T4 denote the four Mandarin lexical tones: high level (T1), low-rising (T2), low-dipping (T3), and high-falling (T4).
Figure 1 Long description
The multi-panel line graph is organized into a two-by-four grid. The vertical axis represents f sub s super open parenthesis i close parenthesis with a bar over the f, ranging from negative 0.4 to 0.8. The horizontal axis represents time, ranging from 0 to 6000.
* The top row shows data for a Mandarin speaker (Mand).
* The bottom row shows data for a Non-Mandarin speaker (Non-Mand).
* Columns from left to right represent tones T 1 (red), T 2 (blue), T 3 (green), and T 4 (purple).
In the top row (Mandarin speaker):
* T 1 shows a dense oscillatory signal centered around 0.0 with moderate amplitude.
* T 2 shows similar oscillations with a slight peak near the 2000 time mark.
* T 3 displays consistent amplitude across the time domain.
* T 4 shows a high-amplitude spike at the beginning (near time 500) followed by lower-amplitude oscillations.
In the bottom row (Non-Mandarin speaker):
* T 1 shows high-frequency oscillations with several peaks reaching 0.4.
* T 2 shows a very regular, wave-like pattern with consistent amplitude.
* T 3 shows a distinct peak near time 3000.
* T 4 shows an increasing amplitude trend toward the end of the time window, with a significant peak near time 5500.
This limited visual separability of FFRs justifies our reliance on quantitative, distance-based comparisons of the entire response rather than waveform inspection. Specifically, we focus on analyzing the Euclidean distances between neural representations of Mandarin tones computed from the subjects’ brain activity measured via scalp-recorded FFRs, yielding for each subject an
$S \times S$
(here,
$S=4$
) symmetric distance matrix with zeros on the diagonal. FFRs were extracted from the dataset described in Llanos et al. (Reference Llanos, Xie and Chandrasekaran2017). They were acquired using three Ag-Cl electrodes placed on the vertex of the scalp (active channel), the left mastoid (ground), and the right mastoid (reference channel). The neural activity captured by these electrodes was amplified and digitized with an actiCHamp system and a dedicated preamplifier (EP-preamp, gain 50
$\! \times $
). The FFR dataset included
$1,000$
artifact-free FFR trials per tone (i.e.,
$1,000$
FFRs in response to
$1,000$
repetitions of each tone). To account for different trial-by-trial noise levels in the neural representation of the brain between the two language groups across subjects and time, some rescaling of the observations is necessary (see, e.g., Feng, Li, et al., Reference Feng, Li, Hsu, Wong, Chou and Chandrasekaran2021; Nili et al., Reference Nili, Wingfield, Walther, Su, Marslen-Wilson and Kriegeskorte2014), which helps eliminate internal variability. Specifically, we compute the average FFR across
$1,000$
repeated measurements (i.e., trials) rescaled by the average (across stimuli) of the within (repeated measurements) standard deviations.
Formally, let
$f^{(i)}_{s,m}(t)$
denote the FFR in the time domain, representing the brain activity of subject i at time t under the exposure to stimulus s in the measurement m. We denote by
$\overline {f}^{(i)}_{s}(t) = \sum _{m=1}^{1,000} f^{(i)}_{s,m}(t)/1,000$
the average over repeated measurements (as shown in Figure 1) and by
$\text {SSE}^{(i)}_{s}(t) = \sum _{m=1}^{1,000} \{f^{(i)}_{s,m}(t)-\overline {f}^{(i)}_{s}(t)\}^{2}$
the sum of squared errors at time t for subject i and stimulus s. Thus, as a measure of noise for the measurements at time t in subject i under the exposure to stimulus s, we use the estimated standard deviation
$\widehat {\sigma }_{s}^{(i)}(t) = \{\text {SSE}^{(i)}_{s}(t)/(1,000-1)\}^{1/2}$
; and as a measure of individual dispersion at time t, the average of these standard deviations across stimuli
$\widehat {\sigma }^{(i)}(t)=\sum _{s=1}^{S} \widehat {\sigma }_{s}^{(i)}(t)/S$
. Finally, we rescale the average empirical distances at time t by these estimated noise variances that cannot be explained by the average FFR, that is,
$\widetilde {f}^{(i)}_{s}(t) = \overline {f}^{(i)}_{s}(t)/\widehat {\sigma }^{(i)}(t)$
. By rescaling the FFRs between average brain activities by the corresponding language group’s inherent noise levels in this manner, we allow for biologically meaningful comparisons of the brain’s ability to distinguish different auditory stimuli across subjects.
Figure 2 shows the average distances between stimuli in the two groups. We can see that, although there are group-specific characteristics due to language, for example, the recorded Mandarin speaker distances are better differentiated across pairs of tones than the ones from the non-Mandarin group, there are also strongly shared patterns across these groups, for example, the pairs
$\{2,3\}$
and
$\{2,4\}$
are well-separated pairs of tones while the pairs
$\{1,4\}$
and
$\{1,2\}$
are the closest in both groups.
Tone neural distance data: Mean (intra-group cross-subject) FFR distances between different Mandarin tones, T1, T2, T3, and T4.
Figure 2 Long description
The figure consists of two side-by-side heatmaps labeled Mandarin on the left and Non-Mandarin on the right. Both heatmaps use a four-by-four grid where the x-axis and y-axis are labeled T 1, T 2, T 3, and T 4. The diagonal cells from top-left to bottom-right are colored gray and contain the value 0.
Mandarin Heatmap (Left):
* Row T 1: T 1 is 0, T 2 is 7.9, T 3 is 8, T 4 is 7.3.
* Row T 2: T 1 is 7.9, T 2 is 0, T 3 is 8.6, T 4 is 8.5.
* Row T 3: T 1 is 8, T 2 is 8.6, T 3 is 0, T 4 is 8.
* Row T 4: T 1 is 7.3, T 2 is 8.5, T 3 is 8, T 4 is 0.
The Mandarin group shows higher distance values indicated by darker red shading, particularly between T 2 and T 3 (8.6).
Non-Mandarin Heatmap (Right):
* Row T 1: T 1 is 0, T 2 is 6.1, T 3 is 6.8, T 4 is 6.5.
* Row T 2: T 1 is 6.1, T 2 is 0, T 3 is 6.9, T 4 is 6.9.
* Row T 3: T 1 is 6.8, T 2 is 6.9, T 3 is 0, T 4 is 7.2.
* Row T 4: T 1 is 6.5, T 2 is 6.9, T 3 is 7.2, T 4 is 0.
The Non-Mandarin group shows lower distance values overall, indicated by lighter pink shading, with the highest distance between T 3 and T 4 (7.2).
The individual tone distances show that there is also substantial individual variability within groups. Note also a higher within-group variability in the non-Mandarin group, as expected. The corresponding data are presented alongside the analysis results later in the manuscript.
We are interested in learning the shared geometry of the lower-dimensional latent space representing these tones in the human brain across language groups and individuals, while also assessing how this space varies across these groups and individuals.
3 Mixed multidimensional scaling
Modeling the distances: In MDS, the observed distances between stimuli s and r, namely,
$d_{s,r}, r<s, r=1,\dots , S-1, s=2,\dots ,S$
, are noisy measurements of latent distances
$\delta _{s,r}$
computed in a lower H-dimensional feature space (
$H <S$
). We therefore treat the observed dissimilarities
$d_{j,s,r}^{(i)}$
as noisy, strictly positive measurements of latent denoised distances
$\delta _{j,s,r}^{(i)}$
.
As in Nguyen and Holmes (Reference Nguyen and Holmes2017) and Yanchenko and Hoff (Reference Yanchenko and Hoff2020), we assume a Gamma likelihood for the observed distances to perform our model-based MDS. However, motivated by the application, we allow for individual differences at the level of the latent features. We also consider different variances in the likelihood terms to take into account the heterogeneity of the noise levels between the two groups in our data.
More precisely, we denote by
$d_{j,s,r}^{(i)}$
the distance between the
$s^{th}$
and the
$r^{th}$
stimuli (
$\! r < s$
) for the
$i^{th}$
individual from the
$j^{th}$
group and we assign a Gamma likelihood for the observed distances
In the context of the experiment described in Section 2,
$j=1$
and
$j=2$
refer to the Mandarin-speaking group and the non-Mandarin-speaking group, respectively. The Gamma likelihood is parameterized in terms of mean and variance to ease the interpretation of the latent parameters, and it is centered around the underlying individual denoised distances
$\delta _{j,s,r}^{(i)}$
, which can vary across subjects within and across the two groups. The variance term
$\sigma _{\epsilon ,j}^{2}$
varies across groups, allowing for different noise levels in different groups. The more standard shape and rate representation can be recovered as shape
$=$
mean
$^{2}$
/variance =
$(\delta _{j,s,r}^{(i)})^{2}/\sigma _{\epsilon ,j}^{2}$
and rate
$=$
mean/variance =
$\delta _{j,s,r}^{(i)}/\sigma _{\epsilon ,j}^{2}$
.
The underlying individual distance
$\delta _{j,s,r}^{(i)}$
is a function of H unobserved individual features
$\eta _{j,s,h}^{(i)}$
, such as
where the number of latent features H is smaller than the number of original stimuli S. Here,
$\boldsymbol {\eta }_{j,s}^{(i)} = (\eta _{j,s,1}^{(i)},\dots ,\eta _{j,s,H}^{(i)})^{\mathrm {T}}$
denotes the latent vector with the feature values for the
$s^{th}$
stimulus for the
$i^{th}$
individual belonging to the
$j^{th}$
group.
On the choice of the Gamma likelihood: Note that in our hierarchical framework, the Gamma likelihood for the observed distances
$d_{j,s,r}^{(i)}$
is conditional on latent distances
$\delta _{j,s,r}^{(i)}$
derived from a lower-dimensional feature space. Under this specification, the Euclidean structure enters through the conditional mean of
$d_{j,s,r}^{(i)}$
, while the sampling variances
$\sigma _{\epsilon ,j}^{2}$
accommodate varying noise levels across groups. The likelihood thus implicitly accounts for an error term,
$\epsilon $
, characterizing the variability of the observed distances
$d_{j,s,r}^{(i)}$
around the latent distances
$\delta _{j,s,r}^{(i)}$
, which is itself a function of the underlying feature parameters
$\eta _{j,s,h}^{(i)}$
, as defined in (2). The model for these features
$\eta _{j,s,h}^{(i)}$
and their associated priors are detailed in Sections 3 and 3.2. By operating strictly within this conditional hierarchical framework, we fulfill our inference goals while avoiding the derivation of the induced marginal distribution, which remains analytically intractable.
The choice of the Gamma distribution for the conditional likelihood provides a good framework for modeling dissimilarities on
$(0,\infty )$
, consistent with established practice in Bayesian MDS (e.g., Nguyen & Holmes, Reference Nguyen and Holmes2017; Yanchenko & Hoff, Reference Yanchenko and Hoff2020). Although parametric, the Gamma family offers some advantages over more rigid alternatives: unlike the log-Normal distribution (e.g., Ramsay, Reference Ramsay1977, Reference Ramsay1982), which is constrained by heavy tails and sparse mass near zero, or the truncated-Normal specification (e.g., Oh & Raftery, Reference Oh and Raftery2001), which is often dominated by its behavior at the origin, the Gamma distribution can smoothly transition from Exponential-like to Gaussian-like shapes (as shape
$=$
mean
$^{2}$
/variance increases) while maintaining a robust intermediate tail behavior.
While more flexible nonparametric approaches, such as mixture models, could be considered when larger sample sizes are available, in the present application, the sample size is relatively modest, and the empirical data exhibit means that substantially exceed their variances (Figure 2). In this high-shape regime (shape
$\gg 1$
), the Gamma, log-Normal, and truncated-Normal distributions all converge toward a symmetric, Gaussian-like form, rendering the specific choice of parametric family largely inconsequential. Overall, the parametric Gamma likelihood thus represents a parsimonious and robust choice for our current modeling requirements.
Modeling the latent features: Our modeling efforts concentrate henceforth on flexibly characterizing the latent features
$\eta _{j,s,h}^{(i)}$
. Specifically, we model
$\eta _{j,s,h}^{(i)}$
as a product of latent features shared across groups and subjects,
$\eta _{s,h}$
, and multiplicative coefficients specific to dimension h, group j, and subject i within group j. That is, we let
The random coefficients
$\widetilde {w}_{j,h}^{(i)}$
vary across features h, allowing variation, and thus importance, of the dimensions across individuals between and within groups.
Finally, in order to learn the subjects’ and the two groups’ variations, we specify an inverse gamma distribution with hyperparameters
$a_{w}$
and
$b_{w}$
such that the multiplicative terms have a mean equal to
$1$
and a large variance equal to
$10$
. That is,
Here, with some abuse of notation,
$w_{j,h}^{(i)}$
’s represent the individual-level weights that do not vary by group; however, the index j is still retained to distinguish individuals belonging to different groups (e.g.,
$w_{1,h}^{(1)}$
vs.
$w_{2,h}^{(1)}$
denote the first individual in groups 1 and 2, respectively). Also, we parameterize the gamma distribution in terms of shape and rate hyperparameters to facilitate updating these hyperparameters in the full conditionals. If the multiplicative terms
$w_{j}$
and
$w_{j,h}^{(i)}$
are equal to
$1$
, the MDS model ignores group and individual differences, entailing the same latent feature values for the stimulus s for all subjects within and across the two groups.
On the choice of shared features: Studies on Mandarin lexical tones in both native and non-native listeners have shown their perception to rely on core pitch-related dimensions across language backgrounds (Gandour, Reference Gandour1978, Reference Gandour1983; Gandour & Harshman, Reference Gandour and Harshman1978). Neurophysiological work further demonstrates that these dimensions are robustly represented in early auditory responses, such as FFRs, in both Mandarin and non-Mandarin listeners (Chandrasekaran, Krishnan, et al., Reference Chandrasekaran, Krishnan and Gandour2007).
This motivates our shared-feature MDS framework, where group and individual differences are characterized as shifts in representational geometry rather than the use of distinct feature sets. This approach aligns with the general consensus that these dimensions are universally relevant to auditory processing, while language experience modulates their relative weighting. While relaxing the shared-feature assumption may suit other domains, our model’s use of shared latent features with group and individual-specific weighting provides an interpretable framework, grounded in neurobiological theory, for comparing neural encoding across different language backgrounds.
On the model for the weights: In model (3), the weights
$\widetilde {w}_{j,h}^{(i)}$
are not identifiable from the
$\eta _{s,h}$
parameters without the imposition of scale constraints. We discuss these constraints as part of our post-processing scheme in Section 4.3. Notwithstanding this specific form of non-identifiability, one may consider various strategies to impose further structure and interpretability on the weights. Specifically, we can consider:
-
(a) unrestricted: $\widetilde {w}_{j,h}^{(i)}$
, or -
(b) group-specific decomposition: $\widetilde {w}_{j,h}^{(i)} = w_{j}w_{j,h}^{(i)}$
, or -
(c) group-and-dimension-specific decomposition: $\widetilde {w}_{j,h}^{(i)} = w_{j,h} w_{j,h}^{(i)}$
, or -
(d) group-and-(dimension-within-group)-specific decomposition: $\widetilde {w}_{j,h}^{(i)} = w_{j}w_{j,h}w_{j,h}^{(i)}$
.
As before, in models (b)–(d),
$w_{j,h}^{(i)}$
’s represent individual-level weights that do not vary directly by group j, but the index j is still retained to distinguish individuals from different groups.
Because the likelihood depends solely on latent distances, these specifications are stochastically equivalent, leaving internal structures unidentifiable without further constraints. However, by imposing consistently defined identifying constraints across these parameterizations, the individual weight components can be uniquely recovered from one another. For instance, the marginal effects remain well-defined in terms of the
$\widetilde {w}_{j,h}^{(i)}$
’s across all models: Group-and-dimension-specific effects can be defined as
$\overline {w}_{j,h} = \sum _{i=1}^{n_{j}}\widetilde {w}_{j,h}^{(i)}/n_{j}$
. Under model (c), this simplifies to
$w_{j,h}$
given the constraint
$\sum _{i=1}^{n_{j}}w_{j,h}^{(i)}/n_{j} = 1$
. Likewise, overall group-specific effects can be defined as
$\overline {w}_{j} = \sum _{h=1}^{H}\sum _{i=1}^{n_{j}}\widetilde {w}_{j,h}^{(i)}/(H n_{j})$
. Under model (b), this equals
$w_{j}$
given the constraint
$\sum _{h=1}^{H}\sum _{i=1}^{n_{j}}w_{j,h}^{(i)}/(H n_{j}) = 1$
.
In a Bayesian hierarchical framework, these distinct parameterizations do necessitate tailored prior specifications, which in turn dictate the mechanisms for information borrowing and shrinkage, and the efficiency of the MCMC sampling; generally, more granular decompositions will yield higher posterior uncertainty and require more robust implementation to manage the resulting computational complexity. However, as we will see in Section 4.3, using the constraints described above, the posterior distributions for weight components under alternative specifications of
$\widetilde {w}_{j,h}^{(i)}$
can be readily derived (under implicit priors) from the MCMC samples of a primary model of choice. Consequently, irrespective of the initial specification, the underlying representational geometry remains accessible via straightforward reparameterization of the posterior samples.
Nevertheless, the practical choice between these models can be guided by considerations of meaningful borrowing of information and computational ease and efficiency. In this article, we adopted model (b) because it keeps the MCMC and the post-processing schemes simple, while inducing a higher correlation between individual distances
$\delta $
within the same group. As discussed above, the
$w_{j,h}$
’s can still be obtained using post-processing (see Section 4.3 and Figures 1 and 2 in the Supplementary Material).
3.1 Model identifiability
Our multiplicative individual effects make the latent feature space identifiable under rotations, enabling stronger interpretations of the latent features. To see this, note that reweighting the latent features is equivalent to multiplying the feature matrix by a diagonal matrix of weights. Therefore, if we rotate the axes (aside from the special case of permutations of the axes and the subclass of signed permutations that can be represented as rotation), the weight matrix would not be diagonal anymore, which is not permissible in our model. Therefore, it allows the axes of the latent features to be uniquely identified and unaffected by rotation invariance.
More precisely, let
${\mathbf H}^{(i)}_{j}=((\eta _{j,s,h}^{(i)})) \in \mathbb {R}^{S \times H}$
,
$i=1,\ldots ,n_{j}, \, j=1,\ldots ,J$
be the individual latent feature matrices,
${\mathbf W}_{j}^{(i)}=\text { diag}(\widetilde {w}_{j,1}^{(i)}, \ldots , \widetilde {w}_{j,H}^{(i)}) \in \mathbb {R}^{H\times H}$
,
$i=1,\ldots ,n_{j}, \, j=1,\ldots ,J$
be the corresponding diagonal weight matrices, and
${\mathbf H}=((\eta _{s,h}))\in \mathbb {R}^{S \times H}$
be the shared latent feature matrix. Then, we can rewrite
${\mathbf H}^{(i)}_{j} = {\mathbf H} \mathbf {W}_{j}^{(i)}$
. Now, let
$\mathbf {R} \in \mathbb {R}^{H\times H}$
be a rotation matrix. If we rotate the individual features as
${\mathbf H}^{(i)}_{j} \mathbf {R}^{\mathrm {T}} = {\mathbf H} \mathbf {W}_{j}^{(i)} \mathbf {R}^{\mathrm {T}}$
, it will not change the implied latent distances
$\delta _{j,s,r}^{(i)}$
, but the new individual weight matrices
$\mathbf {W}_{j}^{(i),\text {new}} = \mathbf {W}_{j}^{(i)} \mathbf {R}^{\mathrm {T}}$
are in general not diagonal anymore. This is why the individual weights, which are shared across stimuli, make the latent feature space identifiable in INDSCAL, contrary to standard MDS analysis. See Carroll and Chang (Reference Carroll and Chang1970) for further discussion. Notably, the use of individual-specific dimension weights is not exclusive to INDSCAL; other prominent MDS frameworks, such as ALSCAL and PROXSCAL, can similarly accommodate individual differences by weighting different latent dimensions according to each subject’s unique perceptual profile (see, e.g., Commandeur & Heiser, Reference Commandeur and Heiser1993; Takane et al., Reference Takane, Young and De Leeuw1977).
That said, some non-identifiability persists since the features
$\eta _{j,s,h}^{(i)}$
are still invariant to signed permutations and translations that preserve the latent distances
$\delta $
(even if the axes are identifiable and cannot be arbitrarily rotated as discussed above). We address these issues separately via a post-processing scheme discussed later in Section 4.3.
Finally, we note that the latent distances
$\delta _{j,s,r}^{(i)}$
and group-specific variance terms
$\sigma _{\epsilon ,j}^{2}$
are identifiable in the traditional sense of Basu (Reference Basu2004) as stated in Swartz et al. (Reference Swartz, Haitovsky, Vexler and Yang2004). Specifically, the identifiability of
$\boldsymbol {\theta } =\{\eta _{j,s,h}^{(i)}, \, \sigma _{\epsilon ,j}^{2}, \, i =1,\ldots , \, n_{j},j=1,\ldots , J, \, s=1, \ldots , S, \, h=1, \ldots , H\}$
up to the aforementioned transformations is equivalent to the identifiability of
$\widetilde {{\boldsymbol {\theta }}}=\{\delta _{j,s,r}^{(i)}, \, \sigma _{\epsilon ,j}^{2}, \, i =1,\ldots , \, n_{j},j=1,\ldots , J, \, s=1, \ldots , S, \, r=1, \ldots , S\}$
. And the latter follows from the moment equations
3.2 Prior specification
To learn the shared latent features
$\eta _{s,h}$
, we adapt the multiplicative gamma process (MGP) prior (Bhattacharya & Dunson, Reference Bhattacharya and Dunson2011) for ordinary factor analysis that allows a stochastic ordering of the latent features in diminishing order of their importance.
Specifically, we assume a Gaussian distribution for the shared feature component
$\eta _{s,h}$
whose variance varies across stimuli s and dimension h
The precision parameter of the shared component
$\eta _{s,h}$
follows the MGP
where
$a_{l}=a_{2}$
for
$l \ge 2$
. The MGP prior shrinks the latent features increasingly toward zero as the dimension h increases. This implies that the prior favors a few relevant features, shrinking the rest to zero, while also inducing a probabilistic ranking of the relevant features such that, on average, the first feature explains more variability than the second one, and so on, coherent with the current understanding of the brain’s behavior.
We follow Durante (Reference Durante2017) recommendations for the choice of the hyperparameters of the MGP. In particular, we set
$a_{1}=2$
and
$a_{2}=3$
since they induce the desired prior shrinkage. More precisely, they induce an order in the probability of the dimensions in a neighborhood around zero (Lemma 2 in Durante, Reference Durante2017). Section 4 of the Supplementary Material provides further discussion on the choice of the hyperparameters.
Finally, we put an inverse-gamma prior to learn the group-specific variance terms
3.3 Model reparameterization
We can reparameterize the model in terms of the individual latent features
$\eta _{j,s,h}^{(i)}$
. Integrating out the shared dimensions
$\eta _{s,h}$
, we obtain
This parameterization allows restating the probabilistic statements directly in terms of one of the main quantities of interest, namely, the individual-specific latent features. Note that, by the definition in (3), the
$\eta _{j,s,h}^{(i)}$
’s in (9) are dependent on each other since they share the same latent
$\eta _{s,h}$
’s. That is,
$\big (\eta _{j,s,h}^{(i)} \big )_{h},\, i=1,\ldots ,n_{j},j=1,\ldots ,J,s=1,\ldots ,S$
belong to the span of the vectors
$\eta _{s,h}, s=1,\ldots ,S$
. This is consistent with our goal to infer a common latent feature’s space across groups as well as subjects.
4 Posterior inference
Posterior inference is performed by sampling from an MCMC algorithm that exploits conditional conjugacy of the parameters when available. When sampling from the full conditionals is not straightforward, we exploit adaptive Metropolis schemes (Roberts & Rosenthal, Reference Roberts and Rosenthal2009). The details are discussed in the next section.
4.1 MCMC algorithm
(1) Sample
$\delta _{1}$
from
$p(\delta _{1} \mid \cdots )$
(2) Sample
$\delta _{h}$
from
$p(\delta _{h} \mid \cdots )$
where
$\tau _{l}^{(h)} = \prod _{t=1,\, t\ne h}^{l} \delta _{t}$
, for
$l=1,\ldots H$
.
(3) Sample
$\phi _{s,h}$
from
(4) Sample
$\eta _{s,h}$
from
We perform an adaptive MH step (Roberts & Rosenthal, Reference Roberts and Rosenthal2009) with a random-walk univariate Gaussian proposal.
(5) Sample
$\sigma ^{2}_{\epsilon ,j}$
from
We perform an adaptive MH step with a random-walk Gaussian proposal for
$\text {log}(\sigma ^{2}_{\epsilon ,j})$
.
(6) Sample
$w_{j,h}^{(i)}$
from
We perform an adaptive MH step with a random-walk Gaussian proposal for
$\text {log}(w_{j,h}^{(i)})$
.
(7) Sample
$w_{j}$
from
We perform an adaptive MH step with a random-walk Gaussian proposal for
$\text {log}(w_{j})$
. More precisely, in all adaptive MH steps, we use an adaptive random walk Metropolis kernel with Gaussian proposals, after transforming parameters to have support on
$\mathbb {R}$
when needed. Every
$50$
iterations, we compute, for each parameter, the empirical acceptance rate over the last
$50$
iterations. If this acceptance rate is below the target value
$0.44$
, we decrease the corresponding proposal standard deviation by a factor
$e^{-\delta _m}$
; if it is above
$0.44$
, we increase it by a factor
$e^{\delta _m}$
. Here,
$\delta _m = \min \{0.01, 1/\sqrt {m}\}$
, with m the current MCMC iteration, so the adaptations diminish over time. This scheme targets an acceptance rate of
$0.44$
for all scalar random walk proposals, which is a standard choice for efficient adaptive random walk Metropolis algorithms (Roberts & Rosenthal, Reference Roberts and Rosenthal2009).
4.2 Selecting the number of features
Our Bayesian mixed MDS model allows performing dimensionality reduction by setting
$H < S$
. If we choose a conservative (i.e., larger than needed) upper bound
$H^{+}$
, the model allows us to give relevant importance to a few latent features
$H<H^{+}$
and small importance to the remaining
$H^{+}-H$
latent features. Here, we think of H as the effective dimensionality of the latent feature space (i.e., the number of latent dimensions) so that the contribution from adding additional dimensions in reconstructing a denoised version of the observed distance matrices is negligible. However, running the MCMC for an upper bound
$H^{+}$
larger than needed can be computationally inefficient. Finding the number of relevant latent features H can therefore be helpful in reducing costs. Finding H can also be of inferential interest itself.
To address this, we use a relative Frobenius error
$D(t)$
, similar in spirit to Kruskal’s STRESS, defined as the average, across subjects i and groups j, of the Frobenius distance between the observed
$S \times S$
distance matrices
$((d_{j,s,r}^{(i)}))$
and the corresponding
$S \times S$
denoised latent distance matrices
$((\delta _{j,s,r}^{(i)}))$
sampled at iteration t, divided by the Frobenius norm of the observed distance matrix. We set, according to the noise level in the measurements, a threshold
$D_{T}$
. In particular, we found that
$D_{T}=0.1$
favors a parsimonious latent space that can reconstruct well all the individual observable distances in our motivating neural tone distance experiment, as shown in Figure 8. More generally,
$D_{T}$
should be chosen to reflect the desired level of reconstruction accuracy for the observed distances and, when possible, calibrated using prior knowledge about the signal-to-noise ratio in the specific application. We perform the following steps with probability
$p(t)=\exp \{-(\alpha _{0}+\alpha _{1}t)\}$
at iteration t, where
$\alpha _{0}\ge 0$
and
$\alpha _{1}>0$
such that the adaptations occur often at the beginning of the chains, but decrease in frequency exponentially fast as the chain settles in. At the
$t^{th}$
iteration, if the current latent features are not sufficient to recover the distances well, that is,
$D(t)>D_{T}$
, we set
$H(t+1)=H(t)+1$
and add a latent feature
$\eta _{s,H(t)+1}$
. Otherwise, if
$D(t)<D_{T}$
, we set
$H(t+1)=H(t)-1$
and delete the feature
$\eta _{s,H(t)}$
. When we add a feature, we sample the parameters from the prior distribution.
The adaptive method allows, in a single run, to perform posterior inference on the individual and group-level latent distances,
$\delta _{j,s,r}^{(i)}$
and
$\delta _{j,s,r}$
, together with selecting the number of features, with the convergence of the chain guaranteed by the diminishing probability condition (Roberts & Rosenthal, Reference Roberts and Rosenthal2007). If, as in our motivating auditory neuroscience application, one is interested in performing inference on the actual latent features, we can fix the number of active dimensions H selected in a preliminary run of the MCMC chain or in an initial set of iterations of the chain after burn-in. In this way, we have the same number of features, H, sampled in the final stages of the chain that we can use to perform posterior inference on the latent features after solving the identifiability issues as described in the next section. Finally, as in Bhattacharya and Dunson (Reference Bhattacharya and Dunson2011), we can perform uncertainty quantification on the number of features around the point estimate of H, in our case, the median of the sampled values
$H(t)$
after burn-in, via credible intervals.
Additional evidence on the empirical performance of this adaptive procedure is reported in Section 5 of the Supplementary Material, where simulation studies show how it recovers the true latent dimension H in application-motivated scenarios.
4.3 Post-processing for feature identifiability
In MDS, inference on the latent features
$\eta _{j,s,h}^{(i)}$
is challenging due to identifiability issues—while the reconstructed distances are identifiable, the features themselves are not. As discussed in Section 3.1, our construction overcomes the rotation invariance issue of the latent dimensions that affects standard MDS methods. However, while the set of the latent feature axes are now identifiable, the values of the latent features are still identifiable only up to translations and signed permutations of the axes (e.g., label switching of the dimensions and reflection). This makes the posterior summaries of the sampled
$\eta _{j,s,h}^{(i)}$
’s, for example, the posterior median, computed from the MCMC samples, not very meaningful.
To address this issue, we adapt recent post-processing ideas developed for Bayesian LFMs by Papastamoulis and Ntzoufras (Reference Papastamoulis and Ntzoufras2022). First, we remove translation invariance by centering the shared coordinates so that
$\sum _{s=1}^{S}\eta _{s,h}=0$
for every h, that is,
Under the parameterization
$\widetilde {w}_{j,h}^{(i)} = w_{j} w_{j,h}^{(i)}$
, we then apply the corresponding translation to the individual product coordinates as
so that the fitted distances remain unchanged. If the goal is inference on the latent distances and on the individual latent dimensions inducing them, together with comparison of the individual configurations across subjects and groups, then this centering step, combined with the signed permutation alignment described below, is sufficient.
Indeed, after centering, we resolve the remaining label switching and sign ambiguity by applying signed permutation transformations that align the latent dimensions across posterior draws. These transformations are applied consistently to the shared coordinates
$\eta _{s,h}$
and to the individual latent coordinates
$\eta _{j,s,h}^{(i)}$
, while the positive weights are affected only through the corresponding permutation of the dimension index h. Based on these adjusted samples, we can compute meaningful posterior summaries for the latent geometry, such as posterior medians and credible intervals for the individual latent coordinates.
As discussed in Section 3, for an additional summary of the contribution of each latent dimension within and across groups, it may also be useful to examine the weight structure directly. To this end, we further normalize the fitted weights
$\widetilde {w}_{j,h}^{(i)}$
. First, we remove the per-dimension scale indeterminacy induced by the transformation
which enforces the normalization
$\frac {1}{n}\sum _{j=1}^{J}\sum _{i=1}^{n_{j}}\widetilde {w}_{j,h}^{(i)} = 1$
via the rescaling
$\eta _{\cdot ,h} \leftarrow c_{h}\eta _{\cdot ,h}$
and
$\widetilde {w}_{j,h}^{(i)} \leftarrow \widetilde {w}_{j,h}^{(i)}/c_{h}$
. Then, we reconstruct group-and-dimension-specific summaries as
which implies the identifying constraint
$\frac {1}{n_{j}}\sum _{i=1}^{n_{j}} w_{j,h}^{(i)} = 1$
. Finally, we define the overall group-specific summary as the average of these dimension-specific effects as
In this way, posterior draws obtained under the original model can be re-expressed post hoc in terms of group-and-dimension-specific weights
$w_{j,h}$
and their group-level average
$w_{j}$
. These weights can approach zero if so warranted by the data, providing a practical way to assess the relative and overall utilization of different latent dimensions within and across groups.
Based on these adjusted samples, we can compute meaningful posterior summaries for inference for each subject that are comparable across subjects and groups, for example, median posterior values and credible intervals. Trace plots of the latent features before and after post-processing that solves the identifiability issues for synthetic data are shown in Figure 6. Similar results for the real data application, along with the estimated posterior distributions of the weights
$w_{j,h}$
and
$w_{j}$
, are included in Section 1 of the Supplementary Material.
4.4 MCMC diagnostics
The results reported in this article are all based on
$10^5$
MCMC iterations with the initial 4,000 iterations discarded as burn-in. The remaining samples were further thinned by an interval of
$10$
. We programmed everything in R (R Core Team, 2024). The analyses were performed on a MacBook Pro with Apple M2 Pro chip, 16 GB RAM, using R version 4.5.1. For the real data set, the MCMC algorithm took approximately 17 minutes, and the post-processing algorithm took approximately 11 minutes. The Supplementary Material provides further details and diagnostics. The code is available in the Supplementary Material.
5 Simulation studies
In this section, we discuss the results of some synthetic numerical experiments. In designing the simulation scenarios, we have tried to closely mimic our motivating “brain activities distances between tones” data set. We thus consider distances between
$S=4$
stimuli from
$n=28$
subjects recorded in two groups with cardinalities
$n_{1}=14$
and
$n_{2}=14$
. We set the observed distances
$d_{j,s,r}^{(i)}$
close to values that correspond to the estimated quantities (i.e., the posterior median of
$\delta _{j,s,r}^{(i)}$
) for the real data set, thus simulated distances can be reconstructed from a three-dimensional latent space.
Figure 3 shows the recovery of the distances at the group level. One can, in principle, define corresponding group-level distances based on the shared features
$\eta _{j,s,h}$
. However, due to the nonlinearity of the distance function, these would not correspond to the average of the individual-level distances within the group. For this reason, we believe it is more meaningful to define group-level distances directly as the average of the corresponding individual-level distances, integrating out the random effects from equation (2). Unfortunately, such integrals are not mathematically tractable in closed form, nor empirically calculable since the full populations are not accessible. To circumvent this, in our illustrations, we approximate the group-level distance parameters
$\delta _{j,s,r}$
using the corresponding sample averages of the individual parameters, that is,
$\delta _{j,s,r} = \sum _{i=1}^{n_{j}} \delta _{j,s,r}^{(i)}/n_{j}$
, which serve as good proxies for the true
$\delta _{j,s,r}$
. Figure 3 suggests that the three shared latent features recover these
$\delta _{j,s,r}$
quite well. Moreover, Figure 4 suggests that we recover the denoised version of the empirical distances quite well at the individual level. Figure 5 shows that our method allows us to recover the underlying latent features very efficiently. Here, we align the signs of the true and the estimated latent features for identifiability. When it comes to permuting the dimension labels, however, we leave them undisturbed. We can do this since our method assigns decreasing importance to the dimensions, consistent with the data-generating truth, and we are able to identify the order of the axes after post-processing. Importantly, contrary to standard MDS methods, we do not need to arbitrarily rotate the latent space, since the directions are identifiable.
Results for synthetic data. Posterior medians and 90% credible intervals of the group-specific latent distances
$\delta _{j,s,r}$
. Black squares represent the median observed distances in each group.
Figure 3 Long description
The figure consists of two side-by-side panels labeled Group 1 and Group 2. The y-axis is labeled delta sub j, s, r and ranges from 6 to 8. The x-axis is labeled s r and contains six categories: T 1 dash T 2, T 1 dash T 3, T 1 dash T 4, T 2 dash T 3, T 2 dash T 4, and T 3 dash T 4. Each data point includes a colored circle representing the posterior median, a vertical line representing the 90 percent credible interval, and a black square representing the median observed distance.
* Group 1 panel: Values are generally higher, ranging between 7 and 8.2. The highest distance is at T 2 dash T 3 (purple dot) and the lowest at T 1 dash T 2 (red dot). The trend shows an increase from T 1 dash T 2 to T 1 dash T 3, a slight dip at T 1 dash T 4, a peak at T 2 dash T 3, followed by a gradual decrease through T 3 dash T 4.
* Group 2 panel: Values are lower, ranging between 5.5 and 6.5. The highest distance is at T 2 dash T 3 (purple dot) and the lowest at T 1 dash T 4 (green dot). The trend shows a slight increase from T 1 dash T 2 to T 1 dash T 3, a drop to the minimum at T 1 dash T 4, a peak at T 2 dash T 3, and a decrease through T 3 dash T 4.
In both groups, the black squares (observed medians) align closely with the colored circles (posterior medians).
Results for synthetic data. Observed distances
$d_{j,s,r}^{(i)}$
versus posterior medians of the denoised distances,
$\delta _{j,s,r}^{(i)}$
, reconstructed from three shared dimensions.
Figure 4 Long description
A scatter plot on a light gray grid. The horizontal x axis is labeled d sub j, s, r super i and the vertical y axis is labeled delta sub j, s, r super i. Both axes range from 4 to 16 with major tick marks at intervals of 4. A solid black diagonal line represents the identity line y equals x, starting from the bottom-left corner and extending to the top-right corner. Data points for two groups are plotted directly on this line. Group 1 is represented by red circular dots and Group 2 is represented by blue circular dots. The points are densely clustered between the values of 4 and 10, with more sparse distribution as values increase toward 16. A legend to the right of the plot area identifies the red dots as Group 1 and the blue dots as Group 2.
Results for synthetic data. True latent features
$\eta _{j,s,h}^{(i)}$
versus their estimated posterior medians
$\widehat {\eta }_{j,s,h}^{(i)}$
.
Figure 5 Long description
A multi-panel figure containing three scatter plots arranged horizontally.
* Axes: The horizontal x-axis is labeled eta sub j, s, h super (i) representing true latent features. The vertical y-axis is labeled hat eta sub j, s, h super (i) representing estimated posterior medians. Both axes use a linear scale.
* Data Flow: In all three panels, the data points fall exactly along a solid black diagonal line with a 45-degree slope, indicating a perfect correlation between true and estimated values.
* Panel 1 (Left): Labeled h equals 1. The data range extends from approximately -5 to 11 on both axes.
* Panel 2 (Center): Labeled h equals 2. The data range extends from approximately -8 to 8 on both axes.
* Panel 3 (Right): Labeled h equals 3. The data range extends from approximately -5 to 8 on both axes.
* Legend: Located to the right of the third panel, it identifies two categories of data points: Group 1 represented by red dots and Group 2 represented by blue dots. Both groups are interspersed along the diagonal line in all panels.
Diagnostics for synthetic data. Trace plots of the individual features
$\eta _{j,s,h}^{(i)}$
sampled in the first subject of group 2 pre- and post-processing.
Figure 6 Long description
The figure consists of two side-by-side line graphs with a shared Y-axis and X-axis. The X-axis is labeled iteration, ranging from 0 to 10,000. The Y-axis is labeled eta sub j, s, h super open parenthesis i close parenthesis, ranging from -6 to 6. To the right of the panels is a legend titled sh, listing 12 color-coded categories from T 1-Dim 1 to T 4-Dim 3.
* The left panel, titled pre, shows high volatility and overlapping paths for all 12 features. The lines exhibit significant fluctuations and non-stationary behavior, particularly between iterations 0 and 5,000, where several lines show sharp upward or downward drifts before somewhat stabilizing at wide intervals.
* The right panel, titled post, shows the same 12 features after processing. The lines are now highly stable, horizontal, and tightly constrained within narrow bands. Each color-coded feature occupies a distinct, non-overlapping vertical position on the Y-axis, indicating successful convergence and separation of the sampled features.
MCMC diagnostics for the simulation experiments are similar to those shown for the real data analysis presented in the next section and hence are omitted. Figure 6 shows the trace plots of the individual latent features before and after applying the algorithm to solve the identifiability issue in the first subject. We can see how the procedure described in Section 4.3 identifies the posterior samples of the different latent features, allowing meaningful inference and comparison at the level of the latent features
$\eta _{j,s,h}^{(i)}$
.
Section 5 of the Supplementary Material provides additional simulations to systematically evaluate latent structure recovery in higher dimensions. Specifically, we consider scenarios with
$S \in \{4,8\}$
stimuli and
$H_{\mathrm {true}} \in \{2,3\}$
. Our model demonstrates robust performance, accurately recovering both the true latent distances and the correct latent dimensions across replications.
6 Application to neural tone distances
In this section, we discuss the results produced by our method applied to the two groups’ tone distance data described in Section 2. Our inference goals, we recall, include understanding similarities in terms of stimulus distances between subjects within and across Mandarin and non-Mandarin listeners as well as comparing individual and group latent feature values. The selected number of dimensions H is
$3$
since, after the burn-in, the adaptive criterion discussed in Section 4.2 suggested sampling from
$H=2$
in
$5\%$
of the MCMC iterations and from
$H=3$
in
$95\%$
of the iterations.
Figure 7 shows the posterior point estimates and
$90\%$
credible intervals of the group distances
$\delta _{j,s,r}$
between each of the
$\binom {4}{2}=6$
pairs of tones. We can see a similar ranking (across groups j) of the pairs of tones according to
$\delta _{j,s,r}$
. Indeed, after discarding individual differences, in both groups, the pairs of Mandarin tones
$\{2,3\}$
and
$\{2,4\}$
exhibited the strongest degree of neural dissimilarity, whereas pairs
$\{1,2\}$
and
$\{1,4\}$
exhibited the strongest degree of neural similarity. These findings are corroborated by empirical evidence from the preliminary analysis of the data in Section 2. Note that although the rankings between the point estimates are similar, they are not exactly the same across the two groups, and that there is more uncertainty, quantified via the posterior credible intervals, in the non-Mandarin group due to greater variability across these individuals. Moreover, consistent with recent neuroscience work (Llanos et al., Reference Llanos, Xie and Chandrasekaran2017; Reetzke et al., Reference Reetzke, Xie, Llanos and Chandrasekaran2018), we observe a better separation between neural representations of tones in the group of native listeners of Mandarin Chinese, relative to the group of non-native listeners. Figure 8 shows the posterior point estimates and
$90\%$
credible intervals of the individual distances
$\delta _{j,s,r}^{(i)}$
as well as the corresponding observed distances. We can see that three latent common features reconstruct the observed distances very well. This indicates that the neural encoding of differences between multidimensional neural representations can be captured using three latent dimensions.
Results for real data. Posterior medians and 90% credible intervals of the group latent distances
$\delta _{j,s,r}$
between stimuli.
Figure 7 Long description
A dot plot with a vertical y-axis labeled delta sub j, s, r ranging from 6 to 9 and a horizontal x-axis labeled s r with six categories. A legend on the right indicates that circles represent Mand and triangles represent Non-Mand.
Data points are color-coded by stimulus pair and show a median value with a vertical line representing the 90 percent credible interval. In all cases, Mand values are higher than Non-Mand values.
* T 1-T 2 (red): Mand is at approximately 7.7 and Non-Mand is at 6.3.
* T 1-T 3 (blue): Mand is at 8.1 and Non-Mand is at 6.6.
* T 1-T 4 (green): Mand is at 7.4 and Non-Mand is at 6.4.
* T 2-T 3 (purple): Mand is at 8.6 and Non-Mand is at 6.8.
* T 2-T 4 (orange): Mand is at 8.4 and Non-Mand is at 7.1.
* T 3-T 4 (magenta): Mand is at 8.0 and Non-Mand is at 6.8.
Results for real data. Posterior medians and 90% credible intervals of the individual latent distances
$\delta _{j,s,r}^{(i)}$
between stimuli. Black points represent the observed distances
$d_{j,s,r}^{(i)}$
.
Figure 8 Long description
The grid is organized into two columns labeled Mand and Non-Mand, and six rows representing time-point comparisons: T 1 minus T 2, T 1 minus T 3, T 1 minus T 4, T 2 minus T 3, T 2 minus T 4, and T 3 minus T 4. The x-axis for all plots is labeled observation, numbered 1 through 14. The y-axis represents distance values ranging from 0 to 15, labeled with the symbols delta sub j, s, r super i and d sub j, s, r super i.
In each plot, black dots represent observed distances, while colored dots with vertical error bars represent posterior medians and 90 percent credible intervals for latent distances. The colors vary by row: red for T 1 minus T 2, blue for T 1 minus T 3, green for T 1 minus T 4, purple for T 2 minus T 3, orange for T 2 minus T 4, and magenta for T 3 minus T 4.
Data Trends:
* In the Mand column, observation 10 consistently shows the highest distance peak across all rows, with values reaching near 15. The latent distance intervals generally overlap with the black observed data points.
* In the Non-Mand column, the peaks are more variable. For T 1 minus T 2 and T 1 minus T 3, observation 14 shows a high peak. For T 1 minus T 4, T 2 minus T 4, and T 3 minus T 4, observation 9 shows the most prominent peak. Observation 4 also shows a significant peak in the T 2 minus T 3 comparison.
* Across most observations, the latent distance medians closely track the observed black points, indicating a strong fit, though some slight deviations occur at the highest peaks where the black point may sit at the upper edge of the credible interval.
We also note greater uncertainty in the estimates of the denoised individual distances for the non-Mandarin subjects. Note that our Bayesian model allows us to quantify such uncertainty via posterior credible intervals while also allowing us to perform proper individual and group comparisons. Figure 8 also highlights the presence of some non-Mandarin speakers (e.g., subjects 7 and 11) who struggle to distinguish (i.e., smaller distances) between neural responses to specific pairs of tones (e.g., {T1, T2} and {T2, T3}) more than other subjects in the same group. Importantly, by addressing the identifiability issues, we are also able to perform inference on the group and individual latent feature values that reconstruct the distances. Notably, the model can also be used to infer structural differences at the level of the evoking stimuli with a high degree of specificity. These differences are depicted in Figure 9. In this figure, the first dimension encodes the difference between the neural representations of T4 (high-falling pitch) and the other tones: T1 (high-level pitch), T2 (low-rising pitch), and T3 (low-dipping pitch). Because T4 is the high-falling tone, this dimension reflects the neural encoding of differences in pitch direction (falling vs. level, rising, and dipping pitch). In contrast, the second dimension encodes the difference between the neural representations of T3 and T2. Because T3 and T2 are the low-dipping and low-rising tones, respectively, this dimension captures neural differences in pitch direction for low-onset tones (dipping vs. rising pitch). Lastly, the third dimension encodes the difference between the neural representations of T1 and T3. Because these tones are the ones with higher (T1) and lower (T3) pitch on average, this dimension may reflect the neural encoding of differences in pitch range (high vs. low range). Figure 9 also highlights the fact that the three latent features play the same role in the two groups and that the tones are more distinguished (i.e.,
$\eta _{j,s,h}$
more distant from
$0$
) in all the features in Mandarin native speakers. See also the Supplementary Material for additional plots summarizing these results. Figure 10 shows the positions of different subjects on the feature space within and across the two groups and tones and the associated
$90\%$
credible intervals, presented separately for clarity. We see this at the individual level as well: the first feature is mainly useful for distinguishing tone
$4$
. Aside from a few outlying large distances, the overall reconstruction is robust, confirming that the three-dimensional latent space effectively denoises and summarizes the individual dissimilarity matrices. The second and third features are mainly useful for representing other tones at both the group and individual levels. Moreover, although we note some similarities between Mandarin and non-Mandarin speakers, as expected, there is more individual variability in the latent features for non-native speakers. The individual latent features in Figure 10 also identify a few non-Mandarin speakers who struggle to distinguish the Mandarin tones, especially in the third dimension. Detailed convergence and identifiability diagnostics for the MCMC algorithm for the real data set are provided in Section 3 of the Supplementary Material.
Results for real data. Upper panel: Three-dimensional scatter plot of posterior medians of the group latent feature values
$\eta _{j,s,h}$
. Lower panel: Two-dimensional representation of the posterior medians and 90% credible intervals of the group latent feature values
$\eta _{j,s,h}$
in the different groups.
Figure 9 Long description
The upper panel is a three-dimensional scatter plot. The x-axis is eta sub j, s, 1 ranging from -3 to 5. The y-axis is eta sub j, s, 2 ranging from -6 to 4. The z-axis is eta sub j, s, 3 ranging from -3 to 5. Data points are categorized by Tone (T 1 red, T 2 blue, T 3 green, T 4 purple) and Group (Mand as circles, Non-Mand as triangles). T 1 points are high on the z-axis. T 2 points are low on all axes. T 3 points are near the origin. T 4 points are high on the y-axis.
The lower panel contains three two-dimensional scatter plots arranged horizontally.
1. Left plot: x-axis is eta sub j, s, 1 and y-axis is eta sub j, s, 2. T 3 is in the upper-left, T 2 in the lower-left, T 1 near the center-left, and T 4 on the far right.
2. Middle plot: x-axis is eta sub j, s, 1 and y-axis is eta sub j, s, 3. T 1 is at the top-left, T 2 and T 3 are at the bottom-left, and T 4 is on the far right.
3. Right plot: x-axis is eta sub j, s, 2 and y-axis is eta sub j, s, 3. T 1 is at the top-center, T 2 is at the center-left, T 4 is at the center, and T 3 is at the bottom-right.
All 2D plots include vertical and horizontal error bars representing 90 percent credible intervals.
Results for real data. Upper panel: Three-dimensional scatter plot of posterior medians of the individual latent features
$\eta _{j,s,h}^{(i)}$
between stimuli in the two groups. Lower panels: Two-dimensional representation of the posterior medians and 90% credible intervals of the individual latent feature values
$\eta _{j,s,h}^{(i)}$
in the different groups.
Figure 10 Long description
The visualization consists of three main sections.
1. Top Panel: A 3D scatter plot with axes labeled eta sub j, s, 1 (x-axis, -5 to 15), eta sub j, s, 2 (y-axis, -10 to 10), and eta sub j, s, 3 (z-axis, -6 to 10). Data points are categorized by Tone (T 1 red, T 2 blue, T 3 green, T 4 purple) and Group (Mand as circles, Non-Mand as triangles). T 1 points cluster high on the z-axis. T 4 points cluster high on the x-axis. T 2 and T 3 cluster near the origin.
2. Middle Section: A 2 by 2 grid of 2D scatter plots. The top row plots eta sub j, s, 2 against eta sub j, s, 1. The bottom row plots eta sub j, s, 3 against eta sub j, s, 1. The left column represents the Mand group and the right column represents the Non-Mand group. Each data point includes a cross-shaped 90% credible interval. T 4 (purple) consistently extends furthest right along the x-axis (eta sub j, s, 1).
3. Bottom Section: Two side-by-side 2D scatter plots showing eta sub j, s, 3 on the y-axis versus eta sub j, s, 2 on the x-axis. The left plot is for Mand and the right is for Non-Mand. In both, T 1 (red) is positioned at the top center, T 2 (blue) at the bottom left, and T 3 (green) at the bottom right, forming a triangular distribution around T 4 (purple) at the center.
7 Discussion
In this article, we proposed a novel MDS model for multi-group and multi-subject distances, motivated by auditory neuroscience research into the latent representation of speech sounds. Applying this model to Mandarin tone perception in native speakers and tone-naive English speakers, we characterized latent features that reconstruct denoised empirical distances across both groups and individuals. Our framework identified structural differences in neural coding linked to long-term language exposure, while also revealing shared behaviors across groups.
Our proposed approach advances the MDS literature to enable neuroscientists to perform statistical comparisons across multi-subject and multi-group brain distances by mapping them to a biologically interpretable lower-dimensional common feature space. Our proposal allows borrowing of information across groups and subjects, providing the means for a more comprehensive understanding of how the human brain processes speech signals. Moving beyond its role as a visualization tool, our proposed MDS framework incorporates formal uncertainty quantification and a data-driven approach for determining the optimal number of latent features.
Although our development is motivated by this auditory neuroscience application, the proposed model can be used to model structural differences between neural representations of other sensory stimuli (e.g., visual signals) captured by neuroimaging technologies, such as EEG, fMRI, and MEG. Prior neuroscience work (Kluender et al., Reference Kluender, Stilp and Llanos2019; Lewicki, Reference Lewicki2002; Lotto & Holt, Reference Lotto and Holt2016; Stilp & Kluender, Reference Stilp and Kluender2012) suggests that the neural coding of complex multidimensional sensory stimuli is informed by principles of dimensionality reduction. Our model could also be used to infer the latent neural dimensions needed by the human brain to encode differences between sensory signals and the impact of different kinds of experiences on the number of dimensions. Additionally, our approach can be used to address scientific problems in other fields, such as in bioinformatics to find lower-dimensional, biologically relevant features that summarize the recorded single-cell RNA gene expressions in several biomarkers, for example, to extend classical MDS used in a standard dimensionality reduction pipeline in bioinformatics research (Perraudeau et al., Reference Perraudeau, Risso, Street, Purdom and Dudoit2017). Such extensions might, however, require additional work to tailor our model to the specific scientific problem. Additionally, computations may need to be scaled up when datasets with a much larger number of stimuli are available. Efficient approximations to the posterior, for example, the variational Bayes algorithm developed for the simpler MDS model in Nguyen and Holmes (Reference Nguyen and Holmes2017) or the Hamiltonian Monte Carlo-based sampler using representative subsamples for likelihood and gradient evaluations instead of the full distance matrix proposed recently by Sheth et al. (Reference Sheth, Smith and Holbrook2026), may offer viable computational templates for such high-dimensional extensions of our work.
Some topics of our ongoing exploration include extensions to nonparametric frameworks utilizing more flexible likelihood functions, clustering listeners into latent subgroups, incorporating covariates, accommodating asymmetric distance matrices, and adaptations to longitudinal experiments capturing the dynamic evolution of the latent features.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/psy.2026.10131.
Data availability statement
The code and dissimilarity matrices necessary to reproduce the analyses are provided in the Supplementary Material.
Acknowledgements
We thank an anonymous associate editor and three anonymous reviewers for comments on an earlier version leading to significant improvements.
Funding statement
This work was supported in part by the National Science Foundation grant NSF DMS-1953712 and the National Institute on Deafness and Other Communication Disorders grants R01DC013315 and R01DC015504. G.R. has been partially supported by the Italian grant MUR, PRIN project 2022CLTYP4.
Competing interests
The authors declare none.
Open access
