2.1.1 High dimensionality

The first challenge to address fixation sequences is their high dimensionality and nonstandard format. Consider the fixation sequence data from 70 participants answering 10 questions in test block 1. Treating each participant’s sequence for one question as a single observation gives 700 observations, each with a varying number of variables. As shown in Figure 3, the length of eye-tracking sequences varies widely. For instance, the average sequence length is 430, but some exceed 1,000, reaching up to 2,000, far surpassing the sample size. This high-dimensional sequence data presents challenges for traditional statistical methods.

Figure 3

The distribution of fixation sequence lengths for test block 1.

Figure 3 Long description

The histogram plots Frequency on the vertical Y-axis, ranging from 0 to 70 in increments of 10, against fixation sequence length on the horizontal X-axis, ranging from 0 to 2000 in increments of 500.

* The distribution is positively skewed, with the highest concentration of data points occurring at the lower end of the scale.

* A primary peak reaches a frequency of 70 at a sequence length of approximately 250.

* A secondary, smaller peak occurs around the 500 mark with a frequency of approximately 40.

* Beyond 500, the frequency steadily declines, showing a long tail that extends toward 2000.

* Between 1000 and 2000, frequencies remain low, mostly under 10, with several small gaps where no data is recorded.

Figure 4

Distribution of fixation duration for each action (test block 1).

Figure 4 Long description

A violin plot displays the distribution of fixation durations across six categories on the X-axis. The Y-axis is labeled Duration Time in m s, ranging from 0 to 600 with increments of 200. Each category features a colored violin shape representing density, with an internal white box plot showing the median, interquartile range, and whiskers extending to outliers.

Moving from left to right along the X-axis:

* look: Light green violin. The distribution is concentrated below 200 m s but has a very long thin tail reaching up to 700 m s. The median is approximately 100 m s.

* look_1: Pale yellow violin. Similar distribution to look, with a tail reaching nearly 600 m s and a median around 100 m s.

* look_2: Light purple violin. Distribution concentrated at the bottom with a tail reaching 500 m s. Median is roughly 100 m s.

* option_correct: Light red violin. A wider base at the bottom with a shorter tail reaching approximately 350 m s. The median is lower, around 60 m s.

* option_wrong: Light blue violin. Similar shape to option_correct, with a tail reaching nearly 400 m s and a median around 70 m s.

* question: Light orange violin. A wider distribution in the lower range with a tail reaching 500 m s. The median is approximately 80 m s.

A legend on the right side maps the colors to the action names: look, look_1, look_2, option_correct, option_wrong, and question.

2.1.2 Motivation of considering fixation duration

To preliminarily examine whether integrating fixation duration with fixation location provides a more informative basis for distinguishing testing and learning patterns, we analyzed the distribution of fixation durations in Test Block 1 (Figure 4). The results show that fixation duration distributions vary notably across locations. For instance, locations, such as “look” and “question,” display wider ranges and longer tails, indicating greater variability, whereas “option_correct” and “option_wrong” exhibit more compact distributions. These findings suggest that fixation duration captures meaningful differences across action types and should be examined in conjunction with fixation location sequences when analyzing testing and learning behaviors.

2.1.3 Motivation of considering elapsed time

To examine how elapsed time between consecutive fixations influences transitions, we divided all fixation pairs into five subsets based on quartiles of the elapsed time distribution: Q1 (elapsed times below the 1st quartile), Q2 (between the 1st quartile and median), Q3 (median to 3rd quartile), Q4 (above the 3rd quartile but below 3000 ms), and a final subset for extremely long elapsed times (>3,000 ms). For each subset, we constructed a $6 \times 6$ transition matrix to represent the probabilities of transitioning between locations. Each element $(i,j)$ in the matrix represents the probability of transitioning from action i to action j, estimated using the proportion of observed transitions in this cell relative to all transitions. Figure 5 displays these matrices from test block 1. In Q1, participants primarily remain within the same location, indicating that shorter elapsed times strengthen the influence of the current fixation. As elapsed time increases (up to Q4), location transitions become more evenly distributed, reflecting a reduced fixation influence. This trend is most noticeable in the final matrix, where transition probabilities approximate a uniform distribution for very long elapsed times (>3,000 ms). These findings highlight the importance of considering elapsed time when analyzing fixation transitions.

Figure 5

Five transition patterns of consecutive fixations for test block 1 sequences.

