ssMousetrack: Analysing computerized tracking data via Bayesian state-space models in {R}

TL;DR

ssMousetrack is an R package that employs Bayesian state-space models to analyze computerized tracking data, offering researchers a comprehensive tool for data preparation, modeling, and diagnostics in behavioral studies.

Contribution

The paper introduces ssMousetrack, a novel R package that applies Bayesian state-space modeling to computerized tracking data, facilitating dynamic response process analysis.

Findings

Effective data modeling and analysis demonstrated in case study

Provides diagnostic tools for model assessment

Enhances understanding of response dynamics in behavioral research

Abstract

Recent technological advances have provided new settings to enhance individual-based data collection and computerized-tracking data have became common in many behavioral and social research. By adopting instantaneous tracking devices such as computer-mouse, wii, and joysticks, such data provide new insights for analysing the dynamic unfolding of response process. ssMousetrack is a R package for modeling and analysing computerized-tracking data by means of a Bayesian state-space approach. The package provides a set of functions to prepare data, fit the model, and assess results via simple diagnostic checks. This paper describes the package and illustrates how it can be used to model and analyse computerized-tracking data. A case study is also included to show the use of the package in empirical case studies.

Figures

Tables

Table 1. Table 1: Overview of the contents of ssMousetrack

function	type	description
generate_data()	main	simulate data according to a user-defined experimental design.
run_ssm()	main	run state-space model on a given mouse-tracking dataset.
check_prior()	main	allows users to define a list of priors for $f\left(𝜸\right)$ prior running run_ssm().
prepare_data()	main	pre-process raw tracking data prior running run_ssm().
evaluate_ssm()	main	run model evaluation given an output of run_ssm(). The function can plot results if requested by users.
compute_D()	internal	compute the matrix of distances $𝐃$ given the observed data $𝐘$ (see Equation 5).
generate_Z()	internal	generate the Boolean trial-by-variable (design) matrix $𝐙$ (see Equation 4).
generate_design()	internal	allows users to specify an experimental design in terms individuals, trials, variables, and design matrix $𝐙$.
congruency	dataset	subset of data from [8].
language	dataset	subset of data from [2].

Table 2. Table 2: Illustrative example: Posterior quantities and Gelman-Rubin indices (Rhat) for the model parameters.

	mean	sd	25%	50%	75%	n_eff	Rhat
gamma1	-0.05	0.19	-0.18	-0.05	0.08	3047.00	1.00
gamma2	-0.02	0.06	-0.06	-0.02	0.02	2764.00	1.00
gamma3	0.16	0.06	0.12	0.16	0.20	2782.00	1.00
gamma4	0.05	0.06	0.01	0.05	0.09	2680.00	1.00

$\left(n=0\right)$		${\widehat{𝐱}}_{I\times 1}^{\left(n\right)}={\mathrm{𝟎}}_{I\times 1}$
		${\widehat{𝝀}}_{I\times 1}^{\left(n\right)}={\mathrm{𝟏}}_{I\times 1}$
$\left(n>0\right)$		${\overline{𝐱}}_{I\times 1}^{\left(n\right)}={\widehat{𝐱}}_{I\times 1}^{\left(n-1\right)}$
		${\overline{𝝀}}_{I\times 1}^{\left(n\right)}={\widehat{𝝀}}_{I\times 1}^{\left(n-1\right)}$
		${\widehat{𝐲}}_{JI\times 1}=𝒢({\overline{𝐱}}_{I\times 1}^{\left(n\right)},{𝐙}_{J\times K}\cdot {𝜸}_{K\times 1}^{†})$
		${𝝈}_{JI\times 1}=({𝐈}_{J\times J}\otimes {\overline{𝝀}}_{I\times 1}^{\left(n\right)}){\mathrm{𝟏}}_{J\times 1}+1\oslash \sqrt{{\exp }^{†}({𝐝}_{JI\times 1}^{\left(n\right)}})$
		${𝐊}_{JI\times 1}=\left({𝐈}_{J\times J}\otimes {\overline{𝝀}}_{I\times 1}^{\left(n\right)}\right){\mathrm{𝟏}}_{J\times 1}\oslash {𝝈}_{JI\times 1}$
		${\widehat{𝐱}}_{I\times 1}^{\left(n\right)}={\overline{𝐱}}_{I\times 1}^{\left(n\right)}+{\left({\left(\left({𝐲}_{JI\times 1}-{\widehat{𝐲}}_{JI\times 1}\right)\odot {𝐊}_{JI\times 1}\right)}^T{𝐀}_{JI\times I}\right)}^T$
		${\widehat{𝝀}}_{I\times 1}^{\left(n\right)}={\overline{𝝀}}_{I\times 1}^{\left(n\right)}+{\left({\left({𝐊}_{JI\times 1}\odot {𝝈}_{JI\times 1}\odot {𝐊}_{JI\times 1}\right)}^T{𝐀}_{JI\times I}\right)}^T$

Equations

$$\displaystyle\mathbf{y}^{(n)}_{JI\times 1}\sim\text{vonMises}\Big{(}\bm{\mu}^{(n)}_{JI\times 1},\bm{\kappa}^{(n)}_{JI\times 1}\Big{)}$$

$$\displaystyle\bm{\mu}^{(n)}_{JI\times 1}=\mathcal{G}\Big{(}\mathbf{x}_{I\times 1}^{(n)},\bm{\beta}_{J\times 1}\Big{)}$$

