Analyzing pupil dilation data with Generalized Additive Models - - PowerPoint PPT Presentation

Jacolien van Rij & Martijn Wieling

Analyzing pupil dilation data with Generalized Additive Models (GAMs)

Thu, June 27, 2013 Ÿ LOT School 2013, Groningen

Pupil dilation

• Primarily determined by light and accommodation reflexes
• Also ‘psychophysiological index of dynamic brain activity in human

cognition.’ (Beatty & Lucero-Wagoner, 2000)

• slow movements
• small scale movements, about .5 mm in pupillary diameter
• First reports from 1875 in Germany, rediscovered in 1960s
• not precisely clear what it means: attention, effort, processing load,

working memory load, emotion, …

2

Characteristics

• Hoeks & Levelt (1993): Estimation of pupil dilation response by

mathematical function

• peak around 930 ms after trigger

pupil dilation from baseline

1000 2000 3000 4000 5000 0.05 0.1

Time (ms)

pupil dilation from baseline

1000 2000 3000 4000 5000 0.05 0.1

Time (ms)

pupil dilation from baseline

1000 2000 3000 4000 5000 0.05 0.1

Time (ms)

pupil dilation from baseline

1000 2000 3000 4000 5000 0.05 0.1

Time (ms)

pupil dilation from baseline

1000 2000 3000 4000 5000 0.05 0.1

Time (ms)

peak latency peak dilation

Pupillometry in language processing research

• Pupil dilation is a precise and consistent measure of processing activity

during on-line language processing

• Sensitive to (among others):
• grammatical complexity (Just & Carpenter, 1993)
• integration of prosodic information and visual context (Engelhardt et

al., 2010)

• focus prosody (Zellin et al., 2011)

4

Overview

• Experimental study
• About the data
• Analyses using GAMs
• random effects structure
• testing fixed effects
• model criticism
• autocorrelation
• Discussion

5

Van Rij, 2012; Van Rij, Van Rijn, & Hendriks, in prep.

Experimental study

6

Experimental study

• Question: is the processing of object pronouns influenced by context?
• Example: The penguin is hitting him with a pan.
• Object pronouns are mainly guided by grammatical principles

(e.g. Binding Theory; Chomsky, 1981)

• In contrast to subject pronouns (he, she)

7

Why is this interesting?

• Investigates the underlying mechanism of sentence processing
• Debate:
• Initial-filter accounts: grammar first filters ungrammatical referents

(a.o., Nicol & Swinney, 1989; Clifton et al., 1997; Lewis et al., 2012) § Hypothesis: no early effect of context, late effect may be possible

• Competing-constraints accounts: grammar competes with other factors,

such as context (a.o., Badecker & Straub, 2002; Clackson et al., 2012)

§ Hypothesis: early effect of context

8

2x2 within-subjects design

9

• Test sentence: The penguin is hitting him with a pan.
• Manipulating linguistic context:
• AP (agent-patient): Here you see a penguin and a sheep.
• PA (patient-agent): Here you see a sheep and a penguin.
• Manipulating visual context:

✓ Congruent: picture matches sentence ✗ Incongruent: picture does not match sentence

Experimental study

• Question: is the processing of object pronouns influenced by context?
• Example: The penguin is hitting him with a pan.
• Task: Is the picture congruent with

this sentence?

• Participants see possible answers on

screen 2500 ms after sentence offset

10

About the data

11

Raw pupil data

12

15871914 480.5 327.1 1201.0 ... 15871918 480.2 327.0 1197.0 ... 15871922 479.9 327.3 1194.0 ... 15871926 479.9 327.8 1194.0 ... 15871930 479.7 328.5 1191.0 ... 15871934 479.5 330.0 1188.0 ... 15871938 479.5 330.1 1188.0 ... 15871942 479.3 330.0 1188.0 ... 15871946 479.6 330.8 1187.0 ... 15871950 480.6 331.1 1186.0 ... 15871954 481.5 331.3 1183.0 ... 15871958 481.2 329.9 1180.0 ... MSG MSG 15871960 15871960 >TrialStart TrialStart 50000 50000 3 3 15871962 480.6 329.1 1178.0 ... 15871966 479.5 328.8 1176.0 ... 15871970 479.1 329.4 1175.0 ... 15871974 480.1 328.4 1178.0 ...

