Applying the pair correlation function to eye movements

Here, we show how to compute the (PCF) for eye movement data in three steps. The PCF reveals whether the distribution of fixation locations during a single scanpath can be explained by the overall inhomogeneity observed across all observers or whether fixation locations of a single scanpath contain additional spatial correlations. For details we refer to our manuscript (Trukenbrod, Barthelmé, Wichmann, & Engbert, 2018).

library(tidyverse)
library(spatstat)
library(ggplot2)
library(parallel)

Data Set

For our analyses we use an experiment, where participants viewed two types of images twice. The repeated presentation of images allows to investigate the influence of visual long-term memory on eye movements. From previous work we expect that the repeated presentation leads to similar fixation densities but shortens saccade amplitudes during the second inspection. The resulting point patterns are similar but differ slightly in the overall inhomogeneity, which makes a direct comparison of the eye movement behavior difficult. The same problem is true for the comparison of different image types. The PCF takes differences in the underlying inhomogeneity into account and allows a direct comparison of the spatial correlations under different viewing conditions (first vs. second presentation) and for different image types (natural vs. texture images).

Our data set contains scanpaths from 35 participants viewing 30 images twice. Images were either natural scenes or images of texture images. The initial fixation, fixations containing blinks, and fixations outside the image boundaries were removed from further analyses. Overall 57344 fixations remained.

# image coordinates
xrange <- c(1.03585769369965,32.1115885046892)
yrange <- c(0.828686154959722,25.6892708037514)
dat <- read.table('./data/SpSt1.dat',header=TRUE)
d <- dat %>%
  mutate(.,xpos=(xR+xL)/2,ypos=(yR+yL)/2) %>% 
  filter(., nth>1 & blinkFix==0 & blinkSac==0 &
           xpos>xrange[1] & xpos<xrange[2] &
           ypos>yrange[1] & ypos<yrange[2]) %>%
  select(.,id,type,nthPres,image,xpos,ypos)
d$id <- factor(d$id)
d$type <- factor(d$type,levels=c(1,2),labels=c('Natural scene','Texture'))
d$nthPres <- factor(d$nthPres,levels=c(1,2),labels=c('First','Second'))
d$image <- factor(d$image)
summary(d)

       id                   type         nthPres          image            xpos             ypos       
 3      : 2545   Natural scene:29061   First :32201   10     : 2125   Min.   : 1.095   Min.   : 0.840  
 21     : 2077   Texture      :28283   Second:25143   8      : 2083   1st Qu.:11.350   1st Qu.: 9.815  
 28     : 2017                                        24     : 2012   Median :16.055   Median :14.005  
 12     : 2010                                        26     : 2007   Mean   :16.561   Mean   :13.705  
 23     : 1992                                        3      : 2001   3rd Qu.:21.995   3rd Qu.:17.765  
 8      : 1814                                        6      : 1997   Max.   :32.110   Max.   :25.670  
 (Other):44889                                        (Other):45119                                    

Step 1. Simulation of inhomogeneous and homogeneous control processes.

For our PCF computations, we need to simulate two control point processes, namely a homogeneous and an inhomogeneous point process. Points (fixation locations) are sampled independently from each other in both control processes and due to the independence of points, we do not expect to observe any spatial correlations between points at distance \(r\). Any observed correlations would be spurious and depend on the data structure (e.g., length of fixation sequences) or a wrong parameterization of the method. Hence, both control processes ensure that correlations in the PCF arise from the empirical data and not by the method itself. In addition, the inhomogeneous point process is used in the second step to estimate an optimal bandwidth for the intensity estimation of the PCF in Step 3.

For simulation of the control processes we need to pick a bandwidth for the estimation of the fixation density. Here, we use Scott’s rule of thumb (\(\tt bw.scott()\)). For each empirical scanpath we simulated one scanpath of equal length (same number of fixations as observed in the experiment) for the inhomogeneous point process and for the homogeneous point process. As a result we now have fixations from three point processes (PP) with the same number of fixations each.

