A state space approach to dynamic modeling of mouse-tracking data

TL;DR

This paper introduces a novel state-space modeling framework for mouse-tracking data that captures movement dynamics and experimental variables, offering a more integrated analysis method than traditional two-step procedures.

Contribution

It presents a unified stochastic state-space approach for analyzing mouse-tracking trajectories, improving the interpretation of real-time cognitive processes.

Findings

Enhanced analysis of mouse-tracking data using the proposed model

Improved understanding of movement dynamics in decision tasks

Case study demonstrating the model's effectiveness

Abstract

Mouse-tracking recording techniques are becoming very attractive in experimental psychology. They provide an effective means of enhancing the measurement of some real-time cognitive processes involved in categorization, decision-making, and lexical decision tasks. Mouse-tracking data are commonly analysed using a two-step procedure which first summarizes individuals' hand trajectories with independent measures, and then applies standard statistical models on them. However, this approach can be problematic in many cases. In particular, it does not provide a direct way to capitalize the richness of hand movement variability within a consistent and unified representation. In this article we present a novel, unified framework for mouse-tracking data. Unlike standard approaches to mouse-tracking, our proposal uses stochastic state-space modeling to represent the observed trajectories in terms of both individual movement dynamics and experimental variables. The model is estimated via a Metropolis-Hastings algorithm coupled with a non-linear recursive filter. The characteristics and potentials of the proposed approach are illustrated using a lexical decision case study. The results highlighted how dynamic modeling of mouse-tracking data can considerably improve the analysis of mouse-tracking tasks and the conclusions researchers can draw from them.

Tables

Table 1. Table 1: Application: Adequacy of the model to reproduce the observed matrices 𝐘 𝐘 \mathbf{Y} ( overall fit ) and 𝐘 i subscript 𝐘 𝑖 \mathbf{Y}_{i} ( by-subject fit ). All the measures are normalized in the range 0% (bad fit) - 100% (otimal fit). See Supplementary Materials for technical details.

	overall	by-subject
Model 1	84%	78%
Model 2	73%	70%
Model 3	75%	71%

Table 2. Table 2: Application: Posterior means ( μ 𝜇 \mu ), 95% posterior intervals ([ q 0.05 , q 0.975 subscript 𝑞 0.05 subscript 𝑞 0.975 q_{0.05},q_{0.975} ]), and Gelman-Rubin R ^ ^ 𝑅 \hat{R} index for the estimated parameters of Models 1-3.

			${\hat{\gamma }}_1$	${\hat{\gamma }}_2$	${\hat{\gamma }}_3$	$\hat{\eta }$	${\hat{\delta }}_1$	${\hat{\delta }}_2$
Model 1		$q_{0.05}$	1.224	1.234	1.211
		$\mu$	1.323	1.337	1.310
		$q_{0.975}$	1.443	1.457	1.432
		$\hat{R}$	1.003	1.002	1.003
Model 2		$q_{0.05}$				0.063
		$\mu$				0.078
		$q_{0.975}$				0.091
		$\hat{R}$				1.001
Model 3		$q_{0.05}$	0.083	1.130	1.217	0.305	-0.505	-0.437
		$\mu$	0.341	1.300	1.314	0.402	-0.385	-0.336
		$q_{0.975}$	0.605	1.468	1.411	0.500	-0.269	-0.235
		$\hat{R}$	1.008	1.013	1.012	1.011	1.013	1.010
All the models		${\hat{\kappa }}_1=22.31$
		${\hat{\kappa }}_2=44.96$

Table 3. Table A.1: Filtering algorithm on the interval { 0 , 1 , … , N } 0 1 … 𝑁 \{0,1,\ldots,N\} . The algorithm takes as input the parameters θ 𝜃 \theta and the data { 𝐲 i ​ j , 0 : N \{\mathbf{y}_{ij,0:N} , 𝐮 i ​ j , 0 : N } \mathbf{u}_{ij,0:N}\} whereas returns as output the filtered states 𝐳 i , 0 : N subscript 𝐳 : 𝑖 0 𝑁 \mathbf{z}_{i,0:N} . Note that the notation n | n − 1 conditional 𝑛 𝑛 1 n|n-1 indicates the prediction of the current n 𝑛 n -th state given the previous n − 1 𝑛 1 n-1 whereas n | n conditional 𝑛 𝑛 n|n denotes the update of the predicted state.

Algorithm 1	Gaussian Approximation filter algorithm
$𝐧=\mathrm{𝟎}:$	$z_{i,0}=0$	initialization
	${\lambda }_{i,0}^2=1$
$𝐧>\mathrm{𝟎}:$	$z_{i,n|n-1}=z_{i,n-1|n-1}$	prediction
	${\lambda }_{i,n|n-1}^2={\lambda }_{i,n-1|n-1}^2+{\sigma }_i^2$
	$z_{i,n|n}=\text{solve}\left(\partial \log f\left(z_{i,n}|{𝐲}_{ij,0:n},\theta \right)/\partial z_{i,n|n}\right)$	update
	${\lambda }_{i,n|n}^2=-\text{inv}\left({\partial }^2\log f\left(z_{i,n}|{𝐲}_{ij,0:n},\theta \right)/\partial z_{i,n|n}^2\right)$

Table 4. Table A.2: Backward smoothing algorithm over the interval { N , N − 1 , … , 0 } 𝑁 𝑁 1 … 0 \{N,N-1,\ldots,0\} . The algorithm takes as input the filtering solutions whereas returns as output the smoothed states 𝐳 i , 0 : N subscript 𝐳 : 𝑖 0 𝑁 \mathbf{z}_{i,0:N} conditional on the whole set of data 𝐲 i ​ j , 0 : N subscript 𝐲 : 𝑖 𝑗 0 𝑁 \mathbf{y}_{ij,0:N} . Note that the the notation . | n .|n refers to the filtering solutions whereas . | N .|N indicates the smoothing results.

Algorithm 2	Fixed-interval smoothing algorithm
$𝐧=𝐍:$	$z_{i,N|N}\sim N(z_{i,N|n},{\lambda }_{i,N|n}^2)$	initialization
	${\lambda }_{i,N|N}^2=1$
	$z_{i,n|N}=z_{i,n|n}+\left(\frac{{\lambda }_{i,n|n}^2}{{\lambda }_{i,n+1|n}^2}\right)\left(z_{i,n+1|N}-z_{i,n+1|n}\right)$	backward update
	${\lambda }_{i,n|N}^2={\lambda }_{i,n|n}^2+{\left(\frac{{\lambda }_{i,n|n}^2}{{\lambda }_{i,n+1|n}^2}\right)}^2\left({\lambda }_{i,n+1|N}^2-{\lambda }_{i,n+1|n}^2\right)$

Table 5. Table A.3: Marginal Metropolis-Hastings to estimate 𝜷 J × 1 subscript 𝜷 𝐽 1 \bm{\beta}_{J\times 1} . Note that U ​ ( 0 , 1 ) U 0 1 \mathcall{U}(0,1) indicates the Uniform distribution over the interval [ 0 , 1 ] 0 1 [0,1] whereas the adaptive phase can be performed at each iteration t 𝑡 t of the chain or rather after a fixed interval t + H 𝑡 𝐻 t+H (with H > 1 𝐻 1 H>1 ).

Algorithm 3	Metropolis-Hasting algorithm
$𝐭=\mathrm{𝟎}:$	Set  ${𝜽}^{\left(t\right)}={𝜽}^{\left(0\right)}$, ${𝚺}^{\left(t\right)}={𝚺}^{\left(0\right)}$	initialization
	Run Algorithms 1-2 to get  ${𝐙}^{\left(0\right)}\sim f\left(𝐙|𝐘,{𝜽}^{\left(\mathrm{𝟎}\right)}\right)$
	Compute  $f\left(𝐘|{𝜽}^{\left(\mathrm{𝟎}\right)}\right)$
$𝐭>\mathrm{𝟎}:$	${𝜽}^{\star }\sim N({𝜽}^{\left(t-1\right)},{𝚺}^{\left(t\right)})$	m-h loop
	Run Algorithms 1-2 to get  ${𝐙}^{\star }\sim f\left(𝐙|𝐘,{𝜽}^{\star }\right)$
	Compute  $f\left(𝐘|{𝜽}^{\star }\right)$
	Compute  ${\alpha }^{\left(t\right)}$ from Eq. (A.16)
	Get  $r\sim U(0,1)$	accept/reject
	Set  ${𝜽}^{\left(t\right)}={𝜽}^{\star }$, ${𝐙}^{\left(t\right)}={𝐙}^{\star }$   if ${\alpha }^{\left(t\right)}\le r$
	Set  ${𝜽}^{\left(t\right)}={𝜽}^{\left(t-1\right)}$, ${𝐙}^{\left(t\right)}={𝐙}^{\left(t-1\right)}$   if ${\alpha }^{\left(t\right)}>r$
	Run  ${𝚺}^{\left(t+1\right)}\leftarrow \text{adapt}\left({𝚺}^{\left(t\right)}\right)$   see [30]	adapting phase

