Shopping Game ET models

Median pupil change from pre-trial baseline, per trip

First, we just model the pupil change from the pre-trial baseline.

ETdata <- read.csv("em6_eyetracking.csv", header = TRUE)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.2.3
## Warning: package 'ggplot2' was built under R version 4.2.3
## Warning: package 'readr' was built under R version 4.2.3
## Warning: package 'stringr' was built under R version 4.2.3
## Warning: package 'forcats' was built under R version 4.2.3
## Warning: package 'lubridate' was built under R version 4.2.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ readr     2.1.4
## ✔ ggplot2   3.4.2     ✔ stringr   1.5.0
## ✔ lubridate 1.9.2     ✔ tibble    3.1.8
## ✔ purrr     0.3.4     ✔ tidyr     1.2.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
library(ggpubr)
library(ggforce)
## Warning: package 'ggforce' was built under R version 4.2.3
library(effectsize)
library(lme4)
## Warning: package 'lme4' was built under R version 4.2.3
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
library(lmerTest)
## Warning: package 'lmerTest' was built under R version 4.2.3
## 
## Attaching package: 'lmerTest'
## 
## The following object is masked from 'package:lme4':
## 
##     lmer
## 
## The following object is masked from 'package:stats':
## 
##     step
library(emmeans)
## Warning: package 'emmeans' was built under R version 4.2.3
library(MuMIn)
## Warning: package 'MuMIn' was built under R version 4.2.3
library(MCMCglmm)
## Warning: package 'MCMCglmm' was built under R version 4.2.3
## Loading required package: coda
## Warning: package 'coda' was built under R version 4.2.3
## Loading required package: ape
## Warning: package 'ape' was built under R version 4.2.3
## 
## Attaching package: 'ape'
## 
## The following object is masked from 'package:ggpubr':
## 
##     rotate
library(mixedpower)
library(DHARMa)
## Warning: package 'DHARMa' was built under R version 4.2.3
## This is DHARMa 0.4.6. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
# Disable scientific notation
options(scipen = 999)

Summarized_ETdata1 <- ETdata %>% 
  group_by(participant,condsFile,trial, adjustTrip) %>% 
  summarise(medianPupilChange=median(change_from_baseline_classic,na.rm = T),
            baseline=first(baseline)
  )
## `summarise()` has grouped output by 'participant', 'condsFile', 'trial'. You
## can override using the `.groups` argument.
# Now we want to just look at regular trials
Summarized_ETdata1 <- Summarized_ETdata1 %>%
  filter(!condsFile %in% c("delayAfterListClosing", "delayBeforeListClosing", "no_costClosing"))

Now we run the model.

Summarized_ETdata1$condsFile <- as.factor(Summarized_ETdata1$condsFile)

# pretty normal
hist(Summarized_ETdata1$medianPupilChange)

pupil_model1 <- glmer(medianPupilChange ~ condsFile + (1|participant) + (1|adjustTrip) + (1|trial), 
                     data = Summarized_ETdata1, 
                     family = gaussian(link = "identity"))
## Warning in glmer(medianPupilChange ~ condsFile + (1 | participant) + (1 | :
## calling glmer() with family=gaussian (identity link) as a shortcut to lmer() is
## deprecated; please call lmer() directly
plot(density(Summarized_ETdata1$medianPupilChange, na.rm = T), xlim=c(-50, 50), ylim=c(0, 0.6), main='M_1 y_hat')
lines(density(predict(pupil_model1, type='response') ), col='red')

# looks great
plot(pupil_model1)

# OK
residuals <- residuals(pupil_model1)
qqnorm(residuals)
qqline(residuals, col = "red")