Figure 5 Long description

A multi-panel figure containing five heatmaps arranged in two rows. The top row contains three heatmaps labeled Q 1, Q 2, and Q 3. The bottom row contains two heatmaps labeled Q 4 and greater than 3000. A vertical color scale on the far right indicates probability from 0.00 in white to 1.00 in dark red.

Each heatmap shares the same axes. The vertical y-axis is labeled From and the horizontal x-axis is labeled To. Both axes contain six categories in the following order from bottom to top and left to right: look, look underscore 1, look underscore 2, option underscore correct, option underscore wrong, and question.

* Transition Heatmap Q 1: Shows high concentration along the diagonal and top-right corner. Key values include look underscore 1 to look underscore 1 at 0.94, look underscore 2 to look underscore 2 at 0.93, and question to question at 0.93.

* Transition Heatmap Q 2: Shows a shift toward the center and right. High probabilities include look underscore 1 to look underscore 1 at 0.58 and option underscore wrong to option underscore wrong at 0.66.

* Transition Heatmap Q 3: Similar to Q 2, with look underscore 1 to look underscore 1 at 0.57 and option underscore wrong to option underscore wrong at 0.63.

* Transition Heatmap Q 4: Shows more distributed probabilities. Highest values are look underscore 2 to look underscore 1 at 0.56 and question to question at 0.57.

* Transition Heatmap greater than 3000: Displays the most dispersed probability distribution. The highest single value is option underscore wrong to option underscore wrong at 0.49, followed by look underscore 1 to look underscore 1 at 0.44.

3.2.1 Input

The categorical actions in $\mathbf {s}_i$ pose a challenge for neural networks due to the absence of inherent numerical relationships. To overcome this, we convert each action into a numerical embedding vector of size k using a learnable embedding layer (Bengio et al., Reference Bengio, Ducharme, Vincent and Jauvin2003; Hancock & Khoshgoftaar, Reference Hancock and Khoshgoftaar2020; Mikolov, Reference Mikolov2013). Mathematically, the embedding layer is parameterized by a learnable matrix $\mathbf {W}_e \in \mathbb {R}^{6 \times k}$ , where $6$ is the number of unique actions and k is the embedding size. For each action $s \in \mathcal {A}$ , the embedding vector $\mathbf {e} \in \mathbb {R}^k$ :

(1) $$ \begin{align}\mathbf{e} = \mathbf{W}_e^{T} \mathbf{u}_s, \end{align} $$

where $\mathbf {u}_s \in \mathbb {R}^6$ is a one-hot encoded vector representing the action s, such that the s-th entry is 1 and all others are 0. For each sequence $\mathbf {s}_i$ , the embedding is represented as a matrix $\mathbf {E}_i=(\mathbf {e}_{i,1}, \dots , \mathbf {e}_{i,T_i})^T$ with shape $T_i\times k$ , where $\mathbf {e}_{i,t}=(e_{i,t,1}\dots ,e_{i,t,k})^T$ . By embedding the actions, we map the discrete categories into a continuous vector space, enabling the neural network to better capture semantic relationships between actions.

To capture how time influences these behaviors, we include not only the action sequence $\mathbf {s}_i$ but also the corresponding time duration $\mathbf {d}_i$ and elapsed time $\mathbf {l}_i$ , that is, feeding all of them into the Autoencoder model for feature extraction. Mathematically, we concatenate $\mathbf {E}_i, \mathbf {l}_i$ , and $\mathbf {d}_i$ along the column dimension to form $\mathbf {X}_i=\left [\mathbf {E}_i, \mathbf {l}_i, \mathbf {d}_i\right ]$ , which serves as the input to the Autoencoder.

3.2.2 Encoder and decoder with T-LSTM framework

The Autoencoder model includes two primary parts: an encoder and a decoder. The encoder compresses the input $\mathbf {X}_i$ into a latent representation $\mathbf {z}_i$ , formulated as

(2) $$ \begin{align}\mathbf{z}_i=f_{\textit{encoder}}(\mathbf{X}_i), \end{align} $$

where $\mathbf {z}_i\in \mathbb {R}^h$ is intended to capture high-level features of the input sequence i with dimension h. That is, $\mathbf {z}_i$ represents the key feature of the original input $\mathbf {X}_i$ and serves as the compressed feature representation used for subsequent analysis. Ideally, if $\mathbf {z}_i$ fully captures the information hidden in $\mathbf {X}_i$ , then we should be able to reconstruct $\mathbf {X}_i$ based solely on $\mathbf {z}_i$ . Thus, the decoder reconstructs the original $\mathbf {X}_i$ from $\mathbf {z}_i$ through

