Pixel Wise Modeling and non parametric statistics

Pixel Wise Modeling

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 ∈ D of the search space

The complete procedure as implemented in iMap4 is explained in the figure below. Eye movement data for each participant is concatenated into one input data matrix. iMap4 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 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. iMap4 then smooths the fixation duration map by convoluting it with a two-dimension Gaussian Kernel function:

Kernel ~ Ν (0 ,σ²Ι) , where I is a two by two identity matrix and the full width at half maximum (FWHM) of the Kernel is 1° visual angle as 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 summarized 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 ,σ²_Ι_). 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](https://github.com/iBMLab/iMap4/wiki/LMMmap). iMap4 offers many parametric and non-parametric methods for hypothesis testing and multiple comparison correction._

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, iMap4 provides robust estimation option by applying Winsorization in order to limit extreme values in the smoothed fixation matrix. The goal here is to reduce the effect of potentially outliers. Additional options include: spatial normalization (z-scored map or probability map), spatial down-sampling (linear transformation using imresize in Matlab) to optimize computing speed, and mask creation to exclude irrelevant pixels.

The resulting 3D fixation matrix is then modeled 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 modeled pixel, iMap4 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 modeling 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 afterward 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 modeling 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_fixed(s) = Xβ(s) + ε(s) = Y(s) - Zb(s)

for s ∈ D of the search space

For any permutation test, iMap4 performs the following algorithms on Y_fixed for each pixel:

Algorithm 1:

For a given hypothesis or linear contrast c, iMap4

1. Performs a linear transformation on the design matrix X to get a new design matrix M so that the partitioning of M = [ M1 , M2 ]. Then iMap4 computes the new coefficients by projecting Y_fixed to the pseudoinverse of M. The design matrix M is created so that the original hypothesis testing is equivalent to the hypothesis regarding M1 coefficients. The matrix transformation and partition are the same as the algorithm described in Winkler et al (2014, appendix A)
2. Computes the residuals related to the hypothesis by subtracting the variance accounted by M2 from Y_fixed to get Y_rr
3. Fits Y_rr to M by solving Y_rr = Mβ_m + ε, and get the statistics value F_rr of M1. Note that to replicate the original hypothesis testing on the fixed effect, the new contrast c1 is just to partition M into M1 and M2
4. Permutes the rows of the design matrix M to obtain the new design matrix M ^*
5. Fits Y_rr to M ^* and gets the F_rr ^* of M1 ^*
6. Repeats step (4) and (5) for a large number of times (k resamplings/repetitions), the p-value is then defined as p =(# F_rr ^* ≥ F_rr ) / k. Importantly, the FWER corrected p-value is computed by comparing the largest F_rr ^* across all tested pixels in one resampling with the original 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 balance 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 in order 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 iMap4, 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 modeling algorithm. Validation and simulation data set indeed showed that the sensitivity and the false alarm rate of the purposed algorithm are not compromised.

iMap4 performs the following algorithm on Y_fixed for each pixel as the bootstrap clustering approach:

Algorithm 2:

1. For each unique categorical variable, iMap4 removes the conditional expectations from Y_fixed for each pixel. A random shuffling is then performed on the centered data to acquire Y_c , 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 in Y_c have expected values equal to zero.
2. Randomly draws with replacement from { X , Z , Y_c } an equal number of subjects { X ^* , Z ^* , Y_c ^* }
3. Fits Y_c ^* to X ^* by solving Y_c ^* = X ^* β ^* + ε. For a given hypothesis or linear contrast c , iMap4 computes the statistics value F ^* and its parametric p-value under the GLM framework.
4. Thresholds statistical maps F ^* at p^* ≤ .05 and records the desired maximum cluster characteristics across all significant clusters. Cluster characteristics considered are: cluster mass (summed F value within a significant cluster), cluster extent (size of the cluster), or cluster density (mean F value within clustering).
5. Repeats step (2), (3) and (4) a large number of times to get the cluster characteristic distribution under the null hypothesis H0.
6. Thresholds the original statistics map F at p ≤ .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 2 extents 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 iMap4 a high-performance algorithm to minimize the computational demands for the 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 optimized by constructing a sparse matrix of the inversion of the covariance matrix. More details of these algribra simplifications could be found in the imapLMMresample function in iMap4.

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 on the statistical (F- value) maps as an option after the permutation and bootstrap clustering procedure (Smith & Nichols, 2009).