Code for generating figures

Packages

Figure 1: Properties of pupil dilation

Figure 2: Baseline correction and normalization

Figure 4: Averages of data

Figure 5: Partial effects (fixed effects only) of the initial GAMM model

Figure 7. Random factor smooths for participants and items estimated by model model1

Figure 8. Estimates of the initial GAMM model model1

Figure 9. Residuals of the initial GAMM model model1

Figure 10 (Autocorrelation in simulation data) is being discussed in the section on Simulations.

Figure 11. Improvement in model fit by adding random smooths for unique time series

Figure 12 (Determining the optimal value of \(\rho\)) is being discussed in the section on Data & Analysis.

Figure 13. Distribution of the data

Figure 14. Evaluation of the scaled-t model model6

Packages and data

## For running GAMM models:
library(mgcv)
packageVersion("mgcv")

## [1] '1.8.24'

## For visualizing and interpretating GAMM models:
library(itsadug)
packageVersion("itsadug")

## [1] '2.3'

## For aggregation:
library(plyr)
packageVersion("plyr")

## [1] '1.8.4'

## Functions for plotting:
library(plotfunctions)
packageVersion("plotfunctions")

## [1] '1.3'

library(colorspace)
packageVersion("colorspace")

## [1] '1.3.2'

## Define colors:
col1 <- 'pink1'
col2 <- 'black'
col3 <- 'indianred'

The data is available for download here.

## load data of Sentences 1 and 2:
load("Data/pupil_S1-S2_alignPronoun.rda")
head(dat0, 3)

##    Event Subject TrialCounter Item  Time Xgaze Ygaze Pupil baseline
## 1 s01.16     s01            3   16 -4010 515.1 389.8  1197     1143
## 2 s01.16     s01            3   16 -3990 515.2 390.0  1195     1143
## 3 s01.16     s01            3   16 -3970 514.7 388.1  1194     1143
##      Condition ImgType Congruency Context ImageMirror StartS1 NP1  NP2
## 1 A1.congruent   other  congruent      A1         Yes     500 995 1875
## 2 A1.congruent   other  congruent      A1         Yes     500 995 1875
## 3 A1.congruent   other  congruent      A1         Yes     500 995 1875
##   EndS1 StartS2 Pronoun EndS2 StartS1-alignPronoun EndS1-alignPronoun
## 1  2490    2960    4005  4980                -3505              -1515
## 2  2490    2960    4005  4980                -3505              -1515
## 3  2490    2960    4005  4980                -3505              -1515
##   StartS2-alignPronoun Pronoun-alignPronoun EndS2-alignPronoun
## 1                -1045                    0                975
## 2                -1045                    0                975
## 3                -1045                    0                975

# load data of Sentence 2 (re-baselined on onset of pronoun):
load("Data/pupil_S2_alignPronoun.rda")
head(dat,3)

##    Event Subject TrialCounter Item  Time Xgaze Ygaze Pupil baseline
## 1 s01.16     s01            3   16 -1510 502.1 400.1  1330     1143
## 2 s01.16     s01            3   16 -1490 501.7 400.2  1327     1143
## 3 s01.16     s01            3   16 -1470 502.0 400.8  1322     1143
##   p.baseline p.length    Condition ImgType Congruency Context ImageMirror
## 1     1195.5       10 A1.congruent   other  congruent      A1         Yes
## 2     1195.5       10 A1.congruent   other  congruent      A1         Yes
## 3     1195.5       10 A1.congruent   other  congruent      A1         Yes
##   StartS1 NP1  NP2 EndS1 StartS2 Pronoun EndS2 StartS1-alignPronoun
## 1     500 995 1875  2490    2960    4005  4980                -3505
## 2     500 995 1875  2490    2960    4005  4980                -3505
## 3     500 995 1875  2490    2960    4005  4980                -3505
##   EndS1-alignPronoun StartS2-alignPronoun Pronoun-alignPronoun
## 1              -1515                -1045                    0
## 2              -1515                -1045                    0
## 3              -1515                -1045                    0
##   EndS2-alignPronoun
## 1                975
## 2                975
## 3                975

Columns:

Event : unique combination of participants and items.
Subject : Participant number.
TrialCounter : Trial in the experiment, ranging between 3 and 66 (the first two trials were practice trials). Note that only the pronoun trials are included in the data, but the practice items and fillers were excluded.
Item : Combination of soundfiles and pictures. Each item has 4 variants, which are described by the columns Congruency and Context.
Time : Time with the trial from onset pronoun.
Xgaze, Ygaze : X and Y coordinates of the gaze position. Coordinates (0,0) indicate the topleft corner of the screen.
Pupil : pupil size measures, baseline corrected.
baseline : baselines (before onset sound files).
Condition : 4 conditions in the experiment, described by Context and Congruency.
ImgType : Picture with other-oriented or self-oriented action.
Congruency : Congruency between picture and sentence. If the image type was a self-oriented picture, the picture is incongruent with the pronoun sentence.
Context : Order or introduction of the two characters. A1 indicates that the actor was introduced first, A2 indicates that the actor was introduced second.
ImageMirror : Whether or not the image was flipped horizontally.
StartS1, NP1, NP2,EndS1 : Onsets of words in the context sentence and offset of the sentence, relative to the onset of the sound file.
StartS2, Pronoun, EndS2 : Onsets of words in the test sentence and offset of the sentence, relative to the onset of the sound file.
Column names ending with -alignPronoun: Sentence information, but then relative to the onset of the pronoun.