$$\displaystyle\mathbf{x}^{(n)}_{I\times 1}\sim\text{Normal}\Big{(}\mathbf{x}^{(n-1)}_{I\times 1},\sigma_{x}\mathbf{I}_{I\times I}\Big{)}$$

$${\mathbf{\beta }}_{J\times 1}={\mathbf{Z}}_{J\times K}\cdot {\mathbf{\gamma }}_{K\times 1}$$

$$\displaystyle\bm{\kappa}_{JI\times 1}^{(n)}=\exp^{\dagger}\Big{(}\lambda\mathbf{d}_{JI\times 1}^{(n)}\Big{)}$$

$$\mathcal{G}=\text{vec}\Bigg{(}\pi^{-1}\Big{[}1+\exp\Big{(}\bm{\beta}_{J\times 1}\mathbf{1}_{1\times I}-\mathbf{1}_{J\times 1}\mathbf{x}_{1\times I}^{(n)}\Big{)}\Big{)}\Big{]}\Bigg{)}$$

$$\mathcal{G}=\text{vec}\Bigg{(}\pi\Big{[}\exp\Big{(}-\bm{\beta}_{J\times 1}\mathbf{1}_{1\times I}\exp\Big{(}\mathbf{1}_{J\times 1}\mathbf{x}_{1\times I}^{(n)}\Big{)}\Big{)}\Big{]}\Bigg{)}$$

$$\mathbf{d}^{(n)}=\begin{cases}\big{|}\mathbf{y}^{(n)}-\frac{3\pi}{4}\big{|},&\leavevmode\nobreak\ \text{if}\leavevmode\nobreak\ \mathbf{y}^{(n)}<\frac{\pi}{2}\\[5.69046pt] \big{|}\mathbf{y}^{(n)}-\frac{\pi}{4}\big{|},&\leavevmode\nobreak\ \text{if}\leavevmode\nobreak\ \mathbf{y}^{(n)}\geq\frac{\pi}{2}\\ \end{cases}$$

$$f(\mathbf{\Theta }|\mathbf{Y})\propto f(\mathbf{\gamma })f(\mathbf{\gamma }|\mathbf{Y})f(\mathbf{X}|\mathbf{Y})$$

$$\displaystyle\text{PA}_{\text{overall}}=1-\Big{(}||\mathbf{Y}_{m}-\mathbf{Y}||^{2}\big{/}||\mathbf{Y}||^{2}\Big{)}$$

$$\displaystyle\text{PA}_{\text{sbj}}=1-\Big{(}||\mathbf{Y}^{(i)}_{m}-\mathbf{Y}^{(i)}||^{2}\big{/}||\mathbf{Y}^{(i)}||^{2}\Big{)}$$

$$i=1,\dots ,I$$

$$\mathbf{P}^{(i)}_{M\times K}=\Big{[}1+\exp\Big{\{}\mathbf{\hat{x}}^{(i)}_{M\times 1}\mathbf{1}_{1\times K}-\mathbf{1}_{M\times 1}\bm{\hat{\gamma}}_{1\times K}\Big{\}}\Big{]}^{-1}$$

$$\mathbf{r}^{(i)}_{K\times 1}=\log\Big{(}{\mathbf{{p}}^{(i)}_{K\times 1}}\Big{/}{1-\mathbf{{p}}^{(i)}_{K\times 1}}\Big{)}$$

Full text

Keywords: state space models, mouse-tracking, dynamic data, bayesian data analysis

1 Introduction

Recent technological advances allow the collection of detailed data on ratings, attitudes, and choices during behavioral tasks. Unlike standard surveys and questionnaires, these tools provide a rich source of data as they adopt tracking devices that collect subject-based information about the dynamics involved during the data collection task [10, 24]. Examples of such devices include eye-tracking, body movement-tracking, computer mouse-tracking, and electrodermal activity. Among these, computer mouse-tracking has become an important and widely used tool in behavioral sciences, as it provides a valid and cost-effective way to measure the usually unknown processes underlying human ratings and decisions [9]. Mouse-tracking data are often collected by means of standard computer mouse, wii instruments, and joystick devices and consist of collections of real-time trajectories recorded during the behavioral task. In a typical mouse-tracking task, individuals are presented with a computer-based interface showing the stimulus at the bottom of the screen (e.g., the image of a “dolphin”) and two labels on the left and right top corners (e.g., the labels “mammal” vs. “fish”). They are asked to decide which of the two labels is appropriate given the task instruction and stimulus (e.g., to decide whether dolphin is mammal or fish). In the meanwhile, the x-y trajectories of the computer device are instantaneously recorded. The real-time trajectories offer an effective way to study the decision process underlying the hand movement behavior by revealing, for instance, the presence of some levels of decisional uncertainty. The applications of mouse-tracking tools spread across many research area, including cognitive sciences [8], neuroscience [33], neurology [28], and forensic studies [25].

Several tools for running mouse-tracking analyses are available in open-source specialized software like MouseTracker [10], EMOT [4], and MouseTrap [20]. In the R environment, only the packages mousetrack [7] and mousetrap [19] are devoted to mouse-tracking data. More generally, there are other packages developed to handle with tracking data such as trajectories [27], trackeR [11], adehabtatLT [5], and move [21]. Similarly, there are many packages developed for state-space models like, for instance, KFAS [17] and bssm [16].