summary(pupil_model1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## medianPupilChange ~ condsFile + (1 | participant) + (1 | adjustTrip) +  
##     (1 | trial)
##    Data: Summarized_ETdata1
## 
## REML criterion at convergence: -110.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8727 -0.5396 -0.0056  0.5987  3.2238 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  participant (Intercept) 0.017823 0.13350 
##  adjustTrip  (Intercept) 0.002077 0.04558 
##  trial       (Intercept) 0.006673 0.08169 
##  Residual                0.041579 0.20391 
## Number of obs: 564, groups:  participant, 24; adjustTrip, 10; trial, 4
## 
## Fixed effects:
##                          Estimate Std. Error t value
## (Intercept)              -0.07024    0.05671  -1.239
## condsFiledelayBeforeList -0.10147    0.02287  -4.436
## condsFileno_cost          0.06385    0.02144   2.978
## 
## Correlation of Fixed Effects:
##             (Intr) cndFBL
## cndsFldlyBL -0.197       
## cndsFln_cst -0.259  0.520
# Get the p value, and R squared
pupil_null <- glmer(medianPupilChange ~ 1 + (1|participant) + (1|adjustTrip) + (1|trial), 
                     data = Summarized_ETdata1, family = gaussian(link = "identity"))
## Warning in glmer(medianPupilChange ~ 1 + (1 | participant) + (1 | adjustTrip) +
## : calling glmer() with family=gaussian (identity link) as a shortcut to lmer()
## is deprecated; please call lmer() directly
pupilAov <- cbind(car::Anova(pupil_model1),
                   r.squaredGLMM(pupil_model1,pupil_null))
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
# full model fits better
anova(pupil_null,pupil_model1)
## refitting model(s) with ML (instead of REML)
## Data: Summarized_ETdata1
## Models:
## pupil_null: medianPupilChange ~ 1 + (1 | participant) + (1 | adjustTrip) + (1 | trial)
## pupil_model1: medianPupilChange ~ condsFile + (1 | participant) + (1 | adjustTrip) + (1 | trial)
##              npar      AIC     BIC logLik deviance Chisq Df        Pr(>Chisq)
## pupil_null      5  -61.982 -40.307 35.991  -71.982                           
## pupil_model1    7 -112.492 -82.147 63.246 -126.492 54.51  2 0.000000000001456
##                 
## pupil_null      
## pupil_model1 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pupilAov
##              Chisq Df            Pr(>Chisq)        R2m      R2c
## condsFile 57.90899  2 0.0000000000002662084 0.06332882 0.428544
# Pairwise comparisons
library(emmeans)
emmeans(pupil_model1, list(pairwise ~ condsFile), adjust = "tukey")
## $`emmeans of condsFile`
##  condsFile         emmean     SE   df lower.CL upper.CL
##  delayAfterList  -0.07024 0.0576 9.34   -0.200   0.0592
##  delayBeforeList -0.17171 0.0577 9.43   -0.301  -0.0422
##  no_cost         -0.00639 0.0557 8.50   -0.134   0.1208
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $`pairwise differences of condsFile`
##  1                                estimate     SE  df t.ratio p.value
##  delayAfterList - delayBeforeList   0.1015 0.0229 534   4.434  <.0001
##  delayAfterList - no_cost          -0.0638 0.0216 534  -2.957  0.0091
##  delayBeforeList - no_cost         -0.1653 0.0219 539  -7.555  <.0001
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates

The above model shows that if we get a median of the pupil for the entire list phase and use the pre-trial baseline, then the median pupil change from baseline in the delay-before-list condition is significantly less than the no-delay and delay-after-list, and that the no-delay condition has the largest pupil (but that it’s still overall a constriction compared to the pre-trial baseline).

Median pupil change from first 500ms of each list trip, per trip

Now we do the same, but instead look at the median pupil change from baseline over the whole list, but where the baseline is the first 500ms of looking at the list (so there is a unique baseline value per trip, but its still not ideal as participants are orienting themselves to pay attention to the list)

Summarized_ETdata2 <- ETdata %>% 
  group_by(participant,condsFile,trial, adjustTrip) %>% 
  summarise(medianPupilChange=median(change_from_baseline2,na.rm = T),
            baseline=first(baseline2)
  )
## `summarise()` has grouped output by 'participant', 'condsFile', 'trial'. You
## can override using the `.groups` argument.
# Now we want to just look at regular trials
Summarized_ETdata2 <- Summarized_ETdata2 %>%
  filter(!condsFile %in% c("delayAfterListClosing", "delayBeforeListClosing", "no_costClosing"))

Summarized_ETdata2$condsFile <- as.factor(Summarized_ETdata2$condsFile)

# pretty normal
hist(Summarized_ETdata2$medianPupilChange)

pupil_model2 <- glmer(medianPupilChange ~ condsFile + (1|participant) + (1|adjustTrip) + (1|trial), 
                     data = Summarized_ETdata2, 
                     family = gaussian(link = "identity"))
## Warning in glmer(medianPupilChange ~ condsFile + (1 | participant) + (1 | :
## calling glmer() with family=gaussian (identity link) as a shortcut to lmer() is
## deprecated; please call lmer() directly
## boundary (singular) fit: see help('isSingular')
plot(density(Summarized_ETdata2$medianPupilChange, na.rm = T), xlim=c(-50, 50), ylim=c(0, 0.6), main='M_1 y_hat')
lines(density(predict(pupil_model2, type='response') ), col='red')

# looks great
plot(pupil_model2)

# OK but less nice
residuals <- residuals(pupil_model2)
qqnorm(residuals)
qqline(residuals, col = "red")

summary(pupil_model2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## medianPupilChange ~ condsFile + (1 | participant) + (1 | adjustTrip) +  
##     (1 | trial)
##    Data: Summarized_ETdata2
## 
## REML criterion at convergence: -512.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1943 -0.4848  0.0194  0.5904  2.7870 
## 
## Random effects:
##  Groups      Name        Variance  Std.Dev.
##  participant (Intercept) 0.0048967 0.06998 
##  adjustTrip  (Intercept) 0.0001084 0.01041 
##  trial       (Intercept) 0.0000000 0.00000 
##  Residual                0.0213691 0.14618 
## Number of obs: 566, groups:  participant, 24; adjustTrip, 10; trial, 4
## 
## Fixed effects:
##                           Estimate Std. Error t value
## (Intercept)              -0.090487   0.019896  -4.548
## condsFiledelayBeforeList  0.070042   0.016275   4.304
## condsFileno_cost         -0.005953   0.015122  -0.394
## 
## Correlation of Fixed Effects:
##             (Intr) cndFBL
## cndsFldlyBL -0.397       
## cndsFln_cst -0.454  0.528
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# Get the p value, and R squared
pupil_null <- glmer(medianPupilChange ~ 1 + (1|participant) + (1|adjustTrip) + (1|trial), 
                     data = Summarized_ETdata2, family = gaussian(link = "identity"))
## Warning in glmer(medianPupilChange ~ 1 + (1 | participant) + (1 | adjustTrip) +
## : calling glmer() with family=gaussian (identity link) as a shortcut to lmer()
## is deprecated; please call lmer() directly
## boundary (singular) fit: see help('isSingular')
pupilAov <- cbind(car::Anova(pupil_model2),
                   r.squaredGLMM(pupil_model2,pupil_null))

# full model fits better
anova(pupil_null,pupil_model2)
## refitting model(s) with ML (instead of REML)
## Data: Summarized_ETdata2
## Models:
## pupil_null: medianPupilChange ~ 1 + (1 | participant) + (1 | adjustTrip) + (1 | trial)
## pupil_model2: medianPupilChange ~ condsFile + (1 | participant) + (1 | adjustTrip) + (1 | trial)
##              npar     AIC     BIC logLik deviance  Chisq Df   Pr(>Chisq)    
## pupil_null      5 -494.31 -472.62 252.16  -504.31                           
## pupil_model2    7 -518.05 -487.68 266.02  -532.05 27.741  2 0.0000009466 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pupilAov
##             Chisq Df      Pr(>Chisq)       R2m       R2c
## condsFile 28.3841  2 0.0000006862308 0.0402094 0.2223507
# Pairwise comparisons
library(emmeans)
emmeans(pupil_model2, list(pairwise ~ condsFile), adjust = "tukey")
## $`emmeans of condsFile`
##  condsFile        emmean     SE   df lower.CL upper.CL
##  delayAfterList  -0.0905 0.0205 24.2  -0.1328  -0.0481
##  delayBeforeList -0.0204 0.0207 25.4  -0.0631   0.0222
##  no_cost         -0.0964 0.0191 23.4  -0.1359  -0.0570
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $`pairwise differences of condsFile`
##  1                                estimate     SE  df t.ratio p.value
##  delayAfterList - delayBeforeList -0.07004 0.0163 541  -4.296  0.0001
##  delayAfterList - no_cost          0.00595 0.0152 531   0.390  0.9194
##  delayBeforeList - no_cost         0.07599 0.0154 537   4.936  <.0001
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates

The above model is saying that if we use these values as the pupil and baseline, then the delay-before-list condition has the largest pupil compared to the no-delay and the delay-after-list. Still, its all a general constriction compared to baseline.

Median pupil change during the LAST 500ms, using pre-trial baseline

The next model will use the LAST 500ms of the list phase as the pupil area of interest (theoretically this is when participants have “loaded” up their memory) and the baseline will be the pre-trial baseline.

Summarized_ETdata3 <- ETdata %>% 
  group_by(participant,condsFile,trial, adjustTrip) %>% 
  summarise(medianPupilChange=median(last500_minus_baseline,na.rm = T),
            baseline=first(baseline)
  )
## `summarise()` has grouped output by 'participant', 'condsFile', 'trial'. You
## can override using the `.groups` argument.
# Now we want to just look at regular trials
Summarized_ETdata3 <- Summarized_ETdata3 %>%
  filter(!condsFile %in% c("delayAfterListClosing", "delayBeforeListClosing", "no_costClosing"))

Summarized_ETdata3$condsFile <- as.factor(Summarized_ETdata3$condsFile)

# pretty normal
hist(Summarized_ETdata3$medianPupilChange)

pupil_model3 <- glmer(medianPupilChange ~ condsFile + (1|participant) + (1|adjustTrip) + (1|trial), 
                     data = Summarized_ETdata3, 
                     family = gaussian(link = "identity"))
## Warning in glmer(medianPupilChange ~ condsFile + (1 | participant) + (1 | :
## calling glmer() with family=gaussian (identity link) as a shortcut to lmer() is
## deprecated; please call lmer() directly
plot(density(Summarized_ETdata3$medianPupilChange, na.rm = T), xlim=c(-50, 50), ylim=c(0, 0.6), main='M_1 y_hat')
lines(density(predict(pupil_model3, type='response') ), col='red')

# looks great
plot(pupil_model3)

# OK but less nice
residuals <- residuals(pupil_model3)
qqnorm(residuals)
qqline(residuals, col = "red")

summary(pupil_model3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## medianPupilChange ~ condsFile + (1 | participant) + (1 | adjustTrip) +  
##     (1 | trial)
##    Data: Summarized_ETdata3
## 
## REML criterion at convergence: -45.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3179 -0.6151 -0.0191  0.5455  2.5696 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  participant (Intercept) 0.027063 0.16451 
##  adjustTrip  (Intercept) 0.000237 0.01540 
##  trial       (Intercept) 0.007274 0.08529 
##  Residual                0.046616 0.21591 
## Number of obs: 560, groups:  participant, 24; adjustTrip, 10; trial, 4
## 
## Fixed effects:
##                          Estimate Std. Error t value
## (Intercept)              -0.08273    0.05832  -1.418
## condsFiledelayBeforeList -0.03449    0.02439  -1.414
## condsFileno_cost          0.04153    0.02250   1.846
## 
## Correlation of Fixed Effects:
##             (Intr) cndFBL
## cndsFldlyBL -0.200       
## cndsFln_cst -0.231  0.525
# Get the p value, and R squared
pupil_null <- glmer(medianPupilChange ~ 1 + (1|participant) + (1|adjustTrip) + (1|trial), 
                     data = Summarized_ETdata3, family = gaussian(link = "identity"))
## Warning in glmer(medianPupilChange ~ 1 + (1 | participant) + (1 | adjustTrip) +
## : calling glmer() with family=gaussian (identity link) as a shortcut to lmer()
## is deprecated; please call lmer() directly
pupilAov <- cbind(car::Anova(pupil_model3),
                   r.squaredGLMM(pupil_model3,pupil_null))

# full model fits better
anova(pupil_null,pupil_model3)
## refitting model(s) with ML (instead of REML)
## Data: Summarized_ETdata3
## Models:
## pupil_null: medianPupilChange ~ 1 + (1 | participant) + (1 | adjustTrip) + (1 | trial)
## pupil_model3: medianPupilChange ~ condsFile + (1 | participant) + (1 | adjustTrip) + (1 | trial)
##              npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)   
## pupil_null      5 -40.740 -19.101 25.370  -50.740                        
## pupil_model3    7 -47.116 -16.820 30.558  -61.116 10.376  2   0.005585 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pupilAov
##              Chisq Df  Pr(>Chisq)        R2m      R2c
## condsFile 11.23944  2 0.003625653 0.01219289 0.432841
# Pairwise comparisons
library(emmeans)
emmeans(pupil_model3, list(pairwise ~ condsFile), adjust = "tukey")
## $`emmeans of condsFile`
##  condsFile        emmean     SE   df lower.CL upper.CL
##  delayAfterList  -0.0827 0.0588 8.92   -0.216   0.0505
##  delayBeforeList -0.1172 0.0590 9.05   -0.251   0.0162
##  no_cost         -0.0412 0.0577 8.44   -0.173   0.0907
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $`pairwise differences of condsFile`
##  1                                estimate     SE  df t.ratio p.value
##  delayAfterList - delayBeforeList   0.0345 0.0244 533   1.413  0.3348
##  delayAfterList - no_cost          -0.0415 0.0226 522  -1.834  0.1596
##  delayBeforeList - no_cost         -0.0760 0.0231 531  -3.297  0.0030
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates

If we use the above measures for the pupil, then this model is saying that the no-delay condition again has the largest pupil change from baseline, and next comes delay-after-list and last is delay-before-list. All still constriction compared to baseline.

Median pupil change during the LAST 500ms, using first 500ms of list time as baseline

Finally, this model will use the last 500ms of the list phase as the pupil time of interest again, but the baseline will be the first 500 ms of list time, so there is a unique baseline per trip.

Summarized_ETdata4 <- ETdata %>% 
  group_by(participant,condsFile,trial, adjustTrip) %>% 
  summarise(medianPupilChange=median(last500_minus_first500,na.rm = T),
            baseline=first(baseline)
  )
## `summarise()` has grouped output by 'participant', 'condsFile', 'trial'. You
## can override using the `.groups` argument.
# Now we want to just look at regular trials
Summarized_ETdata4 <- Summarized_ETdata4 %>%
  filter(!condsFile %in% c("delayAfterListClosing", "delayBeforeListClosing", "no_costClosing"))

Summarized_ETdata4$condsFile <- as.factor(Summarized_ETdata4$condsFile)

# pretty normal
hist(Summarized_ETdata4$medianPupilChange)

pupil_model4 <- glmer(medianPupilChange ~ condsFile + (1|participant) + (1|adjustTrip) + (1|trial), 
                     data = Summarized_ETdata4, 
                     family = gaussian(link = "identity"))
## Warning in glmer(medianPupilChange ~ condsFile + (1 | participant) + (1 | :
## calling glmer() with family=gaussian (identity link) as a shortcut to lmer() is
## deprecated; please call lmer() directly
## boundary (singular) fit: see help('isSingular')
plot(density(Summarized_ETdata4$medianPupilChange, na.rm = T), xlim=c(-50, 50), ylim=c(0, 0.6), main='M_1 y_hat')
lines(density(predict(pupil_model4, type='response') ), col='red')

# looks great
plot(pupil_model4)

# OK but less nice
residuals <- residuals(pupil_model4)
qqnorm(residuals)
qqline(residuals, col = "red")

summary(pupil_model4)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## medianPupilChange ~ condsFile + (1 | participant) + (1 | adjustTrip) +  
##     (1 | trial)
##    Data: Summarized_ETdata4
## 
## REML criterion at convergence: -179.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2183 -0.5181  0.0420  0.5988  2.7873 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  participant (Intercept) 0.007934 0.08907 
##  adjustTrip  (Intercept) 0.001517 0.03895 
##  trial       (Intercept) 0.000000 0.00000 
##  Residual                0.038302 0.19571 
## Number of obs: 562, groups:  participant, 24; adjustTrip, 10; trial, 4
## 
## Fixed effects:
##                          Estimate Std. Error t value
## (Intercept)              -0.07601    0.03127  -2.431
## condsFiledelayBeforeList  0.13872    0.02194   6.323
## condsFileno_cost         -0.03493    0.02058  -1.697
## 
## Correlation of Fixed Effects:
##             (Intr) cndFBL
## cndsFldlyBL -0.343       
## cndsFln_cst -0.444  0.524
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# Get the p value, and R squared
pupil_null <- glmer(medianPupilChange ~ 1 + (1|participant) + (1|adjustTrip) + (1|trial), 
                     data = Summarized_ETdata4, family = gaussian(link = "identity"))
## Warning in glmer(medianPupilChange ~ 1 + (1 | participant) + (1 | adjustTrip) +
## : calling glmer() with family=gaussian (identity link) as a shortcut to lmer()
## is deprecated; please call lmer() directly
## boundary (singular) fit: see help('isSingular')
pupilAov <- cbind(car::Anova(pupil_model4),
                   r.squaredGLMM(pupil_model4,pupil_null))

# full model fits better
anova(pupil_null,pupil_model4)
## refitting model(s) with ML (instead of REML)
## Data: Summarized_ETdata4
## Models:
## pupil_null: medianPupilChange ~ 1 + (1 | participant) + (1 | adjustTrip) + (1 | trial)
## pupil_model4: medianPupilChange ~ condsFile + (1 | participant) + (1 | adjustTrip) + (1 | trial)
##              npar     AIC      BIC logLik deviance  Chisq Df
## pupil_null      5 -116.77  -95.107 63.382  -126.77          
## pupil_model4    7 -182.84 -152.525 98.423  -196.84 70.081  2
##                         Pr(>Chisq)    
## pupil_null                            
## pupil_model4 0.0000000000000006057 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pupilAov
##             Chisq Df                Pr(>Chisq)       R2m       R2c
## condsFile 74.5987  2 0.00000000000000006325548 0.1008339 0.2787909
# Pairwise comparisons
library(emmeans)
emmeans(pupil_model4, list(pairwise ~ condsFile), adjust = "tukey")
## $`emmeans of condsFile`
##  condsFile        emmean     SE   df lower.CL upper.CL
##  delayAfterList  -0.0760 0.0326 19.1 -0.14417 -0.00785
##  delayBeforeList  0.0627 0.0327 19.9 -0.00553  0.13093
##  no_cost         -0.1109 0.0297 18.0 -0.17329 -0.04859
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $`pairwise differences of condsFile`
##  1                                estimate     SE  df t.ratio p.value
##  delayAfterList - delayBeforeList  -0.1387 0.0220 535  -6.315  <.0001
##  delayAfterList - no_cost           0.0349 0.0208 534   1.683  0.2127
##  delayBeforeList - no_cost          0.1736 0.0209 540   8.300  <.0001
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates

The above model says that these pupil measures show that the delay-before-list condition has the largest pupil change and that its actually a dilation response, and then delay-after-list is behind that, and then the no-delay has the smallest pupil. But that the difference between the delay-after-list and the no-delay is not significant.

What each model says is the largest-to-smallest pupil size

Model 1: 1. No-delay (largest) 2. Delay-after-list 3. Delay-before-list (smallest)

Model 2: 1. Delay-before-list (largest) 2. Delay-after-list 3. No-delay (smallest)

Model 3: 1. No-delay (largest) 2. Delay-after-list 3. Delay-before-list (smallest)

Model 4: 1. Delay-before-list (largest) 2. Delay-after-list 3. No-delay (smallest)

So, when we use the pre-trial baseline, we see patterns like in Model 1 and 3. When we use the first 500ms of the list phase itself per trip, we see patterns like in Model 2 and 4. Clearly, choosing the baseline phase is super important! I am not sure what to do atm.