Figure 1: Properties of pupil dilation.

[up \(\uparrow\)]

The effect of cognitive processing on pupil dilation, as described by the pupil dilation function from Hoeks & Levelt (1993). The pupillary response is scaled to 0.5 mm for comparison (cf. Beatty & Lucero-Wagoner, 2000). The red vertical line T1 represents an event that triggers dilation (black solid line). The dashed line shows the adjusted dilation when a second event T2 shortly follows the first event.

## GAMMA ERLANG FUNCTION ########################
## Estimated by Hoeks & Levelt, 1993 (H&L93)

## x:    time stamps in msec
## tmax: average peak dilation; 
##       average of 930 is measured by H&L93
## n:    determines the shape of the Gamma Erlang function; 
##       estimated 10.1 is average determined by H&L93
## f:    scaling parameter

pupil <- function(x, tmax=930, n=10.1, f=1.27e-26){
  # f = 1/(2*max(y)) to get .5
  out <- rep(0, length(x))
  out[x>0] <- (x[x>0]^n)*exp((-1*n*x[x>0])/tmax)
  return(f*out)
}

## PLOT #########################################

par(cex=1.1)

plot(x <- seq(-200,3500), y1 <- pupil(x), type='n',
  main='Pupil dilation model',
  xlab='Time (ms)', ylab='Pupil dilation',
  ylim=c(0,.7), las=1, bty='n')

lines(x, pupil(x)+pupil(x-700), col=col2, lty=2, xpd=T)
lines(x, y1, lwd=2.5)
abline(v=c(0, 700), col=col3, lwd=c(2,1), lty=c(1,1))

text(x=c(0,700), y=.7, labels=c('T1','T2'), pos=3, col=col3, xpd=T)

Extra plot

Below the effects of luminance and field size, as modeled by Watson & Yellott (2012), are plotted. This graph provides an indication of the size of the pupil in mm.

Bottom panel

Example of two actual recorded pupil dilation time series, recorded from two different participants (solid versus dashed lines) in two different trials (red versus black lines). The data is re-aligned on the onset of the pronoun. The horizontal bars indicate the duration of the auditory stimuli (two spoken sentences) of the two trials.

## Selection data & Align time top onset sentence:
subdat <- droplevels(dat0[dat0$Event %in% c("s02.4", "s02.17", "s03.4", "s03.17"),])
subdat <- droplevels(subdat[,c("Event", "Time", "Pupil","Subject","Item",
                               "StartS1","EndS1","StartS2","Pronoun","EndS2")])
row.names(subdat) <- NULL
subdat$Time <- subdat$Time +subdat$Pronoun

## PLOT:

# select 4 trajectories:
subdat$col <- ifelse(subdat$Subject=='s02', col3, col2)
subdat$lty <- ifelse(subdat$Item=='17', 1, 4)
subdat$lwd <- ifelse(subdat$Item=='17', 1.5, 2)

# measures aligned on pronoun onset:
info <- subdat[,c("Subject","Item", "StartS1","EndS1","StartS2","Pronoun","EndS2")]
info <- info[!duplicated(info),]

# PLOT:

emptyPlot(c(0,7000), c(200,1700),
  main='Data examples', xlab='Time (ms)',ylab='Pupil size',
  las=1)

for(i in unique(subdat$Event)){
  col <- subdat[subdat$Event==i,]$col[1]
  lty <- subdat[subdat$Event==i,]$lty[1]
  with(subdat[subdat$Event==i,], lines(Time, Pupil, col=col, lwd=lwd, lty=lty))
}

addInterval(getCoords(c(.2, .1),side=2),
  lowVals=info[1:2,'StartS1'],
  highVals=info[1:2,"EndS1"], lwd=c(1.5,2), lty=c(1, 4), xpd=TRUE)
addInterval(getCoords(c(.2, .1),side=2),
  info[1:2,'StartS2'], info[1:2,'EndS2'], lwd=c(1.5,2), lty=c(1, 4), xpd=TRUE)

# legends
text(5500, getCoords(c(.2, .1),side=2), labels=c('item 4', 'item 17'), adj=0, cex=.85)
tmp <- with(subdat[(subdat$Time) > 6500 & (subdat$Time) <= 7000,],
  aggregate(list(Pupil=Pupil), list(Subject=Subject, col=col), mean))
text(7200, tmp$Pupil, labels=tmp$Subject, pos=4, col=tmp$col, font=2, cex=.85, xpd=TRUE)

Figure 2: Baseline correction and normalization

[up \(\uparrow\)]

The three plots are based on simulated data.

# FUNCTION
# generates modified sine waves
sineWave <- function(x,a=1,b=1,c=0){
    return(a*sin(b*x)+c*x)
}

# simulation data:
simdat <- expand.grid(Time = seq(0,2*pi, length=100),
  Trial=c(1:3) )
simdat$baseline <- ifelse(simdat$Trial==1, 1050, 650)
simdat$scale <- ifelse(simdat$Trial==2, 400, 250)
simdat$Pupil <- sineWave(simdat$Time/1.5, a=simdat$scale)+simdat$baseline

# baseline corrected pupil:
simdat$Pupil2 <- simdat$Pupil - simdat$baseline
# normalized pupil:
simdat$Pupil3 <- (simdat$Pupil - simdat$baseline) / simdat$baseline

Top panel

