Bayesian Semiparametric Longitudinal
Inverse-Probit Mixed Models
for Category Learning

Minerva<EMAIL_ADDRESS>
Jacie R. McHaney2<EMAIL_ADDRESS>
Bharath<EMAIL_ADDRESS>
Abhra Sarkar3<EMAIL_ADDRESS>

1Department of of Mathematics and Statistics,

Indian Institute of Technology,

Kanpur, UP 208016, India

2Department of Communication Science and Disorders,

University of Pittsburgh,

4028 Forbes Tower, Pittsburgh, PA 15260, USA

3Department of Statistics and Data Sciences,

University of Texas at Austin,

105 East 24th Street D9800, Austin, TX 78712, USA

Abstract

Understanding how adult humans learn to categorize can shed novel insights
into the mechanisms underlying experience-dependent brain plasticity. Drift-
diffusion processes are popular in such contexts for their ability to mimic
underlying neural mechanisms but require data on both category responses and
associated response times for inference. Category response accuracies are,
however, often the only reliable measure recorded by behavioral scientists to
describe human learning. Building carefully on drift-diffusion models with
latent response times, we derive a novel biologically interpretable class of
‘inverse-probit’ categorical probability models for such data. The model,
however, presents significant identifiability and inference challenges. We
address these challenges via a novel projection-based approach with a symmetry
preserving identifiability constraint that allows us to work with conjugate
priors in an unconstrained space. We adapt the model for group and individual
level inference in longitudinal settings. Building again on the model’s latent
variable representation, we design an efficient Markov chain Monte Carlo
algorithm for posterior computation. We evaluate the method’s empirical
performances through simulation experiments. The method’s practical efficacy
is illustrated in applications to longitudinal tone learning studies.

Key Words: Category learning, B-splines, Drift-diffusion models, Functional
models, Inverse Gaussian distributions, Longitudinal mixed models, Speech
learning

Short/Running Title: Longitudinal Inverse-Probit Mixed Models

Corresponding Author: Abhra Sarkar<EMAIL_ADDRESS>
##  1 Introduction

##### Scientific background.

Categorization decisions are important in almost all aspects of our lives \-
whether it is a friend or a foe, edible or non-edible, the word /bat/ or
/hat/, etc. The underlying cognitive dynamics are being actively studied
through extensive ongoing research (Smith and Ratcliff, 2004; Heekeren et al.,
2004; Gold and Shadlen, 2007; Schall, 2001; Purcell, 2013; Glimcher and Fehr,
2013).

In typical multi-category decision tasks, the brain accumulates sensory
evidence in order to make a categorical decision. This accumulation process is
reflected in the increasing firing rates at local neural populations
associated with different decisions. A decision is taken when neural activity
in one of these populations reaches a particular threshold level. The decision
category that is finally chosen is the one whose decision threshold is crossed
first (Gold and Shadlen, 2007; Brody and Hanks, 2016). Changes in evidence
accumulation rates and decision thresholds can be induced by differences in
task difficulty and/or cognitive function (Cavanagh et al., 2011; Ding and
Gold, 2013). Decision making is also regulated by demands on both speed and
accuracy of the task (Bogacz et al., 2010; Milosavljevic et al., 2010).

Understanding the neural activity patterns for different decision
alternatives, including their similarities and differences, is a key
scientific interest in modeling brain mechanisms underlying decision making.
Statistical approaches with biologically interpretable parameters that further
allow probabilistic clustering of the parameters associated with different
competing choices can facilitate such inference, the parameters clustering
together indicating similar neural behavior and difficulty levels.

##### Drift-diffusion models.

A biologically interpretable joint model for decision response accuracies and
associated response times is obtained by imitating the underlying neural
evidence accumulation mechanisms using latent drift-diffusion processes racing
toward their respective boundaries, the process reaching its boundary first
producing the final observed decision and the time taken to reach this
boundary giving the associated response time (Figure 1) (Usher and McClelland,
2001). The drift and the boundary parameters then jointly explain the dynamics
of decision-making, including the speed-accuracy trade-off (quick vs more
accurate decisions, etc.).

(a) Drift-diffusion processes for tone learning when data on both response
categories and response times are available.

(b) Drift-diffusion processes for tone learning with identifiability
restrictions when only data on response categories are available.

Figure 1: Drift-diffusion model for tone learning. Here we illustrate a single
trial with input tone T1 ($s=1$) that was eventually correctly identified.
Panel (a) shows a process whose parameters can be inferred from data on both
response categories and response times. Here, after an initial $\delta_{s}$
amount of time required to encode an input signal $s$ (here T1), the evidences
in favor of different response categories $d$ (here tones T1, T2, T3, T4)
accumulate according to latent Wiener diffusion processes (red, blue, green
and purple) with drifts $\mu_{d,s}$. The decision $d^{\prime}$ (here T1) is
eventually taken if the underlying process (here the red one) is the first to
reach its decision boundary $b_{d,s}$. Panel (b) shows a process with
additional identifiability restrictions (for all $d$ and $s$, $\delta_{s}=0$,
$b_{d,s}=b$ fixed, and $\sum_{d=1}^{d_{0}}\mu_{d,s}=d_{0}$) considered in this
article which can be inferred from data on response categories alone.

The literature on drift-diffusion processes for decision making is vast but is
mostly focused on simple binary decision scenarios with a single latent
diffusion process with two boundaries, one for each of the two decision
alternatives (Ratcliff, 1978; Ratcliff et al., 2016; Smith and Vickers, 1988;
Ratcliff and Rouder, 1998; Ratcliff and McKoon, 2008). Multi-category drift-
diffusion models with separate latent processes, one for each decision
category, are mathematically more easily tractable (Usher and McClelland,
2001; Brown and Heathcote, 2008; Leite and Ratcliff, 2010; Dufau et al., 2012;
Kim et al., 2017) but the literature is sparse and focused only on simple
static designs. Learning to make categorization decisions is, however, a
dynamic process, driven by perceptual adjustments in our brain and behavior
over time. In laboratory settings, category learning is thus often studied in
longitudinal experiments (Reetzke et al., 2018). Paulon et al. (2020)
addressed the need for sophisticated statistical methods for such settings,
developing an inverse Gaussian distribution based multi-category longitudinal
drift-diffusion mixed model.

Crucially, measurements on both the final decision categories and the
associated response times are needed to estimate the drift and the boundary
parameters. Unfortunately, however, researchers often only record the
participants’ decision responses as their go-to measure of categorization
performance, ignoring the response times (Chandrasekaran et al., 2014; Filoteo
et al., 2010). Additionally, eliciting accurate response times can be
methodologically challenging, e.g., in the case of experiments conducted
online (Roark et al., 2021) or when the response times from
participant/patients are unreliable due to motor deficits (Ashby et al.,
2003). Participants may also be asked to delay the reporting of their
decisions so that delayed physiological responses that relate to decision
making can be accurately measured (McHaney et al., 2021). In such cases, the
reported response times may not accurately relate to the actual decision times
and hence cannot be used in the analysis. One such data set that motivated the
research being reported here is described below.

Figure 2: Description of PTC1 data: Proportion of times the response to an
input tone was classified into different tone categories over blocks across
different subjects, each for the four input tones (as indicated in the
header). The thick line represents the median performance and the shaded
region indicates the corresponding middle $30\%$ quantiles across subjects.
High-flat tone responses are shown in red; low-rising in blue; low-dipping in
green; and high-falling in purple.

##### The PTC1 data set.

The PTC1 (pupillometry tone categorization experiment 1) data set is obtained
from a Mandarin tone learning study conducted at the Department of
Communication Science and Disorders, University of Pittsburgh (McHaney et al.,
2021). There are four Mandarin tones that all sound like ‘ma’ but have
different lexical meanings depending on the pitch contrasts at the syllable
level. A group of native English speaking younger adult participants learned
to categorize monosyllabic Mandarin tones in a training task. During a single
trial of training, an input tone was presented and the participants were
instructed to categorize the tone into one of four categories via a button
press on a keyboard. Corrective feedback in the form of “Correct” or “Wrong”
was then provided on screen. A total of $n=$ 28 participants completed across
$T=6$ blocks of training, each block comprising $L=40$ trials. Figure 2 shows
the inter-quartile range of the proportion of times the response to an input
tone was classified into different tone categories over blocks across
different subjects, each for the four input tones.

Pupillometry was also recorded during each trial as a measure of cognitive
processes during learning. The pupil response unfolds slowly over time,
therefore, trial events were delayed to allow time for the pupil to dilate.
For example, rather than allowing participants to make a button press
immediately after presentation of the input tone, there was a four-second
delay from the start of the input tone to the response prompt screen where
participants made their category decision via button press. As the subjects
may have already taken their decisions within this period, the recorded
response times are not meaningful measures of their actual decision times.

Conventional drift-diffusion analysis that require data on response times such
as in Paulon et al. (2020) therefore do not apply here. The focus of this
article is to see if the drift-diffusion parameters can still be meaningfully
recovered from the PTC1 data on input-output tone categories alone.

The analysis we present in this paper excludes the pupillometry measurements.
Relating drift-diffusion parameters to measures of cognitive effort such as
pupillometry is another challenging problem that we plan to pursue separately
elsewhere.

##### The inverse-probit model.

Categorical probability models that build on latent drift-diffusion processes
can be useful in providing biologically interpretable inference in data sets
comprising input-output categories but no response times like the PTC1 data
set. To our knowledge, however, the problem has never been considered in the
literature before. We aim to address this remarkable gap in this article.

By integrating out the latent response times from the joint inverse Gaussian
drift-diffusion model for response categories and associated response times in
Paulon et al. (2020), we can arrive at a natural albeit overparametrized model
for the response categories. We refer to this class of models as the ‘inverse-
probit’ categorical probability model. This inverse-probit model serves as the
starting point for the methodology presented in this article but, as we
describe below, it also comes with significant and unique statistical
challenges not encountered in the original drift-diffusion model.

##### Statistical challenges.

While scientifically desirable, unfortunately it is also mathematically
impossible to infer both the drifts and the boundaries in the inverse-probit
model from data only on the decision accuracies. We must thus have to keep the
values of either the drifts or the boundaries fixed and focus on inferring the
other.

Even when we fix either the drift or the decision boundaries, the problem of
overparametrization persists. In the absence of response times, only the
information on relative frequencies, that is empirical probabilities of taking
a decision is available. As the total probability of observing any of the
competing decisions naturally adds up to one, the identifiability problem
remains for the chosen main parameters of interest, and appropriate remedial
constraints needs to be imposed.

Setting an arbitrarily chosen category as the reference provides a simple
solution widely adopted in categorical probability models but comes with
serious limitations, including breaking the symmetry of the problem,
potentially making posterior inference sensitive to the specific choice of the
reference category (Burgette and Nordheim, 2012; Johndrow et al., 2013).

By breaking the symmetry of the problem, a reference category also
additionally makes it difficult to infer potential clustering of the model
parameters. To see this, consider decision making problem with $d_{0}$
categories, with a logistic model for the probabilities $p_{s,d}={\tt
logistic}(\beta_{s,d}),~{}s,d=1,\ldots,d_{0}$, of choosing the $d^{th}$ output
category for the $s^{th}$ input category. By setting the $s^{th}$ category as
reference, one can then cluster the probabilities of incorrect decision
choices, $p_{s,d}$, $d\neq s$. However, it is not clear how to directly
compare the probabilities across different input categories, e.g., how to test
the equality of $p_{1,1}$ and $p_{2,2}$.

Finally, while coming up with solutions for the aforementioned issues, we must
also take into consideration the complex longitudinal design of the
experiments generating the data. Whatever strategy we devise, it should be
amenable to a longitudinal mixed model analysis that ideally allows us to
estimate the smoothly varying longitudinal trajectories of the parameters as
the participants learn over time, accommodate participant heterogeneity, and
compare the estimates at different time points within and between different
input categories.

##### Our proposed approach.

As a first step toward addressing the identifiability issues and related
modeling challenges, we keep the boundaries fixed but leave the drift
parameters unconstrained. The decision to focus on the drifts is informed by
the existing literature on such models cited above where the drifts are almost
always allowed more flexibility. The analysis of Paulon et al. (2020) also
showed that the boundaries remain relatively stable over time, between input-
response category combinations and between participants; it is primarily the
variations in the drift trajectories that explain learning over time, the
heterogeneity between input-response category combinations as well as the
heterogeneity between participants.

As a next step toward establishing identifiability, we apply a ‘sum to a
constant’ restriction on the drifts so that symmetry is maintained in the
constrained model.

Implementation of this restriction brings in significant challenges. One
possibility is to set a prior on the constraint space. However, posterior
computation for such priors would be much complicated in drift-diffusion
models. Instead, we conduct inference with an unconstrained prior on the drift
parameters and project samples from the induced posterior to a positive
simplex through a minimal distance mapping.

To adapt the categorical probability model so obtained to a longitudinal mixed
model setting, we then assume that the drift parameters comprise input-
response decision specific fixed effects and subject specific random effects,
modeling them flexibly by mixtures of locally supported B-spline bases (de
Boor, 1978; Eilers and Marx, 1996) spanning the length of the longitudinal
experiment. These effects are thus allowed to evolve flexibly as smooth
functions of time (Ramsay and Silverman, 2007; Morris, 2015; Wang et al.,
2016) as the participants get more experience and training in their assigned
decision tasks.

We take a Bayesian route to estimation and inference. Carefully exploiting
conditional prior-posterior conjugacies as well as the latent variable
construction of the model, we design an efficient Markov chain Monte Carlo
(MCMC) based algorithm for approximating the posterior, where sampling the
latent response times for each observed response category greatly simplifies
the computations.

We evaluate the numerical performance of the proposed approach in extensive
simulation studies. We apply our method to the PTC1 data described above.
These applications illustrate the efficacy of our method in providing insights
into how the drift parameters that relate to the rates of accumulation of
evidence in the brain evolve over time, differ between input-output category
combinations, differ between individuals, etc.