In this paper we present ssMousetrack, a novel R package to analyse computerized tracking data as they emerge from typical mouse-tracking data recording. The package implements a non-linear state space model to handle with the dynamics of mouse-tracking data. The model is estimated using approximated Kalman filter coupled with MCMC algorithms via the rstan package [31]. The package includes functions for data pre-processing, data generation, and model assessment. It also provides functionalities to set-up designs for mouse-tracking data recording. Despite other R packages are available in this context, ssMousetrack differs in some aspects. For instance, mousetrack and mousetrap focus on descriptive evaluation of mouse-tracking data static measures (e.g., minima, maxima, flips, curvature). By contrast, our package (i) implements a dynamic evaluation of the trajectory data without resorting the use of summary measures and (ii) evaluates the observed data variability in terms of latent dynamics and external covariates (e.g., experimental variables), which are usually of relevance in mouse-tracking data analysis. With respect to KFAS and bssm, our package offers a more focused implementation for mouse-tracking data. Moreover, although complete and useful in many cases, these packages implement a general class of state-space representations which may not be widely applicable to computerized tracking data. Finally, ssMousetrack differs also from trajectories, trackeR, adehabtatLT, and move as they focus on animal tracking and related problems, such as estimation of habitat choices. This makes them not directly suitable for analysing the various aspects of mouse-tracking studies. In this respect, the advantages of ssMousetrack are as follows: (i) it is easy to use as it requires typing a single function to run the entire procedure, (ii) it takes advantages of rstan package to estimate model parameters via MCMC, (iii) it allows modeling and analysing computerized tracking data as they are usually recorded in typical mouse-tracking tasks, (iv) it provides a user-friendly workflow for all processing steps which can be easily understood by non-expert users, (v) it offers users a way to simulate mouse-tracking designs and data as well. In addition, ssMousetrack can be combined with other R packages, including shinystan [32] and ggmcmc [34] to produce further statistical and graphical representation of the output. The package is available from the Comprehensive R Archive Network at https://CRAN.R-project.org/package=ssMousetrack.

This paper is organized as follows. Section 2 gives a statistical overview of the model implemented in the package, the methods of estimation and inference, and the assessment of the model. Section 3 describes the package’s structure and its utilities. Section 4 illustrates the functioning of the package by means of an illustrative case study. Finally, Section 5 concludes the manuscript with a discussion and future directions.

2 Model

Computerized mouse-tracking data typically consist of arrays $(\mathbf{x},\mathbf{y})_{ij}\in\mathbb{R}^{N_{ij}}\times\mathbb{R}^{N_{ij}}$ containing the streaming of x-y Cartesian coordinates of the computer-mouse pointer, for $i=1,\ldots,I$ subjects, $j=1,\ldots,J$ stimuli, and $n=1,\ldots,N_{ij}$ time steps. To simplify data analysis, raw trajectories are usually pre-processed according to the following steps [15, 4]. First, the trajectories $(\mathbf{x},\mathbf{y})_{ij}$ are normalized on a common sampling scale such that $N$ is the same over subjects and stimuli. Next, the arrays $(\mathbf{x},\mathbf{y})_{ij}\in\mathbb{R}^{N}\times\mathbb{R}^{N}$ are rotated and translated into the quadrant $[-1,1]\times[1,1]$ with $(x_{0},y_{0})_{ij}=(0,0)$ and $(x_{N},y_{N})_{ij}=(1,1)$ by convention. Finally, normalized data are projected onto a (lower) 1-dimension space via atan2 function. In this way, the final ordered data \mathbf{y}_{ij}=\big{(}y^{(1)},\ldots,y^{(n)},\ldots,y^{(N)}\big{)}\in(0,\pi]^{N} lie on the arc defined by the union of two disjoint sets, i.e. the set $\{y\in(0,\pi]:y\geq\frac{\pi}{2}\}$ which represents the right-side section of the screen (usually called target, T) and the set $\{y\in(0,\pi]:y<\frac{\pi}{2}\}$ which instead represents the left-side section (usually called distractor, D). The final data are arranged as an $\mathbf{Y}_{N\times JI}$ column-wise stacked matrix.

The state-space model implemented in ssMousetrack is as follows:

[TABLE]

where Equation (1) is a von Mises measurement equation with $\bm{\mu}^{(n)}_{JI\times 1}\in(0,\pi]^{IJ}$ being the mean for the $n$-th time step and $\bm{\kappa}^{(n)}_{JI\times 1}\in\mathbb{R}^{IJ}$ the concentrations around the $n$-th mean vector, Equation (2) represents the locations on the arc defined in $(0,\pi]$ from which the data vector $\mathbf{y}^{(n)}$ is sampled from and it behaves according to the a real function $\mathcal{G}:\mathbb{R}\to(0,\pi]$, which maps reals into radians. In the current version of ssMousetrack, $\mathcal{G}$ can be defined as:

• (i)

$\pi$-scaled logistic function:

[TABLE]

with $\bm{\beta}_{J\times 1}\in\mathbb{R}^{J}$ representing the contribution of the stimuli on $\mathbf{y}^{(n)}$.

• (i)

$\pi$-scaled Gompertz function:

[TABLE]

where $\bm{\beta}_{J\times 1}\in\mathbb{R}_{+}^{J}$ has the same meaning as before.