(3) $$ \begin{align}\hat{\mathbf{O}}_i = f_{\textit{decoder}}(\mathbf{z}_i), \end{align} $$

confirming that $\mathbf {z}_i$ effectively encodes the relevant information. In the Autoencoder model, the encoder block $f_{\textit {encoder}}$ and decoder block $f_{\textit {decoder}}$ play a critical role in effectively extracting hidden features.

To effectively capture the underlying structure in the high-dimensional sequential eye-tracking data, we design $f_{\textit {encoder}}$ and $f_{\textit {decoder}}$ as extensions of the long short-term memory (LSTM) architecture (Hihi & Bengio, Reference Hihi and Bengio1995; Hochreiter, Reference Hochreiter and Schmidhuber1997), referred to as T-LSTM, which incorporates duration and elapsed time information into the fixation action sequences. The standard LSTM framework is widely used in fields like patient subtyping, anomaly detection, and speech recognition (Baytas et al., Reference Baytas, Xiao, Zhang, Wang, Jain and Zhou2017; Irie et al., Reference Irie, Tüske, Alkhouli, Schlüter and Ney2016; Lindemann et al., Reference Lindemann, Maschler, Sahlab and Weyrich2021; Yu et al., Reference Yu, Si, Hu and Zhang2019), effectively capturing sequential dependencies, preserving long-term memory, and mitigating the vanishing gradient problem (Hihi & Bengio, Reference Hihi and Bengio1995; Pascanu, Reference Pascanu, Mikolov and Bengio2013).

LSTMs consist of a sequence of memory blocks, each containing a cell state and three gates: the input gate, forget gate, and output gate. The cell state carries the information of the current memory block, while the gates regulate the flow of information into and out of the cell state. The forget gate decides which information to discard, allowing the LSTM to forget irrelevant data. The input gate controls the addition of new information to the current cell state, enabling the network to update its memory with new inputs. Lastly, the output gate determines which information to pass to the output, based on the current contents of the cell state and the network’s learned parameters. This architecture enables LSTMs to remember important information over long sequences and to discard irrelevant data.

In the implementation of an LSTM memory block, the input gate $i_t$ , forget gate $f_t$ , candidate value of the memory cell $\tilde {c}_t$ , and output gate $o_t$ at time t are computed with Equations 4–7, respectively,

(4) $$ \begin{align} i_t&=\sigma\left(W_i x_t+U_i h_{t-1}+b_i\right), \end{align} $$

(5) $$ \begin{align} f_t=\sigma\left(W_f x_t+U_f h_{t-1}+b_f\right), \end{align} $$

(6) $$ \begin{align} \widetilde{c_t}=\tanh \left(W_c x_t+U_c h_{t-1}+b_c\right), \end{align} $$

(7) $$ \begin{align} o_t=\sigma\left(W_o x_t+U_o h_{t-1}+b_o\right), \end{align} $$

where W and U are weight matrices, b is the bias vector of each unit, and $\sigma $ and $tanh$ are the logistic sigmoid and hyperbolic tangent activation functions, respectively.

The current memory cell’s state $c_t$ is calculated by modulating the current memory candidate value $\tilde {c}_t$ via the input gate $i_t$ and the previous memory cell state $c_{t-1}$ via the forget gate $f_t$ , and $c_t=i_t \tilde {c}_t+f_t c_{t-1}$ . Through this process, a memory cell decides whether to keep or forget the previous memory state and regulates the candidate of the current memory state via the input gate. The memory cell state $c_t$ is then controlled by the output gate $o_t$ to compute the memory cell output $h_t$ of the LSTM block at time t,

(8) $$ \begin{align}h_t=o_t \tanh \left(c_t\right). \end{align} $$

To incorporate the crucial temporal information, we develop T-LSTM by employing the elapsed time ${\mathbf {l}}_i$ to weight the short-term memory built on the work from Baytas et al. (Reference Baytas, Xiao, Zhang, Wang, Jain and Zhou2017). The observations in Section 2.1 show that the passage of time reduces the influence of a previous action on the current one, thereby diminishing how past memory affects the current output. To address this, we propose adding a weight, determined by the elapsed time and pre-defined function $g(\cdot )$ , to adjust the previous memory cell state $c_{t-1}$ in the LSTM framework. Specifically, the adjusted previous memory cell state $c_{t-1}^{\ast}$ is calculated by