This article thus differs in many obvious ways from Paulon et al. (2020) where
response categories and response times were both observed, and therefore, the
drift and boundary trajectories could be modeled jointly with no
identifiability issues. In contrast, the current work is for data on response
categories alone but no response times, leading us to the inverse-probit
categorical probability model likelihood which, with its many identifiability
issues, brings in daunting unique challenges to performing statistical
inference, confining us, unlike Paulon et al. (2020), only to inferring the
drift parameters on a relative scale, achieved via a novel projection-based
approach. The introduction and analysis of the inverse-probit model,
addressing the significant new statistical challenges posed by it, ranging
across identifiability issues, assessment of intra and inter panel
similarities, extension to complex longitudinal mixed effects settings to
accommodate the motivating application, computational implementation of these
new models, etc. are the novel contributions of this article.

##### Broader implications.

Although this article is focused on understanding auditory category learning,
the insights obtained are generalizable to other domains of behavioral
neuroscience research. In addition, the issues of sensitivity of the results
to the choice of a specific reference category and the problems of clustering
within and between panels that we address here are not specific to the drift-
diffusion framework but are broadly relevant for other classes of generalized
linear models as well, including the extensively studied logit and probit
models. To our knowledge, however, in the vast existing literature on such
models, these problems have not been addressed very rigorously before. The
proposed methodology may thus be of interest to a broader audience interested
in categorical probability models in static and longitudinal settings.

##### Outline of the article.

The rest of this article is organized as follows. Section 2 and Section 3
explain our novel longitudinal inverse-probit mixed model. Section 4 outlines
our computational strategies. Section 5 presents the results of simulation
experiments. Section 6 presents the results of the proposed method applied to
our motivating real data set. Section 7 concludes the main article with a
discussion. Additional details, including Markov chain Monte Carlo (MCMC)
based posterior inference algorithms, are deferred to the supplementary
material.

##  2 Inverse-Probit Model

The inverse-probit categorical probability model follows straightforwardly by
integrating out the (unobserved) response times from the joint model for
response categories and associated response times developed in Paulon et al.
(2020). The derivation is still presented here to keep this article self-
contained.

The basic Wiener diffusion process can be specified as $W(\tau)=\mu\tau+\sigma
B(\tau)$, where $B(\tau)$ is the standard Brownian motion, $\mu$ is the drift
rate, and $\sigma$ is the diffusion coefficient (Cox and Miller, 1965; Ross et
al., 1996). The process has independent normally distributed increments, that
is, $\Delta
W(\tau)=\\{W(\tau+\Delta\tau)-W(\tau)\\}\sim\hbox{Normal}(\mu\Delta\tau,\sigma^{2}\Delta\tau)$,
independently from $W(\tau)$. The first passage time of crossing a threshold
$b$, $\tau=\inf\\{\tau^{\prime}:W(0)=0,W(\tau^{\prime})\geq b\\}$, is then
distributed according to an inverse Gaussian distribution (Whitmore and
Seshadri, 1987; Chhikara, 1988; Lu, 1995) with mean $b/\mu$ and variance
$b\sigma^{2}/\mu^{3}$.

Given perceptual stimuli and a set of decision choices, the neurons in the
brain accumulate evidence in favor of the different alternatives. Modeling
this behavior using latent Wiener processes with unit variances, assuming that
a decision is made when the decision threshold for one of the options is
crossed, as illustrated in Figure 1, Panel (a), a probability model for the
time $\tau_{d}$ to reach the threshold for the $d^{th}$ decision category
under the influence of the $s^{th}$ stimulus is obtained as

$\displaystyle\textstyle
f(\tau_{d}\mid\delta_{s},\mu_{d,s},1,b_{d,s})=\displaystyle\frac{b_{d,s}}{\sqrt{2\pi}}(\tau_{d}-\delta_{s})^{-3/2}\exp\left[-\frac{\\{b_{d,s}-\mu_{d,s}(\tau_{d}-\delta_{s})\\}^{2}}{2(\tau_{d}-\delta_{s})}\right],$
(1)

where $\mu_{d,s}$ denotes the rate of accumulation of evidence, $b_{d,s}$ the
decision boundaries, and $\delta_{s}$ an offset representing the collective
time required to encode the $s^{th}$ signal before evidence accumulation
begins, the time to press a computer key to record a response after a decision
is reached, etc. (Figure 1). We let
$\mbox{\boldmath$\theta$}_{d,s}=(\delta_{s},\mu_{d,s},b_{d,s})^{\rm T}$. Since
a decision $d$ is reached at response time $\tau$ if the corresponding
threshold is crossed first, that is when
$\\{\tau=\tau_{d}\\}\cap_{d^{\prime}\neq d}\\{\tau_{d^{\prime}}>\tau_{d}\\}$,
we have $d=\arg\min\tau_{d^{\prime}}$. Assuming simultaneous accumulation of
evidence for all decision categories, modeled by independent Wiener processes,
and termination when the threshold for the observed decision category $d$ is
reached, the joint distribution of $(d,\tau)$ is thus given by

$\displaystyle f(d,\tau\mid
s,\mbox{\boldmath$\theta$})=g(\tau\mid\mbox{\boldmath$\theta$}_{d,s})\prod_{d^{\prime}\neq
d}\\{1-G(\tau\mid\mbox{\boldmath$\theta$}_{d^{\prime},s})\\}.$ (2)

where, to distinguish from the generic notation $f$, we now use
$g(\cdot\mid\mbox{\boldmath$\theta$})$ and
$G(\cdot\mid\mbox{\boldmath$\theta$})$ to denote, respectively, the
probability density function (pdf) and the cumulative distribution function
(cdf) of an inverse Gaussian distribution, as defined in (1).

When the response times $\tau$ are unobserved, the probability of taking
decision $d$ under the influence of stimulus $s$ is thus obtained as

$\displaystyle P(d\mid
s,\mbox{\boldmath$\theta$})=\int_{\delta_{s}}^{\infty}g(\tau\mid\mbox{\boldmath$\theta$}_{d,s})\prod_{d^{\prime}\neq
d}\left\\{1-G(\tau\mid\mbox{\boldmath$\theta$}_{d^{\prime},s})\right\\}d\tau.$
(3)

Model (3) is conceptually similar to traditional multinomial probit/logit
regression models (Borooah, 2002; Agresti, 2018) except that the latent
variables are now inverse Gaussian distributed as opposed to being normal or
extreme-value distributed, and the observed category is associated with the
minimum of the latent variables in contrast to being identified with the
maximum of the latent variables. We thus refer to this model as a ‘multinomial
inverse-probit model’.

With data on both response categories and response times available, the joint
model (2) was used to construct the likelihood function in Paulon et al.
(2020). In the absence of data on the response times, however, the inverse-
probit model in (3) provides the basic building block for constructing the
likelihood function for the observed response categories. The inverse-probit
model, however, brings in many unique identifiable issues, significantly
complicating the inferential challenges, which motivated the research reported
in this article, as detailed in the following section.

###  2.1 Identifiability Issues and Related Modeling Challenges

To begin with, we note that model (3) in itself cannot be identified from data
on only the response categories. The offset parameters can easily be seen to
not be identifiable since
$P\left(\tau_{d}\leq\wedge_{d^{\prime}}\tau_{d^{\prime}}\right)=P\left\\{(\tau_{d}-\delta)\leq\wedge_{d^{\prime}}\left(\tau_{d^{\prime}}-\delta\right)\right\\}$
for any $\delta$, where $\wedge_{d}\tau_{d}$ denotes the minimum of
$\tau_{d}$, $d=1,\ldots,d_{0}$. As is also well known in the literature, in
categorical probability models, the location and scale of the latent
continuous variables are not also separately identifiable. The following lemma
establishes these points for the inverse-probit model.

######  Lemma 1.

The offset parameters $\delta_{s}$ are not identifiable in model (3). The
drift and the boundary parameters, respectively $\mu_{d,s}$ and $b_{d,s}$, are
also not separately identifiable in model (3).

In the proof of Lemma 1 given in Appendix A, we have specifically shown that
$P\left(d\mid s,\mbox{\boldmath$\theta$}\right)=P\left(d\mid
s,\mbox{\boldmath$\theta$}^{\star}\right),$ where the drift and boundary
parameters in $\mbox{\boldmath$\theta$}=\left(\mu_{d,s},b_{d,s}\right)$ and
$\mbox{\boldmath$\theta$}^{\star}=\left(\mu_{d,s}^{\star},b_{d,s}^{\star}\right)$
satisfy $\mu_{d,s}^{\star}=c\mu_{d,s}$ and $b_{d,s}^{\star}=c^{-1}b_{d,s}$ for
some constant $c>0$. The result follows from noting that the transformation
$\tau_{d,s}^{\star}={c}^{-2}\tau_{d,s}$ does not change the ordering between
the $\tau_{d,s}$’s and hence the probabilities of the resulting decisions
$d=\arg\min\tau_{d^{\prime}}=\arg\min\tau_{d^{\prime}}^{\star}$ also remain
the same. This has the simple implication that if the rate of accumulation of
evidence is faster, then the same decision distribution is obtained if the
corresponding boundaries are accordingly closer and conversely.

In fact, given the information on input and output categories alone, if
$d_{0}$ denotes the number of possible decision categories, at most $d_{0}-1$
parameters are estimable. To see this, consider the probabilities $P(d\mid
s,\mbox{\boldmath$\theta$})$, $d=1,\ldots,d_{0}$, where $\theta$ is the
$m$-dimensional vector of parameters, possibly containing drift parameters and
decision boundaries. Given the perceptual stimuli $s$ as input, the
probabilities satisfy $\sum_{d=1}^{d_{0}}P(d\mid
s,\mbox{\boldmath$\theta$})=1$. Thus, the function
$P_{s}(\mbox{\boldmath$\theta$})=\\{P(1\mid
s,\mbox{\boldmath$\theta$}),\cdots,P(d_{0}\mid
s,\mbox{\boldmath$\theta$})\\}^{\rm T}$ lie on a $d_{0}-1$-dimensional
simplex,
$P_{s}(\mbox{\boldmath$\theta$}):\mbox{\boldmath$\theta$}\rightarrow\Delta^{d_{0}-1}$,
and by the model in (3) the mapping is continuous. Thus, it can be shown by
the Invariance of Domain theorem (see, for e.g., Deo, 2018) that if $P_{s}$ is
injective and continuous, then the domain of $P_{s}$ must belong to
$\mathbb{R}^{m}$, where $m\leq d_{0}-1$. Thus in order to ensure
identifiablilty of
$\left\\{P_{s}(\mbox{\boldmath$\theta$});\mbox{\boldmath$\theta$}\right\\}$,
we must parametrize the probability vector with at most ${d_{0}-1}$
parameters.

The existing literature on drift-diffusion models discussed in the
Introduction have traditionally allowed more flexibility in modeling the
drifts. Previous research on joint models for response tones and associated
response times in Paulon et al. (2020) also suggest that the boundaries remain
stable around a value of $2$ and it is primarily the changes in the drift rate
that explain longitudinal learning. In view of this, we keep the boundaries
fixed at the constant $b=2$ and treat the drifts to be the free parameters
instead. In our simulations and real data applications, it is observed that
the estimates of the drift parameters and the associated cluster
configurations are not very sensitive to small-to-moderate deviations of $b$
around $2$. In our codes implementing our method, available as part of the
supplementary materials, we allow the practitioner to choose a value of $b$ as
they see fit for their specific application. A latent drift-diffusion process
with these constraints, namely $\delta_{s}=0$ and $b_{d,s}=b$, is shown in
Figure 1, Panel (b).

While fixing $\delta_{s}=0$ and $b_{d,s}$ to some known constant $b$ reduces
the size of the parameter space to $d_{0}$, to ensure identifiability, we
still need at least one more constraint on the drift parameters
$\mbox{\boldmath$\mu$}_{1:d_{0},s}$. In categorical probability models the
identifiabily problem of the location parameter is usually addressed by
setting one category as a reference and modeling the probabilities of the
others (Albert and Chib, 1993; Borooah, 2002; Agresti, 2018). However,
posterior predictions from Bayesian categorical probability models with
asymmetric constraints may be sensitive to the choice of reference category
(see Burgette and Nordheim, 2012; Johndrow et al., 2013). Further, as also
discussed in the Introduction, the goal of clustering the drift parameters
$\mbox{\boldmath$\mu$}_{1:d_{0},s}$ across $s$ can not be accomplished by this
apparently simple solution.

The problem can be addressed by imposing a symmetric constraint on
$\mbox{\boldmath$\mu$}_{1:d_{0},s}$ instead. Symmetric identifiability
constraints have been previously proposed by Burgette and Hahn (2010). Their
sum to zero constraint solves the problem of relying on a single base
category, instead it is adaptively chosen by the model. However, the
introduction of a base category, even if adaptively chosen, does not
facilitate the clustering of $\mbox{\boldmath$\mu$}_{1:d_{0},s}$ within and
across the different input categories $s$.

###  2.2 Our Proposed Approach

In coming up with solutions for these challenges, we take into consideration
the complex design of our motivating tone learning experiments, so that our
approach is easily extendable to longitudinal mixed model settings, allowing
us to (a) estimate the smoothly varying trajectories of the parameters as the
participants learn over time, (b) accommodate the heterogeneity between the
participants, and (c) compare between the estimates not just within but also
crucially between the different panels.

Similar to the sum to zero constraint in the multinomial probit model of
Burgette and Nordheim (2012), we impose a symmetric sum to a constant
constraint on the drift parameters $\mbox{\boldmath$\mu$}_{1:d_{0},s}$ to
identify our new class of inverse-probit models, although our implementation
is quite different from theirs. We conduct inference by starting with an
unconstrained prior, then sampling from the corresponding unconstrained
posterior, finally projecting these samples to the constrained space through a
minimal distance mapping.

This approach is significantly advantageous both from a modeling and a
computational perspective. On one hand, the basic building blocks are
relatively easily extended to complex longitudinal mixed model settings, on
the other hand, posterior computation is facilitated as this allows the use of
conjugate priors for the unconstrained parameters. Projection of the drift
parameters onto the same space further makes them directly comparable,
allowing clustering within and across the panels.

####  2.2.1 Minimal Distance Mapping