Although they represent two cases of the general family of $S$-shaped functions, logistic and Gompertz models differ in some respects. For instance, unlike the logistic model, the Gompertz function is not symmetric around its inflection point, with the consequence that its growth rises rapidly to its maximum rate occurring before the fixed inflection point [23]. Moreover, the parameters of the Gompertz function are always positive, a constrain which is often required by applications where the covariates of the model cannot take negative values (e.g., reaction times). These two implementations allow users to choose the type of $\mathcal{G}$ function on the basis of the experimental designs they have used in their studies.

In Equation (3) is a Normal states equation which represents a lag-1 autoregressive process with time-fixed variance parameter $\sigma_{x}$. In the current version of ssMousetrack, the covariance matrix of the latent processes is set to an identity matrix $\mathbf{I}$ without loss of generality ($\sigma_{x}=1$). Equation (4) is the linear term modeling the contribution of the experimental design (e.g., two-by-two design) and variables involved (e.g., categorical variables, continuous covariates). Note that $\mathbf{Z}$ is a design (dummy) matrix of main and high-order effects defined by adopting the dummy coding (e.g., treatment contrasts, sum contrasts) whereas $\bm{\gamma}$ is the associated vector of parameters for the columns of $\mathbf{Z}$, with $\gamma_{1}$ being the usual baseline term for the contrasts. Finally, Equation (5) defines the concentrations around the $n$-th location by using the transformed data:

[TABLE]

with $\exp^{\dagger}:(0,\pi]\to[\text{lb},\text{ub}]\subset\mathbb{R}_{+}$ being the exponential function scaled in the natural range of the concentration parameter (e.g., $lb=10$, $ub=200$). In the current implementation of the package, the parameter $\lambda$ is fixed to unity.

The interpretation of Equations (2)-(4) is as follows. The $n$-th mean vector $\bm{\mu}^{(n)}$ is expressed as function of the stimuli-related component $\bm{\beta}$ and subject-based component $\mathbf{x}^{(n)}$, which are integrated together to form the conditional sampling $\mathbf{y}^{(n)}|\bm{\beta},\mathbf{x}^{(n)}$ through the function $\mathcal{G}$. As a result, Equation (3) can be interpreted as the individual latent dynamics that are unaffected by the experimental stimuli whereas Equation (4) represents the experimental effect regardless to individual dynamics. More generally, Equation (3) conveys information about the hand movement process underlying the tracking device and as such it can be used to analyse how much individuals differ in executing the task. By contrast, Equation (4) collects information on how a certain experimental manipulation has an effect or not on the movement responses. Interestingly, when normalized into $[0,1]$, $\bm{\mu}^{(n)}$ can be interpreted as the probability of the $i$-th individual at the $j$-th stimulus to navigate close the distractor cue in the left-side section of the arc. Finally, Equation (5) follows from the fact that hand movements underlying computerized tracking data tend to be smooth over the experimental task, with small changes being more likely close to left (distractor) or right (target) endpoints [3].

2.1 Estimation and inference

The state-space model in Equations (1)-(5) requires estimating the array of latent trajectories $\mathbf{X}\in\mathbb{R}^{I\times N}$ together with the array of parameters $\bm{\gamma}_{K\times 1}$, with $\gamma_{1}\in\mathbb{R}$ and $\bm{\gamma}_{(K-1)\times 1}\in\mathbb{R}^{K-1}$ (logistic case) or $\bm{\gamma}_{(K-1)\times 1}\in[-\gamma_{1},\infty)^{K-1}$ (Gompertz case). The array of unknown quantities $\bm{\Theta}=\{\mathbf{X},\bm{\gamma}\}$ can be estimated in various way, by adopting both a frequentist and Bayesian perspectives [30]. In the ssMousetrack package, the parameters are recovered in Bayesian way by means of a marginal MCMC algorithm through which $\mathbf{X}$ and $\bm{\gamma}$ are alternately recovered [1, 29]. The reason is twofold: (i) MCMC algorithms, as those implemented in rstan package, provide a more efficient and complete solution for sampling from the probability distribution of the parameters. (ii) The Bayesian approach offers an elegant solution for data analysis and inference [13] by means of which the model is adequately assessed by the analysis of (marginal) posterior distributions of the parameters [22].

More in details, the posterior density $f(\bm{\Theta}|\mathbf{Y})$ after factorization of the joint density $f(\bm{\gamma},\mathbf{X},\mathbf{Y})$, is as follows [1]:

[TABLE]

where $f(\bm{\gamma}|\mathbf{Y})$ is the (marginal) likelihood function, $f(\mathbf{X}|\mathbf{Y})$ is the filtering density, whereas $f(\bm{\gamma})$ is the prior ascribed on the model parameters. In the current version of ssMousetrack, $f(\mathbf{X}|\mathbf{Y})$ is approximated via Kalman filtering/smoothing, with $f(\bm{\gamma}|\mathbf{Y})$ being computed as a byproduct of the Kalman theory (see Appendix A).

2.2 Model assessment

In the Bayesian context of data analysis, ssMousetrack provides a simulation-based procedure to evaluate the adequacy of the model to reproduce the observed data $\mathbf{Y}$ [13]. More technically, given the posteriors of parameters and latent states $f(\bm{\Theta}|\mathbf{Y})$, $M$ new (simulated) datasets $\mathbf{Y}^{*}_{1},\ldots,\mathbf{Y}^{*}_{M}$ are generated according to the estimated model structure and, for each new dataset, two discrepancy measures are considered [18]:

[TABLE]

