Behav Res (2017) 49:559–575
DOI 10.3758/s13428-016-0737-x

BRIEF COMMUNICATION

iMap4: An open source toolbox for the statistical fixation
mapping of eye movement data with linear mixed modeling
Junpeng Lao 1 & Sébastien Miellet 1,2 & Cyril Pernet 3 & Nayla Sokhn 1 & Roberto Caldara 1

Published online: 3 May 2016
# Psychonomic Society, Inc. 2016

Abstract A major challenge in modern eye movement research is to statistically map where observers are looking, by
isolating the significant differences between groups and conditions. As compared to the signals from contemporary neuroscience measures, such as magneto/electroencephalography
and functional magnetic resonance imaging, eye movement
data are sparser, with much larger variations in space acrosstrials and participants. As a result, the implementation of aconventional linear modeling approach on two-dimensional
fixation distributions often returns unstable estimations andunderpowered results, leaving this statistical problem unresolved (Liversedge, Gilchrist, & Everling, 2011). Here, we
present a new version of the iMap toolbox (Caldara &Miellet, 2011) that tackles this issue by implementing a statistical framework comparable to those developed in state-ofthe-art neuroimaging data-processing toolboxes. iMap4 uses
univariate, pixel-wise linear mixed models on smoothed fixation data, with the flexibility of coding for multiple betweenand within-subjects comparisons and performing all possible
linear contrasts for the fixed effects (main effects, interactions,etc.). Importantly, we also introduced novel nonparametrictests based on resampling, to assess statistical significance.

* Junpeng Laojunpeng.lao@unifr.ch* Roberto Caldara
roberto.caldara@unifr.ch
1

Department of Psychology, University of Fribourg, Faucigny 2,1700 Fribourg, Switzerland

2

Department of Psychology, Bournemouth University,
Bournemouth, UK

3Centre for Clinical Brain Sciences and Neuroimaging Sciences,
University of Edinburgh, Edinburgh, UK

Finally, we validated this approach by using both experimental and Monte Carlo simulation data. iMap4 is a freely available MATLAB open source toolbox for the statistical fixation
mapping of eye movement data, with a user-friendly interface
providing straightforward, easy-to-interpret statistical graphical outputs. iMap4 matches the standards of robust statisticalneuroimaging methods and represents an important step in thedata-driven processing of eye movement fixation data, an important field of vision sciences.
Keywords Eye movement analysis . Statistical mapping .
Linear mixed models
Human beings constantly move the eyes to sample visual information of interest from the environment. Eye
fixations deliver inputs with the highest resolution to thehuman visual cortex from the fovea, as well as blurry,low-spatial-frequency information from peripheral vision(Rayner, 1998). Thus, isolating statistically where andhow long fixations are deployed to process visual information is of particular interest to behavioral researchers,
psychologists, and neuroscientists. Moreover, fixation
mapping has a wide range of practical applications indetermining marketing strategies and the understandingof consumer behaviour (Duchowski, 2002).
Conventional eye movement data analyses rely on the estimation of probabilities of occurrence of fixations and saccades (or their characteristics, such as duration or length) within predefined regions of interest (ROIs), which are at best
defined a priori—but often also defined a posteriori, on thebasis of data exploration, which inflates the Type I error rate.
Another issue with ROIs is of course that other important
information not included in the ROI is discarded. In a continuous effort to circumvent the limitations of the ROI approach

560(for a detailed discussion on this point, see Caldara & Miellet,2011), we previously developed an unbiased, data-driven approach to compute statistical fixation maps of eye movements:
the iMap toolbox (Caldara & Miellet, 2011). From the very
first version, the toolbox was developed as a MATLAB open
source toolbox freely available for download online. The previous versions (1 and 2) made use of Gaussian smoothing andthe random field theory as a statistical engine (Caldara &Miellet, 2011), which is one of the standard methods appliedin statistical analyses for functional magnetic resonance imaging (fMRI) data (Penny, Friston, Ashburner, Kiebel, &Nichols, 2011). Version 3 introduced pixel-wise t test andbootstrap clustering in order to generate self-contained statistical maps (Miellet, Lao, & Caldara, 2014). However, all of
the previous versions of iMap still suffered from a major limitation: They could only contrast two conditions at a time.
A major revision of the toolbox was necessary to enable
the analysis of more complex experimental designs routinely
used in the field. One of the most suitable and obvious
statistical solutions to overcome this problem would be to
implement a general linear model, a widespread approach inboth behavioral and neural-imaging data analyses. In fact,many modern procedures for hypothesis testing, such as thet test, analysis of variance (ANOVA), regression, and soforth, belong to the family of general linear models.
However, eye movement data are a sparse production ofvisual perceptual sampling. Unlike neuroimaging data, eyemovement data contain many empty cells with little to nodata points across the tested space (e.g., all of the pixels in
an image). This caveat engenders a statistical problem when
the same statistical inference procedure is applied on eachpixel, regardless or whether or not its data are missing. To
account for the sparseness and the high variation of spatial
eye movement data, we developed a specific novel approachfor smoothed fixation maps, which was inspired by the statistical framework implemented in diverse state-of-the-art
neuroimaging data-processing toolboxes: statistical parametric mapping (SPM; Penny et al., 2011), Fieldtrip(Oostenveld, Fries, Maris, & Schoffelen, 2011), and LIMO
EEG (Pernet, Chauveau, Gaspar, & Rousselet, 2011). In thesimplest case, users can apply a massive univariate, pixelwise linear mixed model (LMM) on the smoothed fixation
data with the subject considered as a random effect, which
offers the flexibility to code for multiple between- andwithin-subjects comparisons. Our approach allows users to
perform all possible linear contrasts for the fixed effects(main effects, interactions, etc.) from the resulting model
coefficients and the estimated covariance. Importantly, we
also introduced a novel nonparametric statistical test basedon resampling (permutation and bootstrap spatial clustering)to assess the statistical significance of the linear contrasts(Pernet, Latinus, Nichols, & Rousselet, 2015; Winkler,
Ridgway, Webster, Smith, & Nichols, 2014).

Behav Res (2017) 49:559–575

In the next section, we briefly describe the key concepts of
the LMM approach. We then introduce the novel nonparametric statistical approach on the fixed effects that we implemented
in iMap4, which uses a resampling procedure and spatial clustering. We also report a validation of the proposed resampling
procedures, and illustrate how iMap4 can be used, with both asubset of data from a previous study and computer-simulated
data. Finally, we give an overview of future development anddiscuss technical insights on eye fixation mapping.

Linear mixed models
In this part, we outline the key elements and concepts of LMMs
in comparison with general linear models (GLM) andhierarchical linear models (HLM). Mixed models represent acomplex topic, and the discussion of many underlying
mathematical details goes beyond the scope of this article. For
general, thoughtful introductions to mixed models, users of thetoolbox should refer to Raudenbush and Bryk (2002) andMcCulloch, Searle, and Neuhaus (2011). Users may also wish
to consult the documentation and help files of the
LinearMixedModel class in the MATLAB Statistics Toolbox
for details about parameter estimation and the available methods(www.mathworks.com/help/stats/linearmixedmodel-class.html).
Statistical hypothesis testing methods that make use of theanalysis of variance (regression, t test, ANOVA, analysis ofcovariance, etc.) are the most popular methods of data analysis
in many fields of research. Commonly used in psychology andneuroimaging studies, these methods could all be written as
particular cases of GLM:
yi ¼ β1 x1i þ β2 x2i þ …þ
 β t xti þ εi ;εi ∼N 0; σ2 ;