(9) $$ \begin{align}\begin{aligned} c_{t-1}^S &= \sigma (W_d c_{t-1}+ b_d)c_{t-1} , \\ \hat{c}_{t-1}^S &= c_{t-1}^S * g(l_t) , \\ c_{t-1}^T &= c_{t-1}-c_{t-1}^S , \\ c_{t-1}^{\ast} &= c_{t-1}^T+\hat{c}_{t-1}^S. \end{aligned} \end{align} $$

Accordingly, the current memory cell’s state $c_t$ will be

(10) $$ \begin{align} c_t=i_t \tilde{c}_t+f_t c_{t-1}^{\ast}. \end{align} $$

An illustration of the T-LSTM memory block is provided in Figure 8. Our T-LSTM autoencoder stacks these blocks to form an encoder–decoder. Given the input sequence $\mathbf {X}_i=\left [\mathbf {E}_i, \mathbf {l}_i, \mathbf {d}_i\right ]$ , the encoder processes it in time order and produces the final hidden states $h_{T_i}$ . We take this last hidden state as the sequence representation, $\mathbf {z}_i:=h_{T_i}$ in Equation (2), summarizing the full fixation trajectory together with its timing (durations and elapsed times). The decoder is initialized with the encoder’s terminal states $\left (h_{T_i}, c_{T_i}\right )$ . To ensure that the representation alone carries the information needed for reconstruction, the decoder receives a sequence of zero vectors as input; at each step, it updates its state solely from its internal memory. This “zero-input” design creates a tight bottleneck: if the decoder accurately reconstructs actions and times, then the encoder has indeed captured the essential temporal–sequential structure in $\mathbf {z}_i$ .

Figure 8

An illustration of a T-LSTM memory block.

Figure 8 Long description

The diagram illustrates the internal architecture of a Time-Aware Long Short-Term Memory cell.

* Input signals enter from the left: c sub t minus 1 at the top, h sub t minus 1 in the middle, and x sub t at the bottom.

* The top path processes the cell state. c sub t minus 1 passes through a sigma activation block. This output is multiplied by a signal from l sub t which passes through a g function.

* A subtraction operation occurs between the sigma output and the original c sub t minus 1 signal, which then feeds into an addition gate.

* The central horizontal line represents the main cell state flow. It includes a multiplication gate controlled by f sub t, followed by an addition gate.

* The bottom path processes h sub t minus 1 and x sub t. This combined signal branches into four parallel paths:

1. A sigma block yielding f sub t.

2. A sigma block yielding i sub t.

3. A tanh block yielding c tilde sub t.

4. A sigma block yielding o sub t.

* i sub t and c tilde sub t are multiplied and then added to the main cell state.

* The final cell state c sub t exits at the top right.

* The final hidden state h sub t exits at the bottom right, formed by multiplying o sub t with the tanh of the current cell state.

3.2.3 Output

Each decoder output $\hat {\mathbf {o}}_{i, t}$ is mapped to class probabilities over the six actions by a single-layer perceptron (Almeida, Reference Almeida2020; Rosenblatt, Reference Rosenblatt1958) followed by softmax:

(11) $$ \begin{align} \mathbf{p}_{i,t} = \operatorname{softmax} (\mathbf{W}_0 \hat{\mathbf{o}}_{i,t} + \mathbf{b}_0) \in \mathbb{R}^{6}, \end{align} $$

where $\mathbf {W}_0$ and $\mathbf {b}_0$ are learnable parameters of the one-layer perception. The resulting probability vector $\mathbf {p}_{i,t}$ contains six elements $(p_{i,t,1}, \ldots , p_{i,t,6})^T$ , each representing the likelihood of the corresponding action class. The predicted action for timestep t in sequence i is, $\hat {s}_{i,t} = \operatorname {argmax} \mathbf {p}_{i,t}$ , that is, the action class with the highest probability. For simplicity in notation and implementation, the six distinct actions in action set $\mathcal {A}$ are represented by integers from 1 to 6. We have the reconstructed fixation sequences $\hat {\mathcal {S}}=\{\hat {{\mathbf {s}}}_i\}_i^N$ with $\hat {{\mathbf {s}}}_i=(\hat {s}_{i,t}, \ldots , \hat {s}_{i, T_i})^T$ .