which measure the total amount of data reconstruction based on the overall $JI\times N$ observed array $\mathbf{Y}$ (Equation 7) and the amount of data reconstruction based on the $J\times N$ observed matrix $\mathbf{Y}^{(i)}$ for each subject $i=1,\ldots,I$ (Equation 8). Both the indices are in the range 0-100%, with 100% indicating optimal fit. Note that the measure $\text{PA}_{\text{sbj}}$ allows for evaluating the adequacy of the model to reconstruct the individual-based set of data. In addition, the dynamic time warp distance (dtw), as implemented in dtw package, is also computed between $\mathbf{Y}^{(i)}_{m}$ and $\mathbf{Y}^{(i)}$. Unlike the $\text{PA}_{\text{sbj}}$ index, the dtw distance measures the similarity among time series by considering their different dynamics [14].

3 The ssMousetrack package

The ssMousetrack is distributed via the Comprehensive R Archive Network (CRAN). It is based on rstan [31], the R interface to the probabilistic programming language Stan [6]. The current version of the package allows for (i) simulating data according to a given experimental design, (ii) analysing mouse-tracking data via state-space modeling, and (iii) evaluating the adequacy of model results. The package consists of five main function (generate_data(), run_ssm(), check_prior(), prepare_data(), evaluate_ssm()), two datasets (language, congruency), and three sub-functions (compute_D(), generate_Z(), generate_design()). The main functions generate_data() and run_ssm() are wrappers to previously-compiled Stan codes which implement the model described in Section 2. Table 1 provides an overview of the functions and datasets provided in the ssMousetrack package whereas a description of the usage of the functions is reported in the next subsections.

3.1 Generate artificial data

To simulate artificial data we use the function generate_data(), which requires as input the experimental template for the data generation process. More generally, the function works by first sampling the parameters $\bm{\gamma}$ from the prior density $f(\bm{\gamma})$ and then generates the latent states $\mathbf{X}$ from Equation (3), computes the matrix $\bm{\mu}$ from Equation (2) and the matrix $\mathbf{D}$, drawns the matrix of data $\mathbf{Y}$ by simply applying Equation (1). For instance, an experiment with one categorical independent variable and two levels, each with three trials, can be generated via the following syntax:

Ψprior_list <- list("normal(-0.25,0.5)","normal(2.7,1)") Ψdatagen1_ssm <- generate_data(N = 61, M = 100,I = 2,J = 6, Ψ+ K = 2,Z.formula = "~Z1", priors = prior_list)

where M = 100 is the number of data to be generated, N = 61 is the number of time step for the mouse-tracking trajectories, K = 2 means that we have just one variable with two levels, J = 6 indicates the total number of trials such that J/K is the number of trials for each level of the variable, I = 2 is the number of subject, Z.formula indicates the formula for the contrast matrix $\mathbf{Z}$ with standard R syntax. Note that selective priors are specified for each level of the experimental design using the Stan syntax (see the help of check_prior() for further details).

The output is a list containing three sublists, as follows:

• •

params, which contains the model parameters generated for the $M$ datasets:

ΨΨ## List of 4 ΨΨ## sigmax: num [1:2] 1 1 ΨΨ## lambda: num [1:12] 1 1 1 1 1 ... ΨΨ## gamma : num [1:75, 1:2] 0.228 -0.378 ... ΨΨ## beta : num [1:75, 1:12] 0.228 -0.378 ... Ψ

• •

data, which contains the matrices of latent states $\mathbf{X}$ and trajectories $\mathbf{Y}$, together with $\bm{\mu}$, $\mathbf{D}$, and $\mathbf{Z}$:

ΨΨ## List of 5 ΨΨ## Y : num [1:75, 1:61, 1:12] 1.54 1.53 ... ΨΨ## X : num [1:75, 1:61, 1:2] 1e-04 1e-04 1e-04 1e-04 1e-04 ... ΨΨ## MU: num [1:75, 1:61, 1:12] 1.57 1.57 ... ΨΨ## D : num [1:75, 1:61, 1:12] 0.785 0.785 ... ΨΨ## $ Z : num [1:12, 1:2] 1 1 1 1 1 ... ΨΨ

• •

design, which contains the experimental design used as template to generate the data:

ΨΨ## sbj trial Z1 ΨΨ## 1 1 1 100 ΨΨ## 2 1 2 100 ΨΨ## 3 1 3 100 ΨΨ

Similarly, artificial datasets can be generated using more complex designs. For instance, a bivariate design with two variables is produced by typing:

ΨΨdatagen2_ssm <- generate_data(I = 2,J = 8,K = c(2,4),Z.formula = "~Z1*Z2", ΨΨ+ Z.type=c("symmetric","random")) ΨΨ

where K = c(2,4) codifies two variables each with two and four levels, Z.formula = " Z1*Z2" indicates that the variables interact whereas Z.type=c("symmetric","random") indicates that trials must be assigned to the first variable using the symmetric method and to the second variable using the random method (see the help of generate_Z() for further details).

Figure 1 shows a sample of mouse-tracking data $\mathbf{Y}$ generated in the univariate design case with I = 2, K = 2 and J = 6. We report the univariate case only for the sake of simplicity but the same graphical representations can be done for the more complex designs as well.

3.2 Run state-space analysis

State-space analysis can be run on both simulated and real data. In the first case, after the data-generation process, the state-space model implemented in the ssMousetrack package can be fit using run_ssm(). For instance, the syntax:

ΨΨdatagen2_ssm <- generate_data(I = 2,J = 8,K = c(2,4),Z.formula = "~Z1*Z2", ΨΨ+ Z.type=c("symmetric","random")) ΨΨiid <- 2 ΨΨdatagen2_fit <- run_ssm(N = datagen2_ssm$N,I = datagen2_ssm$I, ΨΨ+ J = datagen2_ssm$J,Y = datagen2_ssm$data$Y[iid,,], ΨΨ+ D = datagen2_ssm$data$D[iid,,],Z = datagen2_ssm$data$Z, ΨΨ+ niter = 5000,nwarmup = 2000,nchains = 2) Ψ

runs the state-space analysis on the iid = 2 artificial data datagen2_ssm. Note that niter indicates the number of total samples to be drawn, nwarmup the number of warmup/burnin iterations per chain, and nchains the number of chains to be executed in parallel. The function run_ssm() allows for parallel computing via the parallel package when nchains > 1. In this case, since ncores="AUTO" (default), the function will run two parallel chains using two cores.

Unlike for the case of artificial data, the analysis of real datasets requires preparing raw data in a proper format via prepare_data(), the function that implements the steps described in Section 2. Generally, raw datasets need to be organized using the long-format, with information being organized as nested. The dataset language is an example of a typical data structure required by prepare_data():

ΨΨdata("language") ΨΨstr(language,vec.len=2) ΨΨ## ’data.frame’:Ψ6060 obs. of 6 variables: ΨΨ## sbj : int 1 1 1 1 1 ... ΨΨ## trial : int 1 1 1 1 1 ... ΨΨ## condition: Factor w/ 4 levels "HF","LF","PW",..: 1 1 1 1 1 ... ΨΨ## timestep : int 1 2 3 4 5 ... ΨΨ## x : num 0 -0.0098 -0.0098 -0.0098 -0.0098 ... ΨΨ## y : num 0 0.0025 0.0025 0.0025 0.0025 ... Ψ

where condition is the categorical variable involved in the study. The pre-processing of raw data is performed by the call:

ΨΨlanguage_proc <- prepare_data(X = language,N = 61,Z.formula = "~condition") Ψ

where the output language_proc is a data frame containing the pre-processed dataset together with the column-wise stacked matrix $\mathbf{Y}$ of angles, the contrast matrix $\mathbf{Z}$, and the matrix of distances $\mathbf{D}$.

Once raw data have been pre-processed, the state-space analysis is performed as for the case of artificial data:

ΨΨlanguage_fit <- run_ssm(N = language_proc$N,I = language_proc$I, ΨΨ+ J = language_proc$J,Y = language_proc$Y, ΨΨ+ D = language_proc$D,Z = language_proc$Z, ΨΨ+ niter = 5000,nwarmup = 2000,nchains = 2) Ψ

The function returns as output a list composed of three sublists, as follows:

• •

params, which contains the posterior samples for the free parameters $\bm{\gamma}$ and $\bm{\beta}$:

ΨΨΨ## List of 6 ΨΨΨ## sigmax : num 1 ΨΨΨ## lambda : num 1 ΨΨΨ## kappa_bnds: num [1:2] 5 300 ΨΨΨ## gamma :’data.frame’:Ψ4000 obs. of 4 variables: ΨΨΨ## beta : num [1:4000, 1:60] -0.26 -0.146 ... ΨΨΨ## :function (z, ...) ΨΨ

• •

data, which contains the posterior samples for the latent states $\mathbf{X}$ and the moving means $\bm{\mu}$:

ΨΨΨ## List of 6 ΨΨΨ## Y : num [1:101, 1:60] 1.56 1.7 ... ΨΨΨ## X : num [1:4000, 1:101, 1:5] 1e-04 1e-04 1e-04 1e-04 1e-04 ... ΨΨΨ## MU : num [1:4000, 1:101, 1:60] 1.76 1.68 ... ΨΨΨ## D : num [1:101, 1:60] 0.592 0.474 ... ΨΨΨ## Z : num [1:60, 1:4] 1 1 1 1 1 ... ΨΨΨ## X_smooth: num [1:4000, 1:101, 1:5] -0.0878 -0.0635 ... ΨΨ

• •

stan_table, containing the typical Stan output (i.e., point estimates, credibility intervals, and Gelman-Rubin index) for the sampling() method as implemented in the rstan package:

ΨΨΨ## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat ΨΨΨ## gamma[1] -0.05 0 0.19 -0.43 -0.18 -0.05 0.08 0.33 3047 1 ΨΨΨ## gamma[2] -0.02 0 0.06 -0.13 -0.06 -0.02 0.02 0.09 2764 1 ΨΨΨ## gamma[3] 0.16 0 0.06 0.04 0.12 0.16 0.20 0.28 2782 1 ΨΨ

Note that users can also export the stanfit object with all the Stan results by specifying stan_object=TRUE in run_ssm().

The function run_ssm() allows for different priors specification. In particular, users can specify different priors for the model parameters $\bm{\gamma}$ as follows:

ΨΨpriors_list <- list("lognormal(1,0.5)","normal(1,2)","chi_square(2)","normal(0,1)") ΨΨlanguage_fit <- run_ssm(..., priors = priors_list) Ψ

which means that $\gamma_{1}\sim\text{lognormal}(1,0.5)$, $\gamma_{2}\sim\text{normal}(1,2)$, $\gamma_{3}\sim\text{chi\_square}(2)$, $\gamma_{4}\sim\text{normal}(0,1)$. The list of probability distributions accepted by run_ssm() is described in the help of the function check_prior(). Specification of priors for single parameters is also allowed, by using NULL attributes:

ΨΨpriors_list <- list(NULL,"normal(1,2)","chi_square(2)",NULL) ΨΨlanguage_fit <- run_ssm(..., priors = priors_list) Ψ

where predefined priors are used for parameters $\gamma_{1}$ and $\gamma_{4}$. Further examples about run_ssm() are illustrated in the manual of the package.

3.3 Evaluate the model results

The methods described in Section 2.2 for the model evaluation are implemented by the function evaluate_ssm(), which requires as input the output of run_ssm(). For instance, considering the fitted object language_fit, the model evaluation can simply be run via the command:

ΨΨlanguage_eval <- evaluate_ssm(ssmfit = language_fit, M = 1000, plotx = FALSE) Ψ

where M = 1000 is the number of replications to compute the indices. The function returns as ouput a list containing the mean values of the indices $\text{PA}_{\text{overall}}$, $\text{PA}_{\text{sbj}}$, and dtw, as well as the distributions obtained over the M replications. Note that, users can also ask for a graphical representation of the indices by setting plotx = TRUE.

4 An Illustrative example

In this section we provide a full example about the way ssMousetrack can be used for state-space analysis of real computerized tracking data. Note that the application reported here has an illustrative purpose only. To this end, we will make use of the dataset language, a subset of data originally presented in [2]. In this typical computerized tracking task, participants saw a printed stimulus on the screen (e.g., the word water) and were requested to perform a dichotomous choice task where stimuli need to be classified as word or non-word. The experimental variable condition was a categorical variables with four levels (HF: High-frequency word; LF: Low-frequency word; PW: Pseudo-word; NW: Non-word). Participants had to classify each stimulus as word vs. non-word by using a computer-mouse tracking device. The dataset contains $I=5$ participants, $J=12$ trials, one categorical variable with $K=4$ levels, each with $J/K=3$ trials. From the data-analysis viewpoint, we evaluat the extent to which the parameters of the state-space model $\bm{\gamma}$ reflect eventual differences associated with the levels of condition.

The raw computerized tracking trajectories in the dataset consist of Cartesian coordinates with $N=101$ ($i=1,\ldots,I$; $j=1,\ldots,J$). The dataset is partially pre-processed as raw trajectories have the same length ($N=101$). However, we need to run prepare_data() in order to rotate/translate the raw data into the quadrant $[-1,1]\times[1,1]$ and compute the atan2 projections. The pre-processing step is called by the command:

Ψdata("language") Ψlanguage_proc <- prepare_data(X = language, N = 101, Z.formula = "~condition")

Figure 2 shows the trajectories $\mathbf{Y}$ associated with the task for all participants and trials.

Next, the state-space model is fit to the pre-processed data by the following call:

Ψpriors_list <- list("normal(0,1)","normal(1,1)","normal(-2,1)","normal(2,1)") Ψlanguage_fit <- run_ssm(N = language_proc$N,I = language_proc$I, Ψ+ J = language_proc$J,Y = language_proc$Y, Ψ+ D = language_proc$D,Z = language_proc$Z, Ψ+ niter = 6000,nwarmup = 2000,nchains=4, Ψ+ priors = priors_list, Ψ+ gfunction = "logistic")

where, in this case, the prior for $\bm{\gamma}$ have been choosen to codify a priori expectations about the effect of the variable condition [2]. Figure 3 shows some MCMC graphical diagnostics for the model parameters $\bm{\gamma}$ computed using bayesplot [12] whereas Table 2 reports the posterior quantities for the model parameters. In the Bayesian context of data-analysis, we evaluate the effects of the variable condition by computing the degree of overlapping among marginal posterior densities for each level of the experimental variable (i.e., the more the overlapping, the weaker the evidence supporting the experimental manipulation). Figure 4 shows the results graphically. Overall, the variable condition showed no strong effect, as the densities of the levels are overlapped. In particular, stimuli in HF, LF, and NW conditions showed no activation of the distractor section of the tracking space as $\hat{\gamma}_{\text{HF}}$, $\hat{\gamma}_{\text{NW}}$, and $\hat{\gamma}_{\text{LF}}$ approach zero. By contrast, stimuli in PW condition showed a small effect on activating the target section ($\hat{\gamma}_{\text{HF}}>0$), possibly due to the fact that PW stimuli require less cognitive workload [2].

Figure 5-A reports the filtered latent states $\mathbf{\hat{X}}$ for the subjects in the dataset. To further investigate how individual dynamics differ over the levels of condition, we can make use of $\mathbf{\hat{X}}$ and ask whether HF, LF, NW, and PW stimuli differ in terms of evidence of mouse-tracking competition. The idea is that the higher the evidence, the larger the difficulty in categorizing stimuli as word (target) or non-word (distractor).

To do this, we follow the findings of [2] and divide the entire respose process $1,\ldots,N$ into three disjoint windows $W_{1}=10-35\%$, $W_{2}=45-65\%$, and $W_{3}=70-85\%$. Usually it is expected that a higher competition would be observed in $W_{1}$ and $W_{2}$ rather than $W_{3}$. More formally, let $\mathbf{\hat{x}}^{(i)}_{M\times 1}=(\hat{x}_{1}^{(i)},\ldots,\hat{x}_{M}^{(i)})$ be the sequence of filtered states for the $i$-th subject and the generic time window $W$, with $M$ being equals to the cardinality of $W$. Next, the probability to select non-word (distractor) responses are computed by normalizing the $\mathcal{G}$ function into the domain $[0,1]$, as follows:

[TABLE]

where $\bm{\hat{\gamma}}$ is the array of posterior means of the model parameters. Note that in this example we use the logistic function because we set gfunction="logistic" in run_ssm(). Finally, the evidence measures can be defined in terms of log-odd ratio using the probability matrix $\mathbf{P}^{(i)}$:

[TABLE]

where \mathbf{p}_{K\times 1}=\frac{1}{M}\Big{(}\mathbf{1}_{1\times M}\mathbf{P}^{(i)}_{M\times K}\Big{)}^{T} is the profile probability for HF, LF, NW, and PW. The interpretation of $\mathbf{r}^{(i)}$ is as follows. For $\mathbf{r}^{(i)}>0$ there is a higher competition in categorizing the stimulus as word (target) vs. non-word (distractor). By contrast, for $\mathbf{r}^{(i)}<0$ there is a lower competition in the response process, as stimuli are easily categorized as word (target). Finally, the case $\mathbf{r}^{(i)}=0$ indicates that there is no difference in terms of evidence between word (target) and non-word (distractor) responses. Figure 5-B shows the results for the four levels of condition. As expected, the competition in the third phase of the response process $W_{3}$ is low, as the probability to select the target is higher. The same applies to $W_{2}$. On the contrary, in the first stage of the process $W_{1}$ the competition is higher although the evidence ratio for all the levels of condition approximate zero. Interestingly, the second phase $W_{2}$ shows a higher whithin-subject variability of competition, which probably indicates that subjects differ in the categorization process just in the middle phase of the response process.

Finally, we assess the adequacy of the model with regards to the observed data by means of evaluate_ssm(), as follows:

Ψlanguage_fit_eval <- evaluate_ssm(ssmfit = language_fit,M = 500,plotx = FALSE)

where language_fit is the fitted object returned by run_ssm() whereas M = 500 is the number of replications used to compute the three fit indices. The output of the function consists of a list containing means and distributions of the fit indices:

Ψ## List of 2 Ψ## dist :List of 3 Ψ## .. PA_ov : num [1:500] 0.944 0.938 ... Ψ## .. PA_sbj: num [1:500, 1:5] 0.991 0.991 ... Ψ## .. DTW : num [1:500, 1:60] 0.0945 0.1105 ... Ψ## indices:List of 3 Ψ## .. PA_ov : num 0.936 Ψ## .. PA_sbj: num 0.99 Ψ## .. DTW : num 0.119

Overall, in this example the fitted model is adequate to reproduce the observed trajectory data as supported by high values of the indices $\text{PA}_{\text{ov}}$, $\text{PA}_{\text{sbj}}$, and dtw. Figure 6 shows the results graphically.

5 Conclusions

In this paper we introduced the R package ssMousetrack that analyses computerized-tracking data using Bayesian state-space modeling. The package provides a set of functions to facilitate the preparation and analysis of tracking data and offers a simple way to assess model fit. The package can be of particular interest to researchers needing tools to analyse computerized-tracking experiments using a complete statistical modeling environment instead of descriptive statistics only. In addition, the package ssMousetrack allows for individual-based analysis of trajectories where latent dynamics are used to obtain richer information which can pave the way to further analyses (e.g., profile analysis). The current version of the package can be extended in several ways. For instance, the inclusion of other state-space representations beyond the simple Gaussian AR(1) model can be a further generalization of the package. Still, model parameters like $\sigma_{x}$ and $\lambda$ can be free to allow for multi-group analysis. Similarly, more comprehensive model diagnostics could also be considered in future releases of the package.

Finally, we believe our package may be a useful tool supporting researchers and practitioners who want to make analysis of computerized-tracking experiments using a statistical modeling environment. This will surely help them to improve the interpretability of data analysis as well as the reliability of conclusions they can draw from their studies.

Appendix A Appendix

Given a candidate sample $\bm{\gamma}^{\dagger}$, the mean $\mathbf{x}$ and variance $\bm{\lambda}$ of the density $f(\mathbf{X}|\mathbf{Y})$ are approximated via the following recursion:

where $\otimes$ is the Kronecker product, $\odot$ the (element-wise) Hadamard product, $\oslash$ the (element-wise) Hadamard division, whereas \mathbf{A}=\mathbf{I}_{I\times I}\otimes\big{(}n\mathbf{1}_{J\times 1}\big{)} is a scaling matrix with $n=1/J$. As a byproduct of the Kalman filter, the marginal likelihood $f(\bm{\gamma}^{\dagger}|\mathbf{Y})$ is multivariate Gaussian with mean $\hat{\mathbf{y}}$ and variance $\text{diag}(\bm{\sigma})$, with $\text{diag}()$ being the linear operator that transforms a vector into a diagonal matrix. Finally, the array $\hat{\mathbf{X}}_{I\times N}$ contains the filtered latent states implied by the model whereas $\hat{\bm{\Lambda}}_{I\times N}$ is the array of variances associated with the filtered states. The smoothing part of the algorithm is implemented using the fixed-interval Kalman smoother [29] where the filtered arrays $\hat{\mathbf{X}}_{I\times N}$ and $\hat{\bm{\Lambda}}_{I\times N}$ are used as input of the backward recursion.