ð1Þwhere yi is the ith experiment measure and β1, β2 ,…, βt are themodel coefficients. The error term εi is normally distributedwith mean 0 and variance σ2. Alternatively, the GLMs (Eq. 1)
could be expressed in matrix form:
Y ¼ Xβ þ ε; ε Ν 0; σ2 I ;eð2Þwhere matrix X = [x1, x2, … , xt] is the design matrix, and I isan n-by-n identity matrix (n being the total number of observation). Usually, one of the columns in X is 1, so that themodel includes a constant or intercept coefficient that represents the overall mean. It is worth noting that the design matrix
could be parameterized in a different way. In conventional
psychology or behavioral researches, a sigma-restricted parameterization is often applied. In a sigma-restricted design
matrix, X is full rank and invertible, and the degrees of freedom are equal to the number of columns. In comparison,many types of neuroimaging analysis software prefer a cell-Behav Res (2017) 49:559–575

561

mean model or an overparameterized design matrix at thesingle-subject level (Penny et al., 2011; Pernet et al., 2011).
Such software uses an overparameterized design matrix, andits solution to Eq. 2 is given by projecting the response vector
Y to the pseudo-inverse of the design matrix X. The form of
the design matrix is important, since it codes different experiment designs and the intended statistical testing. In iMap4,
the design matrix of the fixed effect can be a cell-mean model,
sigma-restricted model (for the Type III ANOVA), or the offset from a reference model (for Type I ANOVA).
^ could be found easily by orThe coefficient estimations β
dinary least squares or other, more robust methods. Finally,
statistical inferences on the model estimations could be
expressed in different forms, depending on the types of design
matrix. In the case of sigma-restricted parameterization, we canseparate the design matrix Xsr and the vector of parameters βsrinto two parts (Kherad-Pajouh & Renaud, 2010):
 
β1
Xsr ¼ ½X1 X2  ; β sr ¼β2where X1 and β1 are the components of interest, with the corresponding hypotheses:
H0 : β1 ¼ 0 versus H1 : β1 ≠ 0:

ð3Þ

Given the Gaussian distribution of the error ε and the existence of the inverse or general inverse of design matrix Xsr,
we can get the statistics for the F test by means of ANOVA
with the following equations (for simplicity, in Eqs. 4–8 wedenote X = Xsr):‐ð4ÞH ¼ X XT X XT
 T ‐ T 
Xresid ¼ I‐X2 X2 X2 X2 X1
ð5Þ
 T
‐ T
Hresid ¼ Xresid Xresid Xresid Xresid
ð6Þ
d f e ¼ Number of observations−rankðXÞ

ð7Þ

YT Hresid Y=rankðX1 Þ;
YT ðI − HÞ Y = ðd f eÞ

ð8Þ

F¼

where H represents the hat matrix of the linear model in Eq. 2;
it projects the response vector Y onto the column space of X.
Hresid is the hat matrix of the hypothesis in Eq. 3, and dfe is themodel’s degrees of freedom. F has a Fisher–Snedecor distribution F ðrank ðX1 Þ ; d f eÞ.
As a comparison, in an overparameterized design matrix orcell-mean model design matrix with design matrix Xcm and thevector of parameters βcm, the statistics of various effects are
performed by linear combinations of the coefficient βcm. For
example, the equivalent Hypothesis (3) could be expressed as:
H0 : c*βcm ¼ 0 versus H1 : c*βcm ≠ 0;

ð9Þ

where rank(c) = rank(X1) and c*βcm= β1 in the sigmarestricted parameterization model. The related F test is thengiven by the quartic form of the linear contrast matrix c andthe inverse of the covariant matrix of βcm (for simplicity, in
Eqs. 10 and 11 we denote X = Xcm):
MSE ¼ YT ðI − HÞ Y=ðd f eÞ

− −

ðc*βcm ÞT MSE*c XT X cT ðc*βcm Þ
F¼rankðcÞ

ð10Þ
;

ð11Þ

where H and dfe are computed using Eqs. 4 and 7, respectively. Moreover, it could be proved that Eqs. 8 and 11 are equivalent. The related details and mathematical proofs could befound in many textbooks (e.g., Christensen, 2011).
The GLM Y = Xβ + ε could be easily extended into a generalized form with ε ~ Ν(0, σ2V) where V is some known
positive definite matrix. Moreover, if a more specific structureof the error ε is available, the GLM (Eq. 2), which has one
random-effect term (the error ε), could be further extended intoa mixed model. Mixed models include additional randomeffect terms that can represent the clusters or classes. In a typical neuroimaging study, this could be the subjects or groups. In
the following example, we consider a simplified case in which
only the subject is considered as the additional random effect.
This type of model is one of the most widely used models inboth fMRI and electroencephalography (EEG).
As a demonstration, here we consider a random intercept
and slope model, with both the intercept (i.e., the overall meanof each subject) and the slope (i.e., the differences among
conditions within each subject) varying independently. This
type of HLM, or so-called two-level linear model, takes theform of an expansion of Eq. 1 into:
yi j ¼ β 1 j x1i j þ β 2 j x2i j þ … þ β t j xti j þ εi j ; εi j e Ν 0; σ2
β 1 j ¼ β 10 þ b1 j ; b1 j e Ν 0; σ21⋮
β t j ¼ β t0 þ bt j ;bt j Ν 0; σ2t ;eð12Þwhere j stands for the jth subject. After substituting the subjectlevel parameters in the first-level model, Eq. 12 becomes

562

Behav Res (2017) 49:559–575

If we express the subject-level predictor xtij in the random
effects by the term ztij, we get the LMM

which corresponds to the standard form of LMMs:
Y ¼ Xβ þ Zb þ εð13Þ