As the drift parameters are positive, the sum to a constant $k$ constraint
leads to constrained space ${\cal S}_{k}=\\{\mbox{\boldmath$\mu$}:{\bf
1}^{T}\mbox{\boldmath$\mu$}=k,~{}\mu_{j}>0,~{}j=1,\ldots,d_{0}\\}$, in which
$\mbox{\boldmath$\mu$}_{1:d_{0},s}$ should be projected. The space ${\cal
S}_{k}$ is semi-closed, and therefore, the projection of any point $\mu$ onto
${\cal S}_{k}$ may not exist. As a simple one dimensional example, let $x=-1$
and $\mathcal{S}=(0,1]$, then
$\arg\min_{y\in\mathcal{S}}|y-x|=0\notin\mathcal{S}$. Further, from a
practical perspective, a drift parameter infinitesimally close to zero makes
the corresponding distribution of response time nearly non-informative.
Therefore, we choose a small $\varepsilon>0$ and project $\mu$ onto ${\cal
S}_{\varepsilon,k}=\\{\mbox{\boldmath$\mu$}:{\bf
1}^{T}\mbox{\boldmath$\mu$}=k,~{}\mu_{j}\geq\varepsilon,~{}j=1,\ldots,d_{0}\\}.$
We then define the projection of a point $\mu$ onto ${\cal S}_{\varepsilon,k}$
through minimal distance mapping as

$\displaystyle\mbox{\boldmath$\mu$}^{\star}=\mathrm{Proj}_{{\cal
S}_{\varepsilon,k}}(\mbox{\boldmath$\mu$}):=\left\\{\mathrm{argmin}_{\mbox{\boldmath$\nu$}}\|\mbox{\boldmath$\mu$}-\mbox{\boldmath$\nu$}\|:\mbox{\boldmath$\nu$}\in{\cal
S}_{\varepsilon,k}\right\\}.$

Note that for appropriate choices of $(k,\varepsilon)$, ${\cal
S}_{\varepsilon,k}$ is non-empty, closed and convex. Therefore,
$\mbox{\boldmath$\mu$}^{\star}$ exists and is unique by the Hilbert projection
theorem (Rudin, 1991). The solution to this projection problem is given in the
following result.

######  Lemma 2.

Let ${\cal S}_{\varepsilon,k}$ be as defined above, and ${\cal
S}_{\varepsilon}=\left\\{\mbox{\boldmath$\mu$}:\mu_{j}\geq\varepsilon,~{}j=1,\ldots,d_{0}\right\\}$.
Then, $\mathrm{Proj}_{{\cal
S}_{\varepsilon,k}}(\mbox{\boldmath$\mu$})=\mathrm{Proj}_{{\cal
S}_{\varepsilon}}(\mbox{\boldmath$\mu$}-u^{\star}{\bf 1})$, where $u^{\star}$
is a solution to the equation ${\bf 1}^{T}\mathrm{Proj}_{{\cal
S}_{\varepsilon}}(\mbox{\boldmath$\mu$}-u^{\star}{\bf 1})=k$.

A proof of Lemma 2 is given in Beck (2017). Although the analytical form of
the solution is not available, as is evident from the above result, the
solution mainly relies on finding a root $u^{\star}$ of the non-increasing
function $\phi(u^{\star})={\bf 1}^{T}\mathrm{Proj}_{{\cal
S}_{\varepsilon}}(\mbox{\boldmath$\mu$}-u^{\star}{\bf 1})-k$. We apply an
algorithm by Duchi et al. (2008) to reach the solution. The algorithm is
described in Appendix C.

####  2.2.2 Identifiability Restrictions

The projection approach solves the problem of identifiability and maps the
probability vector corresponding to an input tone $s$ to the constraint space
of $\mbox{\boldmath$\mu$}_{1:d_{0},s}$, ${\cal S}_{\varepsilon,k}$. The
following theorem shows that the mapping from the constrained space of
$\mbox{\boldmath$\mu$}_{1:d_{0},s}$ to the probability vector
${\mathbf{P}}(\mbox{\boldmath$\mu$}_{1:d_{0},s})=\\{p_{1}(\mbox{\boldmath$\mu$}_{1:d_{0},s}),\dots,p_{d_{0}}(\mbox{\boldmath$\mu$}_{1:d_{0},s})\\}^{\rm
T}$ is injective. To keep the ideas simple, we consider the domain of the
function to be ${\cal S}_{0,k}$ (i.e., $\varepsilon=0$) instead of ${\cal
S}_{\varepsilon,k}$ although a very similar proof would follow if ${\cal
S}_{\varepsilon,k}$ were considered.

######  Theorem 1.

Let $p_{d}(\mbox{\boldmath$\mu$}_{1:d_{0},s})$ be the probability of observing
the output tone $d$ given the input tone $s$ and the drift parameters
$\mbox{\boldmath$\mu$}_{1:d_{0},s}$, as given in (3), for each $d=1:d_{0}$.
Suppose $\mbox{\boldmath$\mu$}_{1:d_{0},s}$ lies on the space ${\cal
S}_{0,k}$. Then, the function from ${\cal S}_{0,k}$ to the space of
probabilities
$\left\\{p_{d}(\mbox{\boldmath$\mu$}_{1:d_{0},s});d=1:d_{0}\right\\}$ is
injective.

A proof is presented in Appendix B.

####  2.2.3 Conjugate Priors for the Unconstrained Drifts

From (3), given $\tau_{1},\ldots,\tau_{d_{0}}$, such that
$\tau_{d}\leq\min\left\\{\tau_{1},\ldots,\tau_{d_{0}}\right\\}$, the posterior
full conditional of $\mbox{\boldmath$\mu$}_{1:d_{0},s}$ is proportional to
$\pi_{\mbox{\boldmath$\mu$}}^{(n)}\propto\pi(\mbox{\boldmath$\mu$}_{1:d_{0},s})\times\prod_{d=1}^{d_{0}}g\left(\tau_{d}~{}\mid\mu_{d,s}\right),$
where $\pi(\cdot)$ is the prior of $\mbox{\boldmath$\mu$}_{1:d_{0},s}$.
Observe that $\prod_{d=1}^{d_{0}}g\left(\tau_{d}~{}\mid\mu_{d,s}\right)$ is
Gaussian in $\mbox{\boldmath$\mu$}_{1:d_{0},s}$. A Gaussian prior on
$\mbox{\boldmath$\mu$}_{1:d_{0},s}$ thus induces a conditional posterior for
$\mbox{\boldmath$\mu$}_{1:d_{0},s}$ that is also Gaussian and hence very easy
to sample from. Importantly, these benefits also extend naturally to
multivariate Gaussian priors for any parameter vector
$\mbox{\boldmath$\beta$}_{d,s}$ that relates to $\mu_{d,s}$ linearly. This
will be crucial in allowing us to extend the basic building block to
longitudinal functional mixed model settings in Section 3 next, where we will
be modeling time varying $\mu_{d,s}(t)$ as flexible mixtures of B-splines with
associated coefficients $\mbox{\boldmath$\beta$}_{d,s}$.

####  2.2.4 Justification as a Proper Bayesian Procedure

Let us define the constrained conditional posterior distribution,
$\tilde{\pi}_{\tilde{\mbox{\boldmath$\mu$}}}^{(n)}$, as

$\displaystyle\tilde{\pi}_{\tilde{\mbox{\boldmath$\mu$}}}^{(n)}\left(B\mid\mbox{\boldmath$\zeta$}\right)=\pi_{\mbox{\boldmath$\mu$}}^{(n)}\left(\left\\{\mbox{\boldmath$\mu$}:\mathrm{Proj}(\mbox{\boldmath$\mu$})\in
B\right\\}\mid\mbox{\boldmath$\zeta$}\right),\quad B\subseteq{\cal
S}_{\varepsilon,k},$

where $\pi_{\mbox{\boldmath$\mu$}}^{(n)}$ is the unconstrained conditional
posterior of $\mbox{\boldmath$\mu$}_{1:d_{0},s}$, given the other variables
$\zeta$. The analytic form of the constrained conditional posterior is not
available.

Sen et al. (2021) established a proper Bayesian justification for the
posterior projection approach by showing the existence of a prior
$\tilde{\pi}\left(\mbox{\boldmath$\mu$}_{1:d_{0},s}\right)$ on the constrained
space ${\cal S}_{\varepsilon,k}$ such that the resulting posterior is same as
the projected posterior $\tilde{\pi}_{\tilde{\mbox{\boldmath$\mu$}}}^{(n)}$.
When ${\cal S}_{\varepsilon,k}$ is non-empty, closed and convex, i.e., the
projection operator is measurable, such a prior exists if the unconstrained
posterior is absolutely continuous with respect to the unconstrained prior
(Sen et al., 2021, Corollary 1). As the unconstrained induced prior and
posterior of the drift parameters are both Gaussian, this result holds in our
case as well.

##  3 Extension to Longitudinal Mixed Models

In this section we adapt the inverse probit model discussed in Section 2 to
complex longitudinal design of the PTC1 or similar data sets described in
Introduction. Let $s_{i,\ell,t}$ denote the input tone for the $i^{th}$
individual in the $\ell^{th}$ trial of block $t$. Likewise, let $d_{i,\ell,t}$
denote, respectively, the output tone selected by the $i^{th}$ individual in
the $\ell^{th}$ trial of block $t$. Setting the offsets at zero, and boundary
parameters to a fixed constant $b$, we now have

$\displaystyle\begin{split}P\\{d_{i,\ell,t}\mid
s_{i,\ell,t}=s,\mbox{\boldmath$\mu$}_{1:d_{0},s}^{(i)}(t)\\}=\int_{0}^{\infty}g\\{\tau\mid\mu_{d,s}^{(i)}(t)\\}\prod_{d^{\prime}\neq
d}\left[1-G\\{\tau\mid\mu_{d^{\prime},s}^{(i)}(t)\\}\right]d\tau,\\\
\mbox{where}~{}~{}g\\{\tau\mid
s_{i,\ell,t}=s,\mu_{d,s}^{(i)}(t)=\mu\\}\textstyle=\displaystyle\frac{b}{\sqrt{2\pi}\tau^{3/2}}\exp\left[-\frac{\\{b-\mu\tau\\}^{2}}{2\tau}\right].\hskip
72.26999pt\end{split}$ (4)

The drift rates $\mu_{d,s}^{(i)}(t)$ now vary with the blocks $t$. In
addition, we accommodate random effects by allowing $\mu_{d,s}^{(i)}(t)$ to
also depend on the subject index $i$. We let
${\mathbf{d}}=\\{d_{i,\ell,t}\\}_{i,\ell,t}$, and $d_{0}$ be the number of
possible decision categories (T1, T2, $\cdots$, T${d_{0}}$). The likelihood
function of our longitudinal inverse-probit mixed model thus takes the form

$\displaystyle
L({\mathbf{d}}\mid{\mathbf{s}},\mbox{\boldmath$\theta$})=\prod_{d=1}^{d_{0}}\prod_{s=1}^{d_{0}}\prod_{t=1}^{T}\prod_{i=1}^{n}\prod_{\ell=1}^{L}\left[P\\{d_{i,\ell,t}\mid
s_{i,\ell,t},\mbox{\boldmath$\mu$}_{1:d_{0},s}^{(i)}(t)\\}\right]^{1\\{d_{i,\ell,t}=d,s_{i,\ell,t}=s\\}}.$

Toward modeling mixed drifts, with $x=(d,s)\in{\cal
X}=\\{(1,1),(1,2),\dots,(d_{0},d_{0})\\}\equiv\\{1,2,\dots,x_{\max}\\}$
representing the input-response category combinations, we let

$\displaystyle\begin{split}&\mu_{x}^{(i)}(t)=f_{x}(t)+u_{x}^{(i)}(t),~{}~{}~{}~{}~{}u_{x}^{(i)}(t)\sim
f_{u}\\{u_{x}^{(i)}(t)\\}.\end{split}$ (5)

Here $f_{x}(t)$ and $u_{x}^{(i)}(t)$ denote, respectively, fixed and random
effects components in the exponential scale; $f_{u}$ denotes the random
effects distribution. To simplify notation, generic data recording time stamps
in $\\{1,\dots,T\\}$ as well as other generic time points in $[1,T]$ will both
be denoted by $t$. Likewise, generic input-response category combinations as
well as their particular values will both be denoted by $x$ and so forth.

In the following subsection, we discuss our model for the drift parameters
without imposing any constraints. The imposition of the constraint via
projection as well as the clustering of the parameters to allow direct
comparisons are discussed thereafter.

###  3.1 Functional Fixed Effects

We model the fixed effects functions $f_{x}(t)$ using flexible mixtures of
B-spline bases (de Boor, 1978) that allow them to smoothly vary with time $t$
while also depending locally on the indexing variable $x$ as

$\displaystyle\textstyle
f_{x}(t)=\sum_{k=1}^{K}\beta_{x,k}B_{k}(t)={\mathbf{B}}(t)\mbox{\boldmath$\beta$}_{x}.$
(6)

Here ${\mathbf{B}}(t)=\\{B_{1}(t),\dots,B_{K}(t)\\}$ are a set of known
locally supported basis functions spanning $[1,T]$,
$\mbox{\boldmath$\beta$}_{x}=(\beta_{x,1},\dots,\beta_{x,K})^{\rm T}$ are
associated unknown coefficients to be estimated from the data. Allowing the
$\mbox{\boldmath$\beta$}_{x}$’s to flexibly vary with $x$ can generate widely
different shapes for different input-response category combinations.

In the absence of data on response times, we have limited statistical ability
to model subtle variations in the drift curves over time. We therefore
consider a model for clustering of the fixed effects curves globally on the
entire time domain, unlike Paulon et al. (2020), who focused on clustering
locally at each time point instead.

To that end, introducing a set of latent variables $z_{x}$ for each input-
response category combination $x$ with a shared state space
$\\{1,\dots,z_{\max}\\}$ and associated coefficient atoms
$\mbox{\boldmath$\beta$}_{z}^{\star}=(\beta_{z,1}^{\star},\dots,\beta_{z,K}^{\star})^{\rm
T}$, we let