• Eye tracker:

Eye Link 1000, monocular 250 Hz

• Markers / messages to

align trials and stimuli

¤time ¤x ¤y ¤pupil

Raw pupil data

• Pupil dilation is continuous

signal

• Slow changing signal, in

contrast with gaze data

• Preprocessing
• remove blinks
• (filter + down sampling)
• baseline + normalization

13

x-gaze

Timestamp

5233282 5242778 250 500 750 1000

y-gaze

Timestamp

5233282 5242778 250 500 750

pupil area

Timestamp

5233282 5242778 100 200 300 400 500

Timestamp

Linear interpolation of blinks

14

pupil area

Timestamp

5233282 5242778 250 500

Filtering, down sampling

15

pupil area

Timestamp

5233282 5242778 250 500 5235500 5236500 380 405 430

unfiltered filtered

• To avoid aliasing, low pass

filter of 25 Hz

• After filtering, down

sampled to 50 Hz (original frequency 250 Hz)

Baseline & normalization

• Pupil size can be measured in different units
• area, diameter, or arbitrary unit - depends on eye tracker / camera
• Converted to pupil dilation:
• pupil <- (pupil_size - baseline) / baseline
• baseline: average pupil size 50 ms before to 50 ms after onset pronoun
• Note: R package for preprocessing of pupil dilation data

16

Align data on pronoun onset

17

• Averaged subject means per condition (±1SE): lot of variation

Is object pronoun interpretation influenced by context?

Analyses using GAMs

18

Data structure

• 17 participants, 21-32 pronoun

trials each

• 32 items, 4 variants

(AP-congr, AP-incongr, PA-congr, PA-incongr)

• Time window: -200 ms before

pronoun onset to 3000 ms after pronoun onset

19

> str(dat) 'data.frame': 78907 obs. of 19 variables: $ Subject : Factor w/ 17 levels $ Item : Factor w/ 32 levels $ trial : int (3-66) $ Time : int $ pTime : int (-199 - 2983) $ trialTime : int $ Context : num (0=PA, 1=AP) $ Congruency : num (0=match, 1=mismatch) $ Interaction : num (Context * Congruency) $ Cond : Factor w/ 4 levels $ ActorRight : num (0=no, 1=yes) $ StartS2 : int $ StartVerb : int $ Anafoor : int $ EndAnafoor : int $ EndS2 : int $ dilationFiltered: num $ baseline : num $ Pupil : num (-.64 - .64) $ prev : num (first NA!)

Baseline model

> summary( m0 <- gam(Pupil ~ s(pTime), data=dat) ) ... Parametric coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.0136331 0.0003368 40.48 <2e-16 ***

• Approximate significance of smooth terms:

edf Ref.df F p-value s(pTime) 8.123 8.793 666.2 <2e-16 ***

• R-sq.(adj) = 0.0691 Deviance explained = 6.92%

fREML score = -74097 Scale est. = 0.008945 n = 78907

20

Baseline model

> m0 <- gam(Pupil ~ s(pTime), data=dat)

§ Sanity check: plot model smooths

> plot(m0, rug=F) > abline(h=0) > abline(v=c(0,1000), lty=3) > abline(h=-1*coef(m0)[1], col='red', lty=2)

21

500 1500 2500

• 0.04

0.00 0.04 pTime s(pTime,8.12) 500 1500 2500

• 0.04

0.00 0.04 pTime s(pTime,8.12)

¡intercept adjustment

Variation between subjects

22

500 1500 2500

• 0.05

0.00 0.05 0.10 500 1500 2500

• 0.05

0.00 0.05 0.10 500 1500 2500

• 0.05

0.00 0.05 0.10 500 1500 2500

• 0.05

0.00 0.05 0.10 500 1500 2500

• 0.05

0.00 0.05 0.10 500 1500 2500

• 0.05

0.00 0.05 0.10 500 1500 2500

• 0.05

0.00 0.05 0.10 500 1500 2500

• 0.05

0.00 0.05 0.10 500 1500 2500

• 0.05

0.00 0.05 0.10 500 1500 2500

• 0.05

0.00 0.05 0.10 500 1500 2500

• 0.05