b ~ Ν(0, σ2D), ε ~ Ν(0, σ2I), b and ε are independent fromeach other,where σ2D is the covariance matrix for the random
effects. In the example here, D would be a j-by-j identity
matrix. An alternative form of Eq. 12, as applied in LIMO
EEG or SPM, can be found in Friston, Stephan, Lund,
Morcom, and Kiebel (2005, Eq. 1).
HLMs are specific cases of LMMs. In a mixed model,
factors are not necessarily hierarchical. Moreover, crossed factors between fixed effects and random effects are much easier
to model in mixed than in hierarchical models. In additions,the fixed effects and random effects are estimated simultaneously in mixed models, which is not always the case inhierarchical models.
Parameter estimation in mixed models is much more complicated than in GLM or HLM. Assuming that the model in
Eq. 13 has the error covariance matrix R: var ðYjbÞ ¼ R, thismodel is equivalent to yeN ðX β; V Þ ; V ¼ ZDZ T þ R. The
estimation of the fixed effects β requires prior knowledge of V,
which is usually unavailable. In practice, the variance compo^ based on one
nent V is commonly replaced by an estimation Vof several approaches, such as ANOVA, maximum likelihood(ML) estimation, restricted maximum likelihood (ReML) estimation, or Monte Carlo approximation (McCulloch, Searle,& Neuhaus, 2011; Pinheiro & Bates, 2000). In general, themodel-fitting procedure of LMM is implemented in major
statistical packages (e.g., R and Stata) by solving
Henderson’s mixed model equation. iMap4 calls the
MATLAB class LinearMixedModel from Statistics Toolbox(versions R2013b or above) to estimate the coefficients (fixed
effect β and random effect b) and the covariance matrix V
with various options (key concepts with regard to parameter
estimations can be found in the MATLAB documentation:www.mathworks.com/help/stats/estimating-parameters-inlinear-mixed-effects-models.html). In brief, model
coefficients are estimated by ML or ReML, and the patternof the covariance matrix of the random effects (D) could takethe form of a full covariance matrix, a diagonal covariancematrix, or other symmetry structure.
Statistical inferences in LMM are also much more complex
than in a GLM. In a balanced design, or with the variancecomponent V known, hypothesis testing of the fixed effectfollows Eq. 8 or 11 as an exact test. However, in an unbalanced design with random effects, no exact F statistics areavailable, since biases in the estimation usually result in an
unknown distribution of F (Kherad-Pajouh & Renaud, 2015).
Although F and t values are available as approximate tests inmost statistical packages, Baayen, Davidson, and Bates(2008) discouraged the usage of t or F statistics, and especiallyreport of the p value, in mixed models. Other approaches havealso been proposed. For example, likelihood ratio tests could
be performed to test composite hypotheses by comparing thedesired model with the reduced model. However, there aremany constraints on the application of likelihood ratio tests(e.g., the method of model fitting and selection of the reduced
model). Moreover, running multiple copies of similar LMMs
is computationally expensive, especially in the context of
pixel-wise testing, such as in iMap4.
Besides the practical problem of statistical inferences with
LMM, another main challenge in the application of LMM tospatial eye movement data is the Type I error from multiple
comparisons. To resolve these issues, we adopted resamplingtechniques for null-hypothesis statistical testing, as is suggested in neuroimaging analysis with GLM or HLM (Pernet
et al., 2015; Winkler et al., 2014). Nonparametric statistics
using Monte Carlo simulation are ideal for both parameter
estimation and hypothesis testing (Baayen et al., 2008;
Kherad-Pajouh & Renaud, 2015). In iMap4, we adapted asimplified version of the permutation test suggested by
Winkler et al. (2014) and a bootstrap clustering method similar to the one applied in LIMO EEG (Pernet et al., 2011).
Details of the proposed algorithm and preliminary validation
result are described in the following section.

Pixel-wise modeling and spatial clustering
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 toconsider the statistical analysis of fixations, especially theirlocation distribution. For example, Barthelmé, Trukenbrod,
Engbert, and Wichmann (2013) recommended using the point
process framework to infer how fixations are distributed inspace. Although we endorse this fruitful approach and its
Bayesian nature, here we aimed to resolve this problem fromthe opposite perspective. Instead of inferring from the spatial

Behav Res (2017) 49:559–575

563

distribution of the fixations, we inferred on each location in
the search space (i.e., each pixel within the eyetracker’s recordable range or each pixel in the visible stimuli). In otherwords, we addressed the question: BHow long is this pixel
being fixated (or what is the probability of this pixel being
fixated) in the function of the experimental conditions?^, byformally applying mixed models independently on each pixel,
we haveYðsÞ ¼ XβðsÞ þ ZbðsÞ þ εðsÞ

ð14ÞFor s ∈ D of the search space.
The complete procedure as implemented in iMap4 is explained in Fig. 1. The eye movement data for each participant
are concatenated into one input data matrix. iMap4 first partitions the data matrix into a fixation characteristic matrix (redbox) and an experiment condition information matrix (green

Fig. 1 Illustration of the basic processing steps implemented in
iMap4. The input data matrix is partitioned into an eye movement
matrix and predictor matrix. Fixation durations are projected into thetwo-dimensional space according to their x- and y-coordinates at thesingle-trial level for each participant. The experimental information ofeach trial is also summarized in a predictor table. Subsequently, the sparse
representation of the fixation duration map is smoothed by convoluting itwith a two-dimensional Gaussian kernel function, kerneleΝ ð0 ; σ2 Ι Þ.

box). The fixation characteristic matrix contains a fixation’s
spatial location (x and y), the fixation duration, and an orderindex of each fixation. The experiment condition matrix contains an index of each subject, an index of each trial/item, andthe different levels of each experimental condition. Fixation
durations are then projected into the two-dimensional spaceaccording 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 ∼N 0; σ2 I ;where I is a two-by-two identity matrix and the full width at
half maximum (FWHM) of the kernel is 1° of visual angle as
the default setting.
This step is essential to account for the spatial uncertainty
of eye movement recordings (both mechanical andAfter estimating the fixation bias of each condition independently for
all observers (by taking the expected values across trials within the same
condition), iMap4 models the 3-D smoothed fixation map(item*xSize*ySize) independently for each pixel using an LMM. The
result is saved as a MATLAB structure in LMMmap. iMap4 offers many
parametric and nonparametric methods for hypothesis testing andmultiple-comparison correction

564

Behav Res (2017) 49:559–575

physiological) and the sparseness of the fixation locations.
The Gaussian kernel could also be replaced by other 2-D
spatial filters to best suit the research question.
The resulting smoothed fixation map is a 3-D matrix. The
last two dimensions of the fixation matrix are the sizes of thestimuli or search space. The information of each entry in thefirst dimension is stored in a predictor table, which is generated from the experiment condition matrix. Each experimentcondition 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 a robust estimation option by
applying Winsorization in order to limit extreme values in thesmoothed fixation matrix. The goal here is to reduce the effectof any potential 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 3-D 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 ofLMMmap are nearly identical to the output from theLinearMixedModel 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 for the LMM.
Additional modeling specifications, as well as other model parameters, including the LMM’s formula, design matrix for fixed
and random effect, and residual degrees of freedom, are also
saved in LMMmap. Linear contrasts and other analyses basedon variance or covariance can be performed afterward from themodel-fitting information. Any other computation on the
LinearMixedModel output can also be replicated with LMMmap.
One of the crucial assumptions of pixel-wise modeling isthat 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 a large 2-D pixel space, we applied nonparametricstatistics to resolve the biases in parameter estimation andproblems 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 β operate on the fixed-effect-related variances. To doso, we simply removed the variance associated with the random effects from the response matrix:
Yfixed ðsÞ ¼ XβðsÞ þ εðsÞ ¼ YðsÞ−ZbðsÞ;ð15Þ

For s ∈ D of the search space.
For any permutation test, iMap4 performs the following
algorithms on Yfixed for each pixel.

Algorithm 1
For a given hypothesis or linear contrast c (as in Eq. 9), iMap4
–

–
–