Three pupil dilation trials (simulated data) with Trial 1 and 3 differing in baseline, but showing the exact same pupillary response. Trials 2 and 3 share the same baseline, but Trial 2 shows a higher peak amplitude.

library(plotfunctions)

emptyPlot(range(simdat$Time)*1000, range(simdat$Pupil), h0=0,
  main="Data (simulated)", ylab="Pupil size", xlab="Time", las=1)
for(i in unique(simdat$Trial)){
  with(simdat[simdat$Trial==i,], lines(Time*1000, Pupil, lwd=2, col=1, lty=i))
  text(max(x)*1000, simdat[simdat$Trial==i & simdat$Time==max(x),'Pupil'], labels=i, 
    pos=3, font=2, cex=.85, xpd=TRUE)
}
points(c(0,0), c(1050,650), pch=16)

Bottomleft panel

The same three trials after the baseline is subtracted. The baselined data for Trials 1 and 3 overlap.

emptyPlot(range(simdat$Time)*1000, range(simdat$Pupil2), h0=0,
  main="Baselined data", ylab="Absolute pupil size change", xlab="Time", las=1)
for(i in unique(simdat$Trial)){
  with(simdat[simdat$Trial==i,], 
    lines(Time*1000, Pupil2, lwd=2, col=ifelse(i==1,alpha(1),1), lty=i))
}
points(c(0,0), c(0,0), pch=16)
legend('topright', legend=1:3, lwd=2, lty=1:3, col=1, 
bty="n", cex=.85, text.font=2)

Bottomright panel

The proportion pupil dilation change with respect to the baseline for same three trials. As the baseline of Trial 1 is much higher than of Trial 3, the pupil dilation change is much lower for Trial 1 than for Trial 3, although the measured pupillary response was exactly the same.

emptyPlot(range(simdat$Time)*1000, range(simdat$Pupil3), h0=0,
  main="Normalized data", ylab="% Pupil dilation", xlab="Time", las=1)
for(i in unique(simdat$Trial)){
  with(simdat[simdat$Trial==i,], lines(Time*1000, Pupil3, lwd=2, col=1, lty=i))
}
points(c(0,0), c(0,0), pch=16)
points(c(0,0), c(0,0), pch=16)
legend('topright', legend=1:3, lwd=2, lty=1:3, col=1, 
bty="n", cex=.85, text.font=2)

Figure 4: Averages of data

[up \(\uparrow\)]

Top panel

Example of two pupil dilation time series, recorded from two different participants (solid versus dashed lines) in two different trials (red versus black lines). The data is aligned on the onset of the pronoun. The horizontal bars indicate the duration of the auditory stimuli of the two trials, which consisted of two sentences. The baseline for the averages in this graph was calculated from a 250 ms time window before sound onset.

# define timebins of 100 ms
dat0$Timebin <- timeBins(dat0$Time, 100)

# calculate averages for each subject per timebin per condition
avg.subj <- ddply(dat0, c("Subject", "Timebin", "Context", "Congruency"), 
                  summarise,
                  Pupil = mean(Pupil-baseline, na.rm=TRUE))
# calculate averages and SE (capturing the variation between subject averages):
avg <- ddply(avg.subj, c("Timebin", "Context", "Congruency"), summarise,
             mean = mean(Pupil, na.rm=TRUE),
             se   =se(Pupil,na.rm = TRUE))


par(cex=1.1, mar=c(7.1,4.1,4.1,2.1))

emptyPlot(c(-3700,3000), c(0,250),
  main='Grand averages', xlab='Time aligned on pronoun onset (ms)', 
  ylab='Pupil size', v0=0, las=1)

avg$Condition <- with(avg, interaction(Context, Congruency))
avg$col <- ifelse(avg$Context=='A1', col2, col3)
avg$lty <- ifelse(avg$Congruency=='congruent', 1, 6)
avg$lwd <- ifelse(avg$Congruency=='congruent', 1, 2)
for(i in unique(avg$Condition)){
  col <- avg[avg$Condition==i,]$col[1]
  lty <- avg[avg$Condition==i,]$lty[1]
  with(avg[avg$Condition==i & avg$Timebin > -3800,], lines(Timebin, mean, col=col, lwd=lwd, lty=lty))
}


# calculate median sentence lengths
info <- dat0[,c("Item", "StartS1-alignPronoun", "EndS1-alignPronoun", "StartS2-alignPronoun", "EndS2-alignPronoun")]
info <- info[!duplicated(info),]
info2 <- apply(info[2:5], 2, median)

addInterval(lowVals=info2[1], highVals = info2[2], pos=250, lwd=2)
text(info2[1], 275, labels="S1", adj=0, xpd=TRUE)
addInterval(lowVals=info2[3], highVals = info2[4], pos=250, lwd=2)
text(info2[3], 275, labels="S2", adj=0, xpd=TRUE)

# mark analysis window:
lines(c(-500,2500), rep(getCoords(0, side=2),2), lwd=6, col=alpha(1), lend=1, xpd=TRUE)
rug(x=c(-500,2500), side=1)

# add legend:
legend_margin('bottomright',
  legend=c('a1 congr.', "A1 incon.", "A2 congr.", "A2 incon."),
  lty=rep(c(1,6), 2), lwd=rep(c(1,2), 2),
  col=rep(c(col2, col3), each=2), seg.len=1.5,
  bty='n', cex=.85, ncol=2)