Alongside action reconstruction, we also predict the time sequences $\mathbf {l}_i$ and $\mathbf {d}_i$ . Two separate regression layers operate on the decoder outputs $\hat {\mathbf {O}}_i$ to estimate $\hat {l}_{i, t}$ and $\hat {d}_{i, t}$ for each time t:

(12) $$ \begin{align}\hat{l}_{i, t}=\operatorname{ReLU}\left(\mathbf{w}_1^{\top} \hat{\mathbf{o}}_{i,t}+b_1\right), \hat{d}_{i,t}=\operatorname{ReLU}\left(\mathbf{w}_2^{\top} \hat{\mathbf{o}}_{i,t}+b_2\right), \end{align} $$

where $\mathbf {w}_1$ and $\mathbf {w}_2$ are learnable weight vectors, and $b_1$ and $b_2$ are scalar biases.

Thus, at each step, the decoder yields (1) action probabilities, (2) a non-negative elapsed time, and (3) a non-negative duration, enabling joint reconstruction of sequence content and timing from $\mathbf {z}_i$ . The accurate reconstructions of these three components from the extracted feature/representation $h_T$ can validate that we successfully extract the essential information from the original sequences.

To achieve the accurate reconstruction, we optimize the model’s learning parameters using the traditional reconstruction loss, that is, the discrepancy between the original and reconstructed data through gradient descent (Ruder, Reference Ruder2016).

The learnable parameters include $\mathbf W_e^T$ in the embedding layer, the parameters in the encoder block $f_{\textit {encoder}}$ and decoder layer $f_{\textit {decoder}}$ that we will discuss in the following section, as well as $\mathbf W_0$ , $\mathbf b_0$ , $\mathbf w_1$ , $\mathbf w_2$ , $b_1$ , and $b_2$ in the classifier and regressors added after the decoder layer. Collectively, these parameters are denoted by $\boldsymbol {\theta }$ .

For the fixation action sequences, the discrepancy between $\mathcal {S}$ and the reconstructed $\hat {\mathcal {S}}$ is measured using the cross-entropy loss

(13) $$ \begin{align}L_{\mathrm{CE}}(\boldsymbol{\theta} \mid \mathcal{S})=-\frac{1}{N} \sum_{i=1}^N \frac{1}{T_i} \sum_{l=1}^{T_i} \sum_{k=1}^6 \mathbb{I}\left(s_{i,t}=k\right) \log \left(p_{i, t, k}\right), \end{align} $$

where $\mathbb {I}(\cdot )$ is the indicator function. This loss function is a widely used criterion for measuring the difference between the true action and the predicted probability. It ensures that the model is penalized more heavily when it assigns a low probability to the correct action. Considering the reconstruction performance on temporal information, we apply the widely used mean squared error (MSE) loss, defined as follows:

(14) $$ \begin{align}L_{\mathrm{MSE}}(\boldsymbol{\theta} \mid \mathcal{T}, \mathcal{D}) = \frac{1}{N} \sum_{i=1}^N \frac{1}{T_i} \sum_{t=1}^{T_i}\left[\left(l_{i, t}-\hat{l}_{i, t}\right)^2+\left(d_{i,t}-\hat{d}_{i,t}\right)^2\right]. \end{align} $$

We combine the cross-entropy loss $L_{CE}$ for action prediction and the MSE loss $L_{MSE}$ for time sequence reconstruction into a unified objective function:

(15) $$ \begin{align} L=L_{\mathrm{CE}}+\lambda L_{\mathrm{MSE}}, \end{align} $$

where $\lambda $ is a hyperparameter balancing the contributions of two parts. By minimizing this combined loss, the Autoencoder effectively learns to extract meaningful hidden representations from the input, enabling accurate reconstruction of both the action sequence and its associated temporal information from the hidden space.

3.3.1 Selection of $g(\cdot )$