$\displaystyle\begin{split}&\textstyle(\mbox{\boldmath$\beta$}_{x}\mid
z_{x}=z)=\mbox{\boldmath$\beta$}_{z}^{\star},~{}~{}\text{implying}~{}~{}\\{f_{x}(t)\mid
z_{x}=z\\}=f_{z}^{\star}(t)=\sum_{k=1}^{K}\beta_{z,k}^{\star}B_{k}(t),\\\
\end{split}$ (7)

To probabilistically cluster the $\mbox{\boldmath$\beta$}_{x}$’s, we next let

$\displaystyle\begin{split}&z_{x}\sim\hbox{Mult}(\mbox{\boldmath$\pi$}_{z})=\hbox{Mult}(\pi_{1},\dots,\pi_{z_{\max}}),\\\
&\mbox{\boldmath$\pi$}_{z}\sim\hbox{Dir}(\alpha/z_{\max},\dots,\alpha/z_{\max}).\end{split}$
(8)

We next consider priors for the atoms $\mbox{\boldmath$\beta$}_{z}^{\star}$.
We let

$\displaystyle\begin{split}&\mbox{\boldmath$\beta$}_{z}^{\star}\sim\hbox{MVN}_{K}\\{\mbox{\boldmath$\mu$}_{\beta,0},(\sigma_{a}^{-2}{\mathbf{I}}_{K}+\sigma_{s}^{-2}{\mathbf{P}})^{-1}\\},~{}~{}~{}~{}~{}\\\
\end{split}$ (9)

where $\hbox{MVN}_{K}(\mbox{\boldmath$\mu$},\mbox{\boldmath$\Sigma$})$ denotes
a $K$ dimensional multivariate normal distribution with mean $\mu$ and
covariance $\Sigma$ and ${\mathbf{P}}={\mathbf{D}}^{\rm T}{\mathbf{D}}$, where
the $(K-1)\times K$ matrix ${\mathbf{D}}$ is such that
${\mathbf{D}}\mbox{\boldmath$\beta$}$ computes the first order differences in
$\beta$. The model thus penalizes
$\sum_{k=1}^{K}(\nabla\beta_{z,k}^{\star})^{2}=\mbox{\boldmath$\beta$}^{\rm
T}{\mathbf{P}}\mbox{\boldmath$\beta$}$, the sum of squares of first order
differences in $\mbox{\boldmath$\beta$}_{u}^{(i)}$ (Eilers and Marx, 1996).
The variance parameter $\sigma_{s}^{2}$ models the smoothness of the
functional atoms, smaller $\sigma_{s}^{2}$ inducing smoother
$f_{z}^{\star}(t)$’s. Additional departures from
$\mbox{\boldmath$\mu$}_{\beta,0}$ are explained by the other variance
component $\sigma_{a}^{2}$. We assign half Cauchy priors on the variance
parameters as

$\displaystyle\sigma_{s}^{2}\sim\hbox{C}^{+}(0,1),~{}~{}~{}~{}~{}\sigma_{a}^{2}\sim\hbox{C}^{+}(0,1).$

###  3.2 Functional Random Effects

We allow different random effects $u_{C}^{(i)}(t)$ and $u_{I}^{(i)}(t)$ for
correct (C) (when $d=s$) and incorrect (I) (when $d\neq s$) identifications,
respectively, as

$\displaystyle
u_{d,s}^{(i)}(t)=u_{C}^{(i)}(t)~{}~{}~{}~{}\text{when}~{}d=s,~{}~{}~{}~{}~{}u_{d,s}^{(i)}(t)=u_{I}^{(i)}(t)~{}~{}~{}~{}\text{when}~{}d\neq
s.$

Suppressing the subscripts to simplify notation, we model the time-varying
random effects components $u^{(i)}(t)$ as

$\displaystyle\begin{split}&\textstyle
u^{(i)}(t)=\sum_{k=1}^{K}\beta_{u,k}^{(i)}B_{k}(t)={\mathbf{B}}(t)\mbox{\boldmath$\beta$}_{u}^{(i)},~{}~{}~{}~{}~{}\\\
&\mbox{\boldmath$\beta$}_{u}^{(i)}\sim\hbox{MVN}_{K}\\{{\mathbf{0}},(\sigma_{u,a}^{-2}{\mathbf{I}}_{K}+\sigma_{u,s}^{-2}{\mathbf{P}})^{-1}\\},~{}~{}~{}~{}~{}\\\
\end{split}$ (10)

where
$\mbox{\boldmath$\beta$}_{u}^{(i)}=(\beta_{1,u}^{(i)},\dots,\beta_{K,u}^{(i)})^{\rm
T}$ are subject-specific spline coefficients. We assign non-informative half
Cauchy priors on the variance parameters as

$\displaystyle\sigma_{u,s}^{2}\sim\hbox{C}^{+}(0,1),~{}~{}~{}~{}~{}\sigma_{u,a}^{2}\sim\hbox{C}^{+}(0,1).$

Given the updated coefficients $\mbox{\boldmath$\beta$}_{x}$ and
$\mbox{\boldmath$\beta$}_{u}^{(i)}$ from the full conditional posterior
distribution, an unconstrained update of $\mu_{x}^{(i)}(t)$ is obtained. Given
an input tone $s$, the updated $\mbox{\boldmath$\mu$}_{1:d_{0},s}^{(i)}(t)$ is
then projected to the space ${\cal S}_{\varepsilon,k}$ as described in Section
2.2.1.

##  4 Posterior Inference

Posterior inference for our proposed inverse-probit mixed model is carried out
using samples drawn from the posterior using MCMC algorithm. The algorithm
carefully exploits the conditional independence relationships encoded in the
model as well as the latent variable construction of the model.

Inference can be greatly simplified by sampling the the passage times
$\tau_{1:d_{0}}$ and then conditioning on them. However, it is not possible to
generate $\tau_{1:d_{0}}$ sequentially, for example, by generating the passage
time of the $d$-th decision choice $\tau_{d}$ independently, and that of the
other decision choices from a truncated inverse-Gaussian distribution, left
truncated at $\tau_{d}$.

We instead implement a relatively simple accept-reject sampler which generates
values from the joint distribution of $\tau_{1:d_{0}}$ and accepts the sample
if $\tau_{d}\leq\tau_{1:d_{0}}$. It is fast and produces a sample from the
desired target conditional distribution. We formalize this result in the
following lemma.

######  Lemma 3.

Let $g\left(\tau_{1:d_{0}}\mid\mu_{1:d_{0}}\right)$ be the joint distribution
of $\tau_{1:d_{0}}$. Consider the following accept-reject algorithm:

Algorithm 1 Generating the passage times $\tau_{1:d_{0}}$ given
$\tau_{d}\leq\tau_{1:d_{0}}$

1:Generate $\tau_{1:d_{0}}$ from the joint distribution
$g\left(\tau_{1:d_{0}}|\mu_{1:d_{0}}\right)$.

2:Accept $\tau_{1:d_{0}}$ if $\tau_{d}\leq\tau_{1:d_{0}}$.

3:Return to Step 1 otherwise.

Algorithm 1 generates samples from the conditional joint distribution of
$\tau_{1:d_{0}}$, conditioned on the event $\tau_{d}\leq\tau_{1:d_{0}}$.

Proof of Lemma 3 is provided in Appendix D.

It can be verified that the acceptance ratio of Algorithm 1 is
$M^{-1}=P\left(\tau_{d}\leq\tau_{1:d_{0}}\right)$ (see Robert and Casella,
2004) which depends on the drift parameters only. If the drift parameters are
ordered accordingly, so as to satisfy $\mu_{d}\geq\mu_{1:d_{0}}$, the
acceptance ratios become close to one. The implication is that the algorithm
becomes faster as the sampler converges.

As noted earlier, sampling the latent inverse-gaussian distributed response
times $\tau_{1:d_{0}}$ greatly simplifies computation. Most of the chosen
priors, including the priors on the coefficients $\beta$ in the fixed and
random effects, are conjugate. Due to space constraints, the details are
deferred to Section S.2 in the supplementary material.

##  5 Simulation Studies

In this section, we discuss the results of a synthetic numerical experiment.
We are not aware of any method for longitudinal categorical learning data
which can perform intra and inter panel clustering simultaneously, while
making inference on the drift parameters in the absence of the response times.
Therefore, we restrict our focus mostly to evaluating the performances of the
proposed model by its proximity to the true data generating process but are
unable to present any comparisons with any other model.

##### Design.

In designing the simulation scenarios, we have tried to mimic our motivating
category learning data sets. We chose $n=20$ as the number of participants
being trained over $T=10$ blocks to identify $d_{0}=4$ tones. For each input
tone, and each block, there are $L=40$ trials. We set the true $\mu_{d,s}(t)$
values in such a way that they are far from satisfying the constraint
$\sum_{1:d_{0}}\mu_{d,s}=k$, and the decision boundary is set to $b=2$ for all
$(d,s)$. The true drift parameters and the true probabilities, averaged over
the participants of each input-response category combinations are shown in
Figure LABEL:fig:_simulated_data_description.

There are four true clusters in total, two in correct categorizations,
$S_{1},~{}S_{2}$, and two in incorrect categorizations, $M_{1},~{}M_{2}$, as
follows: $S_{1}=\\{(1,1),(2,2)\\}$, $S_{2}=\\{(3,3),(4,4)\\}$,
$M_{1}=\\{(1,2),(1,3),(2,1),(2,3),(3,4),(4,3)\\}$,
$M_{2}=\\{(1,4),(2,4),(3,1),(3,2),(4,1),(4,2)\\}$. We may interpret $M_{1}$ as
the cluster of difficult alternatives, and $M_{2}$ as the cluster of easy
alternatives. Thus, there are similarities in overall trajectories of
$\\{T_{1},T_{2}\\}$ and $\\{T_{3},T_{4}\\}$, differentiating between easy and
hard category recognition problems. We experimented with $50$ synthetic data
sets generated according to this design.

##### Results.

As the true drift parameters themselves do not satisfy the constraint, and the
estimated drift parameters are on the constrained space, we cannot validate
our method by its predictive performance of the drift parameters. Instead the
proposed method is validated in terms of the estimated probabilities.

Figure LABEL:fig:_simulation_estimated_prob shows the estimated posterior
probability trajectories along with the $95\%$ credible interval and the
underlying true probability curves for every combination $(d,s)$ in a typical
scenario. The credible interval fails to capture the truth in two situations,
when the true probability is very close to zero, or it is very close to one.
The former case corresponds to classes with very low success probability,
resulting very few observations to estimate. The latter is underestimated as a
consequence of the former, since the probabilities add up to one.

The results produced by our method are mostly stable and consistent across all
synthetic data sets. There are, however, a few cases of incorrect cluster
assignments, resulting in some outliers in each boxplot. Note that if an
incorrect cluster assignment takes places, probabilities of all input-response
combinations are affected by that. For example, if a component of $M_{1}$ is
wrongly assigned to $M_{2}$, then not only the probabilities of input-output
combinations in $M_{1}$ and $M_{2}$ are affected, since the probabilities sum
up to one, those of $S_{1}$ and $S_{2}$ are also affected.

In estimating the probabilities, the overall mean squared error, i.e., the
mean squared difference of the estimated and the true probabilities taking all
combinations of $(d,s,i,t)$ into account, came out to be $0.0028$. Figure 3
provides the detailed description of estimation of the probabilities for two
input categories (one from each similarity group). As described for the
individual simulation results, there are cases of under-estimation of the
probabilities which are close to one, and consequently, over-estimation of the
probabilities close to zero. However, the amount of departure from the true
probability in each case is very small which can also be seen in the small
overall MSE.

Figure 3: Results of the synthetic data: Figure showing boxplots of the
estimated probabilities over $50$ simulations, and true probabilities (in red
dot) of each block and for two panels, one from each similarity group (panel
$T_{1}$ in top and $T_{3}$ in bottom).

Further, the overall efficiency in identifying the true clustering structure
is validated using Rand (Rand, 1971) and adjusted Rand (Hubert and Arabie,
1985) indices. The definitions of Rand and adjusted Rand indices are provided
in Section S.5 in the supplementary material. The average Rand and adjusted
Rand indices for our proposed method over $50$ simulations $0.9105$ and
$0.8277$, respectively, indicating high overall efficacy in correctly
clustering the probability curves.

##  6 Applications

##### Analysis of PTC1 Data Set.

We present the analysis of the PTC1 data set described in the Introduction
using our proposed longitudinal inverse-probit mixed model here. We first
demonstrate the performance of the proposed method in estimating the
probabilities associated with different $(d,s)$ pairs. Figure 4 shows the
$95\%$ credible intervals for the estimated probabilities for different input
tones, along with average proportions of times an input tone was classified
into different tone categories across subjects. The latter serves as the
empirical estimate of the probabilities.

Figure 4: Results for PTC1 data: Estimated probability trajectories compared
with average proportions of times an input tone was classified into different
tone categories across subjects (in dashed line). The means across subjects
are indicated by thick lines and the shaded regions indicate corresponding
$95\%$ coverage regions. Estimates related to high-flat tone responses are
shown in red; those for low-rising tone responses in blue; low-dipping tone
responses in green; and high-falling tone responses in purple.

We observe that, except for input-response combination (1,1) for block 3, and
some situations with low number of data points, the $95\%$ credible intervals
include the corresponding empirical probabilities. Further, the estimated
credible regions are often quite narrow, implying high precision of the
estimates.

Next we examine the clusters identified by the proposed model. Apart from the
two clusters obtained for the success combinations $(d=s)$, three are
additionally identified in the incorrect input-response combinations $(d\neq
s)$. The clusters of success combinations are $S_{1}=\\{(1,1),(2,2),(4,4)\\}$
and $S_{2}=\\{(3,3)\\}$, and of wrong allocations are
$M_{1}=\\{(1,2),(1,4),(2,1),(2,4),(3,2),(4,1),(4,2)\\}$,
$M_{2}=\\{(1,3),(2,3),(3,4),(4,3)\\}$, and $M_{3}=\\{(3,1)\\}$. Figure 5 shows
the input-response tone combinations color codes as per cluster identity, and
the proportion of times each pair of input-response tone combinations appeared
in the same cluster after burnin. Figure 5 indicates that, while the clusters
$S_{1},S_{2},M_{1}$ are stable, there are some instability among the other two
clusters, namely $M_{2}$ and $M_{3}$.