Bottom panel

The grand averages for the four conditions. In contrast with the plot above, the baseline window of this analysis data was at the pronoun onset, indicated with the vertical line.

# define timebins of 100 ms
dat$Timebin <- timeBins(dat$Time, 100)

# calculate averages for each subject per timebin per condition
avg.subj <- ddply(dat, c("Subject", "Timebin", "Context", "Congruency"), 
                  summarise,
                  Pupil = mean(Pupil-p.baseline, na.rm=TRUE))
# calculate averages and SE (capturing the variation between subject averages):
avg <- ddply(avg.subj, c("Timebin", "Context", "Congruency"), summarise,
             mean = mean(Pupil, na.rm=TRUE),
             se   =se(Pupil,na.rm = TRUE))


par(cex=1.1, mar=c(7.1,4.1,4.1,2.1))

emptyPlot(c(-1000,3000), c(-50,125),
  main='Grand averages', xlab='Time aligned on pronoun onset (ms)', 
  ylab='Pupil size', v0=0, h=0, las=1)

avg$Condition <- with(avg, interaction(Context, Congruency))
avg$col <- ifelse(avg$Context=='A1', col2, col3)
avg$lty <- ifelse(avg$Congruency=='congruent', 1, 6)
avg$lwd <- ifelse(avg$Congruency=='congruent', 1, 2)
for(i in unique(avg$Condition)){
  col <- avg[avg$Condition==i,]$col[1]
  lty <- avg[avg$Condition==i,]$lty[1]
  with(avg[avg$Condition==i & avg$Timebin > -1000,], lines(Timebin, mean, col=col, lwd=lwd, lty=lty))
}

# calculate median sentence lengths
info <- dat[,c("Item", "StartS2-alignPronoun", "EndS2-alignPronoun")]
info <- info[!duplicated(info),]
info2 <- apply(info[2:3], 2, median)

addInterval(lowVals=info2[1], highVals = info2[2], pos=100, lwd=2)
text(info2[1], 125, labels="S2", adj=0, xpd=TRUE)

# mark analysis window:
lines(c(-500,2500), rep(getCoords(0, side=2),2), lwd=6, col=alpha(1), lend=1, xpd=TRUE)
rug(x=c(-500,2500), side=1)

# legend:
legend_margin('bottomright',
  legend=c('A1 congr.', "A1 incon.", "A2 congr.", "A2 incon."),
  lty=rep(c(1,6), 2), lwd=rep(c(1,2), 2),
  col=rep(c(col2, col3), each=2), seg.len=1.5,
  bty='n', cex=.85, ncol=2)

Figure 5: Partial effects (fixed effects only) of the initial GAMM model

[up \(\uparrow\)]

Note that this inital GAMM model does not provide a good fit of the data, as explained in the paper. For more information on the model, see the Data & Analysis section.