–
–
–

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 Yfixed to the pseudo-inverse of M.
The design matrix M is created so that the original hypothesis testing is equivalent to the hypothesis regarding
the M1 coefficients. The matrix transformation and partition are the same as the algorithm described in Winkler
et al. (2014, Appx. A).
Computes the residuals related to the hypothesis bysubtracting the variance accounted for by M2 from
Yfixed, to get Yrr.
Fits Yrr to M by solving Yrr = Mβm + ε, and gets thestatistical value Frr of M1 according to Eqs. 10 and 11.
Note that to replicate the original hypothesis testing on
the fixed effect, the new contrast c’ is just used to partition
M into M1 and M2.
Permutes the rows of the design matrix M to obtain the
new design matrix M*.
Fits Yrr to M* and gets the Frr* of M1*.
Repeats the previous two steps a large number of times (kresamplings/repetitions), and the p value is then defined
as in Eq. 16. Importantly, the family-wise error rate
(FWER) corrected p value is computed by comparing
the largest Frr* across all tested pixels in one resamplingwith the original Frr:
# Frr * ≥ Frr
p¼
:
ð16Þ
k

Algorithm 1 is a simplified version of Winkler et al. (2014,
Algorithm 1): The resampling table includes permutation but
not sign-flipping, which assumes the errors to be independentand 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 values and only subjects as a random effect.
As was previously shown in Kherad-Pajouh and Renaud(2015), a general and exact permutation approach for mixedmodel designs should be performed on modified residuals thathave up to second-moment exchangeability. This is done tosatisfy the important assumptions for repeated measures
ANOVA: normality and the sphericity of errors. However,there are strict requirements to achieve this goal: careful transformation and partition of both the fixed- and random-effects
design matrices, and removal of the random effects related toM2 (Kherad-Pajouh & Renaud, 2015). In iMap4, we perform
an approximation version by removing all random effects toincrease the efficiency and speed of the huge amount ofBehav Res (2017) 49:559–575

resampling computation in our pixel-wise modeling algorithm. Validation and simulation data set indeed showed thatthe sensitivity and the false alarm rate of the proposed algorithm were not compromised.
Algorithm 2
iMap4 performs the following algorithm on Yfixed for eachpixel as the bootstrap clustering approach.
–

–
–

–

–
–

For each unique categorical variable, iMap4 removes theconditional expectations from Yfixed for each pixel. A
random shuffling is then performed on the centered datato acquire Yc, so that any potential covariance is also
disrupted. This is done to construct the true empiricalnull-hypothesis distribution in which all elements andtheir linear combinations in Yc have expected values
equal to 0.
Randomly draws with replacement from {X, Z, Yc} equal
numbers of subjects {X*, Z*, Yc*}.
Fits Yc* to X* by solving Yc* = X*β* + ε. For a given
hypothesis or linear contrast c (as in Eq. 9), iMap4 computes the statistics value F* according to Eqs. 10 and 11,
and their parametric p value under the GLM framework.
Thresholds the statistical maps F* at p*≤.05 and records
the desired maximum cluster characteristics across all
significant clusters. The cluster characteristics considered
are cluster mass (summed F value within a cluster), cluster extent (size of the cluster), and cluster density (mean F
value).
The previous three steps are repeated a large number oftimes, to get the cluster characteristic distribution under
the null hypothesis.
Thresholds the original statistical map F at p≤ :05 andcompares the selected cluster characteristic with the valueof the null distribution corresponding to the 95th 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; Pernet et al.,2015) if only a subject intercept is considered as the random
effect. In addition, Algorithm 2 extents the philosophy andapproach presented by Pernet et al. (2011; Pernet et al.,2015) to nonhierarchical mixed-effect models.
It is worth noting that we implemented in iMap4 a highperformance algorithm to minimize the computational demands of the large amount of resampling. The model fitting
in both resampling approaches makes use of ordinary leastsquares. The inversion of the covariance matrices (requiredfor Eq. 11) is computed on the upper triangular factor of theCholesky decomposition. Calculation of the quartic form (asin Eq. 11) for all pixels is optimized by constructing a sparse565

matrix of the inverse of the covariance matrix. More details ofthese algrebraic simplifications can be found in the
imapLMMresample function in iMap4.
Other multiple-comparison correction methods, such as
Bonferroni correction, false discovery rate, or random field
theory (RFT), could also be applied. A threshold-free cluster
enhancement algorithm could also be applied on the statistical(F-value) maps as an option after the permutation and bootstrap clustering procedures (Smith & Nichols, 2009).
We performed a validation study to assess the Type I error
rate when applying the permutation and bootstrap clustering
approach for hypothesis testing. We used a balanced repeatedmeasures ANOVA design with a two-level between-groupfactor and a three-level within-group factor. A total populationof 134 observers (67 in each group) was drawn from previous
face-viewing eye movement studies. We centered the cellmeans for the whole dataset to obtain the validation dataset
under the null hypothesis (similar to Step 1 in Algorithm 2).
Thus, we used real data to warrant realistic distributions andcentered them to ensure that H0 was confirmed. Any significant output from iMap4 performed on this dataset would beconsidered as a false alarm (Type I error).
The validation procedure followed the steps below: We
first randomly sampled without replacement a balanced number of subjects from both groups. We then ran iMap4 under thedefault settings and performed hypothesis testing on the twomain effects and the interaction. To estimate the FWER, we
computed the frequency of significant output under different
statistics and MCC settings. Preliminary results based on 1,
000 randomizations with a sample size of n ∈ [8, 16, 32, 64]
showed that with an alpha of .05, the FWERs were indeed allunder. 05 using nonparametric statistics (see Fig. 2b for the
permutation test, and Fig. 2c and d for the bootstrap clusteringtest). More simulations considering a wider range of scenarios
will be required to understand fully the behavior of the proposed approaches, although the cluster stats are likely to behave as in Pernet et al. (2015).

Graphical user interface (GUI) and command line
handling
iMap4 runs on MATLAB 2013b and above, since it requires
some essential functions and classes from the Image
Processing Toolbox and Statistics Toolbox in these versions.
iMap4 will execute in parallel on multicores or distributed
workers, when available.
We recommend that users install iMap4 as a MATLAB
application. The users can call iMap4 directly in the
MATLAB command window after installation. A general
GUI will open upon >>iMAP, called in the command windowor when launching the app (Fig. 3a). The users can then import
the fixation data, load a preprocessed data matrix for LMM, or566

Behav Res (2017) 49:559–575

a)b)

BetweenWithinInteraction

.2

FWER

.05

.1

8

16

32
Number of Subjects

64

c)

8

16

32

64

8

16

32

64

d)

.05

.05

8

16

32

64

Fig. 2 Validation results of the proposed resampling procedure forstatistical inference. (a) Family-wise error rates (FWERs) using theuncorrected parametric p values. All FWERs are significantly above.05. (b) FWERs using the permutation approach (Algorithm 1). (c)FWERs using the proposed bootstrap-clustering approach (Algorithm2) thresholds on cluster mass. (d) FWERs using the proposedbootstrap-clustering approach (Algorithm 2) thresholds on cluster
extent. Notice that the FWERs of panels a and b are computed at the
pixel level (i.e., the proportions of false-positive pixels acrosssimulations), and the FWERs of panels c and d are calculated at the testlevel (i.e., the percentages of any false positives per test for the 1,000
simulations). Error bar shows the standard errors

display the modeling results and perform statistical hypothesis
testing. These main steps have their own independent GUIs:
Create Fixation Matrix (Fig. 3b), Linear Mixed Model(Fig. 3c), and Display Results (Fig. 3d). Although most features of iMap4 could be obtained via these GUIs, we encourage advanced users to use command lines, especially for the
additional options specification of the LinearMixedModelclass. A short example of the command-line handling of themain functions is shown in Fig. 3e. A user guidebook containing the instructions for each step can be accessed via the Help
button. We have also provided datasets with tutorial files to
explain practically how to use iMap4. As a demonstration,
two examples based on real and simulation data are given in
the next section. MATLAB scripts of the examples are part ofthe iMap4 installation package.