sigma <- NULL
for (img in unique(d$image)){
  for (pres in unique(d$nthPres)){
    idx <- d$nthPres==pres & d$image==img
    dImg <- d[idx,]
    ppImg <- ppp(dImg$xpos,dImg$ypos,window=owin(xrange,yrange))
    bwImg <- mean(bw.scott(ppImg))
    sigma <- rbind(sigma,data.frame(img,pres,bwImg))
    denImg <- density(ppImg,sigma=bwImg)
    sim <- rpoint(nrow(dImg),denImg)
    d$xposInhom[idx] <- round(sim$x,digits=3)
    d$yposInhom[idx] <- round(sim$y,digits=3)
  }
}
sim <- rpoint(nrow(d))
d$xposHom <- round(sim$x*diff(xrange)+xrange[1],digits=3)
d$yposHom <- round(sim$y*diff(yrange)+yrange[1],digits=3)
dx <- gather(d,'pp','x',c(5,7,9))
dy <- gather(d,'pp','y',c(6,8,10))
d <- cbind(dx[,-c(5,6,7)],y=dy$y)
d$pp <- factor(d$pp,levels=c('xpos','xposInhom','xposHom'),
               labels=c('Experiment','Inhomogeneous','Homogeneous'))
table(d$pp)  

   Experiment Inhomogeneous   Homogeneous 
        57344         57344         57344 

The next figure shows all fixations of the three point processes (experimental data, inhomogeneous PP, homogeneous PP) during the first inspection of an image. As expected, fixation locations are not uniformly distributed in the experimental data and indicate inhomogeneity. The estimated intensity of all fixation locations is depicted by gray shading where darker areas represent higher intensities. The intensity of the experimental data was used for the simulation of the inhomogeneous PP.

dImg <- filter(d,image=='1' & nthPres=='First')
bw <- filter(sigma,img=='1' & pres=='First')$bw
ggplot(data=dImg,aes(x=x,y=y,col=pp)) +
  stat_density2d(aes(alpha=..density..),show.legend=FALSE,
                 geom="raster",contour=FALSE,h=bw*c(1,1),fill="black") +
  geom_point(size=.01) +
  facet_grid(.~pp) +
  labs(x='x-Coordinate [°]',y='y-Coordinate [°]',colour='Point Process') +
  coord_fixed(xlim=xrange,ylim=yrange,expand = FALSE) +
  theme_bw(base_size = 16) +
  theme(legend.position="none") 

Examples of an empirical and of the corresponding simulated scanpaths are visualized in the next figure. The estimated intensity of each point process is displayed with gray shading. Fixations are likely to be located in areas of high average intensity. Due to the independence of fixations, scanpaths of the inhomogeneous and homogeneous point processes reveal less systematic exploration behavior than the empirical data. Hence, saccade amplitudes increase considerably for these PPs.

dTrial <- filter(d,image=='1' & nthPres=='First' & id=='1')
dImg <- filter(d,image=='1' & nthPres=='First')
bw <- filter(sigma,img=='1' & pres=='First')$bw
ggplot(data=dTrial,aes(x=x,y=y,col=pp)) +
  stat_density2d(data=dImg,aes(alpha=..density..),show.legend=FALSE,
                 geom="raster",contour=FALSE,h=bw*c(1,1),fill="black") +
  geom_point(size=2) +  geom_path(size=.5) +
  labs(x='x-Coordinate [°]',y='y-Coordinate [°]',colour='Point Proces') +
  coord_fixed(xlim=xrange,ylim=yrange,expand=FALSE) +
  theme_bw(base_size=16) +
  facet_grid(.~pp,drop=TRUE) +
  theme(legend.position="none") 

Step 2. Choose optimal bandwidth for intensity estimation of PCF