Figure 5: Results for PTC1 data: Network plot of similarity groups showing the
intra and inter cluster similarities of tone recognition problems. Each node
is associated with a pair indicating the input-response tone category,
$(s,d)$. The number associated with each edge indicates the proportion of
times the pair in the two connecting nodes appeared in the same cluster after
burnin.

Generally speaking, the trajectories of the low-rising tone $T_{2}$ and the
high-falling tone $T_{4}$ are similar with two wrong allocations from $M_{1}$
and one from $M_{2}$. Although the final clustering results show that the
high-level tone $T_{1}$ has a similar trajectory as $T_{2},T_{4}$, it is clear
(Figure 5) that the input-response tone category $(1,3)$ is less associated
with $M_{2}$ than the other members of $M_{2}$. In fact in $24\%$ cases, it
formed an individual group of its own. The trajectory of low-dipping tone
$T_{3}$, on the other hand, is very different from the others (Figure 5).

Next we consider the results concerning the estimation of the drift parameters
$\mu_{x}^{(i)}(t)$. As discussed throughout, the estimates of
$\mu_{x}^{(i)}(t)$ can only be interpreted on a relative scale. Figure 6 shows
the posterior mean trajectories and associated $95\%$ credible intervals for
the projected drift rates estimated by our method for different combinations
of $x=(d,s)$. The relative similarities and differences of the estimated
drifts are consistent with the knowledge of the domain experts.

Figure 6: Results for PTC1 data: Estimated posterior mean trajectories of the
population level drifts $\mu_{x}(t)$ for the proposed model. The shaded areas
represent the corresponding $95\%$ pointwise credible intervals. Estimated
related to high-flat tone responses are shown in red; those for low-rising
tone responses in blue; low-dipping tone responses in green; and high-falling
tone responses in purple.

Importantly, our proposed mixed model also allows us to assess individual
specific parameter trajectories. Figure 7 shows the posterior mean
trajectories and the associated $95\%$ credible intervals for the drift rates
$\mu_{x}^{(i)}$ estimated by our method for the different success combinations
of $x=(d,s)$ for two participants \- one with the best accuracy averaged
across all blocks, and the other with the worst accuracy averaged across all
blocks. These results suggest significant individual specific heterogeneity.
For the well performing participant, the drift parameters are much higher than
those for the poorly performing individual. These differences persisted over
all blocks with small gradual increase over time.

Figure 7: Results for PTC1 data: Estimated posterior mean trajectories for
individual specific drifts
$\mu_{d,s}^{(i)}(t)=\exp\\{f_{d,s}(t)+u_{C}^{(i)}(t)\\}$ for successful
identification $(d=s)$ for two different participants \- one performing well
(dashed line) and one performing poorly (dotted line). The shaded areas
represent the corresponding $95\%$ point wise credible intervals. Estimates
related to high-flat tone responses are shown in red; low-rising tone
responses in blue; low-dipping tone responses in green; and high-falling tone
responses in purple.

##### Analysis of Benchmark Data.

To validate the proposed method, we also analyzed a tone learning data which,
in addition to response accuracies, included accurate measurements of the
response times. It was previously analyzed in Paulon et al. (2020) using the
drift-diffusion model (2) which allowed inference on both the drift and the
boundary parameters. In our analysis with the method proposed here, we ignored
the response times. We observed that the estimates of the drifts produced by
our proposed methodology match well with the estimates obtained by Paulon et
al. (2020). A description of this ‘benchmark’ data set and other details of
our analyses are provided in Section S.4 of the supplementary material.

##  7 Discussion

In this article, we proposed a novel longitudinal inverse-probit mixed
categorical probability model. Our research was motivated by category learning
experiments where scientists are interested in understanding how the decision
making mechanisms evolve as the participants get more training and experience
and have only recorded data on response categories but have not measured
response times. Biologically interpretable latent drift-diffusion process
based categorical probability models for such scenarios, useful for drawing
inferences about the underlying decision making cognitive dynamics, had never
been considered in the literature before. We addressed the need, building
carefully on our previous work on longitudinal drift-diffusion mixed models,
resulting in a new class of category probability model we referred to as the
inverse-probit model. We explored parameter recoverability in such models,
showing, in particular, that the offset parameters can not be recovered and
drifts and boundaries both can not be recovered from data only on response
categories. In our analyses, we thus focused on estimating the biologically
more important drift parameters but kept the offsets and the boundaries fixed.
These results show that if a researcher is interested in understanding the
detailed interplay between drifts and boundaries behind category learning,
additional data on response times must also be collected. On the positive
side, our results also suggest that, with carefully designed models, the
general longitudinal trends in the drifts can still be estimated well from
data only on response categories.

While we focused here on studying auditory category learning, the method
proposed here is applicable to other domains of behavioral neuroscience
research studying categorical decision making when the response times are
either not measured or are not reliable. The projection based solutions
proposed to address sensitivity and clustering challenges here are not
restricted to inverse-probit models but can be easily adapted to other classes
of generalized linear models such as the widely popular logit and probit
models and hence may also be of interest to a much broader audience.

## Appendix

##  Appendix Appendix A Proof of Lemma 1

###### Proof.

It is easy to check that the offset parameters $\delta_{s}$ are not
identifiable since

$\displaystyle P(d\mid
s,\delta_{s},\mbox{\boldmath$\mu$}_{1:d_{0},s},{\mathbf{b}}_{1:d_{0},s})=\int_{\delta_{s}}^{\infty}g(\tau\mid\delta_{s},\mu_{d,s},b_{d,s})\prod_{d^{\prime}\neq
d}\left\\{1-G(\tau\mid\delta_{s},\mu_{d^{\prime},s},b_{d^{\prime},s})\right\\}d\tau$
$\displaystyle=\int_{0}^{\infty}g(\tau\mid
0,\mu_{d,s},b_{d,s})\prod_{d^{\prime}\neq d}\left\\{1-G(\tau\mid
0,\mu_{d^{\prime},s},b_{d^{\prime},s})\right\\}d\tau=P(d\mid
s,0,\mbox{\boldmath$\mu$}_{1:d_{0},s},{\mathbf{b}}_{1:d_{0},s}).$

Next we will show that the drift parameters and decision boundaries are not
separately identifiable, even if we fix offset parameters to a constant.

First note that equation (3) can also be represented as

$\displaystyle\int_{\delta_{s}}^{\infty}\ldots\int_{\delta_{s}}^{\infty}\prod_{d^{\prime}\neq
d}g(\tau_{d^{\prime}}\mid\mbox{\boldmath$\theta$}_{d^{\prime},s})\int_{\delta_{s}}^{\displaystyle\wedge_{d\neq
d^{\prime}}\tau_{d^{\prime}}}g(\tau_{d}\mid\mbox{\boldmath$\theta$}_{d,s})d\tau_{d}\prod_{d^{\prime}\neq
d}d\tau_{d^{\prime}}.$ (A.1)

First observe that $\tau^{\star}=\wedge_{d^{\prime}\neq
d}\tau_{d^{\prime}}=\tau_{-1}^{\star}\wedge\tau_{1}$, where
$\tau_{-1}^{\star}=\wedge_{d^{\prime}\neq\\{1,d\\}}\tau_{d^{\prime}}$. Thus
the integral above can be written as

$\displaystyle\int_{\delta_{s}}^{\infty}\ldots\int_{\delta_{s}}^{\infty}\prod_{d^{\prime}\neq\\{1,d\\}}g(\tau_{d^{\prime}}\mid\mbox{\boldmath$\theta$}_{d^{\prime},s})\left\\{\int_{\delta_{s}}^{\infty}g(\tau_{1}\mid\mbox{\boldmath$\theta$}_{1,s})\int_{\delta_{s}}^{\tau_{-1}^{\star}\wedge\tau_{1}}g(\tau_{d}\mid\mbox{\boldmath$\theta$}_{d,s})d\tau_{d}\right\\}\prod_{d^{\prime}\neq\\{1,d\\}}d\tau_{d^{\prime}}$
$\displaystyle=\int_{\delta_{s}}^{\infty}\ldots\int_{\delta_{s}}^{\infty}\prod_{d^{\prime}\neq\\{1,d\\}}g(\tau_{d^{\prime}}\mid\mbox{\boldmath$\theta$}_{d^{\prime},s})\left\\{\int_{\delta_{s}}^{\tau_{-1}^{\star}}g(\tau_{d}\mid\mbox{\boldmath$\theta$}_{d,s})\int_{\tau_{d}}^{\infty}g(\tau_{1}\mid\mbox{\boldmath$\theta$}_{1,s})d\tau_{1}d\tau_{d}\right\\}\prod_{d^{\prime}\neq\\{1,d\\}}d\tau_{d^{\prime}}.$

Proceeding sequentially one can show that the integral above is same as the
integral in (3).

Using the above we express the probability in (3) as in (A.1). As the offset
parameter $\delta_{s}$ is already shown to be not identifiable, we need to fix
the same. Without loss of generality, we fix the offset parameter at $0$. The
probability density function of inverse Gaussian distribution, with parameters
$\mbox{\boldmath$\theta$}_{d,s}=(\mu_{d,s},b_{d,s})$ evaluated at $\tau_{d}$,
$g(\tau_{d}\mid\mbox{\boldmath$\theta$}_{d,s})$ can be obtained from (1) by
replacing $\delta_{s}=0$.

Consider the transformation of $\tau_{d}$ to $\tau_{d}^{\star}$ as
$\tau_{d}=c^{2}\tau_{d}^{\star}$, for some constant $c>0$, and for all $d$.
Further, define $b_{d,s}^{\star}=b_{d,s}/c$ and
$\mu_{d,s}^{\star}=c\mu_{d,s}$, for $d$. Then observe that

$\displaystyle
g(\tau_{d}\mid\mbox{\boldmath$\theta$}_{d,s})d\tau_{d}=(2\pi)^{-1/2}b_{d,s}^{\star}(\tau_{d}^{\star})^{-3/2}\exp\left\\{-(2\tau_{d}^{\star})^{-1}\left(b_{d,s}^{\star}-\mu_{d,s}^{\star}\tau_{j}^{\star}\right)^{2}\right\\}d\tau_{d}^{\star}=g(\tau_{d}^{\star}\mid\mbox{\boldmath$\theta$}_{d,s}^{\star})d\tau_{d}^{\star},$

$g(\tau_{d}^{\star}\mid\mbox{\boldmath$\theta$}_{d,s}^{\star})$ where is the
pdf of inverse Gaussian distribution with parameters $\mu_{d,s}^{\star}$ and
$b_{d,s}^{\star}$, evaluated at the point $\tau_{d}^{\star}$.

Applying the transformation on $\tau_{d}$ for all $d$ we get that the integral
in (A.1) with $\delta_{s}=0$ is same as

$\displaystyle\int_{0}^{\infty}\ldots\int_{0}^{\infty}\prod_{d^{\prime}\neq
d}g(\tau_{d^{\prime}}^{\star}\mid\mbox{\boldmath$\theta$}_{d^{\prime},s}^{\star})\int_{0}^{\displaystyle\wedge_{d\neq
d^{\prime}}\tau_{d^{\prime}}^{\star}}g(\tau_{d}^{\star}\mid\mbox{\boldmath$\theta$}_{d,s}^{\star})d\tau_{d}^{\star}\prod_{d^{\prime}\neq
d}d\tau_{d^{\prime}}^{\star}.$

As $c$ is arbitrary, this shows that the drifts and boundaries are not
separately estimable. ∎

##  Appendix Appendix B Proof of Theorem 1

###### Proof.

Let
${\mathbf{P}}(\mbox{\boldmath$\mu$}_{1:d_{0},s})=\\{p_{1}(\mbox{\boldmath$\mu$}_{1:d_{0},s}),\dots,p_{d_{0}}(\mbox{\boldmath$\mu$}_{1:d_{0},s})\\}^{\rm
T}$ be the function, given by (4), from ${\cal S}_{0,k}$ to unit probability
simplex $\Delta^{d_{0}-1}$. For notational simplicity, we write
$\mbox{\boldmath$\mu$}_{1:d_{0},s}=\mbox{\boldmath$\mu$}=(\mu_{1},\dots,\mu_{d_{0}})^{\rm
T}$. We first find the matrix of partial derivative $\nabla{\mathbf{P}}$ with
respect to $\mu$.

For $\mbox{\boldmath$\mu$}\in{\cal S}_{0,k}$, ${\bf
1}^{T}\mbox{\boldmath$\mu$}=k$, and hence the likelihood reduces to

$\displaystyle
p_{d}\left(\mbox{\boldmath$\mu$}\right)=\displaystyle\frac{\left(be^{b}\right)^{d_{0}}}{(2\pi)^{d_{0}/2}}\int_{0}^{\infty}\int_{\tau_{d}}^{\infty}\cdots\int_{\tau_{d}}^{\infty}|\mbox{\boldmath$\tau$}|^{-3/2}\exp\left\\{-\frac{1}{2}\left({\bf
1}^{T}\mbox{\boldmath$\tau$}^{-1}{\bf
1}+\mbox{\boldmath$\mu$}^{T}\mbox{\boldmath$\tau$}{\mbox{\boldmath$\mu$}}\right)\right\\}d\mbox{\boldmath$\tau$}_{-d}d\tau_{d},$

$d=1,\ldots,d_{0}$, where
$\mbox{\boldmath$\tau$}=\hbox{diag}(\tau_{1},\ldots,\tau_{d})$. Next,
differentiating $p_{d}\left(\mbox{\boldmath$\mu$}\right)$ with respect to
$\mu$, we get