Function $g(\cdot )$ essentially captures how elapsed time affects subsequent actions, which can be proposed based on evidence from a data set. As shown in Figure 5, longer elapsed times correspond to a diminished influence of the current fixation on the next one. Thus, a decreasing function is a natural choice in our case. We observed from Figure 5 that as elapsed time increases, the transition pattern changes rapidly at first but stabilizes over longer intervals. For instance, when the elapsed time is around $800$ ms, the transition pattern closely resembles that at very long intervals (e.g., $>3,000$ ms). This empirical pattern supports the use of an exponential decay function, since it decays quickly at the beginning and then gradually levels off, consistent with the observed transition dynamics. The exponential form is also mathematically simple, differentiable, and widely used in modeling memory decay and time-dependent influence in sequential data, making it both interpretable and computationally efficient. Consequently, we set $g(l)=\exp (-\alpha l)$ , with $\alpha =0.003$ , a value carefully selected to match the empirical transition patterns observed in the data (see Figure 9(a)).

Figure 9

(a) The proposed data-driven function. (b)–(d) Three potential alternative functions.

Figure 9 Long description

The figure consists of four panels arranged in a two-by-two grid. All panels share a horizontal x-axis labeled l ranging from 0 to 900 and a vertical y-axis labeled g of l.

* Top-left panel labeled (a) Function 1: The y-axis ranges from 0.0 to 1.0. The function is g of l equals exp(-alpha l). The data shows a decaying exponential curve starting at 1.0 when l is 0 and asymptotically approaching 0 as l increases.

* Top-right panel labeled (b) Function 2: The y-axis ranges from 0.00 to 2.00. The data shows a constant horizontal line at g of l equals 1.00 across the entire range of l.

* Bottom-left panel labeled (a) Function 3: The y-axis ranges from 0.0 to 3.5. The function is g of l equals a times (x minus b) squared plus c. The data shows a parabolic curve with a minimum point near l equals 400, where g of l is approximately 0.1, rising to over 3.0 at l equals 900.

* Bottom-right panel labeled (b) Function 4: The y-axis ranges from 0.0 to 1.8. The function is g of l equals exp(alpha l) plus a. The data shows an increasing exponential curve starting near 0.05 and rising to approximately 1.5 at l equals 900.

For comparison, three additional functions are presented in Figure 9(b)–(d), illustrating common possible patterns in how elapsed time impacts subsequent actions. Specifically, the constant function $g(l)=1$ in Figure 9(b) represents the standard LSTM, where elapsed time has no effect on subsequent actions. In contrast, the other two functions exhibit different trends: one shows an initial decrease in influence followed by an increase, while the other demonstrates a consistently increasing trend.

3.3.2 Weighting parameter $\lambda $

To determine the weighting parameter $\lambda $ in Equation (15), which balances the action and time sequence contributions, we ensure that their losses are on a similar scale during early training. This prevents one part from dominating. In our dataset, $L_{CE}$ is on the scale of 1, while $L_{MSE}$ is on the scale of 0.1. Since our primary goal is to reconstruct the action sequences, with time information as supplementary, we set the maximum $\lambda $ to 10 to align $L_{CE}$ and $\lambda L_{MSE}$ on the same scale. We also test values of 0.1 and 1 to emphasize $L_{CE}$ .

3.3.3 Embedding size k

Theoretical and empirical work suggests that for a large number of unique actions, the optimal embedding size is proportional to the number of actions, with a small scaling factor, such as 0.2 (Guo & Berkhahn, Reference Guo and Berkhahn2016; Yin & Shen, Reference Yin and Shen2018). In our case, with only six unique actions, we consider embedding sizes from the pool ${5,10,15}$ , capping at 15 to prevent over-fitting. A smaller size of 5 is suitable for simple action sequences, while 15 offers greater capacity to capture complex patterns.

3.3.4 Hidden size h

The selection of the hidden size, that is, the dimension of the extracted features, is closely related to the nature of the input data. A smaller hidden size may result in the loss of important information, while a larger one risks the model learning an identity function, thereby failing to extract meaningful features. Given that the average sequence length of our data is around $400$ , and the belief that students’ problem-solving behaviors can be effectively represented by features of relatively low dimensionality, we set the largest hidden size to be $20$ , which is $5\%$ of the average sequence length. The possible pool for the hidden size is defined as $\{5,10,20\}$ , which allows for experimentation to determine the optimal balance between compression and reconstruction fidelity. For more complex datasets or sequences containing richer information, we recommend considering higher dimensions, such as 30 or 50. However, the hidden size should generally remain modest to ensure effective information compression and to help subsequent analyses.

3.3.5 Evaluation