Applications to real and simulation data
In the following examples, we illustrate iMap4’s flexibility
and power with two real data sets and a computer simulation.
All material and codes presented here are available in theiMap4 installation package.
Example 1
We consider first a subset of participants from Bovet, Lao,
Bartholomée, Caldara, & Raymond, (2016), as a demonstration
Fig. 3 The main graphical user interfaces of iMap4 (a–d) and example
command lines handling the core functions (e). For more details, pleaserefer to the online guidebook and demonstration codes

Behav Res (2017) 49:559–575

567

568

Behav Res (2017) 49:559–575

of the analysis procedure in iMap4. A step-by-step demonstration is available in the user guidebook and example code.
In short, the dataset consists of eye movement data from 20
male observers during a gaze-contingent study. Observers
viewed computer-rendered female bodies in different conditions and performed a behavioral task (i.e., subjective ratingsof bodily attractiveness). This was a within-subjects design
with two experimental manipulations: the viewing condition(three levels: 2° spotlight, 4° spotlight, or natural viewing) andbody orientation (two levels: front view or back view). The
aim of the study was to evaluate the use of visual information
for bodily attractiveness evaluation in the male observers.
Other details of the experiment can be found in the article.
Fixation durations were projected into the two-dimensional
space according to their coordinates at the single-trial level.
The fixation duration maps were first smoothed at 1° of visual
angle. We used the Bestimated^ option by taking the expectedvalues across trials within the same condition independently
for each observer. To reduce the computation time, we downsampled the fixation map to 256*205 pixels and applied amask to only model the pixels with average durations longer
than half of the minimum fixation duration input.
Before proceeding to the modeling step, we visualized thepreprocessed fixation maps and the descriptive statistics to geta sense of the data. For each of the categorical conditions,
iMap4 outputs the mean fixation map for each level.
Descriptive statistics for the following eye movement measures are saved in a matrix and will be plotted in a histogram
or boxplot: number of fixations, sum of fixation durations(total viewing time), mean fixation duration, total path length(total eye movement distance in pixels), and mean path length.
See Fig. 4 for an example of the descriptive-results output.
We applied a full model on the fixation duration map without any spatial normalization:
PixelIntensity ðx; yÞ e Viewing condition þ Body orientationþ Viewing condition* Body orientationþ ðfixation duration jsubjectÞ;x; y ∈ fixation map resolution:

Notice that the mean fixation duration for each conditionand subject were treated as random effects to control for thevariation across individuals. The parameters were fitted withrestricted maximum likelihood estimation (ReML).
We encourage users to interpret the result from iMap4 in
the following way. First, check the model fitting by displaying
the model criteria. For example, Fig. 5a shows the R2 values ormultiple-correlation coefficients, which represent the proportions of variability in the fixation matrix explained by thefitted model. Interpretation of the result should be drawn withcaution if the R2 values are too low. The users can then proceed to test their hypotheses, such as through ANOVA or
linear contrast, and perform multiple-comparisons corrections(Fig. 5b and c). A post-hoc analysis is applicable if any interaction is presented, or if any condition contains multiple
levels. The user can select one or more significant area(s) as
data-driven ROI(s) for the post-hoc analysis. iMap4 performst tests between any pairs of categorical conditions within thisROI by using the raw input values from the nonsmootheda)Fixation Number

100

200

Fixation Duration(SUM)

150

Fixation Duration(MEAN)

4° spotlight

Natural Viewing

10

6 × 10 4

0Path Length(TOTAL)

15

500

1000

Path Length(MEAN)

Back view

02° spotlight

Front view

b)

10

8
6

10

4
5
0

2
0

15000

0

0

300

Fig. 4 Descriptive results from iMap4 on the real dataset. (a) Five eye movement measures plotted in histograms. In this case, fixation durations are in
milliseconds and path lengths are in pixels. (b) Mean fixation maps of all levels of the categorical conditions

Behav Res (2017) 49:559–575

569

a)

b)
Ordinary R2Main effect ofViewing Condition

F-value

Variance explained
20%
80%

c)

Main effect ofBody Orientation

n.s

60

d)
Linear contrast

-3.92

NV

_b

ac

NV
_f
ro
nt
k

Tvalue (p < .05 bonferroni corrected)

4°

_f

ro
n

t

5.53

t

-3.49

6.31

6.16

2°

_f

ro
n

3.52

k

4°

_b

ac

k

-4.61

-5.74

-8.87

_b

ac

F-value
60

ac

k

t
ro
n

_b
2°

k

_f

ac
2°

_b
4°

ro
n

t

k

_f

ac
4°

_b
NV

NV
_f
ro

nt

2°

n.s

e)
2° spotlight front view

2° spotlight back view

coefficients (beta)0coefficients (beta)
6000Fig. 5 iMap4 results for Bovet et al. (2016) with different output styles.
(a) Ordinary R2 values for the fitted model. (b) ANOVA results of themain effects and interaction. Here the intensity represents the F values.
iMap4 only displays significant maps. (c) Statistical results of the linear
contrast [2° spotlight–natural viewing] in the back view condition. Herethe F value is represented on a contour map. (d) Post-hoc analysis in theselected mask. The mask is generated from the significant region of the

200

body orientation effect (left panel). The t test results are shown in thematrix in the right panel (labeled conditions on the y-axis minus thoselabeled on the x-axis). Only significant results are shown (p < .05,
Bonferroni corrected). (e) One-tailed t tests against the average over all
fixation intensities for the 2° spotlight front view and 2° spotlight backview conditions. The solid black lines contain the significant regions forall of the panels above

570

Behav Res (2017) 49:559–575

fixation matrix (Fig. 5d). In addition, users can compute theabove-average or above-chance fixation intensity for each categorical predictor (Fig. 5e).
Example 2
As a second demonstration, we reanalyzed the full dataset
from one of our previous studies, Miellet, He, Zhou, Lao,and Caldara (2012).
Previous studies testing Western Caucasian (WC) and East
Asian (EA) observers had shown that people deploy different
eye movement strategies during free viewing of faces. WC
observers fixate systematically toward the eyes and mouth,
following a triangular pattern, whereas EA observers predominantly fixated at the center of the face (Blais, Jack, Scheepers,
Fiset, & Caldara, 2008; Caldara, Zhou, & Miellet, 2010).
Moreover, human observers can flexibly adjust their eye
movement strategies to adapt to environmental constraints,
as has been shown using different gaze-contingent paradigms(Caldara, Zhou, & Miellet, 2010; Miellet et al., 2012). In our2012 study, we tested two groups of observers in a face task inwhich their foveal vision was restricted by a blind spot. This
was a mixed design with the culture of the observers as the
between-subjects factor (WCs or EAs) and the blind spot sizeas the within-subjects factor (four level: natural viewing, 2°
blindspot, 5° blindspot, or 8° blindspot). For more details of
the experiment, please refer to Miellet et al. (2012).
Using iMap4, we created the single-trial 2-D fixation duration map and smoothed at 1° of visual angle. Importantly, tokeep in line with Miellet et al. (2012), spatial normalization
was performed by z-scoring the fixation map across all pixels
independently for each trial (the results are identical without
spatial normalization in this example). We also applied a mask
generated with the default option. No down-sampling was
performed. We then applied a full model on the single-trial
fixation duration map made used of the Bsingle-trial^ option in
iMap4:
PixelIntensity ðx; yÞ e Observer culture þ Blindspot size
þ Observer culture*Blindspot sizeþ ð1jsubjectÞ; x; y ∈ fixation map resolution:

Only the subject predictor was treated as a random effect,
and the model was fitted using ML.
After model fitting, we performed an ANOVA to test thetwo main effects and their interactions. We applied a bootstrap
clustering test using a cluster density of 1,000 as the criterion.
We found a significant interaction and a main effect of blindspot size, but not a main effect of culture (see Fig. 6a). This
result replicates the findings in Miellet et al. (2012).
Moreover, by performing a linear contrast of the model coefficients, we reproduced Fig. 2 from Miellet et al. (2012). The
results using iMap4 are shown in Fig. 6b.

Example 3
We also used simulated data to illustrate the use of iMap4 with
continuous predictor. We created a dataset and manually introduced an effect between the numbers of fixations and thesubjective rating on a single-trial level. Moreover, to maximize the simulation’s efficiency, different linear relationships
were introduced simultaneously. For each subject, we generated a data matrix through the two following steps:––

In a 4*4 grid, we introduced a different linear relationshipin each cell between fixation number and subjective rating. Figure 7a shows the linear relationships we introduced for one subject. We varied the slope and thestrength of the linear association. The correlation wasstrongest on the top row (r = .9), and there was no correlation on the bottom row (r = 0). The slope varied among[1, 0.4, –0.2, –0.8] across the columns. Note that each dot
on a scatterplot represents one trial, and the dots with thesame rating (value on the x-axis) across subplots belong
to the same trial. The resulting matrices after this step
were a one-dimensional array Rating and a twodimensional matrix P (matrix size: 16 * number of trials)
The spatial locations of fixations were generated using
linear Gaussian random fields. For each trial, we createda Gaussian mixture model gm using the gmdistribution
class in MATLAB. The Gaussian mixture model gm contains 16 (4*4) 2-D Gaussian distribution components.
The center of each component aligned with the center of
each grid, and the covariance was an identity matrix with
1° of visual angle on the diagonal. Crucially, the mixing
proportion of each component was decided by the columnof the specific trial in P. A number of random fixations
were then generated from this Gaussian mixture model
gm. See Fig. 7b for a realization of one random trial for
one subject.

The dataset contained 20 subjects performing 100 trials,
each with an average fixation number of 58.02. Figure 7cshows the average map for fixation numbers. We fitted a simple model with ReML:
Pixel Intensityðx; yÞ e 1 þ Ratingþ ð1jsubjectÞ;x; y∈screen resolution:

The significant regression coefficients of Rating are shown in
Fig. 7d. iMap4 accurately rejected the null hypothesis for mostconditions when there was a significant relationship. For the mostrobust effect (r = .9), iMap4 accurately estimated the coefficients.
It also correctly reported a null result for r = 0. Moreover, iMap4
did not report any significant effect for the weakest relationship(slope = –0.2, r = .3), due to the lack of power. Indeed, furtherBehav Res (2017) 49:559–575

571

Main effectof BlindSpot size

Main effectof Culture

Interaction

WC observers

EA observers

WC - EA

N.S.

F-value

35

a)

b)

NV

100
50

5°
0

Fixation duraction (ms)2°8°NV - 8°Linear contrast of theβ (model coefficients)
− 60

0

60

Fig. 6 iMap4 results for Miellet et al. (2012). (a) ANOVA results of the linear mixed model. (b) Replication of the Fig. 2 results for Miellet et al. (2012),
using linear contrasts of the model coefficients. The solid black lines contain the significant regions for all of the panels above

simulations showed that increasing the numbers of fixations,trials, or subjects would lead to significance.

Discussion and future developments
In the present article, we have reported a major update of iMap,
a toolbox for statistical fixation mapping of eye movement data.

While keeping unchanged the general data-driven philosophyof iMap, we significantly improved the underlying statisticalengine, by incorporating pixel-wise LMMs and a variety ofrobust nonparametric statistics. Crucially, the new analysis
pipeline allows for the testing of complex designs while controlling for a wide range of random factors. We also implemented a full GUI to make this approach more accessible to
MATLAB beginners. Examples from empirical and572

computer-simulated datasets showed that this approach has aslightly conservative FWER under H0, while remaining highly
sensitive to actual effects (e.g., Fig. 6d). The present method
represents a significant advance in eye movement data analysis,
particularly for analyzing experimental designs using normalized visual stimuli. In fact, iMap4 uses a statistical inference
method similar to those in fMRI and magnetoencephalography/
EEG analysis. The interpretation of the statistical maps is simply done by looking at which stimulus features/pixels relate to
the significant areas (after multiple-comparison correction).
This procedure is similar to the interpretation of fMRI results:
after a significant region is revealed, we can use its spatial
coordinates to check in which part of the cortex the region
activated above chance level is located.
As a powerful statistical tool, LMMs are gaining popularity
in psychological research and have previously been applied ineye movement studies (e.g., Kliegl, Masson, & Richter, 2010).
Similarly, particular cases of LMM, such as HLM or two-level
models, are now standard data-processing approaches in neuroimaging studies. As a general version of HLMs, LMMs aremuch more flexible and powerful than other multilevel
models. Most importantly, an exact same LMM could be applied to behavior, eye movement, and neuroimaging data,
bridging these different measures to allow drawing more direct and complete conclusions.
However, there are both theoretical and practical challenges in using LMM for the statistical spatial mapping offixation data. First, the fixation locations are too sparse todirectly apply pixel-wise modeling. Similarly to previous versions of iMap, we used spatial smoothing of the fixation locations, a preprocessing step necessary to account for the measurement error of eyetrackers and the imprecision of the physiological system (i.e., the human eye). The second issue isselecting the appropriate hypothesis testing for LMM andthe multiple-comparison problems caused by modeling massive number of pixels in nonbalanced designs. We addressed
this issue by applying nonparametric statistics based on resampling and spatial clustering. Another important challenge
is the constraint of computational resources. Parameter estimations using LMM, the pixel-wise modeling approach, andresampling techniques are very computationally demanding
and time-consuming. To produce a useful but also usable tool,
we adapted many advanced and novel algorithms, such asparallel computing. Preprocessing options such as downsampling and applying a mask also significantly decrease
the computational time of iMap4.
The comparison among ROIs/areas of interest, iMap 2.0,
and the current version
In classical eye movement data analyses, particularly those considering fixation locations, the main challenge for statstially identifying the regions that have been fixated above chance level lies

Behav Res (2017) 49:559–575