Equations

$$\mathbf{s}_{ij}=\big{(}(\tilde{x}_{0},\tilde{y}_{0}),\ldots,(\tilde{x}_{n},\tilde{y}_{n}),\ldots,(\tilde{x}_{N_{ij}},\tilde{y}_{N_{ij}})\big{)}$$

$$\underset{{\mathbf{y}}_{ij}}{\underbrace{(y_0,\dots ,y_n,\dots ,y_N)}}$$

$$f(z_{i,n}|z_{i,n-1},\theta)=\big{(}\sigma^{2}_{i}\sqrt{2\pi}\big{)}^{-1}\cdot\exp\bigl{(}-(z_{i,n}-z_{i,n-1})^{2}/2\sigma_{i}^{2}\bigl{)}$$

$$f(y_{ijn}|{\pi }_{ijn},\theta )=f(y_{ijn}|{\mu }_1,{\kappa }_1){\pi }_{ijn}+f(y_{ijn}|{\mu }_2,{\kappa }_2)(1-{\pi }_{ijn})$$

$$f(y_{ijn}|\mu,\kappa)=\frac{1}{2\pi I_{0}(\kappa)}\exp\biggl{(}\cos(y_{ijn}-\mu)^{\kappa}\biggr{)}$$

$$\pi_{ijn}\leavevmode\nobreak\ \triangleq\leavevmode\nobreak\ \big{(}1+\exp(-\beta_{j}-z_{i,n})\big{)}^{-1}$$

$$\beta_{j}\leavevmode\nobreak\ \triangleq\leavevmode\nobreak\ \sum_{k=1}^{K}d_{jk}\gamma_{k}+x_{j}\bigg{(}\eta+\sum_{k=1}^{K}d_{jk}\delta_{k}\bigg{)}$$

$$\displaystyle\bm{\Theta}=\big{(}(\beta_{1},\ldots,\beta_{j},\ldots,\beta_{J}),(\kappa_{1},\kappa_{2})\big{)}$$

$$\displaystyle\mathbf{Z}=\big{(}(z_{1,0},\ldots,z_{1,N}),\ldots,(z_{i,0},\ldots,z_{i,N}),\ldots,(z_{I,0},\ldots,z_{I,N})\big{)}$$

$$\log f(\mathbf{Z},\mathbf{Y}|\mathbf{\Theta })$$

$$=\sum _{i=1}^I\sum _{j=1}^J\log f({\mathbf{z}}_{i,0:N}|\mathbf{\Theta })+\log f({\mathbf{y}}_{ij,0:N}|{\mathbf{z}}_{i,0:N},\mathbf{\Theta })$$

$$\displaystyle=\sum_{i=1}^{I}\sum_{j=1}^{J}\biggl{(}\log f(z_{i,0}|\theta_{Z_{0}})+\sum_{n=1}^{N}\log f(z_{i,n}|z_{i,n-1},\theta_{Z})+$$

$$\displaystyle+\sum_{n=1}^{N}\log f(y_{ijn}|z_{i,n},\theta_{Y})\biggr{)}$$

$$\log f(\mathbf{Z},\mathbf{\Theta }|\mathbf{Y})$$

$$=\sum _{i=1}^I\sum _{j=1}^J\log f(\mathbf{\Theta }|{\mathbf{y}}_{ij,0:N})+$$

$$+\sum _{i=1}^I\sum _{j=1}^J\log f({\mathbf{z}}_{i,0:N}|{\mathbf{y}}_{ij,0:N})+\log f(\mathbf{\Theta })$$

$${\beta }_j$$

$$\log f(z_{i,n}|{\mathbf{y}}_{ij,0:n},\theta )$$

$$+\log {\int }_{\mathbb{R}}f(z_{i,n}|z_{i,n-1},\theta )f(z_{i,n-1}|{\mathbf{y}}_{ij,0:n-1},\theta )d{z}_{i,n-1}$$

$$\log f(z_{i,n}|{\mathbf{y}}_{ij,0:n},\theta )$$

$$+\log {\int }_{\mathbb{R}}f(z_{i,n}|z_{i,n-1},\theta )f(z_{i,n-1}|{\mathbf{y}}_{ij,0:n-1},\theta )d{z}_{i,n-1}$$

$$y_{ijn}|z_{i,n}\text{leavevmode}\sim \text{leavevmode}mixVM({\mu }_1,{\mu }_2,{\kappa }_1,{\kappa }_2,{\pi }_{ijn})$$

$$z_{i,n}|z_{i,n-1}\text{leavevmode}\sim \text{leavevmode}\text{mathcall}N(z_{i,n-1},{\sigma }_i^2)$$

$$z_{i,n-1}|{\mathbf{y}}_{ij,0:n-1}\sim \text{leavevmode}\text{mathcall}N(z_{i,n-1|n-1},{\lambda }_{i,n-1|n-1}^2)$$

$$z_{i,n}|{\mathbf{y}}_{ij,0:n-1},\theta \sim \text{mathcall}N(z_{i,n-1|n-1},{\lambda }_{i,n-1|n-1}^2+{\sigma }_i^2)$$

$$f(y_{ijn}|{\pi }_{ijn},\theta )$$

$$\displaystyle\cdot\biggl{[}\frac{1-\pi_{ijn}}{2\pi I_{0}(\kappa_{2})}\exp\biggl{(}\cos(y_{ijn}-\mu_{2})^{\kappa_{2}}\biggr{)}\biggr{]}^{1-u_{ijn}}$$

$$\begin{split}\log f(z_{i,n}|\mathbf{y}_{ij,0:n},\theta)&\propto\log f(y_{ijn}|z_{i,n},\theta)+\log f(z_{i,n}|\mathbf{y}_{ij,0:n-1},\theta)\\ &=u_{ijn}\bigl{(}\kappa_{1}\cos(y_{ijn}-\mu_{1})+\log\pi_{ijn}-\log I_{0}(\kappa_{1})\bigr{)}+\\ &+(1-u_{ijn})\bigl{(}\kappa_{2}\cos(y_{ijn}-\mu_{2})+\log(1-\pi_{ijn})-\log I_{0}(\kappa_{2})\bigr{)}+\\ &+\frac{1}{2}\Biggl{(}\log(\lambda^{2}_{i,n-1|n-1}+\sigma^{2}_{i})-\frac{(z_{i,n}-z_{i,n|n-1})^{2}}{\lambda^{2}_{i,n-1|n-1}+\sigma^{2}_{i}}\Biggr{)}\end{split}$$

$$\frac{\partial \log \text{mathcall}F}{\partial z_{i,n|n}}=$$

$$\displaystyle-\frac{\exp\big{(}\beta_{j}-z_{i,n|n}\big{)}(u_{ijn}-1)}{\exp\big{(}2z_{i,n|n}-2\beta_{j}\big{)}-1}\Biggl{)}$$

