## Pixel-Wise Modelling and non-parametric statistics

### Pixel-Wise Modelling

Although the generation mechanism of eye movement data is still largely under debate, recent theories and applications suggest that a spatial model is the most appropriate to consider the statistical analysis of fixation especially its location distribution. For example, Barthelmé, Trukenbrod, Engbert, & Wichmann, (2013) recommend using the point process framework to inference how fixations are distributed in space. While we endorse this fruitful approach and its Bayesian nature, here we aim to resolve this problem from an opposite perspective. Instead of inferring from the spatial distribution of the fixation, we infer on each location in the search space (i.e., each pixel of within the eye tracker recordable range or each pixel in the visible stimuli). In other words, we try to answer the question:

*“How long is this pixel being fixated (or what is the probability of this pixel to be fixated) in the function of the experimental conditions?”*.

Formally, by applying mixed models independently on each pixel, we have:

=Y(s)+Xβ(s)+Zb(s)ε(s)

for s ∈Dof the search space

The complete procedure as implemented in *i*Map4 is explained in the figure below. Eye movement data for each participant is concatenated into one input data matrix. *i*Map4 first partition the data matrix into a fixation characteristic matrix *(red box)* and an experiment condition information matrix *(green box)*. The fixation characteristic matrix contains fixation spatial location (*x* and *y*), fixation duration, and order index of each fixation. The experiment condition matrix contains the index of each subject, index of each trial/item, and different levels of each experimental condition. Fixation durations are then projected into the two-dimensional space according to their *x* and *y* coordinates at the single-trial level. *i*Map4 then smooths the fixation duration map by convoluting it with a two-dimension Gaussian Kernel function:

Kernel~Ν(0,σ^{2}Ι), where I is a two by two identity matrix and the full width at half maximum (FWHM) of the Kernel is1° visual angleas the default setting.

This step is essential to account for the spatial uncertainty of eye movement recording (both mechanical and physiological) and the sparseness of the fixation locations. The Gaussian kernel could also be replaced by other 2D spatial filters to best suit the research question.

*Illustration of the procedure in iMap4. The input data matrix is partitioned into eye movement matrix and predictor matrix. Fixation durations are projected into the two-dimensional space according to their x and y coordinates at the single trial level for each participant. The experimental information of each trial is also summarised in a predictor table. Subsequently, the sparse representation of the fixation duration map is smoothed by convoluting it with a two dimensions Gaussian Kernel function Kernel ~ Ν (0,σ^{2} Ι). After estimating the fixation bias of each condition independently for all the observers (by taking the expected values across trial within the same condition), iMap4 models the 3D smoothed fixation map (item × xSize × ySize) independently for each pixel using a LMM. The result is saved as a Matlab structure in LMMmap. iMap4 offers many parametric and non-parametric methods for hypothesis testing and multiple comparison corrections.*

The resulting smoothed fixation map is a 3D matrix. The last two dimensions of the fixation matrix are the size of the stimuli/search space. Information of each entry in the first dimension is stored in a predictor table, which is generated from the experiment condition matrix. Each experiment condition can be coded at the single trial level in the predictor table, or as one entry by taking the average map across trials.

In addition, *i*Map4 provides robust estimation option by applying Winsorization to limit extreme values in the smoothed fixation matrix. The goal here is to reduce the effect of potential outliers. Additional options include spatial normalisation (z-scored map or probability map), spatial downsampling (linear transformation using *imresize* in Matlab) to optimise computing speed, and mask creation to exclude irrelevant pixels.

The resulting 3D fixation matrix is then modelled in a LMM as the response variable. The results are saved as a Matlab structure (LMMmap as in the examples below). The fields of LMMmap are nearly identical to the output from *LinearMixedModel* class. For each modelled pixel, *i*Map4 saves the model criterion, variances explained, error sum of squares, coefficient estimates and their covariance matrix for both fixed and random effects, and the ANOVA results on the LMM. Additional modelling specifications, as well as other model parameters including LMM formula, design matrix for fixed and random effect, and residual degrees of freedom, are also saved in LMMmap. Linear contrasts and other analyses based on variance or covariance can be performed afterwards from the model fitting information. Any other computation on the *LinearMixedModel* output can also be replicated on LMMmap.

### Non-parametric statistics using permutation and bootstrap spatial clustering

One of the crucial assumptions of pixel-wise modelling is that all pixels are independent and identically distributed. Of course, this assumption is never satisfied neither before nor after smoothing. To ensure valid inferences on activity patterns in large 2D pixel space, we applied non-parametric statistics to resolve the biases in parameter estimation and problems arising from multiple comparisons. We developed two resampling-based statistical hypothesis testing methods for the fixed effect coefficients: a universal permutation test and a universal bootstrap clustering test.

The resampling tests on the model coefficient for fixed effects * β* are operated on the fixed effect related variances. To do so, we simply removed the variance associated with the random effects from the response matrix:

=Y(s)_{fixed}+Xβ(s)=ε(s)-Y(s)Zb(s)

for s ∈Dof the search space

For any permutation test, *i*Map4 performs the following algorithms on * Y_{fixed}* for each pixel:

## Algorithm 1:

For a given hypothesis or linear contrast

,ciMap4

- Performs a linear transformation on the design matrix
to get a new design matrixXso that the partitioning ofM= [M,M1]. ThenM2iMap4 computes the new coefficients by projectingto the pseudoinverse ofY_{fixed}. The design matrixMis created so that the original hypothesis testing is equivalent to the hypothesis regardingMcoefficients. The matrix transformation and partition are the same as the algorithm described in Winkler et al (2014, appendix A)M1- Computes the residuals related to the hypothesis by subtracting the variance accounted by
fromM2to getY_{fixed}Y_{rr}- Fits
toY_{rr}by solvingM=Y_{rr}+Mβ_{m}, and get the statistics valueεofF_{rr}. Note that to replicate the original hypothesis testing on the fixed effect, the new contrastM1is just to partitionc1intoMandM1M2- Permutes the rows of the design matrix
to obtain the new design matrixMM^{ * }- Fits
toY_{rr}and gets theM^{ * }ofF_{rr}^{ * }M1^{ * }- Repeats step (4) and (5) for a large number of times (
k resamplings/repetitions), thep-value is then defined asp=(#≥F_{rr}^{ * }) /F_{rr}k. Importantly, the FWER correctedp-value is computed by comparing the largestacross all tested pixels in one resampling with the originalF_{rr}^{ * }.F_{rr}

**Algorithm 1** is a simplified version of Winkler et al (2014, Algorithm 1): the resampling table includes permutation but not sign-flipping. Sign-flipping assumes the errors to be independent and symmetric. Thus, the underlying assumptions are stronger than with classical permutations, which require only exchangeable errors (Winkler et al, 2014).

Importantly, this test is exact only under a balanced design with no missing value and only subject as the random effect. As previously shown in Kherad-Pajouh and Renaud (2014), a general and exact permutation approach for mixed-model designs should be performed on modified residuals that have up to second-moment exchangeability. This is done to satisfy the important assumptions for repeated measures ANOVA: normality and sphericity of the error. However, there are strict requirements to achieve this goal: careful transformation and partition of both fixed and random effects design matrices, and removal of the random effects related to * M2* (Kherad-Pajouh and Renaud, 2014). In

*i*Map4, we perform an approximation version by removing all random effects to increase the efficiency and the speed of the huge amount of resampling computation in our pixel-wise modelling algorithm. Validation and simulation data set indeed showed that the sensitivity and the false alarm rate of the proposed algorithm are not compromised.

*i*Map4 performs the following algorithm on * Y_{fixed}* for each pixel as the bootstrap clustering approach:

## Algorithm 2:

- For each unique categorical variable,
iMap4 removes the conditional expectations fromfor each pixel. A random shuffling is then performed on the centered data to acquireY_{fixed}, so that any potential covariance is also disrupted. This is done to construct the true empirical null hypothesis distribution in which all elements and their linear combinations inY_{c}have expected values equal to zero.Y_{c}- Randomly draws with replacement from {
,X,Z} an equal number of subjects {Y_{c},X^{ *},Z^{ *}}Y_{c}^{ *}- Fits
toY_{c}^{ *}by solvingX^{ *}=Y_{c}^{ *}X^{ *}+β^{ *}. For a given hypothesis or linear contrastε,ciMap4 computes the statistics valueand its parametricF^{ *}p-value under the GLM framework.- Thresholds statistical maps
atF^{ *}p≤ .05 and records the desired maximum cluster characteristics across all significant clusters. Cluster characteristics considered are: cluster mass (summed^{*}Fvalue within a significant cluster), cluster extent (size of the cluster), or cluster density (meanFvalue within clustering).- Repeats step (2), (3) and (4) a large number of times to get the cluster characteristic distribution under the null hypothesis H0.
- Thresholds the original statistics map
atFp≤ .05 and compares the selected cluster characteristic with the value of the null distribution corresponding to the 95^{th}percentile. Any cluster with the chosen characteristic larger than this threshold is considered significant. The bootstrap clustering approach is identical to the bootstrap procedure described by Pernet et al. (2011; 2014) if only subject intercept is considered as the random effect. In addition,Algorithm 2extents the philosophy and approach presented by Pernet et al. (2011; 2014) to non-hierarchical mixed-effect models.

It is worth noting that we implemented in *i*Map4 a high-performance algorithm to minimise the computational demands for a large amount of resampling. Model fitting in both resampling approaches makes use of Ordinary Least Squares. The inversion of the covariance matrixes is computed on the upper triangular factor of the Cholesky decomposition. Calculation of the quartic form for all pixels is optimised by constructing a sparse matrix of the inversion of the covariance matrix. More details of these algebra simplifications could be found in the imapLMMresample function in *i*Map4.

Other multiple comparison correction methods such as Bonferroni correction, False Discovery Rate (FDR), or Random Field Theory (RFT) could also be applied. A Threshold-Free Cluster Enhancement (TFCE) algorithm could also be applied to the statistical (* F*- value) maps as an option after the permutation and bootstrap clustering procedure (Smith & Nichols, 2009).