Next we need to choose an optimal bandwidth for estimation of the intensity used to calculate the inhomogeneous PCF in Step 3. This is different from the bandwidth estimated in Step 1 to simulate the inhomogeneous point process. Since fixation locations in scanpaths of both control point processes are sampled independent from the preceding fixation history, average PCFs of both point processes are expected to reveal no spatial correlations, i.e. we expect the PCF \(g(r)\approx 1\) at all distances \(r\).

To choose an optimal bandwidth, we computed the deviation from complete spatial randomness (\(\int_a^b(g(r)-1)^2dr\)) of PCFs for a range of bandwidths for each image and presentation. We varied bandwidths from 0° to 10° in steps of 0.1° and computed the deviation from complete spatial randomness of the average PCF on an image. The average deviation at each bandwidth is plotted in the next figure. Lines represent individual images. For all images, the deviation is smallest in the center an increases for smaller and larger bandwidths. For all images the optimal bandwidth is between 1.5° and 5°. The bandwidth yielding the smallest deviation from complete spatial randomness on an image was chosen for the intensity estimation of the PCF in the next step.

# average PCF
mPCF <- function(pcf){
  r <- pcf$r[1,]
  trans <- sapply(data.frame(pcf$pcf),mean)
  mPCF <- data.frame(r,trans,row.names=NULL)
  return(mPCF)
}
# PCF deviation function
devPcfComp <- function(g,rmin,rmax){
  indx <- g$r>=rmin & g$r<=rmax
  dev <- sum((g$trans[indx]-1)^2)*(g$r[2]-g$r[1])
  return(dev)
}
# compute PCF for fixation sequence
fixSeqPCF <- function(fixSeq,den,p){
  
  # start, end, number of trials
  idx <- which(diff(as.numeric(fixSeq$id))!=0)
  i1 <- c(1,idx+1)
  i2 <- c(idx,nrow(fixSeq))
  Nid <- length(i1)
  numA <-  nrow(fixSeq)
  
  alvr <- NULL
  alvt <- NULL
  alvs <- NULL
  hdev <- NULL
  aid  <- NULL
  
  for (t in 1:Nid ) {
    scanpath <- fixSeq[i1[t]:i2[t],]
    
    if ( nrow(scanpath)>p$maxlen )  {
      scanpath <- scanpath[(1:p$maxlen),]
    } # maxlen
    
    if ( nrow(scanpath)>=p$minlen ) {
      numX <- nrow(scanpath)
      id <- unique(scanpath$id)
      X <- ppp(scanpath$x,scanpath$y,window=owin(xrange,yrange))
      g <- pcfinhom(X,den,kernel="epanechnikov",r=p$reval,divisor="r") # divisor="d" | "r"
      g$trans <-  g$trans*numA/numX
      ginteg <- devPcfComp(g,p$rmin,p$rmax)
      hdev <- c(hdev,ginteg)
      alvt <- rbind(alvt,g$trans)
      alvs <- rbind(alvs,rep(t,length(g$trans)))
      alvr <- rbind(alvr,g$r)
      aid <- c(aid,id)
    } # minlen
  } # id
  pcfSummary <- list(r=alvr,pcf=alvt,trial=alvs,dev=hdev,minlen=p$minlen,
                     maxlen=p$maxlen,rmin=p$rmin,rmax=p$rmax,nId=Nid,
                     excluded=Nid-length(hdev),id=aid)
  return(pcfSummary)
} # function
optimalBandwidthPcf <- function(bw,img,pres){
  pExp <- filter(d,pp=='Experiment' & image==img & nthPres==pres)
  ppExp <- ppp(x=pExp$x,y=pExp$y,window=owin(xrange,yrange))
  denExp <- density(ppExp,sigma=bw,edge=FALSE)
  denExp$v[denExp$v<epsilon] <- epsilon
  pInhom <- filter(d,pp=='Inhomogeneous' & image==img & nthPres==pres)
  
  pcfInhom <- fixSeqPCF(pInhom,denExp,p)
  meanPcfInhom <- mPCF(pcfInhom)
  dev <- devPcfComp(meanPcfInhom,p$rmin,p$rmax)
  type <- unique(pExp$type)
  
  return(data.frame(bw,img,pres,type,dev))
}