in the fact that we are facing a high-dimensional data space.
Mathematically, each pixel represents one dimension that could
be potentially important. However, it is trivial to say that many ofthese dimensions are redundant and could be reduced to a particular set of representations or features. In other words, eyefixation data points are embedded in a high-dimensional pixel
space, but they actually occupy only a subspace with much lower
dimensionality (Belkin & Niyogi, 2003). Indeed, in similar highdimensional datasets, a low-dimensional structure is often assumed and is naturally the main focus for investigation. Thus,by arbitrarily choosing one or multiple ROIs, one can represent
the high-dimensional dataset as a low-dimensional manifold.
The fixation map thus projects into this manifold, and all the
pixels within the same ROI are then considered as being in thesame dimension. In this case, each ROI represents one feature.
Such a method is comparable to early neural network and manyother linear dimension reduction methods in the machinelearning literature with hand-coded features (LeCun, Haffner,
Bottou, & Bengio, 1999; Sorzano, Vargas, & Montano, 2014).
The early versions of iMap (1 and 2) adopted a similar
logic, but relied on RFT to isolate data-driven features.
Therefore, the fixation bias in each pixel was projected intoa lower-dimensional subspace, resulting in fixation clusters.
Second-level statistics were then computed at the cluster level
instead of the pixel level to perform statistical inference(Miellet, Lao, & Caldara, 2014).
From iMap 3 onward, we took a very different approach.
We used spatial clustering and multiple-comparison correction to avoid the use of second-level statistics to perform statistical inference. In iMap4, the fixation bias is similarly
modeled on each pixel using a flexible yet powerful statisticalmodel: the LMM. The LMM, in combination with nonparametric statistics and a spatial clustering algorithm, directly
isolates the significant pixels. As a result, the iMap4 outputs
can be interpreted intuitively and straightforwardly at the map
level (i.e., by visualizing the areas reaching significance fromthe tested hypothesis).

Parameter settings and statistical choices
Our aim was and still is the development of a data-driven andfully automatized analysis tool. However, even in iMap4
Fig. 7 iMap4 results on the simulation dataset. (a) Linear relationships
being introduced into the 4*4 grid. The x-axis shows the z-scored rating,and the y-axis shows the expected number of fixations. The slopes
between y and x are the same within each column ([1, 0.4, –0.2, –0.8],
respectively), whereas the correlation rho is the same within each row([0.9, 0.6, 0.3, 0], respectively). (b) One realization of a random trial for
one subject. The left panel shows the raw fixation locations; the rightpanel shows the smoothed fixation number map. (c) The average
fixation map across all trials for the 20 subjects. (d) Estimatedrelationships between rating and fixation number (regression
coefficients). The black circles indicate statistical significanceBehav Res (2017) 49:559–575

a)

573

y=x+c

y = .4 * x + c

y = -.2 * x + c

y = -.8 * x + c

r = 0.9

r = 0.6

Responser = 0.3

r=0

5

0

-2

0

2

Ratingb)c)d)Fixation number map (mean)(beta)

1

0

-0.8

574

some of the parameters in the analysis rely on a user’s expertise and subjective choices, which thus should be considered
carefully before use. These parameters include the kernel sizefor the smoothing procedure, the spatial down-sampling andmasking, the spatial normalization, and the choice of statistics.
The rationale for the determining the kernel size for thesmoothing procedure has been previously discussed (Caldara& Miellet, 2011), and the majority of the arguments we put
forward in this previous article still hold true. Here, we would
remind users that the spatial smoothing procedure mainly resolves the sparseness of fixation data. It also partially accounts
for the spatial covariance, which is ignored in univariate pixelwise modeling. Finally, it accounts for the recording errors
from eyetrackers, such as drift during the calibration, pupilsize variations, and so forth.
We also recommend that users perform down-sampling andapply a mask before modeling their data. This step is importantto reducing computational demands (time, memory, etc.). In
general, we recommend that the down-sampling factor not bebigger than half of the smoothing kernel size. In other words, if
the FWHM of the Gaussian kernel is 10 pixels, the rescale factor
should be less than 5. We are currently running further simulations and validations to investigate the best parameters under
different settings, and hopefully will provide a statistical datadriven solution for this choice in future updates.
Spatial normalization (via a z-scored or probability map) isavailable as an option in iMap4. Spatial normalization used to
be a standard preprocessing procedure in previous versions of
iMap. However, the hypotheses tested on raw fixation
duration/number maps are fundamentally different from their
spatially normalized versions. Importantly, after spatial normalization, the interpretation of results should be drawn on a
spatially relative bias instead of on the absolute differences.
Of course, if the viewing duration in each trial is constantwithin an experiment, spatial normalization will not makeany difference.
For iMap4 we developed two main, nonparametric statistics based on resampling techniques. It is worth noting thatdifferent applicability comes with the choice of permutation
tests versus bootstrap spatial-clustering tests. In our own experience during empirical and simulation studies, permutation
tests are more sensitive for studies with small sample sizes; thebootstrap-clustering approach usually gives more homogeneous results but is biased toward bigger clusters. We suggestthat users adopt a Bwisdom of crowds^ approach and look at
the agreement among different approaches before concluding
on the data analysis (Marbach et al., 2012). Nonconvergent
results should be interpreted carefully.
An alternative to pixel-wise approaches
In recent years, other frameworks have been also developed to model eyetracking data (Boccignone, 2015). One

Behav Res (2017) 49:559–575

such approach is the aforementioned Poisson point process model (Barthelmé et al., 2013). It is a wellestablished statistical model when the point (fixation)
occurrence is the main concern. Under some transformation, the Poisson point processes model of fixation occurrences could be expressed and modeled as a logistic
regression, making it straightforward to apply using conventional statistical software (Barthelmé & Chopin,
2015). For example, Nuthmann and Einhäuser (2015)
made use of logistic mixed models to determine the influence of low- and high- visual properties in scene
images on eye movements. Moreover, smooth effectand spatial covariants could be captured by applying regression splines in a generalized additive model, as demonstrated in Barthelmé and Chopin (2015).
Importantly, the point process model addresses different
questions than does iMap. It is most appropriate when theeffect of spatial location is considered irrelevant, a nuisanceeffect, or a fixed intercept (see, e.g., Barthelmé & Chopin,
2015; Nuthmann & Einhäuser, 2015). As a comparison, in
iMap the parameters of interest are location specific, varying
from pixel to pixel. In other words, the differences or effects
among different conditions are location-specific, forming acomplex pattern in two dimensions. These high dimension
effects are more natural and easy to model using a pixelwise model, as in iMap4.
Conclusion and future development
In conclusion, we have presented an advanced eye movement
analysis approach using LMMs and nonparametric statistics:
iMap4. This method is implemented in MATLAB with a userfriendly interface. We aimed to provide a framework for analyzing spatial eye movement data with the most sophisticatedstatistical modeling to date. The procedure described in thepresent article currently represents our best attempt to conformwith the conventional null-hypothesis testing, while providing
options for robust statistics. We currently are still working on
many improvements, including functions to compare different
fitted models, statistics on the random-effect coefficients, andreplacing LMMs with generalized LMMs for modeling fixation numbers (Bolker, Brooks, Clark, Geange, Poulsen,
Stevens, & White, 2009). In the future, we will also switch
our focus to Bayesian statistics and the generative model (suchas the Gaussian process) in an effort to develop a unified
model of statistical inference for eye movement data (Jaynes& Bretthorst, 2003).
Author note The authors declare no competing financial interests. This
study was supported by the Swiss National Science Foundation (Grant
No. 100014_138627). awarded to R.C. The authors would like to thankDr Simon Barthelmé and an anonymous reviewer for their helpful comments that contributed to improving the final version of the article.