# load model (see analysis section):
load("R-output/m0.rda")
summary(m0)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Pupil_base2 ~ Condition + s(Time, by = Condition, k = 20) + s(Xgaze, 
##     Ygaze) + s(Time, Subject, bs = "fs", m = 1) + s(Time, Item, 
##     bs = "fs", m = 1)
## 
## Parametric coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                1.781      8.889    0.20    0.841    
## ConditionA2.congruent     17.039      1.034   16.48   <2e-16 ***
## ConditionA1.incongruent   18.172      1.069   17.00   <2e-16 ***
## ConditionA2.incongruent   23.059      1.067   21.61   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                     edf  Ref.df      F  p-value    
## s(Time):ConditionA1.congruent     5.149   6.260  4.972 2.61e-05 ***
## s(Time):ConditionA2.congruent     7.159   8.854  8.797 5.50e-13 ***
## s(Time):ConditionA1.incongruent   7.196   8.903 13.477  < 2e-16 ***
## s(Time):ConditionA2.incongruent   6.732   8.315 13.908  < 2e-16 ***
## s(Xgaze,Ygaze)                   27.820  28.900 49.447  < 2e-16 ***
## s(Time,Subject)                 132.951 152.000 56.184  < 2e-16 ***
## s(Time,Item)                    213.731 287.000 19.691  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.236   Deviance explained = 24.1%
## fREML = 4.3251e+05  Scale est. = 9461.2    n = 72030

Top panels

The top four panels show the nonlinear regression lines for each of the four conditions with pointwise 95% confidence intervals, with the value of the parametric estimates for that condition on the right (red numbers).

par(mfrow=c(2,2), cex=1.1)

# determine the parametric effects to be added to the smooths:
add <- round( as.vector( coef(m0)[1:4] ) + c(0, rep(as.vector(coef(m0)[1]),3)) ,1)

for( i in 1:4){
    emptyPlot(c(-500,2500), c(-60,50), h0=0, v0=0, las=1)
    # extract partial effects:
  mt <- get_modelterm(m0, select=i, as.data.frame=TRUE)
  # plot regression lines with CI:
    plot_error(mt$Time, mt$fit, mt$se.fit, shade=TRUE, xpd=TRUE)
    # polish plots:
    mtext(m0$smooth[[i]]$by.level, side=3, line=.5, font=2, cex=1)
    text(2500,0,labels=add[i], pos=4, font=2, col=2, xpd=TRUE)
    mtext("Time", side=1, line=2.5, cex=.85, xpd=TRUE)
}

## Summary:
##  * Time : numeric predictor; with 30 values ranging from -490.000000 to 2490.000000.

## Summary:
##  * Time : numeric predictor; with 30 values ranging from -490.000000 to 2490.000000.

## Summary:
##  * Time : numeric predictor; with 30 values ranging from -490.000000 to 2490.000000.

## Summary:
##  * Time : numeric predictor; with 30 values ranging from -490.000000 to 2490.000000.

The package mgcv provides a simpler way to achieve the same plots (but with less options for graphical optimization):

par(mfrow=c(2,2))
for(i in 1:4){
  plot(m0, select=i, shade=TRUE, ylim=c(-60,50))
  abline(h=0, v=0, lty=c(1,3))
}

The sum of the parametric effects for each condition is added on the right of each plot. These numbers are based on the summary statistics of the model:

round(summary(m0)$p.table,2)

##                         Estimate Std. Error t value Pr(>|t|)
## (Intercept)                 1.78       8.89    0.20     0.84
## ConditionA2.congruent      17.04       1.03   16.48     0.00
## ConditionA1.incongruent    18.17       1.07   17.00     0.00
## ConditionA2.incongruent    23.06       1.07   21.61     0.00

Bottom panel

The bottom panel shows the interaction between Xgaze and Ygaze, and implements the effect of gaze position on the measured pupil size. Note that (0,0) represents the topleft corner of the screen.

# plot partial effects contour plot:
pvisgam(m0, view=c("Xgaze", "Ygaze"), select=5, 
        main="",xlab="",ylab="", xlim=c(0,1024), ylim=c(0, 768),
        color=terrain_hcl(50), nlevels=20, 
        zlim=c(-200,200), las=1, add.color.legend=FALSE)

## [1] "Tensor(s) to be plotted: s(Xgaze,Ygaze)"

# add label:
mtext(m0$smooth[[5]]$label, side=3, line=.5, font=2, cex=1)
# fade areas without supporting data:
fadeRug(m0$model$Xgaze, m0$model$Ygaze)
# add color legend & labels:
gradientLegend(c(-200,200), pos=.875,color=terrain_hcl(50))
mtext("Xgaze", side=1, line=2.5, cex=.85, xpd=TRUE)
mtext("Ygaze", side=2, line=3, cex=.85, xpd=TRUE)

The package mgcv provides a simpler way to achieve the same plot (but with less options for graphical optimization):

plot(m0, select=5, rug=TRUE)

plot(m0, select=5,scheme=2)

Figure 7. Random factor smooths for participants and items estimated by model `model1`.

[up \(\uparrow\)]

Random factor smooths for particpants:

par(cex=1.1, mar=c(3.5,3,2.5,3))

emptyPlot(c(-500,2500), c(-100,100), h0=0, v0=0, las=1)
# plot random effects:
inspect_random(m0, select=6, col=alpha(1), add=TRUE, xpd=TRUE)

## Summary:
##  * Time : numeric predictor; with 30 values ranging from -490.000000 to 2490.000000. 
##  * Subject : factor with 17 values; set to the value(s): s01, s02, s03, s04, s05, s06, s07, s08, s09, s10, ...

mtext(m0$smooth[[6]]$label, side=3, line=.5, font=2, cex=1)
mtext("Time", side=1, line=2.5, cex=.85, xpd=TRUE)

Random factor smooths for items:

par(cex=1.1, mar=c(3.5,3,2.5,3))

emptyPlot(c(-500,2500), c(-100,100), h0=0, v0=0, las=1)
# plot random effects:
inspect_random(m0, select=7, col=alpha(1), add=TRUE, xpd=TRUE)

## Summary:
##  * Time : numeric predictor; with 30 values ranging from -490.000000 to 2490.000000. 
##  * Item : factor with 32 values; set to the value(s): 1, 10, 11, 12, 13, 14, 15, 16, 17, 18, ...

mtext(m0$smooth[[7]]$label, side=3, line=.5, font=2, cex=1)
mtext("Time", side=1, line=2.5, cex=.85, xpd=TRUE)

The package mgcv provides a simpler way to achieve the same plots (but with less options for graphical optimization):

plot(m0, select=6)

plot(m0, select=7)

Figure 8. Estimates of the initial GAMM model `model1`

[up \(\uparrow\)]

Top panel

Summed effects for all conditions, with the random effects set to zero.

plot_smooth(m0, view="Time", cond=list(Condition="A2.congruent"), rm.ranef=TRUE,
  v0=0, col=col2, lwd=2, lty=6, rug=FALSE, se=0,
  main="Estimated effects", ylab="Pupil (baselined)", las=1,
  ylim=c(-40,80))

## Summary:
##  * Condition : factor; set to the value(s): A2.congruent. 
##  * Time : numeric predictor; with 30 values ranging from -490.000000 to 2490.000000. 
##  * Xgaze : numeric predictor; set to the value(s): 517.7. 
##  * Ygaze : numeric predictor; set to the value(s): 328.4. 
##  * Subject : factor; set to the value(s): s03. (Might be canceled as random effect, check below.) 
##  * Item : factor; set to the value(s): 3. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(Time,Subject),s(Time,Item)
##

plot_smooth(m0, view="Time", cond=list(Condition="A2.incongruent"), rm.ranef=TRUE,
  v0=0, col=col3, lwd=2, lty=6, rug=FALSE, add=TRUE, xpd=TRUE, se=0)

## Summary:
##  * Condition : factor; set to the value(s): A2.incongruent. 
##  * Time : numeric predictor; with 30 values ranging from -490.000000 to 2490.000000. 
##  * Xgaze : numeric predictor; set to the value(s): 517.7. 
##  * Ygaze : numeric predictor; set to the value(s): 328.4. 
##  * Subject : factor; set to the value(s): s03. (Might be canceled as random effect, check below.) 
##  * Item : factor; set to the value(s): 3. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(Time,Subject),s(Time,Item)
##

plot_smooth(m0, view="Time", cond=list(Condition="A1.congruent"), rm.ranef=TRUE,
  v0=0, col=col2, lwd=1, lty=1, rug=FALSE, add=TRUE, xpd=TRUE, se=0)

## Summary:
##  * Condition : factor; set to the value(s): A1.congruent. 
##  * Time : numeric predictor; with 30 values ranging from -490.000000 to 2490.000000. 
##  * Xgaze : numeric predictor; set to the value(s): 517.7. 
##  * Ygaze : numeric predictor; set to the value(s): 328.4. 
##  * Subject : factor; set to the value(s): s03. (Might be canceled as random effect, check below.) 
##  * Item : factor; set to the value(s): 3. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(Time,Subject),s(Time,Item)
##

plot_smooth(m0, view="Time", cond=list(Condition="A1.incongruent"), rm.ranef=TRUE,
  v0=0, col=col3, lwd=1, lty=1, rug=FALSE, add=TRUE, xpd=TRUE, se=0)

## Summary:
##  * Condition : factor; set to the value(s): A1.incongruent. 
##  * Time : numeric predictor; with 30 values ranging from -490.000000 to 2490.000000. 
##  * Xgaze : numeric predictor; set to the value(s): 517.7. 
##  * Ygaze : numeric predictor; set to the value(s): 328.4. 
##  * Subject : factor; set to the value(s): s03. (Might be canceled as random effect, check below.) 
##  * Item : factor; set to the value(s): 3. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(Time,Subject),s(Time,Item)
##

# legend
legend('bottomright',
       legend=c('a1, congr.', "a1, incon.", "a2, congr.", "a2, incon."),
       lty=rep(c(1,6), each=2), lwd=rep(c(1,2), each=2),
       col=rep(c(col2, col3), 2), seg.len=1.5,
       bty='n', cex=.85, ncol=2, xpd=TRUE)

Bottom panel

Difference curves, derived from model1. The gray solid line represents the estimated difference (and pointwise 95% confidence intervals) between the incongruent and congruent items when the actor is introduced first (‘a1’), and the dashed red line represents the estimated difference (and pointwise 95% confidence intervals) between the incongruent and congruent items when the actor is introduced second (‘a2’).

plot_diff(m0, view="Time", comp=list(Condition=c("A1.incongruent", "A1.congruent")),
          rm.ranef=TRUE, main="Incongruent - Congruent", 
          las=1, ylab="Est. difference Pupil", 
          col.diff = alpha(1, f=0), # remove difference marking
          v0=0, col=col2)

## Summary:
##  * Time : numeric predictor; with 100 values ranging from -490.000000 to 2490.000000. 
##  * Xgaze : numeric predictor; set to the value(s): 517.7. 
##  * Ygaze : numeric predictor; set to the value(s): 328.4. 
##  * Subject : factor; set to the value(s): s03. (Might be canceled as random effect, check below.) 
##  * Item : factor; set to the value(s): 3. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(Time,Subject),s(Time,Item)
##

## 
## Time window(s) of significant difference(s):
##  21.717172 - 2128.787879

plot_diff(m0, view="Time", comp=list(Condition=c("A2.incongruent", "A2.congruent")),
          col.diff = alpha(1, f=0), # remove difference marking
          rm.ranef=TRUE, add=TRUE, col=col3, lty=2)

## Summary:
##  * Time : numeric predictor; with 100 values ranging from -490.000000 to 2490.000000. 
##  * Xgaze : numeric predictor; set to the value(s): 517.7. 
##  * Ygaze : numeric predictor; set to the value(s): 328.4. 
##  * Subject : factor; set to the value(s): s03. (Might be canceled as random effect, check below.) 
##  * Item : factor; set to the value(s): 3. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(Time,Subject),s(Time,Item)
##  
## 
## Time window(s) of significant difference(s):
##  593.636364 - 2098.686869
##  2429.797980 - 2490.000000

# legend:
legend('topright', legend=c("a1","a2"),
       col=c(col2, col3), lwd=1, lty=2, fill=alpha(c(col2, col3)),
       border=NA, merge=TRUE, bty='n')

Figure 9. Residuals of the initial GAMM model `model1`

[up \(\uparrow\)]

Left panel

acf(resid(m0), bty='n', main="ACF residuals model1")

Center panel

qqnorm(resid(m0), pch=16, col=alpha(1, f=.25), cex=.5, bty='n')
qqline(resid(m0), col=col2)

Right panel

emptyPlot(range(fitted(m0)),range(resid(m0)),  
   h0=0, v0=0,
   main='Model fit', xlab='Fitted effects', ylab='Residuals')
points(fitted(m0), residuals(m0), pch=16, cex=.5, col=alpha(1,f=.25))

Figure 11. Improvement in model fit by adding random smooths for unique time series

[up \(\uparrow\)]

Loading the new GAMM model and selecting the four time series from Figure 1:

# load model (see analysis section):
load("R-output/m2.rda")
# summary(m2, re.test=FALSE)

ev <- c("s02.4", "s02.17", "s03.4", "s03.17")

Left panel

Data (black thick lines) and the model fit of the initial GAMM model (thick red lines) for three random three events in the experiment.

# combine data and fitted estimates
tmp       <- m0$model
tmp$Event <- interaction(tmp$Subject, tmp$Item, drop=TRUE)
tmp$fit   <- fitted(m0)
tmp       <- droplevels(tmp[tmp$Event %in% ev,])

emptyPlot(range(tmp$Time), range(tmp[,1]),
          main="Examples (model 1)", xlab="Time (ms)", ylab="Pupil size", 
          v0=0, h0=0)
for( i in unique(tmp$Event)){
  with(tmp[tmp$Event==i,], lines(Time, Pupil_base2, col=1,  xpd=TRUE))
  with(tmp[tmp$Event==i,], lines(Time, fit, col=alpha('indianred',f=.75), lwd=3, xpd=TRUE))
}

Center panel

Data (black thick lines) and the model fit of the improved GAMM model (thick red lines) for the same three time series.

tmp <- m2$model
tmp$fit <- fitted(m2)
tmp <- droplevels(tmp[tmp$Event %in% ev,])

emptyPlot(range(tmp$Time), range(tmp[,1]),
          main="Examples (model 4)", xlab="Time (ms)", ylab="Pupil size", 
          v0=0, h0=0)
for( i in unique(tmp$Event)){
  with(tmp[tmp$Event==i,], lines(Time, Pupil_base2, col=1,  xpd=TRUE))
  with(tmp[tmp$Event==i,], lines(Time, fit, col=alpha('indianred',f=.75), lwd=3, xpd=TRUE))
}

Right panel

ACF for the improved GAMM model, model4. (The ACF of model1 is presented in Figure 9).

acf(resid(m2), bty='n', main="ACF residuals model4")

Figure 13. Distribution of the data

[up \(\uparrow\)]

load("R-output/m5rho.rda")

Left panel

## same as: 
# qqnorm(m5rho$model[,1])
# qqline(m5rho$model[,1])

qq <- qqnorm(m5rho$model[,1], plot=FALSE)
emptyPlot(range(qq$x), c(-750,750))
mtext("QQ-plot", side=3, line=.5, font=2, cex=1)

points(qq$x, qq$y, pch=16, col=alpha(1), cex=.65)
qqline(m5rho$model[,1], col=1, lwd=.5)

Right panel

d1 <- density(m5rho$model[,1])
emptyPlot(c(-750,750), c(0, max(d1$y)), h0=0, las=1, axes=FALSE)
axis(1)
lines(d1$x, d1$y, lwd=1.5)
mtext("Density data", side=3, line=.5, font=2, cex=1)

Figure 14. Evaluation of the scaled-t model `model6`

[up \(\uparrow\)]

Load scaled-t model:

load("R-output/m4rho.rda")

Top row

Estimated effects (Left panel) and estimated differences with pointwise 95% confidence intervals (Center and Right panels).

Left panel: Estimates

plot_smooth(m4rho, view="Time", cond=list(Condition="A2.congruent"), rm.ranef=TRUE,
  v0=0, col=col2, lwd=2, lty=6, rug=FALSE, se=0,
  main="Estimated effects", ylab="Pupil (baselined)", las=1,
  ylim=c(-40,80))

## Summary:
##  * Condition : factor; set to the value(s): A2.congruent. 
##  * Time : numeric predictor; with 30 values ranging from -490.000000 to 2490.000000. 
##  * Xgaze : numeric predictor; set to the value(s): 517.7. 
##  * Ygaze : numeric predictor; set to the value(s): 328.4. 
##  * Event : factor; set to the value(s): s01.3. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(Time,Event)
##

plot_smooth(m4rho, view="Time", cond=list(Condition="A2.incongruent"), rm.ranef=TRUE,
  v0=0, col=col3, lwd=2, lty=6, rug=FALSE, add=TRUE, xpd=TRUE, se=0)

## Summary:
##  * Condition : factor; set to the value(s): A2.incongruent. 
##  * Time : numeric predictor; with 30 values ranging from -490.000000 to 2490.000000. 
##  * Xgaze : numeric predictor; set to the value(s): 517.7. 
##  * Ygaze : numeric predictor; set to the value(s): 328.4. 
##  * Event : factor; set to the value(s): s01.3. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(Time,Event)
##

plot_smooth(m4rho, view="Time", cond=list(Condition="A1.congruent"), rm.ranef=TRUE,
  v0=0, col=col2, lwd=1, lty=1, rug=FALSE, add=TRUE, xpd=TRUE, se=0)

## Summary:
##  * Condition : factor; set to the value(s): A1.congruent. 
##  * Time : numeric predictor; with 30 values ranging from -490.000000 to 2490.000000. 
##  * Xgaze : numeric predictor; set to the value(s): 517.7. 
##  * Ygaze : numeric predictor; set to the value(s): 328.4. 
##  * Event : factor; set to the value(s): s01.3. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(Time,Event)
##

plot_smooth(m4rho, view="Time", cond=list(Condition="A1.incongruent"), rm.ranef=TRUE,
  v0=0, col=col3, lwd=1, lty=1, rug=FALSE, add=TRUE, xpd=TRUE, se=0)

## Summary:
##  * Condition : factor; set to the value(s): A1.incongruent. 
##  * Time : numeric predictor; with 30 values ranging from -490.000000 to 2490.000000. 
##  * Xgaze : numeric predictor; set to the value(s): 517.7. 
##  * Ygaze : numeric predictor; set to the value(s): 328.4. 
##  * Event : factor; set to the value(s): s01.3. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(Time,Event)
##

legend('bottomright',
  legend=c('A1, congr.', "A1, incon.", "A2, congr.", "A2, incon."),
  lty=rep(c(1,6), each=2), lwd=rep(c(1,2), each=2),
  col=rep(c(col2, col3), 2), seg.len=1.5,
  bty='n', cex=.85, ncol=2, xpd=TRUE)

Center panel: Difference plot A1

plot_diff(m4rho, view="Time", 
          comp=list(Condition=c("A1.incongruent", "A1.congruent")), rm.ranef=TRUE, 
          main="A1: Incongruent - Congruent", las=1, ylab="Est. difference Pupil",
          col.diff = alpha(1, f=0), # remove difference marking
          v0=0, col=col2)

## Summary:
##  * Time : numeric predictor; with 100 values ranging from -490.000000 to 2490.000000. 
##  * Xgaze : numeric predictor; set to the value(s): 517.7. 
##  * Ygaze : numeric predictor; set to the value(s): 328.4. 
##  * Event : factor; set to the value(s): s01.3. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(Time,Event)
##

## 
## Time window(s) of significant difference(s):
##  533.434343 - 1045.151515

legend('topright', legend=c("a1"),
       col=c(col2), lwd=1, lty=2, fill=alpha(c(col2)),
       border=NA, merge=TRUE, bty='n')

Right panel: Difference plot A2

plot_diff(m4rho, view="Time", 
          comp=list(Condition=c("A2.incongruent", "A2.congruent")), rm.ranef=TRUE, 
          main="A2: Incongruent - Congruent", las=1, ylab="Est. difference Pupil", 
          col.diff = alpha(1, f=0), # remove difference marking
          v0=0, col=col3)

## Summary:
##  * Time : numeric predictor; with 100 values ranging from -490.000000 to 2490.000000. 
##  * Xgaze : numeric predictor; set to the value(s): 517.7. 
##  * Ygaze : numeric predictor; set to the value(s): 328.4. 
##  * Event : factor; set to the value(s): s01.3. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(Time,Event)
##

## 
## Difference is not significant.

legend('bottomright', legend=c("a2"),
       col=c(col3), lwd=1, lty=2, fill=alpha(c(col3)),
       border=NA, merge=TRUE, bty='n')

Bottom row

Residuals of model model6.

Left panel: ACF plot

acf_resid(m4rho, bty='n', main="ACF residuals model1")

Note that this are standardixed residuals, corrected for the AR1 model that was included. Compare this ACF plot with the ACF of the uncorrected residuals:

acf(resid(m4rho), bty='n', main="ACF uncorrected residuals")

Center panel: QQ-plot

qq.gam(m4rho, pch=16, col=alpha(1, f=.25), 
  cex=.5, bty='n', main="QQ-plot")

Right panel: Model fit

emptyPlot(range(fitted(m4rho)),range(residuals(m4rho)),  
   h0=0, v0=0,
   main='Model fit', xlab='Fitted effects', ylab='Residuals')
points(fitted(m4rho), residuals(m4rho), pch=16, cex=.5, col=alpha(1,f=.25))

References

Beatty, J. and Lucero-Wagoner, B. (2000). The pupillary system. In: Cacioppo, J.T., Tassinary, L.G. and Berntson, G.G. (eds.) Handbook of psychophysiology. New York: Cambridge University Press., pp. 142–162.
Hoeks, B. and Levelt, W.J.M. (1993). Pupillary dilation as a measure of attention: A quantitative system analysis. Behavior Research Methods, 25(1): 16–26.
Watson, A.B. and Yellott, J.I. (2012). A unified formula for light- adapted pupil size. Journal of Vision, 12(10): 1–16.

Pupil Size AnalysisSupplementary Materials

Jacolien van Rij

November 2018

Code for generating figures

Packages and data

Figure 1: Properties of pupil dilation.

Extra plot

Bottom panel

Figure 2: Baseline correction and normalization

Top panel

Bottomleft panel

Bottomright panel

Figure 4: Averages of data

Top panel

Bottom panel

Figure 5: Partial effects (fixed effects only) of the initial GAMM model

Top panels

Bottom panel

Figure 7. Random factor smooths for participants and items estimated by model model1.

Figure 8. Estimates of the initial GAMM model model1

Top panel

Bottom panel

Figure 9. Residuals of the initial GAMM model model1

Left panel

Center panel

Right panel

Figure 11. Improvement in model fit by adding random smooths for unique time series

Left panel

Center panel

Right panel

Figure 13. Distribution of the data

Left panel

Right panel

Figure 14. Evaluation of the scaled-t model model6

Top row

Left panel: Estimates

Center panel: Difference plot A1

Right panel: Difference plot A2

Bottom row

Left panel: ACF plot

Center panel: QQ-plot

Right panel: Model fit

References

Pupil Size Analysis
Supplementary Materials

Figure 7. Random factor smooths for participants and items estimated by model `model1`.

Figure 8. Estimates of the initial GAMM model `model1`

Figure 9. Residuals of the initial GAMM model `model1`

Figure 14. Evaluation of the scaled-t model `model6`