0.00 0.05 0.10 500 1500 2500

• 0.05

0.00 0.05 0.10 500 1500 2500

• 0.05

0.00 0.05 0.10 500 1500 2500

• 0.05

0.00 0.05 0.10 500 1500 2500

• 0.05

0.00 0.05 0.10 500 1500 2500

• 0.05

0.00 0.05 0.10 500 1500 2500

• 0.05

0.00 0.05 0.10

Different types of random effects with GAMs

1. Random intercept: s(Item, bs="re") 2. Random intercept + random slope: s(Item, pTime,bs="re") 3. Random wiggly curve: s(pTime, Subject, bs="fs", m=1)

• Additive effects:

23

s(pTime)

pTime 500 1500 2500

• 0.04

0.02

s(pTime,Subject)

pTime 500 1500 2500

• 0.05

0.05

s(pTime,Subject)

pTime 500 1500 2500

• 0.05

0.05 0.1 0.15

Pupil ~ s(Time) + s(Time, Subject, Pupil ~ s(Time) + s(Time, Subject, bs bs=" ="fs fs", m=1) ", m=1) Pupil ~ s(Time, Subject, Pupil ~ s(Time, Subject, bs bs=" ="fs fs", m=1) ", m=1)

Test random effects

> m0 <- bam(Pupil ~ s(pTime,k=10), data=dat) > m0a <- bam(Pupil ~ s(pTime,k=10) + s(pTime, Subject, bs='fs', m=1), data=dat) > AIC(m0)-AIC(m0a) [1] 5483.468 > m0b <- bam(Pupil ~ s(pTime, k=10) + s(pTime, Subject, bs='fs', m=1) + s(Item, bs='re'), data=dat) > AIC(m0a)-AIC(m0b) [1] 2756.84

24

Random effects – use AIC

Random effects

> summary(m0b <- bam(Pupil ~ s(pTime) + s(pTime, Subject, bs='fs', m=1) + s(Item, bs='re'), data=dat)) Parametric coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.034755 0.005342 6.506 7.75e-11 ***

• Approximate significance of smooth terms:

edf Ref.df F p-value s(pTime) 7.111 7.587 22.42 <2e-16 *** s(pTime,Subject) 124.549 152.000 53.26 <2e-16 *** s(Item) 30.668 31.000 91.43 <2e-16 ***

• R-sq.(adj) = 0.163 Deviance explained = 16.5%

fREML score = -78028 Scale est. = 0.0080423 n = 78907

25

Visualize random effects

26

1000 2000 3000

• 0.04

0.02 0.08 pTime s(pTime,Subject,115.5)

random wiggly curves

• 2
• 1

1 2

• 0.04
• 0.01

0.02

s(Item,30.9)

Gaussian quantiles effects

Testing for fixed effects: complete model

Parametric coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.035201 0.005381 6.541 6.14e-11 ***

• Approximate significance of smooth terms:

edf Ref.df F p-value s(pTime) 6.867 7.373 21.014 < 2e-16 *** s(pTime):Context 6.525 7.578 9.354 3.15e-12 *** s(pTime):Congruency 5.069 5.978 12.546 4.10e-14 *** s(pTime):Interaction 7.412 8.459 11.918 < 2e-16 *** s(pTime):ActorRight 5.486 6.489 16.426 < 2e-16 *** ti(pTime,trial) 14.851 15.820 61.259 < 2e-16 *** s(pTime,Subject) 124.510 152.000 52.821 < 2e-16 *** s(Item) 30.670 31.000 92.487 < 2e-16 ***

• R-sq.(adj) = 0.178 Deviance explained = 18%

fREML score = -78654 Scale est. = 0.0078961 n = 78907

27

Different types of smooth functions

§ s():

• for 1-dimensional smooths, e.g. s(Time)
• for more dimensional smooths with predictors on the same scale,

e.g. s(Longitude, Latitude)

• for random effects (do not use 'by='), e.g., s(Item, bs="re"), or

s(Time, Subject, bs="fs", m=1) § ti(): for more dimensional smooths; models the pure interaction effect, does not include main effect, e.g. ti(Time, Trial, d=c(1,1)) § te(): for more dimensional smooths; models the full interaction, and includes main effects, e.g. te(Time, Trial, d=c(1,1)) (More information about this in lab session 4)

28

Inspection of the model

• The main effect reflect the baseline condition

29

main effect pTime

500 1500 2500

• 0.04
• 0.02

0.00 0.02 0.04

Inspection of the model

• Binary predictors reflect difference between condition 1 and 0

30

Context

500 1500 2500

• 0.04

0.00 0.04

Congruency

500 1500 2500

• 0.04

0.00 0.04

Interaction

500 1500 2500

• 0.04

0.00 0.04

AP: Here you see a penguin and a sheep. PA: Here you see a sheep and a penguin. AP: Here you see a penguin and a sheep.

Main effect of Time

• The main effect reflect the baseline condition

31

main effect pTime

500 1500 2500

• 0.04
• 0.02

0.00 0.02 0.04

PA: Here you see a sheep and a penguin.

Inspection of the model

• Binary predictors reflect difference between condition 1 and 0

32

Context

500 1500 2500

• 0.04

0.00 0.04

Congruency

500 1500 2500

• 0.04

0.00 0.04

Interaction

500 1500 2500

• 0.04

0.00 0.04

Context

500 1500 2500

• 0.04

0.00 0.04

Congruency

500 1500 2500

• 0.04

0.00 0.04

Interaction

500 1500 2500

• 0.04

0.00 0.04

AP: Here you see a penguin and a sheep. PA: Here you see a sheep and a penguin. AP: Here you see a penguin and a sheep.

Inspection of the model

• Fixed effects predictions:

33

1000 2000 3000

• 0.04

0.00 0.04

PA-con PA-inc AP-con AP-inc

Inspection of the model

• Partial effect of time and trial:

34

500 1500 2500

• 30
• 10

10 20 30 pTime trial

Remarks

• Alternatively: Use categorical variables (factor levels) instead of binary
• tests whether factor levels are different from zero, not whether and

where factors levels are different

• Include all predictors within the same model
• a step-wise modeling procedure using the residuals of other models will

inflate the error terms, yielding too liberal results

• the direction and shape of effects will not be affected

35

Model criticism: check residuals

> gam.check(m1)

36

Model criticism: gam.check

Method: fREML Optimizer: perf newton full convergence after 11 iterations. Gradient range [-1.858345e-06,1.144776e-06] (score -78653.69 & scale 0.007896092). Hessian positive definite, eigenvalue range [1.072725,39448.08]. Basis dimension (k) checking results. Low p-value (k-index<1) may indicate that k is too low, especially if edf is close to k'. k' edf k-index p-value s(pTime) 9.000 6.867 1.009 0.76 s(pTime):Context 10.000 6.525 1.009 0.72 s(pTime):Congruency 10.000 5.069 1.009 0.72 s(pTime):Interaction 10.000 7.412 1.009 0.70 s(pTime):ActorRight 10.000 5.486 1.009 0.72 ti(pTime,trial) 16.000 14.851 0.854 0.00 s(pTime,Subject) 152.000 124.510 1.009 0.70 s(Item) 32.000 30.670 NA NA

37

Check autocorrelation

> acf(resid(m1), main='ACF of res model1')

38

10 20 30 40 0.0 0.4 0.8 Lag ACF

ACF of res model1

HUGE!

Reduce autocorrelation

> rho = .99 # always order data first > dat <- dat[order(dat$Subject, dat$trial, dat$pTime, decreasing=F),] > m2 <- bam(Pupil ~ s(pTime) + s(pTime, by=Context) + s(pTime, by=Congruency) + s(pTime, by=Interaction) + s(pTime, by=ActorRight) + ti(pTime, trial, k=10) + s(pTime, Subject, bs='fs', m=1)+s(Item, bs='re'), data=dat, rho rho=rho rho)

39

Summary model 2

Parametric coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.06766 0.01025 6.602 4.09e-11 ***

• Approximate significance of smooth terms:

edf Ref.df F p-value s(pTime) 7.271 7.761 14.988 < 2e-16 *** s(pTime):Context 7.237 8.504 20.658 < 2e-16 *** s(pTime):Congruency 4.266 5.282 38.187 < 2e-16 *** s(pTime):Interaction 7.002 8.290 4.401 1.98e-05 *** s(pTime):ActorRight 4.898 6.140 102.027 < 2e-16 *** ti(pTime,trial) 49.405 62.512 57.798 < 2e-16 *** s(pTime,Subject) 118.970 152.000 72.180 < 2e-16 *** s(Item) 30.761 31.000 162.778 < 2e-16 ***

• R-sq.(adj) = -0.123 Deviance explained = -12%

fREML score = -2.502e+05 Scale est. = 0.0051284 n = 78907

40

before: 0.035

Normalized residuals

> dat$res.norm <- dat$res.t1 - rho * dat$res.t0 > acf(dat[!is.na(dat$res.norm),]$res.norm)

41

10 20 30 40 0.0 0.4 0.8 Lag ACF

ACF of res model1

10 20 30 40 0.0 0.4 0.8 Lag ACF

ACF of res model2

Normalized residuals

• Note that gam.check does

not correct for auto- correlation in residuals

• These plots are plotted

without gam.check

42

Model criticism: removing outliers

• Most outliers part of 1 trial (≈.2% of data)
• After removing that trial, fit again the same model:

# always order data first > dat <- dat[order(dat$Subject, dat$trial, dat$pTime, decreasing=F),] > m3 <- bam(Pupil ~ s(pTime) + s(pTime, by=Context) + s(pTime, by=Congruency) + s(pTime, by=Interaction) + s(pTime, by=ActorRight) + ti(pTime, trial, k=10) + s(pTime, Subject, bs='fs', m=1)+s(Item, bs='re'), data=dat, rho=rho)

43

Normalized residuals model 3

• Note: gam.check() was

not used here

44

Summary final model

Parametric coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.03357 0.01031 -3.257 0.00113 **

• Approximate significance of smooth terms:

edf Ref.df F p-value s(pTime) 7.305 7.782 15.341 < 2e-16 *** s(pTime):Context 7.439 8.674 20.854 < 2e-16 *** s(pTime):Congruency 4.445 5.507 37.507 < 2e-16 *** s(pTime):Interaction 7.145 8.410 3.771 0.000155 *** s(pTime):ActorRight 4.928 6.178 110.753 < 2e-16 *** ti(pTime,trial) 50.509 63.599 59.254 < 2e-16 *** s(pTime,Subject) 120.339 152.000 78.188 < 2e-16 *** s(Item) 30.769 31.000 171.292 < 2e-16 ***

• R-sq.(adj) = -0.126 Deviance explained = -12.3%

fREML score = -2.52e+05 Scale est. = 0.0048355 n = 78747

45

Effect of Context x Congruency

• The rho parameter results in more uncertainty in the effects of Context,

Congruency and Interaction

46

Context

500 1500 2500

• 0.04

0.00 0.04

Congruency

500 1500 2500

• 0.04

0.00 0.04

Interaction

500 1500 2500

• 0.04

0.00 0.04

Context

500 1500 2500

• 0.04

0.00 0.04

Congruency

500 1500 2500

• 0.04

0.00 0.04

Interaction

500 1500 2500

• 0.04

0.00 0.04

AP: Here you see a penguin and a sheep. PA: Here you see a sheep and a penguin. AP: Here you see a penguin and a sheep.

??

Summary

• A significant interaction between Context

and Congruency has been found:

• When the agent (the penguin) is

mentioned first, the congruency of the picture makes a difference

• When the patient (the sheep) is

mentioned first, the congruency of the picture does not make a difference

• Difference starts early: within 500 ms after

pronoun onset

• against initial filter accounts

47

1000 2000 3000

• 0.04

0.00 0.04

PA-con PA-inc AP-con AP-inc

Discussion

• Experimental data: variation between Subjects, between Items
• Use random intercepts, random slopes, or random wiggly curves
• Random effects may change your fixed effects
• Model criticism is very very important
• However, no model is perfect
• often with experimental data such as EEG / pupil dilation it is difficult

to get normal residuals

• GAMs can be used to model pupil dilation
• but some problems (bugs?) in mgcv need to be resolved

48

Lab session

1. Some basic examples:

• How to average your data in R
• Plotting averages in R

2. What is autocorrelation, and how to calculate normalized residuals (advanced) 3. Practice in setting up a model to test predictions

49