To evaluate the Autoencoder model’s performance with a set of hyperparameters, we focus on its ability to accurately reconstruct the fixation action sequences $\mathcal {S}$ , emphasizing the reconstruction of action sequences using temporal information. Performance is assessed based on three criteria. The first is Accuracy, which measures the proportion of correctly reconstructed actions compared to the original input sequence and is defined as

(16) $$ \begin{align}\text{ Accuracy }=\frac{1}{N} \sum_{i=1}^N \frac{1}{L_i} \sum_{l=1}^{L_i} \mathbb{I} (s_{i, l}=\hat{s}_{i, l} ). \end{align} $$

The second is BLEU (Bilingual Evaluation Understudy Score), which evaluates the similarity between the original and reconstructed sequences based on overlapping n-grams. For each sequence pair $(\mathbf {s}_i, \hat {\mathbf {s}}_i)$ , the BLEU score is computed as

(17) $$ \begin{align}\text{BLEU}=\exp ( \frac{1}{M} \sum_{m=1}^M \log p_m), \end{align} $$

where $p_m=\frac {\text {Number of overlapping m-grams between } \mathbf {s}_i \ \text {and } \hat {\mathbf {s}}_i}{\text { Total number of m-grams in } \hat {\mathbf {s}}_i }$ is the precision of m-grams of length m. In this study, we set $m=4$ and compute the final BLEU score by averaging the BLEU scores across all sequences.

The last is ROUGE (Recall-Oriented Understudy for Gisting Evaluation), which evaluates the quality of the reconstructed sequence by comparing m-gram recall and the longest common subsequence (LCS). For m-grams, the ROUGE score for a sequence pair $\left (\mathbf {s}_i, \hat {\mathbf {s}}_i\right )$ is defined as

(18) $$ \begin{align} \text{ ROUGE-M }=\frac{\text {Number of overlapping m-grams between } \mathbf{s}_i \ \text{and } \hat{\mathbf{s}}_i }{\text{ Total number of m-grams in } \mathbf{s}_i }. \end{align} $$

For the LCS, ROUGE-L is defined as

(19) $$ \begin{align} \text { ROUGE-L }=\frac{\text { LCS length }}{\text { Length of } \mathbf{s}_i}. \end{align} $$

In this work, we calculate the ROUGE score for each sequence by $\text {ROUGE} = \frac {\text {ROUGE-1{+}ROUGE-2{+}ROUGE-L}}{3}$ , and average it over all sequences. These three metrics provide a comprehensive evaluation of the Autoencoder’s performance by assessing both token-level reconstruction (accuracy) and sequence-level quality (BLEU and ROUGE). Higher values indicate better model performance.

4.2.1 Descriptive statistics and visualization

To address research questions RQ1 and RQ4, we calculate descriptive statistics for each cluster, including the number of sequences, average response accuracy, and the distribution of fixation locations, presented via pie charts. We also adjust the transition matrix for each cluster to highlight deviations in transition probabilities from the overall pattern, aiding our interpretation of testing and learning behaviors based on clusters from the testing and learning blocks.

4.2.2 Contingency table analysis

To address research questions RQ2, RQ5, and RQ6, we perform contingency table analysis between testing clusters across test blocks 1 and 2, and across different item types in testing blocks (RQ2), as well as between learning clusters across the two learning blocks and different item types (RQ5), and between learning and testing clusters RQ6. We use mosaic plots, where the area of each rectangular tile represents the conditional relative frequency for a cell in the contingency table, to summarize the patterns of association. Cramer’s V is used to calculate the association between testing clusters and item types, and between learning clusters and item types. Additionally, $\chi ^2$ tests of independence, using likelihood ratio statistics ( $G^2$ ), assess the significance of associations between the two sets of clusters.

4.2.3 Generalized mixed effect model

To explore the relationship between testing clusters and the test performance (RQ3), we fit the generalized linear mix-effect model (GLMM) (Breslow & Clayton, Reference Breslow and Clayton1993) where the outcome variable is the binary score to each question (1 indicating correct and 0 otherwise), the explanatory variables include the testing cluster ID, the type of items, and the block identifier. Additionally, we include a person-level random effect to account for individual-level variability, capturing the within-person correlation in test performance. Note that we reduce the item types to levels of 2, indicating an item measures one direction rotation or two direction rotations to simplify the interpretation. The GLMM is conducted via $lme4$ R package (Bates, Reference Bates, Mächler, Bolker and Walker2015).