Vowel contrast: GAMM analysis
Holger Mitterer
2026-02-18
Note that this Markdown requires the file “mouse_data.txt”, which is generated by running the preprocessing markdown as well as the mousePreProcess.R script.
TARGET FIXATIONS
rm(list = ls())
library(lmerTest)## Loading required package: lme4## Loading required package: Matrix##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## steplibrary(stringr)
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.0 ✔ purrr 1.2.1
## ✔ forcats 1.0.1 ✔ readr 2.2.0
## ✔ ggplot2 4.0.2 ✔ tibble 3.3.1
## ✔ lubridate 1.9.5 ✔ tidyr 1.3.2## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ tidyr::pack() masks Matrix::pack()
## ✖ tidyr::unpack() masks Matrix::unpack()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(mgcv)## Loading required package: nlme
##
## Attaching package: 'nlme'
##
## The following object is masked from 'package:dplyr':
##
## collapse
##
## The following object is masked from 'package:lme4':
##
## lmList
##
## This is mgcv 1.9-3. For overview type 'help("mgcv-package")'.library(itsadug)## Loading required package: plotfunctions
##
## Attaching package: 'plotfunctions'
##
## The following object is masked from 'package:ggplot2':
##
## alpha
##
## Loaded package itsadug 2.5 (see 'help("itsadug")' ).options(warn = 0)
gazeData = read.table("dataCatalan/gazeTarget.txt", T, sep = "\t", quote = "")
trialData = read.table("dataCatalan/trialInfo.txt", T, sep = "\t", quote = "\"")
expFilter = grepl("_", trialData$condition)
gazeDataE = gazeData %>% filter(expFilter)
trialDataE = trialData %>% filter(expFilter)
timeline = -200 + (0:120) * 10
colnames(gazeDataE) = paste0("T_",timeline)
#contrast coding
temp = str_split_fixed(trialDataE$condition,"_",2)
trialDataE$accent = temp[, 1]
trialDataE$contrast = temp[, 2]
trialDataE$accentC = -0.5 + as.numeric(trialDataE$accent == "SEng")
table(trialDataE$accentC)##
## -0.5 0.5
## 3596 3598trialDataE$contrastC = -0.5 + as.numeric(trialDataE$contrast == "DRESS" |trialDataE$contrast == "TRAP" )
table(trialDataE$contrastC)##
## -0.5 0.5
## 3598 3596useData = cbind(trialDataE,gazeDataE)Analysis of the time course using GAMMs
We first need to reorganize the data.
gamData = useData[ ,c(1,2,4,11:14,15:135)] %>%
pivot_longer(., cols = starts_with("T_"),
names_to = "Time",
values_to = "position") %>% as.data.frame
gamData$position = as.numeric(gamData$position)
gamData$Time = gsub("T_","", gamData$Time) %>% as.numeric()
gamData = gamData[ order(gamData$pp,gamData$trial,gamData$Time), ]
gamData$pp = as.factor(gamData$pp)
gamData$item = as.factor(gamData$item)
gamData$isAustralian = ifelse(gamData$accent == "AusE",1,0)
gamData$isAustralianDress = gamData$isAustralian *
as.numeric(gamData$contrast == "DRESS" |gamData$contrast == "TRAP" )
gamData$vowelContrast = ifelse(gamData$contrast == "DRESS" |gamData$contrast == "TRAP",
"midFront","tensenessI") %>% as.factor()
#restricting the data to 0-900 after vowel onset
gamData = gamData %>% filter(Time > 0 & Time < 900)
gamData <- start_event(gamData, event=c("pp","trial"))
gamDataPlot = gamData %>% group_by(vowelContrast, Time, accent) %>%
summarize(meanFix = mean(position))## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by vowelContrast, Time, and accent.
## ℹ Output is grouped by vowelContrast and Time.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(vowelContrast, Time, accent))` for per-operation
## grouping (`?dplyr::dplyr_by`) instead.ggplot(gamDataPlot, aes(x = Time, y = meanFix, col = vowelContrast)) +
geom_line() + facet_wrap(~accent)
This is the descriptive data
Building an initial GAMM
fullModelFile = "target_overall_fullModels.Rdata"
if (file.exists(fullModelFile))
{
load(fullModelFile)
}else{
interaction.gamm0 <- bam( position ~ vowelContrast+
s(Time, by = vowelContrast) +
s(Time, by = isAustralian) +
s(Time, by = isAustralianDress) +
s(Time, pp, by = vowelContrast, bs="fs") +
s(Time, item, by = vowelContrast, bs="fs"),
family = "binomial",
data=gamData, discrete = T, nThreads = 3)
#dealing with autocorrelation
myAutoCorr =acf_resid(interaction.gamm0)
myRhoVal = myAutoCorr[2]
interaction.gamm1 <- bam( position ~ vowelContrast+
s(Time, by = vowelContrast, k=12) +
s(Time, by = isAustralian, k= 12) +
s(Time, by = isAustralianDress) +
s(Time, pp, by = vowelContrast, bs="fs") +
s(Time, item, by = vowelContrast, bs="fs"),
family = "binomial",
data=gamData, rho = myRhoVal, method = "ML",
AR.start = gamData$start.event, nThreads = 3, discrete = T)
save(list= c("interaction.gamm0","interaction.gamm1"), file = fullModelFile)
}Model criticism
How are we doing regarding auto-correlation?
myAutoCorr =acf_resid(interaction.gamm1)
Is the default for the smooth (k = 10) sufficient?
gam.check(interaction.gamm1)
##
## Method: ML Optimizer: outer newton
## full convergence after 7 iterations.
## Gradient range [-0.01483295,0.2122266]
## (score 907788.6 & scale 1).
## Hessian positive definite, eigenvalue range [0.01440749,64.73374].
## Model rank = 1746 / 1746
##
## 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(Time):vowelContrastmidFront 11.00 8.08 1 0.36
## s(Time):vowelContrasttensenessI 11.00 7.25 1 0.41
## s(Time):isAustralian 12.00 2.00 1 0.37
## s(Time):isAustralianDress 10.00 6.60 1 0.39
## s(Time,pp):vowelContrastmidFront 450.00 268.00 1 0.36
## s(Time,pp):vowelContrasttensenessI 450.00 279.05 1 0.38
## s(Time,item):vowelContrastmidFront 400.00 117.29 1 0.39
## s(Time,item):vowelContrasttensenessI 400.00 122.48 1 0.39Let’s have a look at this model no whether there is an overall effect of contrast
summary(interaction.gamm1)##
## Family: binomial
## Link function: logit
##
## Formula:
## position ~ vowelContrast + s(Time, by = vowelContrast, k = 12) +
## s(Time, by = isAustralian, k = 12) + s(Time, by = isAustralianDress) +
## s(Time, pp, by = vowelContrast, bs = "fs") + s(Time, item,
## by = vowelContrast, bs = "fs")
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.27085 0.07610 3.559 0.000372 ***
## vowelContrasttensenessI 0.01484 0.10870 0.137 0.891380
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(Time):vowelContrastmidFront 8.081 9.177 398.41 < 2e-16 ***
## s(Time):vowelContrasttensenessI 7.251 8.346 308.41 < 2e-16 ***
## s(Time):isAustralian 2.003 2.005 15.34 0.000464 ***
## s(Time):isAustralianDress 6.602 7.697 80.92 < 2e-16 ***
## s(Time,pp):vowelContrastmidFront 268.004 448.000 14510.19 < 2e-16 ***
## s(Time,pp):vowelContrasttensenessI 279.047 448.000 14097.18 < 2e-16 ***
## s(Time,item):vowelContrastmidFront 117.288 198.000 4047.50 < 2e-16 ***
## s(Time,item):vowelContrasttensenessI 122.484 198.000 5722.33 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.258 Deviance explained = 20.8%
## -ML = 9.0779e+05 Scale est. = 1 n = 640266There is a main effect of competitor, meaning that, overall, the mouse position is closer to target when the competitor is dissimilar to the target. Note that the significant effects of the tensors only indicate changes over time, they do not constitute evidence for a significant difference between conditions.
smooth <- plot_smooth(interaction.gamm1, view = "Time",
plot_all = "vowelContrast",
rug = FALSE, rm.ranef = TRUE, lwd = 4)## Summary:
## * vowelContrast : factor; set to the value(s): midFront, tensenessI.
## * Time : numeric predictor; with 30 values ranging from 10.000000 to 890.000000.
## * isAustralian : numeric predictor; set to the value(s): 0.
## * isAustralianDress : numeric predictor; set to the value(s): 0.
## * pp : factor; set to the value(s): CAT01. (Might be canceled as random effect, check below.)
## * item : factor; set to the value(s): baggage. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(Time,pp):vowelContrastmidFront,s(Time,pp):vowelContrasttensenessI,s(Time,item):vowelContrastmidFront,s(Time,item):vowelContrasttensenessI
## 
Now model comparison
We succesively take effects out to prepare model comparisons.
reducedModelFile = "target_overall_reducedModels.Rdata"
if (file.exists(reducedModelFile))
{
load(reducedModelFile)
}else{
minusInteraction.gamm <- bam( position ~ vowelContrast+
s(Time, by = vowelContrast) +
s(Time, by = isAustralian) +
#s(Time, by = isAustralianDress) +
s(Time, pp, by = vowelContrast, bs="fs") +
s(Time, item, by = vowelContrast, bs="fs"), family = "binomial",
data=gamData, rho = myRhoVal, method = "ML",
AR.start = gamData$start.event, nThreads = 3, discrete = T)
minusAccent.gamm <- bam( position ~ vowelContrast+
s(Time, by = vowelContrast) +
#s(Time, by = isAustralian) +
#s(Time, by = isAustralianDress) +
s(Time, pp, by = vowelContrast, bs="fs") +
s(Time, item, by = vowelContrast, bs="fs"), family = "binomial",
data=gamData, rho = myRhoVal, method = "ML",
AR.start = gamData$start.event, nThreads = 3, discrete = T)
timeOnly.gamm <- bam( position ~ 1 +
s(Time) +
#s(Time, by = isAustralian) +
#s(Time, by = isAustralianDress) +
s(Time, pp, by = vowelContrast, bs="fs") +
s(Time, item, by = vowelContrast, bs="fs"), family = "binomial",
data=gamData, rho = myRhoVal, method = "ML",
AR.start = gamData$start.event, nThreads = 3, discrete = T)
save(list= c("minusInteraction.gamm","minusAccent.gamm","timeOnly.gamm"), file = reducedModelFile)
}Model comparisons
Full model versus no interaction
compareML(interaction.gamm1, minusInteraction.gamm)## interaction.gamm1: position ~ vowelContrast + s(Time, by = vowelContrast, k = 12) +
## s(Time, by = isAustralian, k = 12) + s(Time, by = isAustralianDress) +
## s(Time, pp, by = vowelContrast, bs = "fs") + s(Time, item,
## by = vowelContrast, bs = "fs")
##
## minusInteraction.gamm: position ~ vowelContrast + s(Time, by = vowelContrast) + s(Time,
## by = isAustralian) + s(Time, pp, by = vowelContrast, bs = "fs") +
## s(Time, item, by = vowelContrast, bs = "fs")
##
## Model minusInteraction.gamm preferred: lower ML score (18.230), and lower df (3.000).
## -----
## Model Score Edf Difference Df
## 1 interaction.gamm1 907788.6 24
## 2 minusInteraction.gamm 907770.4 21 -18.230 3.000
##
## AIC difference: -47.28, model interaction.gamm1 has lower AIC.Both main effects versus minus accent
compareML(minusInteraction.gamm, minusAccent.gamm)## minusInteraction.gamm: position ~ vowelContrast + s(Time, by = vowelContrast) + s(Time,
## by = isAustralian) + s(Time, pp, by = vowelContrast, bs = "fs") +
## s(Time, item, by = vowelContrast, bs = "fs")
##
## minusAccent.gamm: position ~ vowelContrast + s(Time, by = vowelContrast) + s(Time,
## pp, by = vowelContrast, bs = "fs") + s(Time, item, by = vowelContrast,
## bs = "fs")
##
## Chi-square test of ML scores
## -----
## Model Score Edf Difference Df p.value Sig.
## 1 minusAccent.gamm 907772.4 18
## 2 minusInteraction.gamm 907770.4 21 1.977 3.000 0.267## Warning in compareML(minusInteraction.gamm, minusAccent.gamm): Only small difference in ML...## AIC difference: -114.88, model minusInteraction.gamm has lower AIC.Main effect of Contrast versus time-only
compareML(minusAccent.gamm, timeOnly.gamm)## minusAccent.gamm: position ~ vowelContrast + s(Time, by = vowelContrast) + s(Time,
## pp, by = vowelContrast, bs = "fs") + s(Time, item, by = vowelContrast,
## bs = "fs")
##
## timeOnly.gamm: position ~ 1 + s(Time) + s(Time, pp, by = vowelContrast, bs = "fs") +
## s(Time, item, by = vowelContrast, bs = "fs")
##
## Model timeOnly.gamm preferred: lower ML score (34.865), and lower df (3.000).
## -----
## Model Score Edf Difference Df
## 1 minusAccent.gamm 907772.4 18
## 2 timeOnly.gamm 907737.5 15 -34.865 3.000
##
## AIC difference: 0.12, model timeOnly.gamm has lower AIC.COMPETITOR FIXATIONS
Setting up the environment and then Getting all the mouse positions as a matrix
rm(list = ls())
library(lmerTest)
library(stringr)
library(tidyverse)
library(mgcv)
library(itsadug)
options(warn = 0)
gazeData = read.table("dataCatalan/gazeCompet.txt", T, sep = "\t", quote = "")
trialData = read.table("dataCatalan/trialInfo.txt", T, sep = "\t", quote = "\"")
expFilter = grepl("_", trialData$condition)
gazeDataE = gazeData %>% filter(expFilter)
trialDataE = trialData %>% filter(expFilter)
timeline = -200 + (0:120) * 10
colnames(gazeDataE) = paste0("T_",timeline)
#contrast coding
temp = str_split_fixed(trialDataE$condition,"_",2)
trialDataE$accent = temp[, 1]
trialDataE$contrast = temp[, 2]
trialDataE$accentC = -0.5 + as.numeric(trialDataE$accent == "SEng")
table(trialDataE$accentC)##
## -0.5 0.5
## 3596 3598trialDataE$contrastC = -0.5 + as.numeric(trialDataE$contrast == "DRESS" |trialDataE$contrast == "TRAP" )
table(trialDataE$contrastC)##
## -0.5 0.5
## 3598 3596useData = cbind(trialDataE,gazeDataE)Analysis of the time course using GAMMs
We first need to reorganize the data.
gamData = useData[ ,c(1,2,4,11:14,15:135)] %>%
pivot_longer(., cols = starts_with("T_"),
names_to = "Time",
values_to = "position") %>% as.data.frame
gamData$position = as.numeric(gamData$position)
gamData$Time = gsub("T_","", gamData$Time) %>% as.numeric()
gamData = gamData[ order(gamData$pp,gamData$trial,gamData$Time), ]
gamData$pp = as.factor(gamData$pp)
gamData$item = as.factor(gamData$item)
gamData$isAustralian = ifelse(gamData$accent == "AusE",1,0)
gamData$isAustralianDress = gamData$isAustralian *
as.numeric(gamData$contrast == "DRESS" |gamData$contrast == "TRAP" )
gamData$vowelContrast = ifelse(gamData$contrast == "DRESS" |gamData$contrast == "TRAP",
"midFront","tensenessI") %>% as.factor()
#restricting the data to 0-900 after vowel onset
gamData = gamData %>% filter(Time > 0 & Time < 900)
gamData <- start_event(gamData, event=c("pp","trial"))
gamDataPlot = gamData %>% group_by(vowelContrast, Time, accent) %>%
summarize(meanFix = mean(position))## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by vowelContrast, Time, and accent.
## ℹ Output is grouped by vowelContrast and Time.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(vowelContrast, Time, accent))` for per-operation
## grouping (`?dplyr::dplyr_by`) instead.ggplot(gamDataPlot, aes(x = Time, y = meanFix, col = vowelContrast)) +
geom_line() + facet_wrap(~accent)
This is the descriptive data
Building an initial GAMM
fullModelFile = "competitor_overall_fullModels.Rdata"
if (file.exists(fullModelFile))
{
load(fullModelFile)
}else{
interaction.Comp.gamm0 <- bam( position ~ vowelContrast+
s(Time, by = vowelContrast) +
s(Time, by = isAustralian) +
s(Time, by = isAustralianDress) +
s(Time, pp, by = vowelContrast, bs="fs") +
s(Time, item, by = vowelContrast, bs="fs"),
family = "binomial",
data=gamData, discrete = T, nThreads = 3)
#dealing with autocorrelation
myAutoCorr =acf_resid(interaction.Comp.gamm0)
myRhoVal = myAutoCorr[2]
interaction.Comp.gamm1 <- bam( position ~ vowelContrast+
s(Time, by = vowelContrast, k=12) +
s(Time, by = isAustralian, k= 12) +
s(Time, by = isAustralianDress) +
s(Time, pp, by = vowelContrast, bs="fs") +
s(Time, item, by = vowelContrast, bs="fs"),
family = "binomial",
data=gamData, rho = myRhoVal, method = "ML",
AR.start = gamData$start.event, nThreads = 3, discrete = T)
save(list= c("interaction.Comp.gamm0","interaction.Comp.gamm1"), file = fullModelFile)
}Model criticism
How are we doing regarding auto-correlation?
myAutoCorr =acf_resid(interaction.Comp.gamm1)
Is the default for the smooth (k = 10) sufficient?
gam.check(interaction.Comp.gamm1)
##
## Method: ML Optimizer: outer newton
## full convergence after 5 iterations.
## Gradient range [-0.01541229,0.01586758]
## (score 904971.4 & scale 1).
## Hessian positive definite, eigenvalue range [0.01509869,67.37419].
## Model rank = 1746 / 1746
##
## 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(Time):vowelContrastmidFront 11.00 6.80 0.98 0.80
## s(Time):vowelContrasttensenessI 11.00 8.29 0.98 0.76
## s(Time):isAustralian 12.00 2.00 0.98 0.80
## s(Time):isAustralianDress 10.00 4.03 0.98 0.71
## s(Time,pp):vowelContrastmidFront 450.00 283.62 0.98 0.75
## s(Time,pp):vowelContrasttensenessI 450.00 285.41 0.98 0.78
## s(Time,item):vowelContrastmidFront 400.00 124.26 0.98 0.75
## s(Time,item):vowelContrasttensenessI 400.00 121.27 0.98 0.78Let’s have a look at this model no whether there is an overall effect of contrast
summary(interaction.Comp.gamm1)##
## Family: binomial
## Link function: logit
##
## Formula:
## position ~ vowelContrast + s(Time, by = vowelContrast, k = 12) +
## s(Time, by = isAustralian, k = 12) + s(Time, by = isAustralianDress) +
## s(Time, pp, by = vowelContrast, bs = "fs") + s(Time, item,
## by = vowelContrast, bs = "fs")
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.73400 0.05714 -30.347 <2e-16 ***
## vowelContrasttensenessI 0.05727 0.08820 0.649 0.516
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(Time):vowelContrastmidFront 6.803 7.857 76.06 < 2e-16 ***
## s(Time):vowelContrasttensenessI 8.293 9.301 128.62 < 2e-16 ***
## s(Time):isAustralian 2.002 2.004 10.94 0.004244 **
## s(Time):isAustralianDress 4.027 4.762 23.02 0.000745 ***
## s(Time,pp):vowelContrastmidFront 283.620 448.000 7655.26 < 2e-16 ***
## s(Time,pp):vowelContrasttensenessI 285.413 448.000 7919.94 < 2e-16 ***
## s(Time,item):vowelContrastmidFront 124.255 198.000 3053.46 < 2e-16 ***
## s(Time,item):vowelContrasttensenessI 121.265 198.000 4759.45 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.081 Deviance explained = 9.3%
## -ML = 9.0497e+05 Scale est. = 1 n = 640266There is a main effect of competitor, meaning that, overall, the mouse position is closer to target when the competitor is dissimilar to the target. Note that the significant effects of the tensors only indicate changes over time, they do not constitute evidence for a significant difference between conditions.
smooth <- plot_smooth(interaction.Comp.gamm1, view = "Time",
plot_all = "vowelContrast",
rug = FALSE, rm.ranef = TRUE, lwd = 4)## Summary:
## * vowelContrast : factor; set to the value(s): midFront, tensenessI.
## * Time : numeric predictor; with 30 values ranging from 10.000000 to 890.000000.
## * isAustralian : numeric predictor; set to the value(s): 0.
## * isAustralianDress : numeric predictor; set to the value(s): 0.
## * pp : factor; set to the value(s): CAT01. (Might be canceled as random effect, check below.)
## * item : factor; set to the value(s): baggage. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(Time,pp):vowelContrastmidFront,s(Time,pp):vowelContrasttensenessI,s(Time,item):vowelContrastmidFront,s(Time,item):vowelContrasttensenessI
## 
Now model comparison
We succesively take effects out to prepare model comparisons.
reducedModelFile = "competitor_overall_reducedModels.Rdata"
if (file.exists(reducedModelFile))
{
load(reducedModelFile)
}else{
minusInteraction.Comp.gamm <- bam( position ~ vowelContrast+
s(Time, by = vowelContrast) +
s(Time, by = isAustralian) +
#s(Time, by = isAustralianDress) +
s(Time, pp, by = vowelContrast, bs="fs") +
s(Time, item, by = vowelContrast, bs="fs"),
family = "binomial",
data=gamData, rho = myRhoVal, method = "ML",
AR.start = gamData$start.event, nThreads = 3, discrete = T)
minusAccent.Comp.gamm <- bam( position ~ vowelContrast+
s(Time, by = vowelContrast) +
#s(Time, by = isAustralian) +
#s(Time, by = isAustralianDress) +
s(Time, pp, by = vowelContrast, bs="fs") +
s(Time, item, by = vowelContrast, bs="fs"),
family = "binomial",
data=gamData, rho = myRhoVal, method = "ML",
AR.start = gamData$start.event, nThreads = 3, discrete = T)
timeOnly.Comp.gamm <- bam( position ~ 1 +
s(Time) +
#s(Time, by = isAustralian) +
#s(Time, by = isAustralianDress) +
s(Time, pp, by = vowelContrast, bs="fs") +
s(Time, item, by = vowelContrast, bs="fs"),
family = "binomial",
data=gamData, rho = myRhoVal, method = "ML",
AR.start = gamData$start.event, nThreads = 3, discrete = T)
save(list= c("minusInteraction.Comp.gamm","minusAccent.Comp.gamm","timeOnly.Comp.gamm"), file = reducedModelFile)
}Model comparisons
Full model versus no interaction
compareML(interaction.Comp.gamm1, minusInteraction.Comp.gamm)## interaction.Comp.gamm1: position ~ vowelContrast + s(Time, by = vowelContrast, k = 12) +
## s(Time, by = isAustralian, k = 12) + s(Time, by = isAustralianDress) +
## s(Time, pp, by = vowelContrast, bs = "fs") + s(Time, item,
## by = vowelContrast, bs = "fs")
##
## minusInteraction.Comp.gamm: position ~ vowelContrast + s(Time, by = vowelContrast) + s(Time,
## by = isAustralian) + s(Time, pp, by = vowelContrast, bs = "fs") +
## s(Time, item, by = vowelContrast, bs = "fs")
##
## Model minusInteraction.Comp.gamm preferred: lower ML score (33.743), and lower df (3.000).
## -----
## Model Score Edf Difference Df
## 1 interaction.Comp.gamm1 904971.4 24
## 2 minusInteraction.Comp.gamm 904937.7 21 -33.743 3.000
##
## AIC difference: -8.95, model interaction.Comp.gamm1 has lower AIC.Both main effects versus minus accent
compareML(minusInteraction.Comp.gamm, minusAccent.Comp.gamm)## minusInteraction.Comp.gamm: position ~ vowelContrast + s(Time, by = vowelContrast) + s(Time,
## by = isAustralian) + s(Time, pp, by = vowelContrast, bs = "fs") +
## s(Time, item, by = vowelContrast, bs = "fs")
##
## minusAccent.Comp.gamm: position ~ vowelContrast + s(Time, by = vowelContrast) + s(Time,
## pp, by = vowelContrast, bs = "fs") + s(Time, item, by = vowelContrast,
## bs = "fs")
##
## Chi-square test of ML scores
## -----
## Model Score Edf Difference Df p.value Sig.
## 1 minusAccent.Comp.gamm 905009.8 18
## 2 minusInteraction.Comp.gamm 904937.7 21 72.114 3.000 < 2e-16 ***
##
## AIC difference: -59.24, model minusInteraction.Comp.gamm has lower AIC.Main effect of Contrast versus time-only
compareML(minusAccent.Comp.gamm, timeOnly.Comp.gamm)## minusAccent.Comp.gamm: position ~ vowelContrast + s(Time, by = vowelContrast) + s(Time,
## pp, by = vowelContrast, bs = "fs") + s(Time, item, by = vowelContrast,
## bs = "fs")
##
## timeOnly.Comp.gamm: position ~ 1 + s(Time) + s(Time, pp, by = vowelContrast, bs = "fs") +
## s(Time, item, by = vowelContrast, bs = "fs")
##
## Chi-square test of ML scores
## -----
## Model Score Edf Difference Df p.value Sig.
## 1 timeOnly.Comp.gamm 905031.1 15
## 2 minusAccent.Comp.gamm 905009.8 18 21.300 3.000 2.991e-09 ***
##
## AIC difference: 0.91, model timeOnly.Comp.gamm has lower AIC.VISUALIZATION
#rm(list = ls())
library(tidyverse)
library(mgcv)
library(itsadug)
library(stringr)
options(warn = 0)
# Load data
gazeTarget <- read.table("dataCatalan/gazeTarget.txt",
header = TRUE, sep = "\t", quote = "")
gazeComp <- read.table("dataCatalan/gazeCompet.txt",
header = TRUE, sep = "\t", quote = "")
trialData <- read.table("dataCatalan/trialInfo.txt",
header = TRUE, sep = "\t", quote = "\"")
# Remove fillers
expFilter <- grepl("_", trialData$condition)
gazeTargetE <- gazeTarget %>% filter(expFilter)
gazeCompE <- gazeComp %>% filter(expFilter)
trialDataE <- trialData %>% filter(expFilter)
# Timeline
timeline <- -200 + (0:120) * 10
colnames(gazeTargetE) <- paste0("Tt_", timeline)
colnames(gazeCompE) <- paste0("Tc_", timeline)
# Accent - Contrast
temp <- str_split_fixed(trialDataE$condition, "_", 2)
trialDataE$accent <- temp[,1]
trialDataE$contrast <- temp[,2]
# Combine datasets (same style as original)
useData <- cbind(trialDataE, gazeTargetE, gazeCompE) %>% as.data.frame()gamData <- useData %>%
pivot_longer(
cols = matches("^T[tc]_-?\\d+$"),
names_to = c("which", "Time"),
names_pattern = "^T([tc])_(-?\\d+)$",
values_to = "fix"
) %>%
pivot_wider(
names_from = which,
values_from = fix
) %>%
rename(target = t, competitor = c) %>%
mutate(
target = as.numeric(target),
competitor = as.numeric(competitor),
diffFix = target - competitor,
Time = as.numeric(Time)
) %>%
as.data.frame()
# Sort (like original)
gamData <- gamData[order(gamData$pp, gamData$trial, gamData$Time), ]
# Factors (like original)
gamData$pp <- as.factor(gamData$pp)
gamData$item <- as.factor(gamData$item)
# Accent/contrast coding (same as your current approach)
gamData$isAustralian <- ifelse(gamData$accent == "AusE", 1, 0)
gamData$vowelContrast <- factor(
ifelse(gamData$contrast %in% c("DRESS","TRAP"), "midFront", "tensenessI")
)
gamData$isAustralianDress <- gamData$isAustralian *
as.numeric(gamData$contrast %in% c("DRESS","TRAP"))
# Time window
gamData <- gamData %>% filter(Time > 0 & Time < 900)
gamData <- start_event(gamData, event = c("pp","trial"))library(dplyr)
library(tidyr)
library(ggplot2)
gamData <- gamData %>%
mutate(
accent = factor(accent),
contrast = toupper(contrast)
)
overall_data <- gamData %>%
group_by(Time, accent, vowelContrast) %>%
summarise(
target = mean(target, na.rm = TRUE),
competitor = mean(competitor, na.rm = TRUE),
.groups = "drop"
) %>%
mutate(
diff = target - competitor,
row = "Overall",
condition = as.character(vowelContrast)
)
dress_trap_data <- gamData %>%
filter(contrast %in% c("DRESS", "TRAP")) %>%
group_by(Time, accent, contrast) %>%
summarise(
target = mean(target, na.rm = TRUE),
competitor = mean(competitor, na.rm = TRUE),
.groups = "drop"
) %>%
mutate(
diff = target - competitor,
row = "Dress vs Trap",
condition = contrast
)
kit_fleece_data <- gamData %>%
filter(contrast %in% c("KIT", "FLEECE")) %>%
group_by(Time, accent, contrast) %>%
summarise(
target = mean(target, na.rm = TRUE),
competitor = mean(competitor, na.rm = TRUE),
.groups = "drop"
) %>%
mutate(
diff = target - competitor,
row = "KIT vs Fleece",
condition = contrast
)
plot_data <- bind_rows(overall_data, dress_trap_data, kit_fleece_data)
plot_data_long <- plot_data %>%
pivot_longer(
cols = c(target, competitor, diff),
names_to = "Measure",
values_to = "Value"
) %>%
mutate(
row = factor(row, levels = c("Overall", "Dress vs Trap", "KIT vs Fleece")),
Measure = factor(Measure, levels = c("target", "competitor", "diff")),
Measure = recode(Measure,
target = "Looks to Target",
competitor = "Looks to Competitor",
diff = "Target - Competitor")
)
Fig.Fixations <- ggplot(plot_data_long,
aes(x = Time, y = Value, colour = condition, linetype = accent)) +
geom_line(linewidth = 0.5) +
facet_grid(row ~ Measure) +
labs(
x = "Time (ms)",
y = "Proportion (or difference)",
colour = "Condition",
linetype = "Accent"
) +
theme_grey() +
theme(
strip.background = element_blank(),
legend.position = "bottom"
)
print(Fig.Fixations)
Across all three conditions, looks to the target increase over time while looks to the competitor decrease in a broadly similar pattern. The two accents (solid vs. dotted lines) also show very similar time courses, with only minor divergence between them within each condition.