$$\displaystyle\frac{\partial^{2}\log\mathcall{F}}{\partial z_{i,n|n}^{2}}=\leavevmode\nobreak\ \sum_{j=1}^{J}\Biggl{(}-\frac{1}{4\cosh\bigl{(}\frac{1}{2}(\beta_{j}-z_{i,n|n}\bigr{)}^{2}}-\frac{1}{(\sigma^{2}_{i}+\lambda^{2}_{i,n-1|n-1})}\Biggr{)}$$

$$\displaystyle\hat{\kappa}_{1}=I^{-1}\biggl{(}\frac{\sum_{i=1}^{I}\sum_{j=1}^{J}\sum_{n=0}^{N}u_{ijn}\cos\bigl{(}y_{ijn}-\mu_{1}\bigr{)}}{\sum_{i=1}^{I}\sum_{j=1}^{J}\sum_{n=0}^{N}u_{ijn}}\biggl{)}$$

$$\displaystyle\hat{\kappa}_{2}=I^{-1}\biggl{(}\frac{\sum_{i=1}^{I}\sum_{j=1}^{J}\sum_{n=0}^{N}(1-u_{ijn})\cos\bigl{(}y_{ijn}-\mu_{2}\bigr{)}}{\sum_{i=1}^{I}\sum_{j=1}^{J}\sum_{n=0}^{N}(1-u_{ijn})}\biggl{)}$$

$$q(\{{\mathbf{\theta }}^{(t)},{\mathbf{Z}}^{(t)}\}|\{{\mathbf{\theta }}^{(t-1)},{\mathbf{Z}}^{(t-1)}\})\text{leavevmode}=\text{leavevmode}q({\mathbf{\theta }}^{(t)}|{\mathbf{\theta }}^{(t-1)})f({\mathbf{Z}}^{(t)}|\mathbf{Y})$$

$$q({\mathbf{\theta }}^{(t)}|{\mathbf{\theta }}^{(t-1)})\text{leavevmode}=\text{leavevmode}\text{mathcall}N({\mathbf{\theta }}^{(t-1)},{\mathbf{\Sigma }}^{(t)})$$

Full text

Keywords: mouse tracking, state space modeling, dynamic systems, categorization task, aimed movements, bayesian filtering

1 Introduction

Over the last decades, the study of computer-mouse trajectories has brought to light new perspectives into the investigation of a wide range of cognitive processes (e.g., for a recent review see [21]). Unlike traditional behavioral measures such as reaction times and accuracies, mouse trajectories may offer a valid and cost-effective way to measure the real-time evolution of ongoing cognitive processes during experimental tasks [25]. This has also been supported by recent researches investigating mouse-tracking in association to more consolidated experimental devices, such as eye-tracking and fMRI [55, 68]. In a typical mouse-tracking experiment, participants are presented with a computer-based interface showing the stimulus at the bottom of the screen and two competing categories on the left and right top corners. Participants are asked to select the most appropriate label given the task instruction and stimulus while the x-y trajectories are instantaneously recorded. The main idea is that trajectories of reaching movements can unfold the decision process underlying the hand movement behavior. For instance, the curvature of computer-mouse trajectories might reveal competing processes activated in discriminating the two categories. Mouse-tracking has been successfully applied in several cognitive research studies, including lexical decision [41, 37], social categorization [24, 11], numerical cognition [19, 20], memory [54], moral decision [43], and lie detection [50]. Moreover, the availability of specialized and freely-available software for mouse-tracking experiments have strongly contributed to the wide-spread application of such a methodology in the more general psychological domain [22, 42]. Recently, the debate on the nature of cognitive processes tracked by this type of reaching trajectories have also received attention from the motor control literature [71, 25].

So far, mouse-tracking data have been analysed using simple strategies based on the conversion of x-y trajectories into summary measures, such as maximum deviation, area under the curve, response time, initiation time [35]. Although these steps are still meaningful in case of simple and well-behaved x-y trajectories, they can also provide biased results if applied to more complex and possibly noisy data. To circumvent these problems, other approaches have been proposed more recently [16, 10, 75, 44]. However, also the more recent proposals require modeling empirical trajectories before the data-analysis. Although these approaches potentially provide informative results in many empirical cases, they can also suffer from a number of issues, which revolve around the reduction of x-y data to simple scalar measures. For instance, problems may arise in the case of trajectories showing multiple phases, averaging with non-homogeneous curves, and signal-noise discrimination [10]. As far as we know, a proper framework to simultaneously model and analyse mouse-tracking data in a unified way is still lacking.

In this paper we describe an alternative perspective based on a state-space approach with the aim to simultaneously model and analyse mouse-tracking data. State-space models are very general time-series methods that allow estimating unobserved dynamics which gradually evolve over discrete time. As for diffusion models, which are widely used in modeling the temporal evolution of cognitive decision processes [64], they belong to the general family of stochastic processes and offer optimal discrete approximation to many continuous differential systems used to represent dynamics with autoregressive patterns [15]. In particular, we used a non-linear and discrete-time model that represents mouse trajectories as a function of some typical experimental manipulations. The model is estimated under a Bayesian framework, using a conjunction of a non-linear recursive filter and a Metropolis-Hastings algorithm. Data analyses is then performed using posterior distributions of model parameters [26].

The reminder of this article is organized as follows. In Section 2 we motivate our proposal by reviewing the main issues of a typical mouse-tracking experiment. In Section 3 we present our proposal and describe its main characteristics. In Section 4 we describe the application of our method to a psycholinguistic case study. Section 5 provides a general discussion of the results, comments and suggestions for further investigations. Section 6 concludes the article by summarizing its main findings.

2 A motivating example

To begin with, consider a two-choice semantic categorization task [17], in which participants have to classify semantic stimuli (e.g., name of animals) into their corresponding categories (e.g., mammal, fish). In the most typical implementation of a mouse-tracking task, participants would sit in front of a computer screen showing a resting frame (see Figure 1-a). They start a trial by clicking a starting button at the bottom-center of the screen, after which they are presented with a given stimulus (e.g. hen, see Figure 1-b). To finalize the trial, participants move the cursor on the screen by means of a well tuned computer-mouse in order to reach and select one of the two labels presented on the top-left and top-right corners of the screen (e.g., mammal vs. bird, see Figure 1-c). In the meanwhile, x-y coordinates, initiation time, and final clicking time are registered for each participant and trial. The basic idea is that x-y trajectories reflect the extent to which the real-time categorization response is affected by the experimental manipulation. More precisely, as a result of the assumption that co-activation of competing processes continuously drive the explicit hand response [66], one would suppose to see more curved - or generally irregular - trajectories in association with stimuli showing higher ambiguity. In our case, for instance, it would be expected that atypical exemplars, such as hen, dolphin, and penguin, globally produce more curved or irregular trajectories than typical exemplars like dog, rabbit, and lion (see Figure 1-d/e).

In the mouse-tracking literature, data analysis commonly proceeds summarizing the recorded trajectories by means of few indices, which are then used as input to standard statistical techniques. In the current example, for instance, the typicality manipulation could be tested by assessing whether the distribution of maximum deviations (i.e., the maximum curvature showed by trajectories) over trials and participants is bimodal or not [23]. In a similar way, linear models could be employed to test whether the typicality effect varies as a function of external covariates, such as psycholinguistic variables.

However, the two-step approach does have some issues. For instance, it lacks a way to represent both the experimental variability - that is induced by task manipulations - and individual variability - that is instead produced by individual-specific motor programs. Likewise, in some cases, it might neglect relevant characteristics of x-y data, with the consequence that similar classes of trajectories are treated as if they were different. Still, a two-step approach does ignore the data generation process underlying observed trajectories. This does not allow, for example, making predictions or simulations on new data given the experimental settings.

In the next section, we will present a dynamic probabilistic model that handles mouse-tracking data in a unified way. Our proposal is based upon a Bayesian non-linear state space approach, which offers a good compromise between model flexibility and model simplicity while overcoming many drawbacks of the standard mouse-tracking analyses.

However, the two-step approach does have some issues. For instance, it lacks a way to represent both the experimental variability - that is induced by task manipulations - and individual variability - that is instead produced by individual-specific motor programs. Likewise, in some cases, it might neglect relevant characteristics of x-y data, with the consequence that similar classes of trajectories are treated as if they were different. Still, a two-step approach does ignore the data generation process underlying observed trajectories. This does not allow, for example, making predictions or simulations on new data given the experimental settings.

In the next section, we will present a dynamic probabilistic model that handles mouse-tracking data in a unified way. Our proposal is based upon a Bayesian non-linear state space approach, which offers a good compromise between model flexibility and model simplicity while overcoming many drawbacks of the standard mouse-tracking analyses.

3 State-space modeling of mouse-tracking data

A state-space model is a mathematical description used to represent linear or generally non-linear dynamic models. In their general form, state-space systems consist of (i) a measurement density $f_{y}(\mathbf{y}_{n};\mathbf{z}_{n},\bm{\theta}_{y})$ that describes how the observed vector of data $\mathbf{y}_{n}$ at time step $n$ is linked to a possibly underlying process $\mathbf{z}_{n}$ and (ii) a state density $f_{z}(\mathbf{z}_{n};\bm{\theta}_{z})$ describing the transition dynamics that drive the vector of states $\mathbf{z}_{n}$. Temporal dynamics can be discrete or continuous and, in the latter case, stochastic differential equations are used to model the transition dynamics. By and large, there are two aims of any analysis involving state-space models. The first is to infer the unobserved process $\widetilde{\mathbf{Z}}=(\mathbf{z}_{0},\ldots,\mathbf{z}_{N})$ given the data $\mathbf{Y}=(\mathbf{y}_{0},\ldots,\mathbf{y}_{N})$. This task is usually accomplished by means of filtering and smoothing procedures [38]. The second aim regards estimating the parameters $(\bm{\theta}_{y},\bm{\theta}_{z})$ given the complete set of data $(\widetilde{\mathbf{Z}},\mathbf{Y})$. This is commonly performed using gradient-based methods on the likelihood of the model [61]. Although state space models were originally used in the area of aerospace modeling [39], they are now applied in a wide variety of domains, including control theory, remote sensing, economics, and statistics [32, 62]. Recently, there has also been an increasing interest in psychology, where state-space models have been used to analyse, for example, dyadic interactions [65], affective dynamics [46, 7], facial electromyography data [72], individual differences [31, 12], and path analysis [29].

In line with this, we developed a state-space representation to simultaneously model and analyse mouse-tracking data. In particular, our proposal is to represent the empirical collection of computer-mouse trajectories as a function of two independent sub-models, one representing the experimental manipulations (stimuli equation) and the other capturing the main features of the mouse movement process (states equation). Thus, the goal of our analysis is twofold: (i) to determine the states equation for each participant over a set of experimental trials, (ii) to estimate the parameters governing the stimuli equation. The first goal will provide information on how participants differ from each other in terms of movement dynamics. By contrast, the second goal will find out to what extent the experimental manipulations affect the individual variations in producing mouse-tracking responses.

3.1 Data

Let $\mathbf{S}$ be a $I$ (individuals) $\times$ $J$ (trials) array representing the observed data. The element $\mathbf{s}_{ij}$ of $\mathbf{S}$ defines the sub-array containing the streaming of Cartesian coordinates of the computer mouse movements:

[TABLE]

with [math] and $N_{ij}$ being the first and the last coordinates for the $i$-th participant in the $j$-th trial. The coordinates in $\mathbf{s}_{ij}$ are temporally ordered $(0<\ldots<n<\ldots<N_{ij})$ because they are usually collected while the computer-mouse is moving along its surface with a constant sampling rate. Further, to make the observed data comparable, we rescale and normalize $\mathbf{s}_{ij}$ as a function of a common ordered scale, which indicates the cumulative amount of progressive time from $0\%$ to $N=100\%$ [69, 56]. Thus, the final trajectories $\mathbf{s}_{ij}$ lie on the real plane defined by the hyper-rectangles $[-1,1]^{N}\times[0,1]^{N}$, with the first movement being equal to $(\tilde{x}_{0i},\tilde{y}_{0i})=(0,0)$ by convention. Since we are interested in studying the co-activation of competing processes as reflected in some spatial properties of the response - such as location, direction, and amplitude of the action dynamics [66, 21] - we need to simplify the original data structure so that these properties can easily emerge. Inspired by some of the work on this problem [28, 40, 10], we reduce the dimensionality of the data by projecting $\mathbf{s}_{ij}$ in a proper lower-dimensional subspace of movement via the restricted four-quadrant inverse tangent mapping (atan2, see [9]) from the real coordinates to the interval $[0,\pi]^{N}$ as follows:

[TABLE]

where $y_{0}$ is the angle at the beginning of the process whereas $y_{N}$ is the angle at the end of the process. Figure 1-F shows a graphical example of the atan2 function for a hypothetical movement path. Finally, the array of angles $\mathbf{y}_{ij}$ is the input for our state-space model.

3.2 Model representation

The unobserved states equation of the model is a AR(1) Gaussian model $Z_{i,n}|Z_{i,n-1}$ with transition density equal to:

[TABLE]

which models how the movement process of the $i$-th subject changes from the step $n-1$ to the next step $n$. The stochastic dynamics for the $i$-th subject is constrained by the variance parameter $\sigma^{2}_{i}\in\mathbb{R}^{+}$ that represents the uncertainty about the future location $z_{i,n+1}$ given the current state $z_{i,n}$.

The measurement equation for the observations $\mathbf{y}_{ij}=(y_{0},\ldots,y_{n},\ldots,y_{N})$ is modeled by means of a two-component von-Mises mixture distribution with density equal to:

[TABLE]

where the generic density is the standard von-Mises law:

[TABLE]

In the density formula, the term I_{0}(.)\leavevmode\nobreak\ \triangleq\leavevmode\nobreak\ (2\pi)^{-1}\int_{0}^{2\pi}e^{\kappa\cos x}dx\leavevmode\nobreak\ is the exponentially scaled Bessel function of order zero [1]. The parameters of the mixture density are mapped to the experimental interface of the two-choice categorization task (see Figure 1-d/e). In particular, the means $\{\mu_{1},\mu_{2}\}\in[-3.14,3.14)^{2}$ are mapped to the label categories C1 and C2 whereas the concentrations $\{\kappa_{1},\kappa_{2}\}\in\mathbb{R}^{2}_{+}$ indicate how the observations are spread around the means. In this sense, since $\{\mu_{1},\mu_{2}\}$ are determined by the fixed and known positions of the labels C1 and C2 on the screen, they are not treated as parameters to be estimated. Finally, the terms ${\pi}_{ijn}$ and $(1-{\pi}_{ijn})$ are the probabilities to activate the first and second density of the von-Mises components and are expressed as function of the latent states $\mathbf{z}_{i,0:N}$ and some additional covariates. The model is Markovian, in the sense that the unobserved states $\{Z_{n};n>1\}$ form a Markov sequence and the measurements $\{Y_{n};n>1\}$ are conditionally independent given the unobserved states.

To further characterize our state-space representation, the probability $\pi_{ijn}$ is defined according to a logistic function:

[TABLE]

with $\beta_{j}\in\mathbb{R}$ being the intercept of the model. Equation (3) can be interpreted as the probability for the $i$-th subject at step $n$ to categorize the $j$-th stimuli as belonging to C1 ($\pi_{ijn}$ tends to $1$) or C2 ($\pi_{ijn}$ tends to [math]). In addition, when the categories C1 and C2 are expressed in terms of distractor and target [21], the sequences $\bm{\pi}_{ij,0:N}$ can be interpreted as the attraction probability that the distractor has exerted on the trajectory process $\mathbf{z}_{i,0:N}$.

The state-space representation is completed by linearly expanding the intercept term $\beta_{j}$ as follows:

[TABLE]

where $\{\gamma_{k},\eta,\delta_{k}\}\in\mathbb{R}^{3}$, $x_{j}$ is an element of the array $\mathbf{x}\in\mathbb{R}^{J}$, whereas $d_{jk}$ is an element of the (Boolean) partition matrix $\mathbf{D}_{J\times K}$, with $d_{jk}=1$ indicating whether the $j$-th stimulus belongs to the $k$-th level of the variable $D$. Note that the matrix $\mathbf{D}$ satisfies the property $\sum_{k=1}^{K}d_{jk}=1$, for all $j=1,\ldots,J$. In our model representation, Equation (4) is the stimuli equation and conveys information about the experiment. It consists of three main terms. (i) A categorical term $\sum_{k=1}^{K}d_{jk}\gamma_{k}$ describing how the stimuli $\mathcall{J}=\{1,\ldots,j,\ldots,J\}$ have been arranged into $K<J$ distinct levels of a categorical variable $D$. (ii) A continuous term $x_{j}\eta$ that expresses the stimuli as a function of a continuous variable $X$ weighted by the coefficient $\eta$. (iii) An interaction term $x_{j}(\eta+\sum_{k=1}^{K}d_{jk}\delta_{k})$ between the levels of $D$ and $X$, where $\delta_{k}\in\mathbb{R}$ and $delta_{1}=0$. This definition allows for modeling all the cases implied by an univariate experimental design with at most one covariate variable. Indeed, for $\eta=0$ and $\bm{\delta}_{K}=\mathbf{0}_{K}$ this formulation boils down to the simplest experimental case with a single categorical variable $D$. By contrast, for $\bm{\delta}_{K}=\mathbf{0}_{K}$ and $\bm{\gamma}_{K}=\mathbf{0}_{K}$ it reduces to the case where stimuli are simply paired with a continuous predictor $X$. Finally, when $\mathbf{D}_{J\times K}=\mathbf{I}_{J\times J}$, the stimuli equation reduces to the most simple case where we have as many parameters as trials. Figure 2 shows a graphical representation of state-space model whereas Figure 3 illustrates the inner-working of the model for the simplest design with a two-level experimental factor.

In our model representation, the observed movement angles $\mathbf{y}_{ij,0:N}$ are sampled from C1 (resp. C2) with probabilities equal to $\bm{\pi}_{ij,0:N}$ (resp. $\bm{\pi}^{c}_{ij,0:N}=1-\bm{\pi}_{ij,0:N}$), which in turn are expressed as a function of the AR(1) latent trajectory $\mathbf{z}_{i,0:N}$. Hence, an increase in the latent process $z_{i,n}>0$ will also increase the probability that $y_{ijn}$ is sampled from the hemispace C1 (e.g., $\pi_{ijn}>0.5$). By contrast, a decrease in the latent process $z_{i,n}<0$ will increase the chance to sample $y_{ijn}$ from C2 (e.g., $\pi_{ijn}<0.5$). As a result of Eq. (4) such an increasing (or decreasing) pattern can be modulated by the stimuli component $\bm{\beta}_{J}$. Moreover, as the coefficients $\bm{\beta}_{J}$ are decomposed as a function of $\eta$, $\bm{\gamma}_{K}$, and $\bm{\delta}_{K}$, we can also analyse the effect of $\bm{\beta}_{J}$ on $\bm{\pi}_{ij,0:N}$ in terms of the experimental manipulation $D$, the covariate $X$, or the interaction term $DX$. Figure 3 shows a conceptual representation of the modeling steps involved by our approach. Panel (A) shows an example of the random-walk used to represent the movement process (Eq. 3). Instead, panel (B) shows the logistic function used to form the stimuli equation (Eq. 3) for two typical cases of $\bm{\beta}_{J}$. Panel (C) represents the probability $\bm{\pi}_{ij,0:N}$ to activate the distractor C1 (upper panel) and the probability $\bm{\pi}^{c}_{ij,0:N}$ to activate the target C2 (lower panel) as a function of $\mathbf{z}_{i,0:N}$ and $\bm{\beta}_{J}$. Finally, panel (D) depicts two cases of observed radians that are associated to $\bm{\pi}_{ij,0:N}$ and $\bm{\pi}^{c}_{ij,0:N}$. In particular, the upper panel shows an example of data with a pronounced attraction toward C1, which is in turn reflected by the blue probability curve of the panels (B)-(C). By contrast, the lower panel represents data with little attraction toward C1, as also reflected by the red probability curve of the panels (B)-(C). In this sense, as Equation (3) represents an intercept model, the parameter $\bm{\beta}_{J}$ does not affect the shape of the movement dynamics $\mathbf{z}_{i,0:N}$. On the contrary, it acts by shifting the movement dynamics upward ($\bm{\beta}<0$) or downward ($\bm{\beta}>0$) toward the C1 or C2 hemispaces, respectively.

3.3 Model identification

State-space model identification consists of inferring the unobserved sequence of states by means of filtering and smoothing algorithms and estimating the model’s parameters via Likelihood-based approximations [62, 58]. For instance, in the simplest linear gaussian case, where both the states and measurement equations are linear with additive gaussian noise, inference of latent states is usually performed via Kalman filter whereas parameter’s estimation is realised with the Expectation-Maximization algorithm. In our case, as Eqs.(3)-(4) describe a more complex non-linear model, we adopted a recursive Gaussian approximation filter for the inference problem [63], coupled with a marginal Metropolis-Hastings MCMC for the parameters estimation [2].

To formulate the problem more precisely, let:

[TABLE]

be the arrays representing all the $J\times 2$ unknown parameters and $I\times N$ unobserved states that characterize the model’s behavior. In this context, $\bm{\sigma}^{2}_{I}$ can be set to $\bm{1}_{I}$ without loss of model adequacy.222Indeed, the constraint $\bm{\sigma}^{2}_{I}=\bm{1}_{I}$ still guarantees the mapping $\pi_{ijn}:\mathbb{R}\to[0,1]$ to cover the needed time-to-time variability of the random walk, as Eq. (3) acts as a shrinkage operator on the support of the r.vs $\{Z_{i,0},\ldots,Z_{i,n}\}$. This has also been confirmed by several pilot simulations we ran on our model. Note that this assumption is not overly limiting, since our state-space representation is built under the smoothness assumption on the movement behavior of the hand, according to which large abrupt changes in the small interval $[n,n+1]$ are not allowed [74]. The joint log-density of the complete-data given the array of parameters and the observed data is defined as follows:

[TABLE]

where the state and measurement equations are given as in (1)-(2) whereas the term $f(z_{ij0}|\theta_{Z_{0}})$ is the a-priori density function for the initial state of the process. Note that the factorization (9) is due to the Markovian properties of the model. By adopting the Bayesian perspective, we perform inference conditional on the observed sample of angles $\mathbf{Y}$, with $\bm{\Theta}$ being an unknown term. The posterior density $f(\mathbf{Z},\bm{\Theta}|\mathbf{Y})$ of hidden states and parameters is as follows:

[TABLE]

where $f(\bm{\Theta})$ is a prior density ascribed on the vector of model’s parameters $\bm{\Theta}$. Note that Equation (10) comes from the standard conditional definition f(\mathbf{Z},\bm{\Theta}|\mathbf{Y})=f(\mathbf{Z},\bm{\Theta},\mathbf{Y})\big{/}f(\mathbf{Y}), where the joint density $f(\mathbf{Z},\bm{\Theta},\mathbf{Y})$ is re-arranged by factorization using the Markovian properties of the model [2]. Since our aim is to get samples from the posterior $f(\mathbf{Z},\bm{\Theta}|\mathbf{Y})$, we proceed by jointly updating $\bm{\Theta}$ and $\mathbf{Z}$ using a marginal Metropolis-Hastings. This alternates between proposing a candidate sample $\bm{\Theta}^{(t)}$ given $\bm{\Theta}^{(t-1)}$ and filtering the sequences $\mathbf{Z}^{(t)}$ conditioned on $\bm{\Theta}^{(t)}$. Finally, the candidate couple \bigl{(}\bm{\Theta}^{(t)},\mathbf{Z}^{(t)}\bigr{)} is jointly evaluated by the Metropolis-Hasting ratio.

The evaluation of both the densities $f(\mathbf{Z}|\mathbf{Y})$ and $f(\bm{\Theta}|\mathbf{Y})$ involve computing the expression in Eq. (11). To do so, we derived the first term by means of filtering and smoothing procedures [38] whereas the second term was evaluated by implementing a Metropolis-Hasting algorithm. All the technical steps for the model identification are included in Appendix A-C.

3.4 Model evaluation

The state-space model formulated can be evaluated in many ways under the Bayesian framework of analysis [26, 60]. For instance, adequacy of the algorithm can be assessed via standard diagnostic measures, such as traceplot of the chains, autocorrelation measures, and the Gelman-Rubin statistics whereas the recovery of the true model structure can be done by simulations from the priors ascribed to the model [26]. Similarly, the adequacy of the model to reproduce the observed data can be assessed by means of simulation-based procedures (e.g., posterior predictive checks) where the fitted model is used to generate new simulated datasets that are then compared to the observed data [27, 14]. In our context, the robustness of the model formulation in recovering the true model structure as well as the goodness of fit to the observed data are assessed by adopting a simulation-based approach. Technical details on this procedure are available in Supplementary Materials.

4 Application

In this section, we will present an application of the model to the analysis of an already published lexical decision dataset [6]. The state-space modeling framework will be evaluated via three different instances of model representation with an increasing level of complexity. Note that the application we report here has only an illustrative purpose with the main goal to introduce and highlight the interpretation of the model’s parameters and the flexibility of its representation with dynamic data. All the models were estimated using 20 (chains) $\times$ 10000 (iterations), with a burning-in period of 2500 iterations. Starting values $\bm{\theta}_{0}$ for the MH algorithm were determined by maximizing the observed likelihood of the model in Eq. (2). Similarly, the starting covariance matrix $\bm{\Sigma}^{(0)}$ was computed by using the Hessian of the observed likelihood at $\bm{\theta}_{0}$. The adaptive phase of the MH algorithm was performed at fixed interval $t+H$ (with $H=25$) to prevent the degeneracy of the adaptation. For each model, the prior densities were defined as $f(\bm{\theta})=\mathcall{N}(\bm{\mu}=\mathbf{0},\mathbf{1}\sigma^{2}=25)$, where the variance was sufficiently large to cover the natural range of the model parameters. The adequacy of the model to reproduce the data was evaluated with a simulation-based approach, where a series of $M=5000$ new datasets $(\mathbf{Y}^{*}_{1},\ldots,\mathbf{Y}^{*}_{M})$ were generated through the fitted model and compared with the observed data $\mathbf{Y}$ [14]. The goodness of fit was evaluated overall (i.e., the adequacy of the model to reproduce the complete observed matrix $\mathbf{Y}$) and subject-based (i.e., the adequacy of the model to reproduce for each subject $i=1,\ldots,I$ the observed matrix $\mathbf{Y}_{i}$). Comparisons were computed by means of 0-100% normalized measures, with 0% indicating bad fit and 100% optimal fit. Technical details as well as extended graphical results are included in Supplementary Materials.

4.1 General context and motivation

Lexical decision is one of the most known and widely used task to study visual world recognition and reading in the cognitive psycholinguistic literature [34, 52, 73]. Generally, this task is very simple and versatile and provides an ideal context for applying the state-space modeling framework when lexical decision data are collected via the mouse tracking paradigm. In this application, we evaluated the extent to which the parameters of the state-space model reflect eventual differences associated with the manipulation of a stimulus type factor composed by words (with either high-frequency or low-frequency) and random strings (i.e., random sequence of letters that are phonotactically illegal in the language) in the lexical decision task. Moreover, we will take advantage of this psycholinguistic case study to show how our state-space model can deal with both categorical and (pseudo)quantitative predictive variables considered either individually or in interaction in the model. In particular, the first model instance will illustrate the application of our modeling framework when a simple categorical variable (stimulus type factor) is considered to affect the observed mouse-tracking trajectories collected using the lexical decision task. By contrast, the second model will be based on a simple regression-type model with a single quantitative independent variable (bigram frequency) used to predict the attraction toward the distractor category. Finally, the third model will integrate these two variables (stimulus type factor and bigram frequency) into a unified model including the main effects of the two variables as well as their interaction. In our context, the first two models will be considered as simple toy examples to illustrate the main features of the state-space model representation when applied to real data, whereas the third model will be discussed in more details according to a group analysis evaluation as well as an individual analysis representation.

4.2 Model 1

Data structure and variables. In the original work by [6], the lexical decision experiment was run in Italian and based only on one stimulus type factor with four different levels: Words of high written frequency (HF, e.g., acqua “water"), words of low written frequency (LF, e.g., cervo “deer"), pseudowords (PW, e.g., “dorto"), and strings of letters that are orthographically illegal in Italian (NW, e.g., “btfpr"). In their study, participants saw a total of 96 stimuli, one at the time, and were required to categorize each stimulus as either a word or a nonword by using the mouse-tracking paradigm. Trajectories were recorded using the Mouse Tracker software [22] with sampling rate of approximately 70Hz [6]. As usual, raw trajectories were normalized into $N=101$ time steps using linear interpolation with equal spaces between coordinate samples. However, for our analysis we preferred to select only three of the four levels of the experimental factor (that is to say, HF,LF, and NW) for a total of 72 stimuli equally distributed within each level.333The motivation for this selection was due to some technical reasons regarding the lack of design balance in the original dataset, as the PW level showed a large number of errors when compared with the other three categories. In addition, the three-level representation of the stimulus type factor simplifies the interpretation of the results when we consider the full model with interaction. Finally, the dependent variable $\mathbf{Y}$ of Model 1 consisted of the movement angles array associated with the mouse-movement trajectory recorded for each distinct stimulus in the stimulus set.

Data analysis and results. In this first model the term $\beta_{j}$ in the stimuli equation boils down to the simple expression:

[TABLE]

where the indices $k=1,2,3$ refer to HF, LF, and NW stimuli. The MCMC convergences of the algorithm are reported in Supplementary Materials. The model fitted the data very satisfactorily, with overall fit of 84% and subject-based fit of 74% (see Table 1). The posterior quantiles (5%,50%, and 95%) are reported in Table 2 whereas figure 4-A shows the probability graph, that is to say, the probability to activate the distractor cue for each of the three levels HF, LF, NW as a function of the latent variable $Z$.

The results of this first analysis clearly show that the dynamics of the state-space model were unaffected by the different categories represented in the recoded experimental factor. This pattern finds further support in the post-hoc comparisons between the three experimental conditions (Figure 4-B). In sum, these findings indicate that for a dynamic model represented according to a state-space modeling framework, the three stimulus categories (HF, LF, and NW) were all processed in a very similar way, as the original trajectories were not sufficiently different among the three stimulus categories. In substantive terms, the results of the categorical model showed how the attraction probability toward the distractor was definitively modest in all the three experimental conditions. This is evident from a direct inspection of figure 4-B where the probability activation function (logistic function) is shifted towards right ($Z>0$) which in turn means that the average activation of the distractor category was relatively poor in HF, LF, and NW items. In this respect, the results of our simple spatial model were partially at odds with the outcomes observed using temporal measures (response time variables) [6].

4.3 Model 2

Data structure and variables. Also for the second model, the dependent variable was represented by the movement angles array $\mathbf{Y}$. However, unlike model 1, in model 2 the original independent categorical variable (stymulus type factor) was replaced with a quantitative psycholinguistic variable called bigram frequency. Bigram frequency is defined as the frequency with which adjacent pairs of letters (bigrams) occur in printed texts; for its characteristics, it may be considered as a measure of orthographic typicality [33]. In this second application, only bigram frequency was used as quantitative variable, since it was the only psycholinguistic variable that could be computed for all the 72 stimuli in the stimulus set. This second model instance nicely provides a simple but effective example of application of our state-space model when a continuous variable is considered to predict the attraction toward distractor.

Data analysis and results. In model 2 the term $\beta_{j}$ simply reduces to:

[TABLE]

as the first and third terms in formula (4) cancel out. In this case, the variable $x_{j}$ denotes the value of the bigram frequency for stimulus $j$ in the stimulus set. For the model results, the posterior quantiles are reported in Table 2 whereas MCMC convergences of the algorithm are reported in Supplementary Materials. Also in this case, the model fitted the data very well, with overall fit of 73% and subject-based fit of 70% (see Table 1). Figure 5 shows the probability graph for model 2. This graph represents the probability to activate the distractor hemispace at three representative levels of the variable, i.e. the lowest, the medium, and the highest values of bigram. As evident from the graph, bigram frequency affects the probability to activate the target, with a higher probability for stimuli with low bigram frequency.

In substantive terms, the results of the quantitative model supported the evidence that the attraction probability toward the distractor was slightly affected by the specif value of the quantitative predictor (bigram frequency). In particular, low-level bigram frequencies were characterized by an average larger activation probability ($0.55$) for the distractor, whereas medium or large frequencies were associated with a logistic function slightly shifted toward positive values of the latent space $Z>0$, thus reflecting a lower chance for the distractor category (average activation probability of $0.45$). Moreover, by an inspection of the contingency table for the joint representation of bigram frequency (as a transformed categorical variable) and stimulus type, we noted that low bigram frequency values were mainly characterized by string letters (NW: $94\%$) whereas high bigram frequency values were predominantly associated with high frequency words (HF: $55\%$) or low frequency words (LF: $44\%$).

4.4 Model 3

Data structure and variables. The final and more complex model included both the three-level categorical predictor (stimulus type factor: HF, LF, STR) and the continuous predictor (bigram frequency) as well as the interaction term between these two variables. The dependent variable was the movement angles array $\mathbf{Y}$.

Data analysis and results. The stimuli equation which characterizes the third model is defined as follows:

[TABLE]

The MCMC diagnostics together with the estimated marginal posterior densities for the model’s parameters are reported in Supplementary Materials. The model fit was good, with an overall fit of 75% whereas the subject-based fit was equal to 71% (see Table 1). The posterior quantiles are reported in Table 2. Figure 6 shows the probability graph for model 3. This graph represents the probability to activate the distractor hemispace for each of the three levels HF, LF, NW of the categorical factor as a function of the latent variable $Z$ and three distinct levels (high, medium, and low) for bigram frequency. The inspection of Figure 6 shows a clear interaction between stimulus type factor and bigram frequency indicating that the impact of stimulus category, in particular word frequency, increases with the decrease of stimulus bigram frequency. In other words, at high level of bigram frequency, the probability to activate the distractor is similarly low in all conditions (.17 $\leq$ p-distractor $\leq$ .2). By contrast, when bigram frequency decreases - that is stimuli become orthographically atypical - the probability of distractor activation increases, but only for the more lexically-familiar stimuli, i.e., words of high frequency (e.g., p-distractor raises from 0.17 to 0.70, in low and high bigram frequency condition respectively).

Finally, it is also worth mentioning the emergence of the main effect of stimulus category which was instead missing in model 1. By a quick inspection of Figure 7, one may clearly observe that HF words differ from both LF words and letter strings (NW), whereas LF words and letter strings do not differ with respect to the probability of activation of the distractor hemispace. Interestingly, the addition of the covariate bigram frequency in the model allowed the main effect of stimulus category to show up. Indeed, while at the medium and high levels of bigram frequency the results are in line with those observed at a sample level in the original study (see Figures 1,2,5 in [6]) and in a recent re-analysis (see Table 2 in [10]), in the case of low bigram the probability to activate the distractor increases with respect to high frequency words (HF). This might be somewhat related to a moderate difficulty in the orthographic processing of low frequency bigram words [57] even in the case of stimuli with richer lexical representation.

4.5 Profiles analysis

To further investigate the dynamic characteristics involved in the lexical decision task, we extend here the results of the third model to include also a profiles analysis. Figure 8 shows the estimated latent movement states $\mathbf{Z}_{I\times N}$ for all the participants involved in the study. The profiles appear regular, as they evolve smoothly toward the target cue (T). We grouped the dynamics into four well-separated clusters (Figure 8, smallest panels on the right) according to their functional similarities [56]. Particularly, the first group is characterized by a higher exploration of the distractor’s hemispace, especially in the first 30% of the process. The same applies to the third and fourth groups, although they show a gradual activation of the distractor. Finally, the second group clearly represents those profiles with no uncertainty in the categorization process, as they show no activation of the distractor’s hemispace at all. Although well-separated among them, these clusters still show some level of inner heterogeneity (for example, see group 1 and 4). To study this latter issue in terms of experimental manipulations, we focused on group 1 and considered the low vs. high frequency conditions (HF vs. LF). We also selected the middle phase of the process ($\Delta=30\%-50\%$), where it is expected to observe larger cognitive competitions in the categorization [6]. Figure 9 shows the participants’ profiles in terms of attraction probability $\bm{\pi}_{4\times N}$ for the two lexical conditions. As expected, the profiles differ between these conditions, with LF eliciting higher attraction probability. This is in line with the fact that low frequency words have a weaker lexical representation than high frequency stimuli and consequently they are more difficult to process [6]. Interestingly, the individual profiles also differ in the way they activate the distractor. For instance, the participant 6 had higher probability in both LF ($p_{\Delta}(D)=0.67$) and HF ($p_{\Delta}(D)=0.54$) conditions whereas the participant 7 had a more pronounced activation just in the LF condition ($p_{\Delta}(D)=0.57$) than HF ($p_{\Delta}(D)=0.43$). Similarly, participants 6 and 7 seemed to prolong the competing dynamics up to the 50% of the process, by contrast participants 8 and 15 seemed to resolve the lexical competition earlier as showed by the abrupt decreasing of their curves. We complete our analysis by evaluating how individual profiles are linked to empirical measurements. Figure 10 represents this scenario for two stimuli belonging to HF and LF conditions. As we can notice, the curves present the same dynamics (due to the latent states $\mathbf{z}_{i,0:N}$) although they clearly differ in terms of attraction exerted by the stimulus (due to the $\beta$ component of the model). In this case, the LF stimulus produced larger conflict than HF in the lexical categorization. This is evident when we turn back to the observed data: as expected, the rose diagrams of LF showed larger directions in the distractor’s hemispace.

5 Discussion

We have described a new approach to model and analyse dynamic data coming from mouse-tracking experiments. Our proposal took the advantages of a state-space representation, in which the observed data $\mathbf{Y}$ were thought as being function of two independent sub-models, one representing the movement process and its properties ($\mathbf{Z}$) and the second modeling the two-choice experimental task ($\bm{\beta}$) according to which the data were collected. These sub-models were integrated by means of an inverse-logit function ($\pi$) that expressed how the experimental manipulations acted on the movement processes in selecting the final correct response against the competing one. This formulation was flexible enough to take into account the complexity of some dynamic behaviors showed by the reaching trajectories. Moreover, it allowed for separately accounting for the motor heterogeneity and experimental variability in $\mathbf{Y}$. Indeed, when $\bm{\beta}=\mathbf{1}\beta_{0}$ our state-space representation simply reduced to a model where the experimental manipulations had no relevant effect in reproducing the observed data. This instance has been illustrated in Section 4.3 (Model 1). In this case, as $\mathbf{Z}=\mathbf{0}$ was not allowed in our model formulation, all the variability of $\mathbf{Y}$ can be ascribed to $\mathbf{Z}$. This is relevant in view of the fact that movement variability may reflect only individual motor executions in absence of any experimental manipulations [74]. The movement component $\mathbf{Z}$ was modeled to be Markovian with gaussian transition density.

Although more complex models can be used to represent movement dynamics, simple random walks still allows a great deal of flexibility in modeling reaching trajectories under weak assumptions on the movement behavior [74, 53]. In particular, in the case of mouse-tracking tasks, they allow representations of the following three properties: (i) Each movement is goal-oriented as individuals have to finalize the action by clicking on one of the two categories shown on the screen. (ii) Mouse-tracking trajectories generally start at rest, proceed out in the movement space, and end at rest. (iii) Hand trajectories tend to be smooth during the reaching process, i.e. small changes in the interval $[n,n+1]$ are more likely than large and abrupt changes [8, 67]. The stimuli component $\bm{\beta}$ was defined to be a linear combination of information typically involved in a univariate design, namely a categorical variable $D$ containing the levels of the experimental factor and a continuous covariate $X$. This gave researchers the opportunity to additionally analyse which component of the experimental design is relevant in producing the effect of $\bm{\beta}$ on $\mathbf{Y}$. The data-generation process was defined according to a mixture of two von-Mises distributions representing the categories of a two-choice categorization task. Among others, we chose this distribution because it provides a flexible representation for angular ordered data, especially because it simplifies mathematical computations involved in the model’s derivation [48, 51].

There are other existing methods that offer alternative ways to model mouse-tracking data. For instance, [71] proposed the use of the movement superposition model [36] to model and analyse the typical two-choice lexical decision task. In particular, they modeled mouse-tracking trajectories by representing the complete hand movement as a summation of sub-movements, which were obtained by the solution of the minimum-jerk equation for the standard reaching trajectory (i.e., a movement characterized by a bell-shaped speed profile that minimizes the sum of the squared rates of jerks over the movement duration). Similarly, [25] discussed how an intermittent model of arms movement can be used for reaching trajectories in random-dot experiments. They used both Wiener’s diffusion process and Flash and Hogan’s movement equation to predict reaction times (RTs) and movement data. Their goal was to assess the link between movement trajectories and underlying cognitive processing. Our model differs in some respects from these works. With regards to [71], for instance, we used a stochastic state-space approach to model the movement trajectories instead of deterministic equations. Instead, with respect to [25], we tailor-made our model to a typical two-choice categorization task, making use of few assumptions on the nature of the movement process (as those implied by the Gaussian AR(1) process). By and large, our goal was not to provide a mathematical representation of the cognitive components underpinning mouse-trajectories since the model does not describe the cognitive processing per se. By contrast, we simply provided a statistical model for the analysis of mouse-tracking data, which can offer a good compromise between data modeling and data analysis.

5.1 Model’s advantages and limitations

Our non-linear state-space model has several advantages. For instance, when comparing with the standard approaches, our proposal provides a unified analytic framework to simultaneously model and analyse trajectories data. By modeling movement heterogeneity and task variability together, we can evaluate how experimental variables directly act on the observed series of trajectories, with no need to use any kind of summary measures. An additional advantage of our model concerns the study of individual differences in terms of latent dynamics. While this is impractical in standard two-step approaches, in our proposal researchers can assess individual variations by studying the movement profiles $\mathbf{\tilde{Z}}$ once they are estimated. For instance, they can be analysed in terms of similarity/dissimilarity with regards to external individual covariates (e.g., vocabulary knowledge and bilingualism in psycholinguistic experiments; IQ, risk-taking propensity, or more generally clinical variables in decision-making tasks). Still, individual dynamics can be compared each other qualitatively in terms of chance to activate the distractor or target cues. As the dynamics are normalized on a common cumulative scale, researchers can assess whether the chance to activate the distractor cue at a certain percentage of the process and for an experimental manipulation, is particularly higher in a sub-group of participants (this case, for example, has been described in Section 4.6).

As for any modeling approach, also the current proposal can potentially suffer from some limitations. A first limitation concerns the only-intercept model $\pi(\mathbf{Z},\bm{\beta})$ used to integrate individual dynamics and experimental information. Although this was enough to represent whether or not certain stimuli can increase the probability to select the distractor cue, we may want to known whether some experimental manipulations can modify the individual dynamics as well. However, this would particularly pronounce the computational costs required for the model identification (especially with regards to filtering), as we need to appropriately generalize Eq. (3) to include more parameters. Lastly, in the current study we used univariate non-linear state-space models to represent individual dynamics for the sake of parsimony. However, more complicated situations may require models including further movement characteristics like step-length, velocity, acceleration, and jerk [45], which may be modeled as statistical constraints of the model [13, 10].

5.2 Further extensions

Our non-linear state-space model can be improved in many aspects. For instance, the stimuli equation (4) can be generalized to cope with more complex experimental designs, like those involving multiple factors and covariates together with high-order interaction terms. Likewise, the current model restrictions can be relaxed to allow changes in slopes of $\pi(\mathbf{Z},\bm{\beta})$ as a function of the experimental stimuli. Further, the development of a hierarchical representation of the model, with a random-effect component in the state equation (3), would offer a way to model the inter-individual variations as resulting from an underlying common population. Still, the development of a multivariate state-space model to include other movement components will surely be considered a future extension of the present work. Further studies may lead to generalize the AR(1) process used for the movement dynamics to include former knowledge on the deterministic constraints of the hand movement as those used, for instance, by [71] and [25]. Moreover, further studies may also lead to generalize the AR(1) process used for the movement dynamics to include former knowledge on the deterministic constraints of the hand movement as those used, for instance, by [71] and [25]. Finally, an open issue which deserves greater consideration in future investigations is the need for a formal comparative framework with which we may eventually contrast and compare spatial modeling perspectives (like the one presented in this contribution) and currently used methods for analyzing mouse tracking data based on descriptive statistics [21].

6 Conclusions

In this paper we presented a novel and comprehensive analytic framework for modeling and analyse mouse-tracking trajectories. In particular, a non-linear state-space approach was used to model the observed trajectories as a function of both individual movement dynamics and experimental variables. Model identification was performed under the umbrella of Bayesian methods, in which a Metropolis-Hastings algorithm was coupled with a recursive gaussian approximation filter to get posterior distributions of model parameters. For the sake of illustration, we applied our new approach to a real mouse-tracking dataset concerning a two-choice lexical categorization task. The results indicated how our proposal can provide valuable insights to assess the dynamics involved in the decision task and identify how the experimental variables significantly contributed to the observed movement heterogeneity. Moreover, the analysis of individual profiles allowed for comprehensive and reliable identification of individual and group-based differences in the dynamics of decision making.

In conclusion, we think that this work yielded interesting findings in the development of computational models able to capture the unfolding high-level cognitive processes as reflected by motor executions which are typically involved in mouse-tracking tasks. To our knowledge, this is the first time that mouse-tracking data are fully modeled and analysed within a process-oriented approach. We believe our contribution will offer a novel strategy that may help cognitive researchers to understand the roles of cognition and action in mouse-tracking based experiments.

Acknowledgments. The authors thank Dr. Simone Sulpizio for many helpful discussions on earlier versions of the manuscript and Dr. Laura Barca for providing the dataset used in the case study. This work was supported by CARITRO, Fondazione Cassa di Risparmio di Trento e Rovereto (ref. grant 2016.0189).

Appendix A: Filtering and smoothing

The term $f(\mathbf{Z}|\mathbf{Y})$ in Eq. (11) is recursively computed given all the measurements up to the $n$-th step. Let:

[TABLE]

be the filtering density at step $n$. The first term of the right-side of the equation is the observation equation whereas the second term represents the Chapman-Kolmogorov equation. Given the initialization at $n=0$, the filter proceeds by solving the integral in the Chapman-Kolmogorov equation (prediction step) and compute the log posterior filtering density (update step). In our model, we solve Eq. (A.1) by means of a Gaussian approximation filter [63], which computes a gaussian approximation to the posterior density $f(z_{i,n}|\mathbf{y}_{ij,0:n},\theta)$ and determines its posterior mode $z_{i,n|n}$ and variance $\lambda^{2}_{i,n|n}$ recursively.

More technically, let:

[TABLE]

be the filtering density at step $n$. Consider the following definitions:

[TABLE]

where $z_{i,n-1|n-1}$ and $\lambda^{2}_{i,n-1|n-1}$ represent the mode and the variance of the gaussian approximation in the prediction step. Under definitions (A.4)-(A.5), integrating out for $z_{i,n-1}$ in Eq.(A.2) yields to:

[TABLE]

For the sake of computational simplicity, we rewrite the measurement density in Eq. (2) as follows:

[TABLE]

where $u_{ijn}$ is a (deterministic) indicator variable taking the value 1 when $y_{ijn}$ is in the area of the screen associated to C1, and 0 when $y_{ijn}$ is in the area associated to C2 (see Fig. 1) [5]. Next, using these results together with Eqs. (A.3)-(A.5) we obtain:

[TABLE]

Differentiate $\mathcall{F}\triangleq f(z_{i,n}|\mathbf{y}_{ij,0:n},\theta)$ around $z_{i,n|n}$ gives:

[TABLE]

[TABLE]

where \cosh(x)\triangleq\bigl{(}1+\exp(-2x)\bigr{)}/\bigl{(}2\exp(-x)\bigr{)} is the hyperbolic cosine function.

Finally, the posterior moment $z_{i,n|n}$ is obtained by solving Eq. (A.9) whereas $\lambda^{2}_{i,n|n}$ is computed by the negative inverse of Eq. (A.10) [70]. As Eq. (A.9) is non-linear, it can be solved numerically (e.g., using the Broyden’s method). The complete filtering algorithm is summarized in Table A.1.

Finally, to ensure that the unobserved sequence $\mathbf{z}_{i,0:N}$ is an approximate realization from $f(\mathbf{z}_{i,0:N}|\mathbf{y}_{ij,0:N})$ we need to refine the filtering results conditional on the whole observed data $\mathbf{y}_{ij,0:N}$. This task is achieved by means of a standard fixed-interval smoothing algorithm [4, 49], which uses the posterior filtering moments $z_{i,n|n}$ and $\lambda^{2}_{i,n|n}$ as input. The smoothing algorithm is described in Table A.2.

Appendix B: Posteriors computation and estimation

In what follows, we describe the steps for computing the term $f(\bm{\Theta}|\mathbf{Y})$. First, we note that the array $\bm{\Theta}$ consists of two blocks of parameters associated to the observation equation and the stimuli equation of the model (3)-(4), namely $J$ scalars $\{\beta_{1},\ldots,\beta_{J}\}$ paired with the set of stimuli/trials and two parameters $\{\kappa_{1},\kappa_{2}\}$ for the vonMises concentrations. To simplify the computations in the Metropolis-Hastings algorithm, the terms $\{\kappa_{1},\kappa_{2}\}$ can be distinctively determined prior running the MH algorithm [47]. Moreover, since we are not interested in the posterior distributions of these parameters, as long as they are not involved in the state-space dynamics, we compute them using the following Maximum-Likelihood estimators [5]:

[TABLE]

where $\mu_{1}$ and $\mu_{2}$ are the location parameters fixed by the experimenter, $u_{ijn}$ is defined as in Eq. (Appendix A: Filtering and smoothing), whereas $I^{-1}$ is the inverse of the modified Bessel function which is evaluated numerically [1]. Given the above results and the constrains $\bm{\sigma}_{I\times 1}=\mathbf{1}_{I}$, the array of parameters simply reduces to $\bm{\theta}_{J\times 1}$. This simplifies the inner-working of the MH algorithm as it now works on a smaller and more compact parameter space. To proceed further, the decomposition (10) involves the following definition for the MH proposal density [2]:

[TABLE]

where $f(\mathbf{Z}^{(t)}|\mathbf{Y})$ is evaluated through filtering/smoothing. This is appealing since the posterior density $f(\mathbf{Z},\bm{\Theta}|\mathbf{Y})$ from which it might be difficult to sample from, reduces now to $f(\bm{\Theta}|\mathbf{Y})$ that is conveniently defined on a smaller parameters space [3]. In our case, we can set:

[TABLE]

with $\bm{\Sigma}^{(t)}$ being a suitable $J\times J$ covariance matrix. Consequently, the MH acceptance ratio is as follows:

[TABLE]

where $f(\mathbf{Y}|\bm{\theta})$ is the density for the marginal likelihood computation, $f(\bm{\theta})$ indicates the prior density over the parameters, whereas $q(\bm{\theta}^{(t)}|\bm{\theta}^{(t-1)})$ is the MH proposal density. Note that under (A.14), the terms $q(\bm{\theta}^{(t)}|\bm{\theta}^{(t-1)})$ and $q(\bm{\theta}^{(t-1)}|\bm{\theta}^{(t)})$ in Eq. (A.15) vanish as they refer to the same probability value. This yields to:

[TABLE]

where the MH ratio is now expressed as a function of the marginal likelihood and the priors. The term $f(\mathbf{Y}|\bm{\theta})$ can be easily computed as a byproduct of the filtering calculations (see Appendix C). Finally, the choice of a well-suited covariance matrix $\bm{\Sigma}^{(t)}$ for the proposal distribution is crucial in order to achieve chains’ convergences. In our context, we used the Haario’s adaptive solution where the proposal covariance is iteratively adapted during the chains using the current proposal covariance up to the adaptation step [30]. The complete MH algorithm is summarized in Table A.3.

Appendix C: Marginal Likelihood computation

The marginal likelihood $f(\mathbf{Y}|\bm{\theta})$ is computed as follows. Let the observed-data log likelihood be:

[TABLE]

According to the standard prediction error decomposition, this functional has been factorized as a function of both the measurement model density (2) and the one-step ahead predictive density (A.6) [18]. By substituting the above definitions with our model’s densities, we get the following marginal likelihood:

[TABLE]

where the von-Mises and Gaussian densities have been written by dropping the constant terms. Using the predictive moments $z_{i,n|n-1}$ and $\lambda^{2}_{i,n|n-1}$ from the filter algorithm (see Table A.1), a recursion on $\{0,1,\ldots,N\}$ can be written for the likelihood computation, which consists of solving $N$ integrals over the support of the r.v. $Z_{i,n}$. Since no analytical solutions are available for the functional above, numerical integration methods can be used to solve the integrals [59].