epsilon <- .Machine$double.eps # set smallest possible computable number
p <- NULL
p$minlen <- 10     # minimum length of a scanpath
p$maxlen <- 100    # maximum length of a scanpath, otherwise truncated
p$rmin <- .1       # evaluation window of PCF (minimum)
p$rmax <- 6.5      # evaluation window of PCF (maximum)
p$sigmaMin <- .1   # minimum bandwidth for optimal bandwidth evaluation  (Step 3)
p$sigmaMax <- 10   # maximum bandwidth for optimal bandwidth evaluation (Step 3)
p$sigmaSteps <- .1 # stepsize for optimal bandwidth evaluation (Step 3)
p$reval <- seq(0,7,length.out=513) # compute PCF from to
# evaluate PCFs for all bandwidths on all images and each presentation separately
eval_bw <- seq(p$sigmaMin,p$sigmaMax,by=p$sigmaSteps)
eval_img <- unique(d$image)
eval_pres <- unique(d$nthPres)
parOptBw <- cbind(data.frame(bw=rep(eval_bw,each=length(eval_img)),img=rep(eval_img,times=length(eval_bw))),pres=rep(eval_pres,each=length(eval_bw)*length(eval_img)))
cl <- makeForkCluster(detectCores()-1)
deviation <- parLapply(cl, seq_len(nrow(parOptBw)), function(i) {
  do.call(optimalBandwidthPcf, as.list(parOptBw[i,]))
})
stopCluster(cl)
# deviation of the inhomogeneous point process for each bandwith, image, and presentation
devInhom <- do.call(rbind,lapply(deviation,as.data.frame))
# save result
save(d,devInhom,p,xrange,yrange,file='devInhom.Rdata')
# load(file='devInhom.Rdata')
ggplot(devInhom,aes(x=bw,y=dev,group=img)) +
  geom_line() +
  facet_grid(.~type+pres) +
  labs(x='Bandwidth Sigma',y='Deviation') +
  coord_cartesian(ylim=c(0,5)) +
  theme_bw(base_size=16)

Step 3. Compute PCF for each trial

In the last step we compute the PCF of our empirical data. We use the optimal bandwidth that resulted in the smallest deviation from complete spatial randomness of the PCF of the inhomogeneous PP in Step 2. PCFs of individual scanpaths on an image are plotted in the next figure (left panel, gray lines). Three example scanpaths are highlighted in black. PCFs vary strongly between individual trials for all point processes. The average empirical PCF across all scanpaths on an image (red line) deviates strongly from complete spatial randomness, i.e. \(g(r) \neq 1\) for distances smaller than 4°. At distances beyond 4° the average PCF suggests independence of points, i.e., \(g(r)\approx 1\). Thus, fixations aggregate in close proximity during individual trials. Conversely, areas further away are fixated as predicted by chance, i.e., by the overall inhomogeneity observed across all participants. As expected, inspection of the control point processes demonstrates the absence of spatial correlations. The average PCF of the inhomogeneous and homogeneous point process are constant with \(g(r) \approx 1\). An artifact of the estimation procedure at short distances is present in all estimates.