$\displaystyle\frac{\partial{p_{d}\left(\mbox{\boldmath$\mu$}\right)}}{\partial\mbox{\boldmath$\mu$}}\hskip
373.99747pt$
$\displaystyle~{}=\displaystyle\frac{\left(be^{b}\right)^{d_{0}}}{(2\pi)^{d_{0}/2}}\int_{0}^{\infty}\int_{\tau_{d}}^{\infty}\cdots\int_{\tau_{d}}^{\infty}|\mbox{\boldmath$\tau$}|^{-3/2}\left(-\mbox{\boldmath$\tau$}\mbox{\boldmath$\mu$}\right)\exp\left\\{-\frac{1}{2}\left({\bf
1}^{T}\mbox{\boldmath$\tau$}^{-1}{\bf
1}+\mbox{\boldmath$\mu$}^{T}\mbox{\boldmath$\tau$}{\mbox{\boldmath$\mu$}}\right)\right\\}d\mbox{\boldmath$\tau$}_{-d}d\tau_{d},$
$\displaystyle~{}=\begin{bmatrix}\mu_{1}\eta_{2}~{}&\cdots&~{}\mu_{d-1}\eta_{2}&\mu_{d}\eta_{1}&\mu_{d+1}\eta_{2}~{}&\cdots&~{}\mu_{d_{0}}\eta_{2}\end{bmatrix}^{T},\hskip
139.11996pt$

where
$\eta_{1}=-E\left\\{\tau_{1}I\left(\tau_{2}>\tau_{1},\cdots,\tau_{d_{0}}>\tau_{1}\right)\left|\mbox{\boldmath$\mu$}\right.\right\\}$,
and
$\eta_{2}=-E\left\\{\tau_{2}I\left(\tau_{2}>\tau_{1},\cdots,\tau_{d_{0}}>\tau_{1}\right)\left|\mbox{\boldmath$\mu$}\right.\right\\}$,
and $I(A)$ is the indicator function of the event $A$. Here the expectation is
considered under the joint distribution of
$\left(\tau_{1},\ldots,\tau_{d}\right)$, which is independent inverse
Gaussian. Clearly $\eta_{1}>\eta_{2}>0$.

From the above derivation it is easy to obtain that

$\displaystyle\nabla{\mathbf{P}}\left(\mbox{\boldmath$\mu$}\right)=\begin{bmatrix}\mu_{1}\eta_{1}&\mu_{2}\eta_{2}&\cdots&\mu_{d_{0}}\eta_{2}\\\
\mu_{1}\eta_{2}&\mu_{2}\eta_{1}&\cdots&\mu_{d_{0}}\eta_{2}\\\
\vdots&\vdots&\cdots&\vdots\\\
\mu_{1}\eta_{2}&\mu_{2}\eta_{2}&\cdots&\mu_{d_{0}}\eta_{1}\end{bmatrix}={\mathbf{M}}\left\\{\left(\eta_{1}-\eta_{2}\right)I+\eta_{2}{\bf
1}{\bf 1}^{T}\right\\},$

where ${\mathbf{M}}=\hbox{diag}\left(\mu_{1},\ldots,\mu_{d_{0}}\right)$.

Now, suppose there exists $\mu$ and $\nu$ in ${\cal S}_{k}$ such that
$\mbox{\boldmath$\mu$}\neq\mbox{\boldmath$\nu$}$ and
${\mathbf{P}}\left(\mbox{\boldmath$\mu$}\right)={\mathbf{P}}\left(\mbox{\boldmath$\nu$}\right)$.
Define $\mbox{\boldmath$\gamma$}:[0,1]\rightarrow\mathbb{R}^{d_{0}}$ such that
$\mbox{\boldmath$\gamma$}(t)=\mbox{\boldmath$\mu$}+t\left(\mbox{\boldmath$\nu$}-\mbox{\boldmath$\mu$}\right)$,
$t\in[0,1]$. Further, define
$h(t)=\left<{\mathbf{P}}\left(\mbox{\boldmath$\gamma$}(t)\right)-{\mathbf{P}}\left(\mbox{\boldmath$\mu$}\right),\mbox{\boldmath$\nu$}-\mbox{\boldmath$\mu$}\right>$,
as the cross product of
${\mathbf{P}}\left(\mbox{\boldmath$\gamma$}(t)\right)-{\mathbf{P}}\left(\mbox{\boldmath$\mu$}\right)$
and $\mbox{\boldmath$\nu$}-\mbox{\boldmath$\mu$}$. Then $h(1)=h(0)=0$ under
the proposition that
${\mathbf{P}}\left(\mbox{\boldmath$\mu$}\right)={\mathbf{P}}\left(\mbox{\boldmath$\nu$}\right)$.
Therefore, by the Mean Value Theorem, as
$\mbox{\boldmath$\mu$}\neq\mbox{\boldmath$\nu$}$, there exists some point
$c\in(0,1)$ such that $\left.\partial h(t)/\partial t\right|_{t=c}=0$. Now,

$\displaystyle\frac{\partial h(t)}{\partial t}$ $\displaystyle=$
$\displaystyle\sum_{d=1}^{d_{0}}\left(\nu_{d}-\mu_{d}\right)\frac{\partial}{\partial
t}\left[p_{d}\left\\{\mbox{\boldmath$\gamma$}(t)\right\\}-p_{d}\left(\mbox{\boldmath$\mu$}\right)\right]$
$\displaystyle=$
$\displaystyle\sum_{d=1}^{d_{0}}\left(\nu_{d}-\mu_{d}\right)\left\\{\frac{\partial}{\partial\mbox{\boldmath$\gamma$}}p_{d}\left(\mbox{\boldmath$\gamma$}\right)\right\\}^{T}\frac{\partial\mbox{\boldmath$\gamma$}(t)}{\partial
t}$ $\displaystyle=$
$\displaystyle\left(\mbox{\boldmath$\nu$}-\mbox{\boldmath$\mu$}\right)^{T}\nabla{\mathbf{P}}\\{\mbox{\boldmath$\gamma$}(t)\\}\left(\mbox{\boldmath$\nu$}-\mbox{\boldmath$\mu$}\right)$
$\displaystyle=$
$\displaystyle\left(\eta_{1}-\eta_{2}\right)\left(\mbox{\boldmath$\nu$}-\mbox{\boldmath$\mu$}\right)^{T}\mbox{\boldmath$\Gamma$}(t)\left(\mbox{\boldmath$\nu$}-\mbox{\boldmath$\mu$}\right)+\eta_{2}\left(\mbox{\boldmath$\nu$}-\mbox{\boldmath$\mu$}\right)^{T}M{\bf
1}{\bf 1}^{T}\left(\mbox{\boldmath$\nu$}-\mbox{\boldmath$\mu$}\right)$
$\displaystyle=$
$\displaystyle\left(\eta_{1}-\eta_{2}\right)\left(\mbox{\boldmath$\nu$}-\mbox{\boldmath$\mu$}\right)^{T}\mbox{\boldmath$\Gamma$}(t)\left(\mbox{\boldmath$\nu$}-\mbox{\boldmath$\mu$}\right),$

as ${\bf 1}^{T}\left(\mbox{\boldmath$\nu$}-\mbox{\boldmath$\mu$}\right)=0$,
where
$\mbox{\boldmath$\Gamma$}(t)=\hbox{diag}\\{\mbox{\boldmath$\gamma$}(t)\\}$.

As every component of $\mu$ and $\nu$ is positive, for any $c\in(0,1)$, the
matrix $\mbox{\boldmath$\Gamma$}(c)$ is positive definite. Further, as
$\eta_{1}>\eta_{2}$, $\left.\partial h(t)/\partial t\right|_{t=c}=0$ only if
$\mbox{\boldmath$\mu$}=\mbox{\boldmath$\nu$}$, which contradicts the
proposition. ∎

##  Appendix Appendix C Algorithm for Minimal Distance Mapping

The problem of finding projection of a point $\mu$ onto the space ${\cal
S}_{k,\varepsilon}$ is equivalent to the following non-linear optimization
problem:

$\displaystyle\mathrm{minimize}_{{\mathbf{w}}}\|{{\mathbf{w}}}-\mbox{\boldmath$\mu$}\|^{2}\qquad\mbox{such
that}\quad\sum_{i=1}^{d_{0}}w_{i}=k,~{}w_{i}\geq\varepsilon.$

Duchi et al. (2008, Algorithm 1) provides a solution to the problem of
projection of a given point $\mu$ onto the space ${\cal S}_{k,\varepsilon}$
for $\varepsilon=0$, which is modified for any given $\varepsilon$ below.

Algorithm 2

1:

2:INPUT: A vector $\mu$, and a pair $(k,\varepsilon)$.

3:

4:Sort $\mu$ into $\mbox{\boldmath$\mu$}^{\star}$ such that the elements of
$\mbox{\boldmath$\mu$}^{\star}$ is in descending order.

5:Find
$\rho=\max\left\\{j:\mu^{\star}_{j}-j^{-1}\left(\sum_{l=1}^{j}\mu_{l}^{\star}-k\right)>\varepsilon\right\\}$.

6:Define
$\theta=\rho^{-1}\left\\{\sum_{l=1}^{j}\mu_{l}^{\star}-k+(k-\rho)\varepsilon\right\\}$.

7:

8:OUTPUT: ${{\mathbf{w}}}$ such that
$w_{i}=\max\left\\{\mu_{i}-\theta,\varepsilon\right\\}$.

9:

##  Appendix Appendix D Proof of Lemma 3

###### Proof.

We consider the unconditional distribution of $\tau_{1:d_{0}}$, given the
parameters $\mu_{1:d_{0}}$ as the proposal distribution, $g$. Clearly, the
proposal distribution $g$ and the target conditional joint distribution $f$
satisfies $f(\tau_{1:d_{0}}|\mu_{1:d_{0}})/g(\tau_{1:d_{0}}|\mu_{1:d_{0}})\leq
M$, where $M^{-1}=P\left(\tau_{d}\leq\tau_{1:d_{0}}\right)$. Therefore, for
any random sample $U\sim U(0,1)$, $f(\tau_{1:d_{0}}|\mu_{1:d_{0}})\geq
MUg(\tau_{1:d_{0}}|\mu_{1:d_{0}})$ if the sample satisfies the condition
$\tau_{d}\leq\tau_{1:d_{0}}$, and
$f(\tau_{1:d_{0}}|\mu_{1:d_{0}})<MUg(\tau_{1:d_{0}}|\mu_{1:d_{0}})$ otherwise.
Hence by Lemma 2.3.1 of Robert and Casella (2004), algorithm above produces
samples from the target distribution. ∎

## Supplementary Materials

The supplementary materials detail the choice of the prior hyper-parameters,
the MCMC algorithm used to sample from the posterior and some performance
diagnostics, and the analysis of a real benchmark data set. Separate files
additionally include R programs implementing the longitudinal inverse-probit
mixed model developed in this article and the PTC1 data set analyzed in
Section 6.

## Acknowledgments

This work was supported by the National Science Foundation grant NSF
DMS-1953712 and National Institute on Deafness and Other Communication
Disorders grants R01DC013315 and R01DC015504, awarded to Sarkar and
Chandrasekaran.