Behav Res (2017) 49:559–575

References
Baayen, R. H., Davidson, D. J., & Bates, D. M. (2008). Mixed-effects
modeling with crossed random effects for subjects and items.
Journal of Memory and Language, 59, 390–412. doi:10.1016/j.
jml.2007.12.005
Barthelmé, S., & Chopin, N. (2015). The Poisson transform for
unnormalised statistical models. Statistics and Computing, 25,
767–780.
Barthelmé, S., Trukenbrod, H., Engbert, R., & Wichmann, F. (2013).
Modeling fixation locations using spatial point processes. Journalof Vision, 13(12), 1. doi:10.1167/13.12.17
Belkin, M., & Niyogi, P. (2003). Laplacian eigenmaps for dimensionality
reduction and data representation. Neural Computation, 15, 1373–
1396.
Blais, C., Jack, R. E., Scheepers, C., Fiset, D., & Caldara, R. (2008).
Culture shapes how we look at faces. PLoS ONE, 3, e3022.
Boccignone, G. (2015). Advanced statistical methods for eye movement
analysis and modeling: A gentle introduction. arXiv preprint. arXiv:
1506.07194.
Bolker, B. M., Brooks, M. E., Clark, C. J., Geange, S. W., Poulsen, J. R.,
Stevens, M. H. H., & White, J. S. (2009). Generalized linear mixed
models: A practical guide for ecology and evolution. Trends in
Ecology and Evolution, 24, 127–135. doi:10.1016/j.tree.2008.10.
008
Bovet, J., Lao, J., Bartholomée, O., Caldara, R., & Raymond, M. (2016).
Mapping females’ bodily features of attractiveness. Scientific
Reports, 6, 18551. doi:10.1038/srep18551
Caldara, R., & Miellet, S. (2011). iMap: A novel method for statistical
fixation mapping of eye movement data. Behavior Research
Methods, 43, 864–878. doi:10.3758/s13428-011-0092-x
Caldara, R., Zhou, X., & Miellet, S. (2010). Putting culture under the
Bspotlight^ reveals universal information use for face recognition.
PLoS ONE, 5, e9708. doi:10.1371/journal.pone.0009708
Christensen, R. (2011). Plane answers to complex questions: The theoryof linear models. Berlin: Springer.
Duchowski, A. T. (2002). A breadth-first survey of eye-tracking applications. Behavior Research Methods, Instruments, & Computers, 34,
455–470.
Friston, K. J., Stephan, K. E., Lund, T. E., Morcom, A., & Kiebel, S.
(2005). Mixed-effects and fMRI studies. NeuroImage, 24, 244–252.
Jaynes, E. T., & Bretthorst, G. L. (2003). Probability theory: The logic ofscience. Cambridge: Cambridge University Press.
Kherad-Pajouh, S., & Renaud, O. (2010). An exact permutation method
for testing any effect in balanced and unbalanced fixed effect
ANOVA. Computational Statistics and Data Analysis, 54, 1881–
1893.
Kherad-Pajouh, S., & Renaud, O. (2015). A general permutation approach for analyzing repeated measures ANOVA and mixedmodel designs. Statistical Papers, 56, 947–967. doi:10.1007/
s00362-014-0617-3
Kliegl, R., Masson, M. E., & Richter, E. M. (2010). A linear mixed model
analysis of masked repetition priming. Visual Cognition, 18, 655–
681.
LeCun, Y., Haffner, P., Bottou, L., & Bengio, Y. (1999). Object recognition with gradient-based learning Shape, contour and grouping incomputer vision (pp. 319–345). New York: Springer.
Liversedge, S., Gilchrist, I., & Everling, S. (2011). The Oxford handbookof eye movements. Oxford: Oxford University Press.

575
Marbach, D., Costello, J. C., Küffner, R., Vega, N. M., Prill, R. J.,
Camacho, D. M., . . . Stolovitzky, G. (2012). Wisdom of crowds
for robust gene network inference. Nature Methods, 9, 796–804.
doi:10.1038/nmeth.2016
McCulloch, C. E., Searle, S. R., & Neuhaus, J. M. (2011). Generalized,
linear, and mixed models. New York: Wiley.
Miellet, S., He, L., Zhou, X., Lao, J., & Caldara, R. (2012). When East
meets West: Gaze-contingent blindspots abolish cultural diversity ineye movements for faces. Journal of Eye Movement Research, 5(5),1–12.
Miellet, S., Lao, J., & Caldara, R. (2014). An appropriate use of iMap
produces correct statistical results: A reply to McManus (2013)
BiMAP and iMAP2 produce erroneous statistical maps of eyemovement differences.^. Perception, 43, 451–457.
Nuthmann, A., & Einhäuser, W. (2015). A new approach to modeling theinfluence of image features on fixation selection in scenes. Annals of
the New York Academy of Sciences, 1339, 82–96.
Oostenveld, R., Fries, P., Maris, E., & Schoffelen, J.-M. (2011).
FieldTrip: Open source software for advanced analysis of MEG,
EEG, and invasive electrophysiological data. Computational
Intelligence and Neuroscience, 2011, 9.
Penny, W. D., Friston, K. J., Ashburner, J. T., Kiebel, S. J., & Nichols, T.
E. (2011). Statistical parametric mapping: The analysis of functional brain images: The analysis of functional brain images. San
Diego: Academic Press.
Pernet, C. R., Chauveau, N., Gaspar, C., & Rousselet, G. A. (2011).
LIMO EEG: A toolbox for hierarchical LInear MOdeling ofElectroEncephaloGraphic data. Computational Intelligence and
Neuroscience, 2011, 11.
Pernet, C. R., Latinus, M., Nichols, T. E., & Rousselet, G. A. (2015).
Cluster-based computational methods for mass univariate analysesof event-related brain potentials/fields: A simulation study. Journalof Neuroscience Methods, 250, 85–93. doi:10.1016/j.jneumeth.
2014.08.003
Pinheiro, J. C., & Bates, D. M. (2000). Mixed-effects models in S and SPLUS. Berlin: Springer.
Raudenbush, S. W., & Bryk, A. S. (2002). Hierarchical linear models:Applications and data analysis methods. Thousand Oaks: Sage.
Rayner, K. (1998). Eye movements in reading and information processing: 20 years of research. Psychological Bulletin, 124, 372–422. doi:
10.1037/0033-2909.124.3.372
Smith, S. M., & Nichols, T. E. (2009). Threshold-free cluster enhancement: Addressing problems of smoothing, threshold dependence
and localisation in cluster inference. NeuroImage, 44, 83–98.
Sorzano, C. O. S., Vargas, J., & Montano, A. P. (2014). A survey ofdimensionality reduction techniques. arXiv preprint. arXiv:
1403.2877.
Winkler, A. M., Ridgway, G. R., Webster, M. A., Smith, S. M., &Nichols, T. E. (2014). Permutation inference for the general linear
model. NeuroImage, 92, 381–397.

Supplemental material
The full iMap4 toolbox and the data samples shown in this article arefreely available for download (http://ibmlab.github.io/iMap4/). A fulluser guide can be found separately from https://github.com/iBMLab/
iMap4/wiki