computePcfImage <- function(img,pres,bw){
  pExp <- filter(d,pp=='Experiment' & image==img & nthPres==pres)
  ppExp <- ppp(x=pExp$x,y=pExp$y,window=owin(xrange,yrange))
  denExp <- density(ppExp,sigma=bw,edge=FALSE)
  pInhom <- filter(d,pp=='Inhomogeneous' & image==img & nthPres==pres)
  ppInhom <- ppp(x=pInhom$x,y=pInhom$y,window=owin(xrange,yrange))
  denInhom <- density(ppInhom,sigma=bw,edge=FALSE)
  pHom <- filter(d,pp=='Homogeneous' & image==img & nthPres==pres)
  ppHom <- ppp(x=pHom$x,y=pHom$y,window=owin(xrange,yrange))
  denHom <- density(ppHom,sigma=bw,edge=FALSE)
  
  pcfExp <- fixSeqPCF(pExp,denExp,p)
  pcfInhom <- fixSeqPCF(pInhom,denInhom,p)
  pcfHom <- fixSeqPCF(pHom,denHom,p)
  
  pcfImage <- rbind(data.frame(pp='Experiment',id=c(pcfExp$trial),
      r=c(pcfExp$r),pcf=c(pcfExp$pcf)),
    data.frame(pp='Inhomogeneous',id=c(pcfInhom$trial),
      r=c(pcfInhom$r),pcf=c(pcfInhom$pcf)),
    data.frame(pp='Homogeneous',id=c(pcfHom$trial),
      r=c(pcfHom$r),pcf=c(pcfHom$pcf)))
  
  type <- unique(pExp$type)
  
  return(data.frame(img,pres,type,pcfImage))
}
parPcf <- devInhom %>% 
  group_by(.,img,pres) %>%
  summarize(.,bw=bw[which.min(dev)])
cl <- makeForkCluster(detectCores()-1)
pcfDataSet <- parLapply(cl, seq_len(nrow(parPcf)), function(i) {
  do.call(computePcfImage, as.list(parPcf[i,]))
})
stopCluster(cl)
pcfDataSet <- do.call(rbind,lapply(pcfDataSet,as.data.frame))
pcfTrial <- pcfDataSet %>%
  filter(.,img==1 & pres=='First') %>%
  group_by(.,pp,id,r) %>%
  summarize(.,pcf=mean(pcf))
meanPcf <- pcfDataSet %>%
  group_by(.,pp,pres,type,r) %>%
  summarize(.,pcf=mean(pcf))
meanPcfImg <- pcfDataSet %>%
  group_by(.,pp,pres,type,img,r) %>%
  summarize(.,pcf=mean(pcf))

ggplot(pcfTrial,aes(x=r,y=pcf,col=pp)) +
  geom_line(aes(group=id),col='gray70') +
  geom_line(data=filter(pcfTrial,id%in%c(1,2,3)), aes(lty=as.factor(id)),col='black') +
  geom_line(stat='summary',fun.y='mean') +
  facet_grid(.~pp) +
  coord_cartesian(ylim=c(0,3)) +
  labs(x='Distance r [°]',y='Pair correlation function g(r)',colour='Point Process',lty='Participant') +
  theme_bw(base_size=16)

The same procedure can be repeated for each image. The next Figure shows the PCFs of each image (gray lines) as well as the average across all images (colored lines). While inhomogeneous and homogeneous point processes reveal no spatial correlations, empirical PCFs show spatial aggregation at short distances, \(r < 4\)° in all conditions. A detailed discussion of the results can be found in the manuscript.

ggplot(meanPcfImg,aes(x=r,y=pcf,col=pp)) +
  geom_line(aes(group=img),col='gray70') +
  geom_line(data=meanPcf) +
#  geom_line(data=pcfTrial,col='black',stat='summary',fun.y='mean') +
  facet_grid(type+pres~pp) +
  coord_cartesian(ylim=c(0.5,2)) +
  labs(x='Distance r [°]',y='Pair correlation function g(r)',colour='Point Process') +
  theme_bw(base_size=16)

References

Trukenbrod, H. A., Barthelmé, S., Wichmann, F. A., & Engbert, R. (2018). Spatial statistics for gaze patterns in scene viewing: Effects of repeated viewing. ArXiv. Retrieved from https://arxiv.org/abs/1704.01761