## References

  * Agresti (2018) Agresti, A. (2018).  An introduction to categorical data analysis.  Wiley.
  * Albert and Chib (1993) Albert, J. H. and Chib, S. (1993).  Bayesian analysis of binary and polychotomous response data.  Journal of the American statistical Association, 88, 669–679.
  * Ashby et al. (2003) Ashby, F. G., Noble, S., Filoteo, J. V., Waldron, E. M., and Ell, S. W. (2003).  Category learning deficits in parkinson’s disease.  Neuropsychology, 17, 115.
  * Beck (2017) Beck, A. (2017).  First-order methods in optimization.  SIAM.
  * Bogacz et al. (2010) Bogacz, R., Wagenmakers, E.-J., Forstmann, B. U., and Nieuwenhuis, S. (2010).  The neural basis of the speed-accuracy tradeoff.  Trends in Neurosciences, 33, 10–16.
  * Borooah (2002) Borooah, V. K. (2002).  Logit and probit: ordered and multinomial models.  Sage.
  * Brody and Hanks (2016) Brody, C. D. and Hanks, T. D. (2016).  Neural underpinnings of the evidence accumulator.  Current Opinion in Neurobiology, 37, 149–157.
  * Brown and Heathcote (2008) Brown, S. D. and Heathcote, A. (2008).  The simplest complete model of choice response time: Linear ballistic accumulation.  Cognitive Psychology, 57, 153–178.
  * Burgette and Hahn (2010) Burgette, L. F. and Hahn, P. R. (2010).  Symmetric Bayesian multinomial probit models.  Duke University Statistical Science Technical Report, pages 1–20.
  * Burgette and Nordheim (2012) Burgette, L. F. and Nordheim, E. V. (2012).  The trace restriction: an alternative identification strategy for the Bayesian multinomial probit model.  Journal of Business & Economic Statistics, 30, 404–410.
  * Cavanagh et al. (2011) Cavanagh, J. F., Wiecki, T. V., Cohen, M. X., Figueroa, C. M., Samanta, J., Sherman, S. J., and Frank, M. J. (2011).  Subthalamic nucleus stimulation reverses mediofrontal influence over decision threshold.  Nature Neuroscience, 14, 1462.
  * Chandrasekaran et al. (2014) Chandrasekaran, B., Yi, H.-G., and Maddox, W. T. (2014).  Dual-learning systems during speech category learning.  Psychonomic Bulletin & Review, 21, 488–495.
  * Chhikara (1988) Chhikara, R. (1988).  The inverse Gaussian distribution: Theory, methodology, and applications.  CRC Press.
  * Cox and Miller (1965) Cox, D. R. and Miller, H. D. (1965).  The theory of stochastic processes.  CRC Press.
  * de Boor (1978) de Boor, C. (1978).  A practical guide to splines.  Springer-Verlag.
  * Deo (2018) Deo, S. (2018).  Algebraic topology, volume 27 of Texts and Readings in Mathematics.  Hindustan Book Agency, New Delhi.
  * Ding and Gold (2013) Ding, L. and Gold, J. I. (2013).  The basal ganglia’s contributions to perceptual decision making.  Neuron, 79, 640–649.
  * Duchi et al. (2008) Duchi, J., Shalev-Shwartz, S., Singer, Y., and Chandra, T. (2008).  Efficient projections onto the l 1-ball for learning in high dimensions.  In Proceedings of the 25th international conference on Machine learning, pages 272–279.
  * Dufau et al. (2012) Dufau, S., Grainger, J., and Ziegler, J. C. (2012).  How to say “no” to a nonword: A leaky competing accumulator model of lexical decision.  Journal of Experimental Psychology: Learning, Memory, and Cognition, 38, 1117.
  * Eilers and Marx (1996) Eilers, P. H. and Marx, B. D. (1996).  Flexible smoothing with b-splines and penalties.  Statistical Science, 11, 89–102.
  * Filoteo et al. (2010) Filoteo, J. V., Lauritzen, S., and Maddox, W. T. (2010).  Removing the frontal lobes: The effects of engaging executive functions on perceptual category learning.  Psychological science, 21, 415–423.
  * Glimcher and Fehr (2013) Glimcher, P. W. and Fehr, E. (2013).  Neuroeconomics: Decision making and the brain.  Academic Press.
  * Gold and Shadlen (2007) Gold, J. I. and Shadlen, M. N. (2007).  The neural basis of decision making.  Annual Review of Neuroscience, 30, 535–574.
  * Heekeren et al. (2004) Heekeren, H. R., Marrett, S., Bandettini, P. A., and Ungerleider, L. G. (2004).  A general mechanism for perceptual decision-making in the human brain.  Nature, 431, 859.
  * Hubert and Arabie (1985) Hubert, L. and Arabie, P. (1985).  Comparing partitions.  Journal of Classification, 2, 193–218.
  * Johndrow et al. (2013) Johndrow, J., Dunson, D., and Lum, K. (2013).  Diagonal orthant multinomial probit models.  In Artificial Intelligence and Statistics, pages 29–38.
  * Kim et al. (2017) Kim, S., Potter, K., Craigmile, P. F., Peruggia, M., and Van Zandt, T. (2017).  A Bayesian race model for recognition memory.  Journal of the American Statistical Association, 112, 77–91.
  * Leite and Ratcliff (2010) Leite, F. P. and Ratcliff, R. (2010).  Modeling reaction time and accuracy of multiple-alternative decisions.  Attention, Perception, & Psychophysics, 72, 246–273.
  * Lu (1995) Lu, J. (1995).  Degradation processes and related reliability models.  Ph.D. thesis, McGill University, Montreal, Canada.
  * McHaney et al. (2021) McHaney, J. R., Tessmer, R., Roark, C. L., and Chandrasekaran, B. (2021).  Working memory relates to individual differences in speech category learning: Insights from computational modeling and pupillometry.  Brain and Language, 22, 1–15.
  * Milosavljevic et al. (2010) Milosavljevic, M., Malmaud, J., Huth, A., Koch, C., and Rangel, A. (2010).  The drift diffusion model can account for the accuracy and reaction time of value-based choices under high and low time pressure.  Judgment and Decision Making, 5, 437–449.
  * Morris (2015) Morris, J. S. (2015).  Functional regression.  Annual Review of Statistics and Its Application, 2, 321–359.
  * Paulon et al. (2020) Paulon, G., Llanos, F., Chandrasekaran, B., and Sarkar, A. (2020).  Bayesian semiparametric longitudinal drift-diffusion mixed models for tone learning in adults.  Journal of the American Statistical Association,, 116, 1114–1127.
  * Purcell (2013) Purcell, B. A. (2013).  Neural mechanisms of perceptual decision making.  Vanderbilt University.
  * Ramsay and Silverman (2007) Ramsay, J. O. and Silverman, B. W. (2007).  Applied functional data analysis: methods and case studies.  Springer.
  * Rand (1971) Rand, W. M. (1971).  Objective criteria for the evaluation of clustering methods.  Journal of the American Statistical Association, 66, 846–850.
  * Ratcliff (1978) Ratcliff, R. (1978).  A theory of memory retrieval.  Psychological review, 85, 59.
  * Ratcliff and McKoon (2008) Ratcliff, R. and McKoon, G. (2008).  The diffusion decision model: Theory and data for two-choice decision tasks.  Neural Computation, 20, 873–922.
  * Ratcliff and Rouder (1998) Ratcliff, R. and Rouder, J. N. (1998).  Modeling response times for two-choice decisions.  Psychological Science, 9, 347–356.
  * Ratcliff et al. (2016) Ratcliff, R., Smith, P. L., Brown, S. D., and McKoon, G. (2016).  Diffusion decision model: Current issues and history.  Trends in Cognitive Sciences, 20, 260–281.
  * Reetzke et al. (2018) Reetzke, R., Xie, Z., Llanos, F., and Chandrasekaran, B. (2018).  Tracing the trajectory of sensory plasticity across different stages of speech learning in adulthood.  Current Biology, 28, 1419–1427.
  * Roark et al. (2021) Roark, C. L., Smayda, K. E., and Chandrasekaran, B. (2021).  Auditory and visual category learning in musicians and nonmusicians.  Journal of Experimental Psychology: General.
  * Robert and Casella (2004) Robert, C. P. and Casella, G. (2004).  Monte Carlo statistical methods.  Springer Texts in Statistics. Springer-Verlag, New York, second edition.
  * Ross et al. (1996) Ross, S. M., Kelly, J. J., Sullivan, R. J., Perry, W. J., Mercer, D., Davis, R. M., Washburn, T. D., Sager, E. V., Boyce, J. B., and Bristow, V. L. (1996).  Stochastic processes.  Wiley New York.
  * Rudin (1991) Rudin, W. (1991).  Functional analysis.  International Series in Pure and Applied Mathematics. McGraw-Hill, Inc., New York, second edition.
  * Schall (2001) Schall, J. D. (2001).  Neural basis of deciding, choosing and acting.  Nature Reviews Neuroscience, 2, 33.
  * Sen et al. (2021) Sen, D., Patra, S., and Dunson, D. (2021).  Constrained inference through posterior projections.  arXiv preprint arXiv:1812.05741.
  * Smith and Ratcliff (2004) Smith, P. L. and Ratcliff, R. (2004).  Psychology and neurobiology of simple decisions.  Trends in Neurosciences, 27, 161–168.
  * Smith and Vickers (1988) Smith, P. L. and Vickers, D. (1988).  The accumulator model of two-choice discrimination.  Journal of Mathematical Psychology, 32, 135–168.
  * Usher and McClelland (2001) Usher, M. and McClelland, J. L. (2001).  The time course of perceptual choice: The leaky, competing accumulator model.  Psychological Review, 108, 550.
  * Wang et al. (2016) Wang, J.-L., Chiou, J.-M., and Müller, H.-G. (2016).  Functional data analysis.  Annual Review of Statistics and Its Application, 3, 257–295.
  * Whitmore and Seshadri (1987) Whitmore, G. and Seshadri, V. (1987).  A heuristic derivation of the inverse gaussian distribution.  The American Statistician, 41, 280–281.

Supplementary Materials for
Bayesian Semiparametric Longitudinal
Inverse-Probit Mixed Models
for Category Learning

Minerva<EMAIL_ADDRESS>
Jacie McHaney2<EMAIL_ADDRESS>
Bharath<EMAIL_ADDRESS>
Abhra Sarkar3<EMAIL_ADDRESS>

1Department of of Mathematics and Statistics,

Indian Institute of Technology,

Kanpur, UP 208016, India

2Department of Communication Science and Disorders,

University of Pittsburgh,

4028 Forbes Tower, Pittsburgh, PA 15260, USA

3Department of Statistics and Data Sciences,

University of Texas at Austin,

105 East 24th Street D9800, Austin, TX 78712, USA

The supplementary materials detail the choice of the prior hyper-parameters,
the MCMC algorithm used to sample from the posterior and some performance
diagnostics, and the analysis of a real benchmark data set. Separate files
additionally include R programs implementing the longitudinal inverse-probit
mixed model developed in this article and the PTC1 data set analyzed in
Section 6 in the main paper.

##  Appendix S.1 Prior Hyper-parameters and Initializations

The random effects of the inverse-probit mixed model (5) are all initialized
at zero. The variance and smoothing parameters are initially set to $0.1$
each. The location parameter of the prior on
$\mbox{\boldmath$\beta$}_{z}^{\star}$, $\mbox{\boldmath$\mu$}_{\beta,0}$ is
set to $(1,\ldots,1)$. This choice of $\mbox{\boldmath$\beta$}_{z}^{\star}$
would set the expected value of $\mu_{x}^{(i)}(t)$ to $1$, which is supported
empirically. The value of the parameter $\alpha$ is set to $1$.

##  Appendix S.2 Posterior Inference

Posterior inference for the longitudinal drift-diffusion mixed model,
described in Section 3 in the main paper, is based on samples drawn from the
posterior using an MCMC algorithm. The algorithm carefully exploits the
conditional independence relationships encoded in the model as well as the
latent variable construction of the model. Sampling the latent inverse-
gaussian distributed response times, in particular, greatly simplifies
computation.

In what follows, $\zeta$ denotes a generic variable that collects all other
variables not explicitly mentioned, including the data points. Also, $p_{0}$
will sometimes be used as a generic for a prior distribution without
explicitly mentioning its hyper-parameters. The notation $x$ is used to
abbreviate $\\{d,s\\}$. The sampler for the inverse-probit mixed model of
Section 3 iterates between the following steps.

  1. 1.

Sampling $\tau_{1:d_{0}}^{(i,l)}(t)$: Suppose the $i$-th individual selects
the output tone $d$, in the $t$-th block, $l$-th trial, given the input tone
$s$. Then $\tau_{1}^{(i,l)}(t),\ldots,\tau_{d_{0}}^{(i,l)}(t)$ is generated as
in Algorithm 1 (see Section 4) from the joint distribution of
$\tau_{1}^{(i,l)}(t),\ldots,\tau_{d_{0}}^{(i,l)}(t)$ given
$\mu_{1,s}^{(i)}(t),\ldots,\mu_{d_{0},s}^{(i)}(t)$, followed by an accept
reject step.

  2. 2.

Updating the components of fixed effects $f_{x}(t)$:

    1. (a)

The latent variable ${\mathbf{z}}_{x}$, indicating the group identities of
$\mbox{\boldmath$\beta$}_{x}$, follows multinomial distribution with
$z_{\max}$ labels and probabilities
$P(z_{x}=z|\mbox{\boldmath$\zeta$}),~{}z=1,\ldots,z_{\max}$ a posteriori. The
probability $P(z_{x}=z|\mbox{\boldmath$\zeta$})\propto\pi_{z}\times l_{z}$,
where $l_{z}$ is the likelihood of $\mbox{\boldmath$\beta$}_{x}$ evaluated at
$\mbox{\boldmath$\beta$}_{z}$. Let $L^{\star}$ be the set of all trials
corresponding to input-output tones $x=(s,d)$ for $i$-th individual and $t$-th
block, and $n_{x}^{(i)}(t)$ be the cardinality of $L^{\star}$. Furthermore,
let $\tau_{x}^{(i)}(t)=\sum_{l\in L^{\star}}\tau_{x}^{(i,l)}(t)$,
$\tau_{x}(t)=\sum_{i}\tau_{x}^{(i)}(t)$, and
$n_{x}(t)=\sum_{i}n_{x}^{(i)}(t)$. A little algebra shows that the likelihood
of $\mbox{\boldmath$\beta$}_{x}$ is Gaussian with variance matrix
$\mbox{\boldmath$\Sigma$}_{\beta,x}=\left\\{\sum_{t}\tau_{x}(t){\mathbf{B}}(t)^{T}{\mathbf{B}}(t)\right\\}^{-1}$,
and mean vector
$\mbox{\boldmath$\mu$}_{\beta,x}=\mbox{\boldmath$\Sigma$}_{\beta,x}\left\\{\sum_{t}{\mathbf{B}}(t)M_{x}(t)\right\\}$,
where $M_{x}(t)=2n_{x}(t)-\sum_{i}u_{x}^{(i)}(t)\tau_{x}^{(i)}(t)$. Therefore,
$l_{z}$ is the Gaussian likelihood with mean $\mbox{\boldmath$\mu$}_{\beta,x}$
and variance $\mbox{\boldmath$\Sigma$}_{\beta,x}$, evaluated at
$\mbox{\boldmath$\beta$}_{z}$.

    2. (b)

Let $N_{z}=\sum_{x}1(z_{x}=z)$, $z=1,\ldots,z_{\max}$, where $1(\cdot)$ is the
indicator function. Then the conditional posterior of
$\mbox{\boldmath$\pi$}_{z}$ is Dirichlet with parameters
$\alpha/z_{\max}+N_{1},\ldots,\alpha/z_{\max}+N_{z_{\max}}$.

    3. (c)

The full conditional posterior distribution of the coefficient atoms
$\mbox{\boldmath$\beta$}_{z}^{\star}$ is Gaussian with variance-covariance
matrix $\mbox{\boldmath$\Sigma$}_{\beta,z}^{\star}$ and
$\mbox{\boldmath$\mu$}_{\beta,z}^{\star}$, where
$\mbox{\boldmath$\Sigma$}_{\beta,z}^{\star,-1}=\sum_{x:z_{x}=z}\mbox{\boldmath$\Sigma$}_{\beta,x}^{-1}+\mbox{\boldmath$\Sigma$}_{\beta,0}^{-1}$,
and
$\mbox{\boldmath$\mu$}_{\beta,z}^{\star}=\mbox{\boldmath$\Sigma$}_{\beta,z}^{\star}\left[\sum_{x:z_{x}=z}\mbox{\boldmath$\Sigma$}_{\beta,x}^{-1}\mbox{\boldmath$\mu$}_{\beta,x}+\mbox{\boldmath$\Sigma$}_{\beta,0}^{-1}\mbox{\boldmath$\mu$}_{\beta,0}\right],$
where
$\mbox{\boldmath$\Sigma$}_{\beta,0}^{-1}=\left(\sigma_{a}^{-2}{\mathbf{I}}_{K}+\sigma_{s}^{-2}{\mathbf{P}}\right)$.

  3. 3.

Updating the components of random effects: We use the generic notation $U$ to
indicate the correct ($C$, i.e., $d=s$) or incorrect ($I$, i.e, $d\neq s$)
cases. Define $\tau_{U}^{(i)}(t)=\sum_{x:x\in U}\tau_{x}^{(i)}(t)$,
$n_{U}^{(i)}(t)=\sum_{x:x\in U}n_{x}^{(i)}(t)$,
$f\tau_{U}^{(i)}(t)=\sum_{x:x\in U}\tau_{x}^{(i)}(t)f_{x}(t)$,
$\mbox{\boldmath$\Sigma$}_{U,0}^{-1}=\sigma_{U,a}^{-2}{\mathbf{I}}_{K}+\sigma_{U,s}^{-2}{\mathbf{P}}$,
and
$\mbox{\boldmath$\Sigma$}_{U}^{(t)-1}=\sum_{t}\tau_{U}^{(i)}(t){\mathbf{B}}(t)^{T}{\mathbf{B}}(t)$.
The conditional posterior of $\mbox{\boldmath$\beta$}_{U}^{(i)}$ is Gaussian
with covariance
$\mbox{\boldmath$\Sigma$}_{U,\mathrm{post}}^{(i)}=\left\\{\mbox{\boldmath$\Sigma$}_{U,0}^{-1}+\mbox{\boldmath$\Sigma$}_{U}^{(i)-1}\right\\}^{-1}$,
and location parameter
$\mbox{\boldmath$\mu$}_{C}^{(i)}=\mbox{\boldmath$\Sigma$}_{U,\mathrm{post}}^{(i)}\mbox{\boldmath$\Sigma$}_{U}^{(i)-1}\left[\sum_{t}\left\\{2n_{U}^{(i)}(t)-f\tau_{U}^{(i)}(t)\right\\}{\mathbf{B}}(t)\right]$,
respectively.

  4. 4.

Updating the precision and smoothing parameters: The precision and smoothness
parameters involved in the fixed effects part are $\sigma_{a}^{2}$ and
$\sigma_{s}^{2}$, and those involved in the random effects part are
$\sigma_{U,a}$ and $\sigma_{U,s}^{2}$, $U=C,I$. We update these variance
components using Metropolis Hastings algorithm with log-normal proposal
distributions centered on the previous sample values.

  5. 5.

Estimation of probability: For each $(s,i,t)$, we calculate the probability of
selecting the $d$-th response in the following way: Let
$g\\{\cdot\mid\mu_{d,s}^{(i)}(t)\\}$ be the pdf of inverse Gaussian
distribution of the form (1) with parameters $\delta_{s}=0$, $b_{d,s}=2$ and
$\mu_{d,s}=\mu_{d,s}^{(i)}(t)$. We generate $M=2000$ independent samples
$\mbox{\boldmath$\tau$}_{m}=\left[\tau_{1,m},\ldots,\tau_{d_{0},m}\right]^{T}$,
$m=1,\ldots,M$, where $\tau_{d,m}$ is generated independently from
$g\\{\cdot\mid\mu_{d,s}^{(i)}(t)\\}$. Among these $M$ independent samples, the
proportion of occurrences of
$\\{\tau_{d,m}\leq\wedge_{d^{\prime}=1:d_{0}}\tau_{d^{\prime},m}\\}$ is
considered as the estimated probability of selecting $d$-th response.

The results reported in this article are all based on $5,000$ MCMC iterations
with the initial $2,000$ iterations discarded as burn-in. The remaining
samples were further thinned by an interval of $5$. We programmed in R. The
codes are available as part of the supplementary material. A ‘readme’ file,
providing additional details for a practitioner, is also included in the
supplementary material. In all experiments, the posterior samples produced
very stable estimates of the population and individual level parameters of
interest. MCMC diagnostic checks were not indicative of any convergence or
mixing issues.

##  Appendix S.3 MCMC Diagnostics

This section presents some convergence diagnostics for the MCMC sampler
described in the main manuscript. The results presented here are for the PTC1
data set. Diagnostics for the simulation experiments and the benchmark data
were similar and hence omitted.

Figure S.1: Analysis of PTC1 data: Trace plots (top) and auto-correlation
plots (bottom) of the individual drift rates $\mu_{1,1}^{(1)}(t)$
corresponding to the success categorization of tone $T_{1}$ evaluated at each
of the training blocks. In each panel, the solid red line shows the running
mean. Results for other drift parameters were very similar.

Figure S.1 shows the trace plots and auto-correlation of some individual level
parameters at different training blocks. These results are based on the MCMC
thinned samples. As these figures show, the running means are very stable and
there seems to be no convergence issues. Additionally, the Geweke test
(Geweke, 1992) for stationarity of the chains, which formally compares the
means of the first and last part of a Markov chain, was also performed. If the
samples are drawn from the stationary distribution of the chain, the two means
are equal and Geweke statistic has an asymptotically standard normal
distribution. The results of the test, reported in Table S.1, indicate that
convergence was satisfactory.

| $t=1$ | $t=2$ | $t=2$ | $t=2$ | $t=2$ | $t=2$
---|---|---|---|---|---|---
Geweke statistics | -1.233 | -0.392 | -0.678 | -0.136 | 0.440 | 0.339
$p$-value | 0.217 | 0.695 | 0.498 | 0.892 | 0.660 | 0.734
Table S.1: Geweke statistics and associated p-values assessing convergence of
the of the individual level drift parameters $\mu_{1,1}^{(1)}(t)$
corresponding to the success categorization of tone $T_{1}$ evaluated at each
of the training blocks. Results for other drift parameters were very similar.

##  Appendix S.4 Analysis of Benchmark Data

##### Description of the data.

The data set we consider next is a multi-day longitudinal speech category
training study reported previously in Reetzke et al. (2018) and analyzed
previously in Paulon et al. (2020). In this study, $n=20$ participants were
trained to learn $4$ tones, namely, high-level (T1), low-rising (T2), low-
dipping (T3), or high-falling (T4) tone, respectively. The trials were
administered in blocks, each comprising $40$ categorization trials.
Participants were trained across several days, with five blocks on each day.
On each trial, participants indicated the tone category they heard via button
press on a computer keyboard. Following the button press, they were given
corrective feedback.

The data consist of tone responses and associated response times for different
input tones for the $20$ participants. We focus here on the first two days of
training (10 blocks in total) as they exhibited the steepest improvement in
learning as well as the most striking individual differences relative to any
other collection of blocks.

##### Analysis.

We first demonstrate the performance of the proposed method in estimating the
probabilities associated with different $(d,s)$ pairs. Figure S.2 shows the
$95\%$ credible intervals for the estimated probabilities for different input
tones along with the average proportions of times an input tone was classified
into different tone categories across subjects.

Figure S.2: Results for the benchmark data: Estimated probability trajectories
compared with average proportions of times an input tone was classified into
different tone categories across subjects (in dashed line). High-flat tone
responses are shown in red; low-rising in blue; low-dipping in green; and
high-falling in purple.

Observe that, except in situations with very small number of data points the
$95\%$ credible intervals include the empirical probabilities. Further, the
estimated credible region is narrow enough implying high precision of the
inference.

Next consider the clustering results. We obtained two clusters each in pairs
of success combinations $(d=s)$, and in the wrong allocations $(d\neq s)$. The
clusters of success combinations are $S_{1}=\\{(1,1),(3,3)\\}$ and
$S_{2}=\\{(2,2),(4,4)\\}$, and that in wrong allocations are
$M_{1}=\\{(1,2),(2,1),(2,3),(3,2),(4,1),(4,2)\\}$, and
$M_{2}=\\{(1,3),(1,4),(2,4),(3,1),(3,4),(4,3)\\}$. The network plot in Figure
S.3 shows the stability of the clusters over the MCMC iterations.

From an overall perspective, the trajectory of ‘High-level’ ($T_{1}$) and
‘Low-dipping’ ($T_{3}$) are similar with two wrong allocations from $M_{2}$
and one from $M_{1}$, and that of ‘Low-rising’ ($T_{2}$) and ‘High-failing’
($T_{4}$) are similar with two wrong allocations from $M_{1}$ and one from
$M_{2}$. These similarities in the overall trajectories of $\\{T_{1},T_{3}\\}$
and $\\{T_{2},T_{4}\\}$ were also noted by Paulon et al. (2020).

Figure S.3: Results for the benchmark data: Network plot of similarity groups
showing the intra and inter cluster similarities. Each node is associated with
a pair indicating the input-response tone category, $(s,d)$. The number
associated with each edge indicates the proportion of times the pair in the
two connecting nodes appeared in the same cluster after burnin.

Next we consider the estimation of the underlying drift parameters
$\mu_{x}^{(i)}(t)$. Due to the identifiability constraints, the estimates of
$\mu_{x}^{(i)}(t)$ can only be observed in a relative scale. Figure
LABEL:fig:_population_effect shows the posterior mean trajectories and
associated $95\%$ credible intervals for the projected drift rates estimated
by our method for different combinations of $x=(d,s)$. On comparison with the
previous analysis of Paulon et al. (2020), the trajectories of our estimated
drift rates show significant similarity throughout.

Figure LABEL:fig:_random_effect shows the posterior mean trajectories and
associated $95\%$ credible intervals for the drift rates $\mu_{x}^{(i)}(t)$
for the different success combinations of $x=(d,s)$ for two participants \-
the one with the best accuracy averaged, and the one with the worst accuracy
averaged across all blocks. For the well performing participant, the drift
trajectories increase rapidly and for the poorly performing candidate on the
other hand, the drift trajectories increase very slowly. Once again, in spite
of the limitation of inferring on a relative scale, the relative differences
of the best and worst performing participants across blocks show great
similarity with the inference of Paulon et al. (2020).

##  Appendix S.5 Rand and Adjusted Rand Indices

##### Rand Index.

Given a set of $n$ objects $S=\\{s_{1},\ldots,s_{n}\\}$, let ${\bf
U}=\left\\{U_{1},\ldots,U_{R}\right\\}$ and ${\bf
V}=\left\\{V_{1},\ldots,V_{C}\right\\}$ represent two different partitions of
the objects in $S$ such that $\cup_{i=1}^{R}U_{i}=S=\cup_{j=1}^{C}V_{j}$ and
$U_{i}\cap U_{i^{\prime}}=\emptyset=V_{j}\cap V_{j^{\prime}}$ for $1\leq i\neq
i^{\prime}\leq R$ and $1\leq j\neq j^{\prime}\leq C$. Rand index estimates the
similarity between the allocations of $S$ in ${\bf U}$ and ${\bf V}$.

Let $a$ be the number of pairs of objects that are placed in the same
partition in ${\bf U}$ and same partition in ${\bf V}$, and $b$ be the number
of pairs of objects that are in different partitions of ${\bf U}$, as well as
in different partitions of ${\bf V}$. The quantities $a$ and $b$ can be
interpreted as agreements in ${\bf U}$ and ${\bf V}$, and the the total number
of pairs is $\binom{n}{2}$. The Rand index (Rand, 1971) is

${\tt RI}=(a+b)/\binom{n}{2}.$

The Rand index lies between 0 and 1. When the two partitions agree perfectly,
the ${\tt RI}$ takes the value 1.

##### Adjusted Rand Index.

The expected value of the Rand index of two random partitions does not take a
constant value. The adjusted Rand index (Hubert and Arabie, 1985) assumes
generalized hypergeometric distribution as the model of randomness, and makes
a base and scale change of the quantity $(a+b)$, defined above, so that the
resultant quantity is bounded by $[-1,1]$ and has expected value $0$ under
completely random allocation.

Let $n_{i,j}$ be the number of object that are both in $i^{th}$ partition of
${\bf U}$ and $j^{th}$ partition of ${\bf V}$, $n_{i}$ and $n_{j}$ be the
total number of components in $i^{th}$ partition of ${\bf U}$, and $j^{th}$
partition of ${\bf V}$, respectively.

The expression $a+d$ can be simplified to a linear transformation of
$\sum_{i,j}\binom{n_{i,j}}{2}$. Further, under the generalized hypergeometric
model, it can be shown that

$\mathbb{E}\left[\sum_{i,j}\binom{n_{ij}}{2}\right]=\left[\sum_{i}\binom{n_{i}}{2}\sum_{j}\binom{n_{j}}{2}\right]/\binom{n}{2}.$

Therefore, scaled the difference of linear transformed $(a+b)$ and its
expectation is the adjusted Rand index, defined as:

${\tt
ARI}=\frac{\sum_{i,j}\binom{n_{i,j}}{2}-\left[\sum_{i}\binom{n_{i}}{2}\sum_{j}\binom{n_{i,j}}{2}\right]/\binom{n}{2}}{\frac{1}{2}\left[\sum_{i}\binom{n_{i}}{2}+\sum_{j}\binom{n_{j}}{2}\right]-\left[\sum_{i}\binom{n_{i}}{2}\sum_{j}\binom{n_{j}}{2}\right]/\binom{n}{2}}.$

The expected value of ${\tt ARI}$ index is zero and the range is $[-1,1]$.
Like the RI, the ARI also takes the value $1$, when the two partitions agree
perfectly.

## References

  * Geweke (1992) Geweke, J. (1992).  Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments.  In Bayesian statistics, 4 (Peñíscola, 1991), pages 169–193. Oxford Univ. Press, New York.
  * Hubert and Arabie (1985) Hubert, L. and Arabie, P. (1985).  Comparing partitions.  Journal of Classification, 2, 193–218.
  * Paulon et al. (2020) Paulon, G., Llanos, F., Chandrasekaran, B., and Sarkar, A. (2020).  Bayesian semiparametric longitudinal drift-diffusion mixed models for tone learning in adults.  Journal of the American Statistical Association,, 116, 1114–1127.
  * Rand (1971) Rand, W. M. (1971).  Objective criteria for the evaluation of clustering methods.  Journal of the American Statistical Association, 66, 846–850.
  * Reetzke et al. (2018) Reetzke, R., Xie, Z., Llanos, F., and Chandrasekaran, B. (2018).  Tracing the trajectory of sensory plasticity across different stages of speech learning in adulthood.  Current Biology, 28, 1419–1427.
  * Robert and Casella (2004) Robert, C. P. and Casella, G. (2004).  Monte Carlo statistical methods.  Springer Texts in Statistics. Springer-Verlag, New York, second edition.