Assessing regional accent barriers for L2 spoken word recognition
Holger Mitterer & Gil Verbeke
2026-03-04
0. Overview
1. Participants
2. Stimulus Materials
2.1. (Lexical)
Properties
2.2. Acoustic Properties
3. Eye-Tracking Data
3.1. Visualization
3.2. Analysis of the Time Course
3.2.1.
Overall Effect of Contrast (FLEECE-KIT vs. DRESS-TRAP)
3.2.2.
FLEECE-KIT Contrast
3.2.3. DRESS-TRAP Contrast
1. Participants
Load datasets and remove data from three participants: (i)
CAT14: born in Dakar, Senegal; (ii) CAT17: .asc-file
with gaze patterns corrupted; (iii) CAT43: L1 Urdu speaker.
Questionnaire <- read.csv("WP3_Main_Questionnaire.csv", row.names = NULL, stringsAsFactors = TRUE)
Familiarity <- read.csv("WP3_Main_Familiarity.csv", row.names = NULL, stringsAsFactors = TRUE)
# Merge datasets by Participant
Questionnaire <- Questionnaire[!(Questionnaire$Part_ID %in% c("CAT14", "CAT17", "CAT43")),]
Questionnaire$Part_ID <- droplevels(Questionnaire$Part_ID)
Familiarity <- Familiarity[!(Familiarity$Part_ID %in% c("CAT14", "CAT17", "CAT43")),]
Familiarity$Part_ID <- droplevels(Familiarity$Part_ID)# Sample Size
str(Questionnaire$Part_ID)## Factor w/ 45 levels " CAT21"," CAT36",..: 3 4 5 6 7 8 9 10 11 12 ...# Age
sapply(list(mean = mean, sd = sd, min = min, max = max), function(f) f(Questionnaire$Age))## mean sd min max
## 19.11111 1.41778 18.00000 24.00000# Gender
table(Questionnaire$Gender)##
## Female Male Non-binary / third gender
## 35 9 1# Language Dominance
table(Questionnaire$LangDomParticipant)##
## Both equally Catalan Spanish
## 0 18 18 9# Programme
table(Questionnaire$FieldStudy)##
## Communication Sciences Engineering
## 1 5
## Linguistics Sea Sciences
## 21 1
## Traduction and interpretation Translation
## 1 1
## Translation and interpretation
## 15# Age of Onset of Acquisition
sapply(list(mean = mean, sd = sd, min = min, max = max), function(f) f(Questionnaire$AgeOnsetEnglish))## mean sd min max
## 4.788889 1.447917 2.000000 8.000000# Self-reported proficiency - speaking
sapply(list(mean = mean, sd = sd, min = min, max = max), function(f) f(Questionnaire$Prof_Speaking))## mean sd min max
## 5.1777778 0.9837056 3.0000000 7.0000000# Self-reported proficiency - listening
sapply(list(mean = mean, sd = sd, min = min, max = max), function(f) f(Questionnaire$Prof_Listening))## mean sd min max
## 6.0000000 0.7687061 4.0000000 7.0000000# Self-reported proficiency - reading
sapply(list(mean = mean, sd = sd, min = min, max = max), function(f) f(Questionnaire$Prof_Reading))## mean sd min max
## 5.977778 1.097288 3.000000 7.000000# Self-reported proficiency - writing
sapply(list(mean = mean, sd = sd, min = min, max = max), function(f) f(Questionnaire$Prof_Writing))## mean sd min max
## 5.3111111 0.9728641 3.0000000 7.0000000# Self-reported proficiency - overall
sapply(list(mean = mean, sd = sd, min = min, max = max), function(f) f(Questionnaire$Prof_Overall))## mean sd min max
## 5.4666667 0.8146388 4.0000000 7.0000000# LexTALE
sapply(list(mean = mean, sd = sd, min = min, max = max), function(f) f(Questionnaire$LexTALE))## mean sd min max
## 72.111111 9.511589 55.000000 92.500000# Familiarity
library(dplyr)
library(knitr)
Table1 <- Familiarity %>%
group_by(Accent) %>%
summarise(
Mean_Familiarity = mean(Familiarity, na.rm = TRUE),
SD_Familiarity = sd(Familiarity, na.rm = TRUE)
)
kable(Table1)| Accent | Mean_Familiarity | SD_Familiarity |
|---|---|---|
| A.Eng | 51.333333 | 19.522714 |
| AusE | 3.911111 | 5.455698 |
| CanE | 2.622222 | 4.608731 |
| N.Eng | 3.711111 | 5.092875 |
| Other | 4.955556 | 15.265760 |
| S.Eng | 29.577778 | 15.980797 |
| SA.Eng | 3.888889 | 5.296749 |
2. Stimulus Materials
2.1 (Lexical) Properties
library(dplyr)
library(ggplot2)
library(ggrepel)
library(knitr)
Stimuli <- read.csv("WP3_Stimuli.csv", row.names = NULL, stringsAsFactors = TRUE)
contrast_labels <- c(
"/i/-/I/" = "FLEECE-KIT",
"/e/-/a/" = "DRESS-TRAP",
"/u/-/o/" = "GOOSE-LOT"
)
fmt_m_sd <- function(x, digits = 2) {
sprintf(paste0("%.", digits, "f<br>(%.", digits, "f)"),
mean(x, na.rm = TRUE),
sd(x, na.rm = TRUE))
}Table2 <- Stimuli %>%
mutate(Contrast = recode(as.character(Contrast), !!!contrast_labels, .default = as.character(Contrast))) %>%
group_by(Contrast) %>%
summarise(
Freq1 = fmt_m_sd(LexFreq_1),
Freq2 = fmt_m_sd(LexFreq_2),
`ΔFreq` = fmt_m_sd(Delta_LexFreq),
PND1 = fmt_m_sd(PND_1),
PND2 = fmt_m_sd(PND_2),
`ΔPND` = fmt_m_sd(Delta_PND),
Fam1 = fmt_m_sd(Familiarity_1),
Fam2 = fmt_m_sd(Familiarity_2),
`ΔFam` = fmt_m_sd(Delta_Familiarity),
.groups = "drop"
) %>%
arrange(factor(Contrast, levels = c("FLEECE-KIT", "DRESS-TRAP", "GOOSE-LOT")))
kable(Table2)| Contrast | Freq1 | Freq2 | ΔFreq | PND1 | PND2 | ΔPND | Fam1 | Fam2 | ΔFam |
|---|---|---|---|---|---|---|---|---|---|
| FLEECE-KIT | 3.62 (0.70) |
3.39 (0.56) |
0.68 (0.44) |
2.50 (3.27) |
1.80 (3.12) |
1.70 (1.95) |
27.70 (3.13) |
26.50 (2.64) |
1.40 (1.51) |
| DRESS-TRAP | 3.85 (0.65) |
3.90 (0.50) |
0.56 (0.49) |
2.80 (3.91) |
4.20 (3.97) |
3.80 (3.49) |
28.30 (2.50) |
28.30 (2.11) |
2.40 (2.41) |
| GOOSE-LOT | 3.28 (0.44) |
3.93 (0.45) |
0.78 (0.46) |
2.50 (3.87) |
2.60 (3.06) |
2.50 (2.95) |
28.00 (2.00) |
29.00 (1.05) |
1.80 (1.69) |
Table 1 shows the mean values and SDS for lexical frequency (Freq), phonological neighborhood density (PND) and familiarity ratings (Fam) for the item pairs used across the three vowel contrasts. Subscripts 1-2 denote the first and second member of each pair, respectively. Difference scores (ΔFreq, ΔPND, ΔFam) represent the absolute difference between the two words within each pair, averaged across all pairs in the contrast.
2.2 Acoustic Properties
Create an F1-F2 vowel plot showing FLEECE, KIT, DRESS and TRAP for Southern British English (S.Eng) and Australian English (AusE).
library(dplyr)
library(ggplot2)
library(ggrepel)
Stimuli_Acoustics <- read.csv("WP3_Stimuli_Acoustics.csv", row.names = NULL, stringsAsFactors = TRUE)
stimuli_F1_F2 <- Stimuli_Acoustics %>%
mutate(
Accent = factor(Accent, levels = c("S.Eng", "AusE")),
Vowel = factor(toupper(Vowel))
)
stimuli_F1_F2 <- stimuli_F1_F2 %>%
filter(stimuli_F1_F2$Contrast != "GOOSE-LOT") %>%
droplevels()
means_va <- stimuli_F1_F2 %>%
group_by(Accent, Vowel) %>%
summarise(
F1_mean = mean(F1, na.rm = TRUE),
F2_mean = mean(F2, na.rm = TRUE),
.groups = "drop"
)Fig.Stimuli.F1_F2 <- ggplot(stimuli_F1_F2, aes(x = F2, y = F1, colour = Accent)) +
stat_ellipse(aes(group = interaction(Accent, Vowel)),
type = "norm",
level = 0.95,
linewidth = 0.6,
alpha = 0.5) +
geom_point(data = means_va,
aes(x = F2_mean, y = F1_mean),
size = 2) +
geom_text_repel(data = means_va,
aes(x = F2_mean, y = F1_mean, label = Vowel),
show.legend = FALSE,
size = 3,
point.padding = 0.3,
box.padding = 0.4,
min.segment.length = 0,
segment.alpha = 0.35) +
scale_colour_manual(values = c(
"S.Eng" = "#E41A1C",
"AusE" = "#377EB8"
)) +
scale_x_reverse() +
scale_y_reverse() +
labs(x = "F2 (Hz)", y = "F1 (Hz)", colour = "Accent") +
theme_grey() +
facet_wrap(.~ Accent) +
theme(
strip.background = element_blank(),
strip.text = element_blank(),
legend.position = "bottom",
legend.title = element_text(size = 9, face = "bold"))
plot(Fig.Stimuli.F1_F2)
Create a bar plot showing vowel duration by accent.
library(ggplot2)
library(plyr)
stimuli_duration <- Stimuli_Acoustics %>%
mutate(
Accent = factor(Accent, levels = c("S.Eng", "AusE")),
Vowel = factor(toupper(Vowel))
)
stimuli_duration <- stimuli_duration %>%
filter(stimuli_duration$Contrast != "GOOSE-LOT") %>%
mutate(Vowel_Duration = Vowel_Duration * 1000) %>%
droplevels()
stimuli_duration$Vowel <- factor(
stimuli_duration$Vowel,
levels = c("FLEECE", "KIT", "DRESS", "TRAP")
)
data_summary <- function(data, varname, groupnames){
require(plyr)
summary_func <- function(x, col){
c(mean = mean(x[[col]], na.rm=TRUE),
sd = sd(x[[col]], na.rm=TRUE))
}
data_sum<-ddply(data, groupnames, .fun=summary_func,
varname)
data_sum <- rename(data_sum, c("mean" = varname))
return(data_sum)
}
stimuli_duration <- data_summary(stimuli_duration, varname="Vowel_Duration",
groupnames=c("Vowel", "Accent"))Fig.Stimuli.Duration <- ggplot(stimuli_duration, aes(x = Vowel, y = Vowel_Duration, fill = Accent)) +
geom_bar(stat = "identity", color = "black",
position = position_dodge()) +
geom_errorbar(aes(ymin = Vowel_Duration - sd,
ymax = Vowel_Duration + sd),
width = .2,
position = position_dodge(.9)) +
scale_fill_manual(values = c(
"S.Eng" = "#E41A1C",
"AusE" = "#377EB8"
)) +
labs(y = "Duration (ms)") +
theme_grey() +
theme(
strip.background = element_blank(),
strip.text = element_blank(),
legend.position = "bottom",
legend.title = element_text(size = 9, face = "bold"))
plot(Fig.Stimuli.Duration)
detach("package:plyr", unload=TRUE)3. Eye-Tracking Data
3.1 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
useData <- cbind(trialDataE, gazeTargetE, gazeCompE) %>% as.data.frame()library(dplyr)
library(tidyr)
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
gamData <- gamData[order(gamData$pp, gamData$trial, gamData$Time), ]
# Factors
gamData$pp <- as.factor(gamData$pp)
gamData$item <- as.factor(gamData$item)
# Accent/contrast coding
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.
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 = "")
trialData <- read.table("dataCatalan/trialInfo.txt",
header = TRUE, sep = "\t", quote = "\"")
# Timeline
timeline <- -200 + (0:120) * 10
colnames(gazeTarget) <- paste0("T_", timeline)
# Combine datasets (same style as original)
useData <- cbind(trialData, gazeTarget) %>% as.data.frame()
# Remove practice trials
useData <- useData %>%
filter(useData$condition != "practice") %>%
droplevels()
# Rename levels of condition
accent <- ifelse(grepl("^SEng", useData$wavfile), "SEng",
ifelse(grepl("^AusE", useData$wavfile), "AusE", NA))
vowel <- ifelse(grepl("oo", useData$item, ignore.case = TRUE),
"GOOSE", "LOT")
useData$condition <- ifelse(
useData$condition == "filler" & !is.na(accent),
paste(accent, vowel, sep = "_"),
useData$condition
)
# Add accent and vowel variables
useData$accent <- sub("_.*", "", useData$condition)
useData$vowel <- sub(".*_", "", useData$condition)ggplot(useData, aes(x = vowel, y = rt)) +
geom_boxplot() + facet_wrap(.~ accent)
Reaction times (i.e. time interval between vowel onset and mouse click) appear broadly similar across vowels and accents, with comparable medians and variability, though there are several high outliers, particularly for LOT in AusE.
useData_trimmed <- useData %>%
group_by(pp, condition) %>%
mutate(
mean_rt = mean(rt, na.rm = TRUE),
sd_rt = sd(rt, na.rm = TRUE)
) %>%
filter(rt >= mean_rt - 2.5 * sd_rt,
rt <= mean_rt + 2.5 * sd_rt) %>%
ungroup()
ggplot(useData_trimmed, aes(x = rt)) +
geom_histogram() + facet_wrap(.~ condition)
useData_trimmed2 <- useData %>%
filter(rt <= 3000) %>%
droplevels()
ggplot(useData_trimmed2, aes(x = rt)) +
geom_histogram() + facet_wrap(.~ condition)
3.2 Analysis of the time course
The analysis is organised into three main blocks: (i) an overall analysis including both contrasts (FLEECE–KIT and DRESS–TRAP), (ii) a focused analysis of the FLEECE–KIT contrast, and (iii) a focused analysis of the DRESS–TRAP contrast.
Within each block, target and competitor fixations are analysed separately. For each fixation type, we conduct two complementary analyses. First, we apply a time-window analysis to examine differences in fixation proportions within the predefined time window. Second, we fit Generalised Additive Mixed Models (GAMMs) to capture the continuous time course of fixation behaviour.
3.2.1 Overall Effect of Contrast
3.2.1.1. TARGET FIXATIONS
A. TIME-WINDOW ANALYSIS
rm(list = ls())
library(lmerTest)
library(stringr)
library(tidyverse)
library(mgcv)
library(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("KIT|FLEECE|DRESS|TRAP", 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 == "FLEECE" | trialDataE$contrast == "KIT")
table(trialDataE$contrastC)##
## -0.5 0.5
## 3596 3598useData = cbind(trialDataE,gazeDataE)windowData = useData[ ,c(1,2,4:8,10:134)] %>%
pivot_longer(., cols = starts_with("T_"),
names_to = "Time",
values_to = "position") %>% as.data.frame
windowData$position = as.numeric(windowData$position)
windowData$Time = gsub("T_","", windowData$Time) %>% as.numeric()
windowData = windowData[ order(windowData$pp,windowData$trial,windowData$Time), ]
windowData$pp = as.factor(windowData$pp)
windowData$item = as.factor(windowData$item)
windowData$isAustralian = ifelse(windowData$accent == "AusE",1,0)
windowData$isAustralianFleece = windowData$isAustralian *
as.numeric(windowData$contrast == "FLEECE")
windowData$vowelContrast = ifelse(windowData$contrastC < 0, "dressTrap", "kitFleece") %>% as.factor()
#restricting the data to 0-900 after vowel onset
windowData = windowData %>% filter(Time > 0 & Time < 900)
windowData <- start_event(windowData, event=c("pp","trial"))
windowDataPlot = windowData %>% group_by(vowelContrast, Time, accent) %>%
summarize(meanFix = mean(position))
Fig.Target_Overall <- ggplot(windowDataPlot, aes(x = Time, y = meanFix, col = vowelContrast)) + geom_line() + facet_wrap(~accent)
Fig.Target_Overall
The two contrasts show very similar trajectories overall. In AusE, DRESS–TRAP appears to rise slightly more steeply and reaches a marginally higher asymptote than KIT–FLEECE, though the difference seems negligible. In SEng, the curves almost completely overlap across the time course, indicating hardly any differences between the two contrasts.
Next, we restrict the dataset to the time window between 150 and 600 ms.
lmerData = windowData %>% filter(Time > 150 & Time < 600) %>%
group_by(pp, item, trial, isAustralian, contrastC) %>% summarize(meanLooks = mean(position))In addition, proportion looks were logit-transformed prior to analysis. Because the logit function is undefined for proportions of exactly 0 or 1, a small correction was applied to these values before transformation. The transformation maps bounded proportions (0–1) onto the real line, resulting in an approximately linear relationship in the mid-range and greater separation of values near the boundaries.
lmerData$accentC = lmerData$isAustralian - 0.5
outcomes = sort(unique(lmerData$meanLooks))
if (outcomes[1] == 0){
correction = outcomes[2] * 0.5
}
lmerData$logitLooks = ifelse(lmerData$meanLooks == 0, qlogis(correction),
ifelse(lmerData$meanLooks == 1, qlogis(1- correction),
qlogis(lmerData$meanLooks))
)
plot(lmerData$meanLooks, lmerData$logitLooks)
lmerData$trialScaled = scale(lmerData$trial)
target_Overall.lmer = lmer(logitLooks ~ contrastC*accentC*trialScaled +
(1+contrastC+contrastC||pp) +
(1+accentC*trialScaled||item), data = lmerData,
control=lmerControl(optimizer="bobyqa",
optCtrl=list(maxfun=2e5)))summary(target_Overall.lmer)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logitLooks ~ contrastC * accentC * trialScaled + (1 + contrastC +
## contrastC || pp) + (1 + accentC * trialScaled || item)
## Data: lmerData
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## REML criterion at convergence: 34748.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.28676 -0.60502 -0.00796 0.56105 2.62053
##
## Random effects:
## Groups Name Variance Std.Dev.
## pp (Intercept) 0.5502189 0.74177
## pp.1 contrastC 0.0000000 0.00000
## item (Intercept) 0.1746046 0.41786
## item.1 accentC 0.1046005 0.32342
## item.2 trialScaled 0.0002669 0.01634
## item.3 accentC:trialScaled 0.0000000 0.00000
## Residual 7.1127009 2.66696
## Number of obs: 7194, groups: pp, 45; item, 40
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.444e-01 1.326e-01 6.543e+01 -1.089 0.280
## contrastC -7.680e-03 1.464e-01 3.802e+01 -0.052 0.958
## accentC 4.684e-02 8.108e-02 3.798e+01 0.578 0.567
## trialScaled 4.418e-02 3.167e-02 3.687e+01 1.395 0.171
## contrastC:accentC -3.857e-02 1.622e-01 3.798e+01 -0.238 0.813
## contrastC:trialScaled 8.604e-03 6.359e-02 3.745e+01 0.135 0.893
## accentC:trialScaled -1.931e-02 6.331e-02 7.117e+03 -0.305 0.760
## contrastC:accentC:trialScaled -1.800e-01 1.266e-01 7.116e+03 -1.422 0.155
##
## Correlation of Fixed Effects:
## (Intr) cntrsC accntC trlScl cntC:C cntC:S accC:S
## contrastC 0.000
## accentC 0.000 0.000
## trialScaled 0.000 -0.009 -0.021
## cntrstC:ccC 0.000 0.000 0.000 0.000
## cntrstC:trS -0.005 0.000 0.000 0.001 -0.021
## accntC:trlS -0.006 0.000 0.000 -0.004 -0.016 0.003
## cntrstC:C:S 0.000 -0.012 -0.016 0.003 0.000 -0.004 0.002
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')The linear mixed-effects model revealed no significant main effects or interactions of contrast, accent, or trial, indicating that none of these predictors reliably influenced logit-transformed looks in the analyzed time window.
B. Generalized Additive Mixed Model
gamData <- windowDataBuilding an initial GAMM.
fullModelFile = "target_overall_fullModels.Rdata"
if (file.exists(fullModelFile))
{
load(fullModelFile)
}else{
interaction.target.overall.gamm0 <- bam( position ~ vowelContrast+
s(Time, by = vowelContrast) +
s(Time, by = isAustralian) +
s(Time, by = isAustralianFleece) +
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.target.overall.gamm0)
myRhoVal = myAutoCorr[2]
interaction.target.overall.gamm1 <- bam( position ~ vowelContrast+
s(Time, by = vowelContrast, k=12) +
s(Time, by = isAustralian, k= 12) +
s(Time, by = isAustralianFleece) +
s(Time, pp, by = vowelContrast, bs="fs") +
s(Time, item, by = vowelContrast, bs="fs"),
family = "binomial",
data=gamData, rho = myRhoVal,
AR.start = gamData$start.event, nThreads = 3, discrete = T)
save(list= c("interaction.target.overall.gamm0","interaction.target.overall.gamm1"),
file = fullModelFile)
}Model criticism: How are we doing regarding auto-correlation?
myAutoCorr =acf_resid(interaction.target.overall.gamm1)
Is the default for the smooth (k = 10) sufficient?
gam.check(interaction.target.overall.gamm1)
##
## Method: fREML Optimizer: perf chol
## $grad
## [1] 4.601007e-04 5.483779e-05 -3.345560e-03 -2.430295e-03 1.087673e-03
## [6] 2.775880e-03 1.093321e-03 1.240446e-03 3.668401e-03 1.710888e-03
## [11] -8.197439e-05 1.447076e-03 7.881274e-04 7.402990e-04 1.392357e-03
## [16] 8.132879e-04
##
## $hess
## [,1] [,2] [,3] [,4] [,5]
## [1,] 3.635972e+00 -1.271894e-04 7.633337e-04 -1.749604e-06 1.137318e-01
## [2,] -1.271894e-04 3.879995e+00 6.609466e-04 3.257504e-04 -1.259463e-05
## [3,] 7.633337e-04 6.609466e-04 3.328010e-03 2.759053e-06 2.821343e-05
## [4,] -1.749604e-06 3.257504e-04 2.759053e-06 2.435581e-03 -1.401612e-07
## [5,] 1.137318e-01 -1.259463e-05 2.821343e-05 -1.401612e-07 1.292756e+01
## [6,] -3.698917e-02 -4.491365e-06 -1.487140e-05 4.632376e-08 -9.822703e-01
## [7,] -7.949496e-03 -4.356437e-06 -1.924396e-06 -1.386686e-07 1.262210e-01
## [8,] -7.431995e-06 3.791765e-02 5.568153e-06 -1.961521e-05 -1.805990e-05
## [9,] -4.111156e-06 -1.697391e-02 -3.563526e-06 1.674737e-06 -1.419291e-05
## [10,] -1.948643e-06 -1.272003e-02 -1.870184e-06 -3.357926e-06 -5.370002e-06
## [11,] 2.612137e-02 -1.716056e-06 3.275396e-06 -2.597939e-08 -8.467355e-03
## [12,] -5.583301e-03 -3.097468e-06 -6.767423e-06 -6.884038e-08 -1.873573e-02
## [13,] -7.195226e-03 -3.415163e-06 -1.740821e-06 -1.025153e-07 -1.811374e-02
## [14,] 5.063515e-06 5.559144e-02 3.378979e-05 2.602375e-03 1.846895e-06
## [15,] 4.405906e-05 -6.290424e-03 -1.504397e-06 -2.658043e-04 3.388592e-05
## [16,] 2.828652e-04 -1.760488e-02 -3.870113e-05 -2.025821e-05 8.301354e-04
## [,6] [,7] [,8] [,9] [,10]
## [1,] -3.698917e-02 -7.949496e-03 -7.431995e-06 -4.111156e-06 -1.948643e-06
## [2,] -4.491365e-06 -4.356437e-06 3.791765e-02 -1.697391e-02 -1.272003e-02
## [3,] -1.487140e-05 -1.924396e-06 5.568153e-06 -3.563526e-06 -1.870184e-06
## [4,] 4.632376e-08 -1.386686e-07 -1.961521e-05 1.674737e-06 -3.357926e-06
## [5,] -9.822703e-01 1.262210e-01 -1.805990e-05 -1.419291e-05 -5.370002e-06
## [6,] 1.261420e+01 -9.334767e-01 -1.398970e-05 -1.112082e-05 -4.313497e-06
## [7,] -9.334767e-01 2.658938e+00 -6.812602e-06 -5.220342e-06 -3.346492e-06
## [8,] -1.398970e-05 -6.812602e-06 1.016423e+01 -1.012778e+00 -8.774159e-02
## [9,] -1.112082e-05 -5.220342e-06 -1.012778e+00 1.367902e+01 -1.158049e+00
## [10,] -4.313497e-06 -3.346492e-06 -8.774159e-02 -1.158049e+00 5.061275e+00
## [11,] -9.262492e-03 -2.463990e-03 -2.101602e-06 -1.627035e-06 -8.588874e-07
## [12,] -3.384085e-02 -4.809104e-03 -6.379611e-06 -4.940295e-06 -2.957988e-06
## [13,] -1.634815e-02 -6.577264e-03 -5.572333e-06 -4.283565e-06 -2.637611e-06
## [14,] 4.725344e-06 3.191183e-05 -4.233348e-02 -2.373596e-02 -1.871917e-02
## [15,] 3.409327e-05 1.071020e-04 -1.614355e-02 -3.846921e-02 -1.291031e-02
## [16,] 6.432591e-04 2.982025e-04 -3.685074e-02 -2.476608e-02 -1.605395e-02
## [,11] [,12] [,13] [,14] [,15]
## [1,] 2.612137e-02 -5.583301e-03 -7.195226e-03 5.063515e-06 4.405906e-05
## [2,] -1.716056e-06 -3.097468e-06 -3.415163e-06 5.559144e-02 -6.290424e-03
## [3,] 3.275396e-06 -6.767423e-06 -1.740821e-06 3.378979e-05 -1.504397e-06
## [4,] -2.597939e-08 -6.884038e-08 -1.025153e-07 2.602375e-03 -2.658043e-04
## [5,] -8.467355e-03 -1.873573e-02 -1.811374e-02 1.846895e-06 3.388592e-05
## [6,] -9.262492e-03 -3.384085e-02 -1.634815e-02 4.725344e-06 3.409327e-05
## [7,] -2.463990e-03 -4.809104e-03 -6.577264e-03 3.191183e-05 1.071020e-04
## [8,] -2.101602e-06 -6.379611e-06 -5.572333e-06 -4.233348e-02 -1.614355e-02
## [9,] -1.627035e-06 -4.940295e-06 -4.283565e-06 -2.373596e-02 -3.846921e-02
## [10,] -8.588874e-07 -2.957988e-06 -2.637611e-06 -1.871917e-02 -1.291031e-02
## [11,] 1.099006e+00 -2.708893e-01 5.270369e-02 5.822178e-06 2.079864e-05
## [12,] -2.708893e-01 5.151946e+00 -6.586046e-01 2.575402e-05 8.717375e-05
## [13,] 5.270369e-02 -6.586046e-01 3.283708e+00 2.373937e-05 8.041618e-05
## [14,] 5.822178e-06 2.575402e-05 2.373937e-05 6.554145e+00 -1.045720e-02
## [15,] 2.079864e-05 8.717375e-05 8.041618e-05 -1.045720e-02 4.719529e+00
## [16,] 9.349311e-05 2.810872e-04 2.451063e-04 1.716178e-01 -5.130372e-01
## [,16]
## [1,] 2.828652e-04
## [2,] -1.760488e-02
## [3,] -3.870113e-05
## [4,] -2.025821e-05
## [5,] 8.301354e-04
## [6,] 6.432591e-04
## [7,] 2.982025e-04
## [8,] -3.685074e-02
## [9,] -2.476608e-02
## [10,] -1.605395e-02
## [11,] 9.349311e-05
## [12,] 2.810872e-04
## [13,] 2.451063e-04
## [14,] 1.716178e-01
## [15,] -5.130372e-01
## [16,] 4.890690e+00
##
## 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):vowelContrastdressTrap 11.00 8.11 0.96 0.69
## s(Time):vowelContrastkitFleece 11.00 8.79 0.96 0.67
## s(Time):isAustralian 12.00 2.01 0.96 0.66
## s(Time):isAustralianFleece 10.00 2.01 0.96 0.66
## s(Time,pp):vowelContrastdressTrap 450.00 121.51 0.96 0.68
## s(Time,pp):vowelContrastkitFleece 450.00 120.16 0.96 0.62
## s(Time,item):vowelContrastdressTrap 400.00 35.45 0.96 0.65
## s(Time,item):vowelContrastkitFleece 400.00 60.80 0.96 0.72We now examine this model to determine whether there is an overall effect of contrast.
summary(interaction.target.overall.gamm1)##
## Family: binomial
## Link function: logit
##
## Formula:
## position ~ vowelContrast + s(Time, by = vowelContrast, k = 12) +
## s(Time, by = isAustralian, k = 12) + s(Time, by = isAustralianFleece) +
## 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.60814 0.04973 -32.336 <2e-16 ***
## vowelContrastkitFleece 0.05554 0.07392 0.751 0.452
## ---
## 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):vowelContrastdressTrap 8.115 9.673 172.409 <2e-16 ***
## s(Time):vowelContrastkitFleece 8.790 10.119 232.412 <2e-16 ***
## s(Time):isAustralian 2.012 2.023 4.758 0.103
## s(Time):isAustralianFleece 2.011 2.021 2.775 0.232
## s(Time,pp):vowelContrastdressTrap 121.510 448.000 351.024 <2e-16 ***
## s(Time,pp):vowelContrastkitFleece 120.157 449.000 352.939 <2e-16 ***
## s(Time,item):vowelContrastdressTrap 35.447 198.000 142.210 <2e-16 ***
## s(Time,item):vowelContrastkitFleece 60.798 199.000 196.735 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0731 Deviance explained = 8.17%
## fREML = 1.6335e+05 Scale est. = 1 n = 640266There is no evidence for an overall (parametric) effect of contrast, given that the coefficient for vowelContrastkitFleece was not significant (p = .452).
smooth <- plot_smooth(interaction.target.overall.gamm1, view = "Time",
plot_all = "vowelContrast",
rug = FALSE, rm.ranef = TRUE, lwd = 4)## Summary:
## * vowelContrast : factor; set to the value(s): dressTrap, kitFleece.
## * Time : numeric predictor; with 30 values ranging from 10.000000 to 890.000000.
## * isAustralian : numeric predictor; set to the value(s): 0.
## * isAustralianFleece : 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):vowelContrastdressTrap,s(Time,pp):vowelContrastkitFleece,s(Time,item):vowelContrastdressTrap,s(Time,item):vowelContrastkitFleece
## 
We now proceed to model comparison by successively removing effects to evaluate their contribution to model fit.
reducedModelFile = "target_overall_reducedModels.Rdata"
if (file.exists(reducedModelFile))
{
load(reducedModelFile)
}else{
minusInteraction.target.overall.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,
AR.start = gamData$start.event, nThreads = 3, discrete = T)
minusAccent.target.overall.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,
AR.start = gamData$start.event, nThreads = 3, discrete = T)
timeOnly.target.overall.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,
AR.start = gamData$start.event, nThreads = 3, discrete = T)
save(list= c("minusInteraction.target.overall.gamm",
"minusAccent.target.overall.gamm",
"timeOnly.target.overall.gamm"),
file = reducedModelFile)
}Model comparisons:
# Full model versus no interaction
compareML(interaction.target.overall.gamm1, minusInteraction.target.overall.gamm)## interaction.target.overall.gamm1: position ~ vowelContrast + s(Time, by = vowelContrast, k = 12) +
## s(Time, by = isAustralian, k = 12) + s(Time, by = isAustralianFleece) +
## s(Time, pp, by = vowelContrast, bs = "fs") + s(Time, item,
## by = vowelContrast, bs = "fs")
##
## minusInteraction.target.overall.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.target.overall.gamm preferred: lower fREML score (107.393), and lower df (1.000).
## -----
## Model Score Edf Difference Df
## 1 interaction.target.overall.gamm1 163345.1 16
## 2 minusInteraction.target.overall.gamm 163237.7 15 -107.393 1.000
##
## AIC difference: 228.44, model minusInteraction.target.overall.gamm has lower AIC.# Both main effects versus minus accent
compareML(minusInteraction.target.overall.gamm, minusAccent.target.overall.gamm)## minusInteraction.target.overall.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.target.overall.gamm: position ~ vowelContrast + s(Time, by = vowelContrast) + s(Time,
## pp, by = vowelContrast, bs = "fs") + s(Time, item, by = vowelContrast,
## bs = "fs")
##
## Model minusAccent.target.overall.gamm preferred: lower fREML score (14.330), and lower df (1.000).
## -----
## Model Score Edf Difference Df
## 1 minusInteraction.target.overall.gamm 163237.7 15
## 2 minusAccent.target.overall.gamm 163223.4 14 -14.330 1.000
##
## AIC difference: 374.04, model minusAccent.target.overall.gamm has lower AIC.# Main effect of Contrast versus time-only
compareML(minusAccent.target.overall.gamm, timeOnly.target.overall.gamm)## minusAccent.target.overall.gamm: position ~ vowelContrast + s(Time, by = vowelContrast) + s(Time,
## pp, by = vowelContrast, bs = "fs") + s(Time, item, by = vowelContrast,
## bs = "fs")
##
## timeOnly.target.overall.gamm: position ~ 1 + s(Time) + s(Time, pp, by = vowelContrast, bs = "fs") +
## s(Time, item, by = vowelContrast, bs = "fs")
##
## Model timeOnly.target.overall.gamm preferred: lower fREML score (436.458), and lower df (1.000).
## -----
## Model Score Edf Difference Df
## 1 minusAccent.target.overall.gamm 163223.4 14
## 2 timeOnly.target.overall.gamm 162786.9 13 -436.458 1.000
##
## AIC difference: 1703.02, model timeOnly.target.overall.gamm has lower AIC.Model comparison shows that progressively simpler models fit the data better, with lower fREML scores and substantially lower AIC values at each step. The best-supported model includes only a general smooth of Time (with random smooths for participants and items), providing no evidence that vowel contrast, accent, or their interactions significantly improve model fit.
3.2.1.2. Competitor Fixations
A. TIME-WINDOW ANALYSIS
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("KIT|FLEECE|DRESS|TRAP", 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 == "FLEECE" | trialDataE$contrast == "KIT")
table(trialDataE$contrastC)##
## -0.5 0.5
## 3596 3598useData = cbind(trialDataE,gazeDataE)windowData = useData[ ,c(1,2,4:8,10:134)] %>%
pivot_longer(., cols = starts_with("T_"),
names_to = "Time",
values_to = "position") %>% as.data.frame
windowData$position = as.numeric(windowData$position)
windowData$Time = gsub("T_","", windowData$Time) %>% as.numeric()
windowData = windowData[ order(windowData$pp,windowData$trial,windowData$Time), ]
windowData$pp = as.factor(windowData$pp)
windowData$item = as.factor(windowData$item)
windowData$isAustralian = ifelse(windowData$accent == "AusE",1,0)
windowData$isAustralianFleece = windowData$isAustralian *
as.numeric(windowData$contrast == "FLEECE")
windowData$vowelContrast = ifelse(windowData$contrastC < 0, "dressTrap", "kitFleece") %>% as.factor()
#restricting the data to 0-900 after vowel onset
windowData = windowData %>% filter(Time > 0 & Time < 900)
windowData <- start_event(windowData, event=c("pp","trial"))
windowDataPlot = windowData %>% group_by(vowelContrast, Time, accent) %>%
summarize(meanFix = mean(position))
Fig.Target_Overall <- ggplot(windowDataPlot, aes(x = Time, y = meanFix, col = vowelContrast)) + geom_line() + facet_wrap(~accent)
Fig.Target_Overall
Across both accents, the KIT-FLEECT contrast is consistently slightly higher than the DRESS-TRAP contrast, especially around the peak, but the difference between accents themselves appears minimal.
Next, we restrict the dataset to the time window between 150 and 600 ms.
lmerData = windowData %>% filter(Time > 150 & Time < 600) %>%
group_by(pp, item, trial, isAustralian, contrastC) %>% summarize(meanLooks = mean(position))lmerData$accentC = lmerData$isAustralian - 0.5
outcomes = sort(unique(lmerData$meanLooks))
if (outcomes[1] == 0){
correction = outcomes[2] * 0.5
}
lmerData$logitLooks = ifelse(lmerData$meanLooks == 0, qlogis(correction),
ifelse(lmerData$meanLooks == 1, qlogis(1- correction),
qlogis(lmerData$meanLooks))
)
plot(lmerData$meanLooks, lmerData$logitLooks)
lmerData$trialScaled = scale(lmerData$trial)
compet_overall.lmer = lmer(logitLooks ~ contrastC*accentC*trialScaled +
(1+contrastC+contrastC||pp) +
(1+accentC*trialScaled||item), data = lmerData,
control=lmerControl(optimizer="bobyqa",
optCtrl=list(maxfun=2e5)))summary(compet_overall.lmer)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logitLooks ~ contrastC * accentC * trialScaled + (1 + contrastC +
## contrastC || pp) + (1 + accentC * trialScaled || item)
## Data: lmerData
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## REML criterion at convergence: 32277.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.46379 -0.95676 0.01107 0.88138 3.14977
##
## Random effects:
## Groups Name Variance Std.Dev.
## pp (Intercept) 0.04894 0.2212
## pp.1 contrastC 0.01785 0.1336
## item (Intercept) 0.05627 0.2372
## item.1 accentC 0.06724 0.2593
## item.2 trialScaled 0.01102 0.1050
## item.3 accentC:trialScaled 0.02940 0.1715
## Residual 5.09308 2.2568
## Number of obs: 7194, groups: pp, 45; item, 40
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.20838 0.05661 53.63943 -39.008 <2e-16 ***
## contrastC 0.26640 0.09416 36.95613 2.829 0.0075 **
## accentC -0.06662 0.06725 38.03244 -0.991 0.3282
## trialScaled -0.02683 0.03147 37.44282 -0.853 0.3994
## contrastC:accentC 0.06292 0.13451 38.03270 0.468 0.6426
## contrastC:trialScaled -0.05899 0.06304 37.66334 -0.936 0.3554
## accentC:trialScaled 0.06065 0.06004 35.80464 1.010 0.3192
## contrastC:accentC:trialScaled 0.03724 0.12007 35.78925 0.310 0.7582
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cntrsC accntC trlScl cntC:C cntC:S accC:S
## contrastC 0.000
## accentC 0.000 0.001
## trialScaled 0.000 -0.009 -0.018
## cntrstC:ccC 0.000 0.000 -0.001 0.000
## cntrstC:trS -0.008 0.000 0.000 0.001 -0.018
## accntC:trlS -0.011 0.000 0.000 -0.003 -0.014 0.002
## cntrstC:C:S 0.000 -0.013 -0.014 0.002 0.000 -0.003 0.002The analysis revealed a significant main effect of contrast, indicating that fixation proportions differed reliably between the two contrasts (FLEECE–KIT vs. DRESS–TRAP). Specifically, FLEECE-KIT showed higher fixation proportions than DRESS-TRAP.
There was no significant main effect of accent, suggesting that fixation behaviour did not differ between accents. Likewise, there was no effect of trial, indicating no systematic change in fixation proportions over the course of the experiment.
B. Generalized Additive Mixed Model
gamData <- windowDataBuilding an initial GAMM.
fullModelFile = "compet_overall_fullModels.Rdata"
if (file.exists(fullModelFile))
{
load(fullModelFile)
}else{
interaction.compet.overall.gamm0 <- bam( position ~ vowelContrast+
s(Time, by = vowelContrast) +
s(Time, by = isAustralian) +
s(Time, by = isAustralianFleece) +
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.compet.overall.gamm0)
myRhoVal = myAutoCorr[2]
interaction.compet.overall.gamm1 <- bam( position ~ vowelContrast+
s(Time, by = vowelContrast, k=12) +
s(Time, by = isAustralian, k= 12) +
s(Time, by = isAustralianFleece) +
s(Time, pp, by = vowelContrast, bs="fs") +
s(Time, item, by = vowelContrast, bs="fs"),
family = "binomial",
data=gamData, rho = myRhoVal,
AR.start = gamData$start.event, nThreads = 3, discrete = T)
save(list= c("interaction.compet.overall.gamm0","interaction.compet.overall.gamm1"),
file = fullModelFile)
}Model criticism: How are we doing regarding auto-correlation?
myAutoCorr =acf_resid(interaction.compet.overall.gamm1)
Is the default for the smooth (k = 10) sufficient?
gam.check(interaction.compet.overall.gamm1)
##
## Method: fREML Optimizer: perf chol
## $grad
## [1] 4.601007e-04 5.483779e-05 -3.345560e-03 -2.430295e-03 1.087673e-03
## [6] 2.775880e-03 1.093321e-03 1.240446e-03 3.668401e-03 1.710888e-03
## [11] -8.197439e-05 1.447076e-03 7.881274e-04 7.402990e-04 1.392357e-03
## [16] 8.132879e-04
##
## $hess
## [,1] [,2] [,3] [,4] [,5]
## [1,] 3.635972e+00 -1.271894e-04 7.633337e-04 -1.749604e-06 1.137318e-01
## [2,] -1.271894e-04 3.879995e+00 6.609466e-04 3.257504e-04 -1.259463e-05
## [3,] 7.633337e-04 6.609466e-04 3.328010e-03 2.759053e-06 2.821343e-05
## [4,] -1.749604e-06 3.257504e-04 2.759053e-06 2.435581e-03 -1.401612e-07
## [5,] 1.137318e-01 -1.259463e-05 2.821343e-05 -1.401612e-07 1.292756e+01
## [6,] -3.698917e-02 -4.491365e-06 -1.487140e-05 4.632376e-08 -9.822703e-01
## [7,] -7.949496e-03 -4.356437e-06 -1.924396e-06 -1.386686e-07 1.262210e-01
## [8,] -7.431995e-06 3.791765e-02 5.568153e-06 -1.961521e-05 -1.805990e-05
## [9,] -4.111156e-06 -1.697391e-02 -3.563526e-06 1.674737e-06 -1.419291e-05
## [10,] -1.948643e-06 -1.272003e-02 -1.870184e-06 -3.357926e-06 -5.370002e-06
## [11,] 2.612137e-02 -1.716056e-06 3.275396e-06 -2.597939e-08 -8.467355e-03
## [12,] -5.583301e-03 -3.097468e-06 -6.767423e-06 -6.884038e-08 -1.873573e-02
## [13,] -7.195226e-03 -3.415163e-06 -1.740821e-06 -1.025153e-07 -1.811374e-02
## [14,] 5.063515e-06 5.559144e-02 3.378979e-05 2.602375e-03 1.846895e-06
## [15,] 4.405906e-05 -6.290424e-03 -1.504397e-06 -2.658043e-04 3.388592e-05
## [16,] 2.828652e-04 -1.760488e-02 -3.870113e-05 -2.025821e-05 8.301354e-04
## [,6] [,7] [,8] [,9] [,10]
## [1,] -3.698917e-02 -7.949496e-03 -7.431995e-06 -4.111156e-06 -1.948643e-06
## [2,] -4.491365e-06 -4.356437e-06 3.791765e-02 -1.697391e-02 -1.272003e-02
## [3,] -1.487140e-05 -1.924396e-06 5.568153e-06 -3.563526e-06 -1.870184e-06
## [4,] 4.632376e-08 -1.386686e-07 -1.961521e-05 1.674737e-06 -3.357926e-06
## [5,] -9.822703e-01 1.262210e-01 -1.805990e-05 -1.419291e-05 -5.370002e-06
## [6,] 1.261420e+01 -9.334767e-01 -1.398970e-05 -1.112082e-05 -4.313497e-06
## [7,] -9.334767e-01 2.658938e+00 -6.812602e-06 -5.220342e-06 -3.346492e-06
## [8,] -1.398970e-05 -6.812602e-06 1.016423e+01 -1.012778e+00 -8.774159e-02
## [9,] -1.112082e-05 -5.220342e-06 -1.012778e+00 1.367902e+01 -1.158049e+00
## [10,] -4.313497e-06 -3.346492e-06 -8.774159e-02 -1.158049e+00 5.061275e+00
## [11,] -9.262492e-03 -2.463990e-03 -2.101602e-06 -1.627035e-06 -8.588874e-07
## [12,] -3.384085e-02 -4.809104e-03 -6.379611e-06 -4.940295e-06 -2.957988e-06
## [13,] -1.634815e-02 -6.577264e-03 -5.572333e-06 -4.283565e-06 -2.637611e-06
## [14,] 4.725344e-06 3.191183e-05 -4.233348e-02 -2.373596e-02 -1.871917e-02
## [15,] 3.409327e-05 1.071020e-04 -1.614355e-02 -3.846921e-02 -1.291031e-02
## [16,] 6.432591e-04 2.982025e-04 -3.685074e-02 -2.476608e-02 -1.605395e-02
## [,11] [,12] [,13] [,14] [,15]
## [1,] 2.612137e-02 -5.583301e-03 -7.195226e-03 5.063515e-06 4.405906e-05
## [2,] -1.716056e-06 -3.097468e-06 -3.415163e-06 5.559144e-02 -6.290424e-03
## [3,] 3.275396e-06 -6.767423e-06 -1.740821e-06 3.378979e-05 -1.504397e-06
## [4,] -2.597939e-08 -6.884038e-08 -1.025153e-07 2.602375e-03 -2.658043e-04
## [5,] -8.467355e-03 -1.873573e-02 -1.811374e-02 1.846895e-06 3.388592e-05
## [6,] -9.262492e-03 -3.384085e-02 -1.634815e-02 4.725344e-06 3.409327e-05
## [7,] -2.463990e-03 -4.809104e-03 -6.577264e-03 3.191183e-05 1.071020e-04
## [8,] -2.101602e-06 -6.379611e-06 -5.572333e-06 -4.233348e-02 -1.614355e-02
## [9,] -1.627035e-06 -4.940295e-06 -4.283565e-06 -2.373596e-02 -3.846921e-02
## [10,] -8.588874e-07 -2.957988e-06 -2.637611e-06 -1.871917e-02 -1.291031e-02
## [11,] 1.099006e+00 -2.708893e-01 5.270369e-02 5.822178e-06 2.079864e-05
## [12,] -2.708893e-01 5.151946e+00 -6.586046e-01 2.575402e-05 8.717375e-05
## [13,] 5.270369e-02 -6.586046e-01 3.283708e+00 2.373937e-05 8.041618e-05
## [14,] 5.822178e-06 2.575402e-05 2.373937e-05 6.554145e+00 -1.045720e-02
## [15,] 2.079864e-05 8.717375e-05 8.041618e-05 -1.045720e-02 4.719529e+00
## [16,] 9.349311e-05 2.810872e-04 2.451063e-04 1.716178e-01 -5.130372e-01
## [,16]
## [1,] 2.828652e-04
## [2,] -1.760488e-02
## [3,] -3.870113e-05
## [4,] -2.025821e-05
## [5,] 8.301354e-04
## [6,] 6.432591e-04
## [7,] 2.982025e-04
## [8,] -3.685074e-02
## [9,] -2.476608e-02
## [10,] -1.605395e-02
## [11,] 9.349311e-05
## [12,] 2.810872e-04
## [13,] 2.451063e-04
## [14,] 1.716178e-01
## [15,] -5.130372e-01
## [16,] 4.890690e+00
##
## 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):vowelContrastdressTrap 11.00 8.11 0.93 0.100 .
## s(Time):vowelContrastkitFleece 11.00 8.79 0.93 0.115
## s(Time):isAustralian 12.00 2.01 0.93 0.110
## s(Time):isAustralianFleece 10.00 2.01 0.93 0.105
## s(Time,pp):vowelContrastdressTrap 450.00 121.51 0.93 0.065 .
## s(Time,pp):vowelContrastkitFleece 450.00 120.16 0.93 0.115
## s(Time,item):vowelContrastdressTrap 400.00 35.45 0.93 0.095 .
## s(Time,item):vowelContrastkitFleece 400.00 60.80 0.93 0.090 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We now examine this model to determine whether there is an overall effect of contrast.
summary(interaction.compet.overall.gamm1)##
## Family: binomial
## Link function: logit
##
## Formula:
## position ~ vowelContrast + s(Time, by = vowelContrast, k = 12) +
## s(Time, by = isAustralian, k = 12) + s(Time, by = isAustralianFleece) +
## 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.60814 0.04973 -32.336 <2e-16 ***
## vowelContrastkitFleece 0.05554 0.07392 0.751 0.452
## ---
## 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):vowelContrastdressTrap 8.115 9.673 172.409 <2e-16 ***
## s(Time):vowelContrastkitFleece 8.790 10.119 232.412 <2e-16 ***
## s(Time):isAustralian 2.012 2.023 4.758 0.103
## s(Time):isAustralianFleece 2.011 2.021 2.775 0.232
## s(Time,pp):vowelContrastdressTrap 121.510 448.000 351.024 <2e-16 ***
## s(Time,pp):vowelContrastkitFleece 120.157 449.000 352.939 <2e-16 ***
## s(Time,item):vowelContrastdressTrap 35.447 198.000 142.210 <2e-16 ***
## s(Time,item):vowelContrastkitFleece 60.798 199.000 196.735 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0731 Deviance explained = 8.17%
## fREML = 1.6335e+05 Scale est. = 1 n = 640266The binomial GAMM shows significant nonlinear time effects for both vowel contrasts and substantial participant and item variability over time, but no significant main effect of vowel contrast nor any significant accent-related smooths, indicating that while fixation trajectories change over time, they do not reliably differ by contrast or accent.
smooth <- plot_smooth(interaction.compet.overall.gamm1, view = "Time",
plot_all = "vowelContrast",
rug = FALSE, rm.ranef = TRUE, lwd = 4)## Summary:
## * vowelContrast : factor; set to the value(s): dressTrap, kitFleece.
## * Time : numeric predictor; with 30 values ranging from 10.000000 to 890.000000.
## * isAustralian : numeric predictor; set to the value(s): 0.
## * isAustralianFleece : 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):vowelContrastdressTrap,s(Time,pp):vowelContrastkitFleece,s(Time,item):vowelContrastdressTrap,s(Time,item):vowelContrastkitFleece
##
We now proceed to model comparison by successively removing effects to evaluate their contribution to model fit.
reducedModelFile = "compet_overall_reducedModels.Rdata"
if (file.exists(reducedModelFile))
{
load(reducedModelFile)
}else{
minusInteraction.compet.overall.gamm <- bam( position ~ vowelContrast+
s(Time, by = vowelContrast) +
s(Time, by = isAustralian) +
#s(Time, by = isAustralianFleece) +
s(Time, pp, by = vowelContrast, bs="fs") +
s(Time, item, by = vowelContrast, bs="fs"),
family = "binomial",
data=gamData, rho = myRhoVal,
AR.start = gamData$start.event, nThreads = 3, discrete = T)
minusAccent.compet.overall.gamm <- bam( position ~ vowelContrast+
s(Time, by = vowelContrast) +
#s(Time, by = isAustralian) +
#s(Time, by = isAustralianFleece) +
s(Time, pp, by = vowelContrast, bs="fs") +
s(Time, item, by = vowelContrast, bs="fs"),
family = "binomial",
data=gamData, rho = myRhoVal,
AR.start = gamData$start.event, nThreads = 3, discrete = T)
timeOnly.compet.overall.gamm <- bam( position ~ 1 +
s(Time) +
#s(Time, by = isAustralian) +
#s(Time, by = isAustralianFleece) +
s(Time, pp, by = vowelContrast, bs="fs") +
s(Time, item, by = vowelContrast, bs="fs"), family = "binomial",
data=gamData, rho = myRhoVal,
AR.start = gamData$start.event, nThreads = 3, discrete = T)
save(list= c("minusInteraction.compet.overall.gamm",
"minusAccent.compet.overall.gamm",
"timeOnly.compet.overall.gamm"),
file = reducedModelFile)
}Model comparisons:
# Full model versus no interaction
compareML(interaction.compet.overall.gamm1, minusInteraction.compet.overall.gamm)## interaction.compet.overall.gamm1: position ~ vowelContrast + s(Time, by = vowelContrast, k = 12) +
## s(Time, by = isAustralian, k = 12) + s(Time, by = isAustralianFleece) +
## s(Time, pp, by = vowelContrast, bs = "fs") + s(Time, item,
## by = vowelContrast, bs = "fs")
##
## minusInteraction.compet.overall.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.compet.overall.gamm preferred: lower fREML score (107.393), and lower df (1.000).
## -----
## Model Score Edf Difference Df
## 1 interaction.compet.overall.gamm1 163345.1 16
## 2 minusInteraction.compet.overall.gamm 163237.7 15 -107.393 1.000
##
## AIC difference: 228.44, model minusInteraction.compet.overall.gamm has lower AIC.# Both main effects versus minus accent
compareML(minusInteraction.compet.overall.gamm, minusAccent.compet.overall.gamm)## minusInteraction.compet.overall.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.compet.overall.gamm: position ~ vowelContrast + s(Time, by = vowelContrast) + s(Time,
## pp, by = vowelContrast, bs = "fs") + s(Time, item, by = vowelContrast,
## bs = "fs")
##
## Model minusAccent.compet.overall.gamm preferred: lower fREML score (14.330), and lower df (1.000).
## -----
## Model Score Edf Difference Df
## 1 minusInteraction.compet.overall.gamm 163237.7 15
## 2 minusAccent.compet.overall.gamm 163223.4 14 -14.330 1.000
##
## AIC difference: 374.04, model minusAccent.compet.overall.gamm has lower AIC.# Main effect of Contrast versus time-only
compareML(minusAccent.compet.overall.gamm, timeOnly.compet.overall.gamm)## minusAccent.compet.overall.gamm: position ~ vowelContrast + s(Time, by = vowelContrast) + s(Time,
## pp, by = vowelContrast, bs = "fs") + s(Time, item, by = vowelContrast,
## bs = "fs")
##
## timeOnly.compet.overall.gamm: position ~ 1 + s(Time) + s(Time, pp, by = vowelContrast, bs = "fs") +
## s(Time, item, by = vowelContrast, bs = "fs")
##
## Model timeOnly.compet.overall.gamm preferred: lower fREML score (436.458), and lower df (1.000).
## -----
## Model Score Edf Difference Df
## 1 minusAccent.compet.overall.gamm 163223.4 14
## 2 timeOnly.compet.overall.gamm 162786.9 13 -436.458 1.000
##
## AIC difference: 1703.02, model timeOnly.compet.overall.gamm has lower AIC.Model comparisons again show that progressively simpler models provide a better fit, with lower fREML scores and substantially lower AIC values at each step. The best-supported model is the time-only model (with random smooths for participants and items), indicating no reliable effects of vowel contrast, accent, or their interaction on the trajectory.
3.2.2 FLEECE-KIT contrast
3.2.2.1. Target fixations
A. TIME-WINDOW ANALYSIS
rm(list = ls())
library(lmerTest)
library(stringr)
library(tidyverse)
library(mgcv)
library(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("KIT|FLEECE", 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
## 1799 1799trialDataE$contrastC = -0.5 + as.numeric(trialDataE$contrast == "FLEECE")
table(trialDataE$contrastC)##
## -0.5 0.5
## 1799 1799useData = cbind(trialDataE,gazeDataE)We first need to reorganize the data.
windowData = useData[ ,c(1,2,4:8,10:134)] %>%
pivot_longer(., cols = starts_with("T_"),
names_to = "Time",
values_to = "position") %>% as.data.frame
windowData$position = as.numeric(windowData$position)
windowData$Time = gsub("T_","", windowData$Time) %>% as.numeric()
windowData = windowData[ order(windowData$pp,windowData$trial,windowData$Time), ]
windowData$pp = as.factor(windowData$pp)
windowData$item = as.factor(windowData$item)
windowData$isAustralian = ifelse(windowData$accent == "AusE",1,0)
windowData$isAustralianFleece = windowData$isAustralian *
as.numeric(windowData$contrast == "FLEECE")
windowData$vowelContrast = windowData$contrast %>% as.factor()
#restricting the data to 0-900 after vowel onset
windowData = windowData %>% filter(Time > 0 & Time < 900)
windowData <- start_event(windowData, event=c("pp","trial"))
windowDataPlot = windowData %>% group_by(vowelContrast, Time, accent) %>%
summarize(meanFix = mean(position))
Fig.target_kitFleece <- ggplot(windowDataPlot, aes(x = Time, y = meanFix, col = vowelContrast)) + geom_line() + facet_wrap(~accent)
print(Fig.target_kitFleece)
aggregatedData = windowData %>% filter(Time > 200 & Time < 600) %>%
group_by(pp, item, contrast) %>% summarize(meanLooks = mean(position)) %>%
ungroup()
meanLooksValues = sort(unique(aggregatedData$meanLooks))
correction = meanLooksValues[2]/2
aggregatedData$logitLooks = ifelse( aggregatedData$meanLooks == 0,
qlogis(correction),
ifelse(aggregatedData$meanLooks == 1,
qlogis(1 - correction),
qlogis(aggregatedData$meanLooks)
)
)
plot(aggregatedData$meanLooks, aggregatedData$logitLooks)
library(lmerTest)
target_fleeceKit.lmer = lmer(meanLooks ~ contrast +
(1+contrast|pp)+
(1|item), data = aggregatedData)
summary(target_fleeceKit.lmer)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanLooks ~ contrast + (1 + contrast | pp) + (1 | item)
## Data: aggregatedData
##
## REML criterion at convergence: -401.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.97962 -0.63052 -0.00039 0.66759 2.62561
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## pp (Intercept) 0.008309 0.09115
## contrastKIT 0.002434 0.04934 -0.36
## item (Intercept) 0.003378 0.05812
## Residual 0.032292 0.17970
## Number of obs: 900, groups: pp, 45; item, 20
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.47826 0.02438 32.63326 19.62 <2e-16 ***
## contrastKIT 0.05291 0.02955 19.96635 1.79 0.0886 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## contrastKIT -0.619B. Generalized Additive Mixed Model
gamData <- windowDatafullModelFile = "target_kitFleece_fullModels.Rdata"
if (file.exists(fullModelFile))
{
load(fullModelFile)
}else{
interaction.target_kitFleece.gamm0 <- bam( position ~ vowelContrast+
s(Time, by = vowelContrast) +
s(Time, by = isAustralian) +
s(Time, by = isAustralianFleece) +
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.target_kitFleece.gamm0)
myRhoVal = myAutoCorr[2]
interaction.target_kitFleece.gamm1 <- bam( position ~ vowelContrast+
s(Time, by = vowelContrast, k=12) +
s(Time, by = isAustralian, k= 12) +
s(Time, by = isAustralianFleece) +
s(Time, pp, by = vowelContrast, bs="fs") +
s(Time, item, by = vowelContrast, bs="fs"), family = "binomial",
data=gamData, rho = myRhoVal,
AR.start = gamData$start.event, nThreads = 3, discrete = T)
save(list= c("interaction.target_kitFleece.gamm0","interaction.target_kitFleece.gamm1"), file = fullModelFile)
}Model criticism: How are we doing regarding auto-correlation?
myAutoCorr =acf_resid(interaction.target_kitFleece.gamm1)
Is the default for the smooth (k = 10) sufficient?
gam.check(interaction.target_kitFleece.gamm1)
##
## Method: fREML Optimizer: perf chol
## $grad
## [1] -4.063860e-12 2.150946e-11 -2.189917e-03 -1.377428e-03 -2.582112e-11
## [6] 3.940315e-11 1.943476e-10 -3.286011e-10 1.482423e-10 1.460698e-10
## [11] -1.557332e-11 1.327738e-11 1.055689e-10 -1.530118e-10 5.630341e-11
## [16] -1.641086e-09
##
## $hess
## [,1] [,2] [,3] [,4] [,5]
## [1,] 2.263304e+00 -1.837073e-04 7.009594e-04 3.679420e-04 9.012070e-02
## [2,] -1.837073e-04 3.124999e+00 9.377747e-04 -1.371116e-06 -1.397019e-05
## [3,] 7.009594e-04 9.377747e-04 2.210810e-03 2.028191e-05 6.695191e-05
## [4,] 3.679420e-04 -1.371116e-06 2.028191e-05 1.388854e-03 2.845913e-05
## [5,] 9.012070e-02 -1.397019e-05 6.695191e-05 2.845913e-05 1.111587e+01
## [6,] 5.478076e-02 -4.345524e-06 3.422655e-05 2.682132e-05 -3.568181e-02
## [7,] -6.503518e-02 4.901514e-06 -1.496909e-05 -7.396163e-06 -5.960655e-01
## [8,] -4.758950e-06 1.556852e-01 1.382991e-05 5.263263e-08 -4.360170e-07
## [9,] -1.300686e-06 2.573962e-02 -2.151230e-06 1.725119e-08 3.967844e-08
## [10,] 1.602039e-06 -5.446143e-02 -2.712552e-05 1.061477e-07 1.651173e-07
## [11,] 1.460777e-01 -1.339621e-05 5.475124e-05 2.903803e-05 -4.583996e-02
## [12,] 8.416323e-03 -5.912161e-07 8.215272e-06 7.494480e-06 -2.994154e-02
## [13,] -2.978886e-02 2.181263e-06 -8.809437e-06 -5.027699e-06 -3.515644e-02
## [14,] -2.561218e-06 4.459439e-02 6.223590e-06 2.416267e-08 -1.886760e-07
## [15,] -4.419886e-07 6.805948e-03 4.275629e-06 -1.616029e-08 -1.943994e-08
## [16,] 5.099160e-08 -1.223311e-02 -4.769035e-06 2.563815e-08 8.134387e-09
## [,6] [,7] [,8] [,9] [,10]
## [1,] 5.478076e-02 -6.503518e-02 -4.758950e-06 -1.300686e-06 1.602039e-06
## [2,] -4.345524e-06 4.901514e-06 1.556852e-01 2.573962e-02 -5.446143e-02
## [3,] 3.422655e-05 -1.496909e-05 1.382991e-05 -2.151230e-06 -2.712552e-05
## [4,] 2.682132e-05 -7.396163e-06 5.263263e-08 1.725119e-08 1.061477e-07
## [5,] -3.568181e-02 -5.960655e-01 -4.360170e-07 3.967844e-08 1.651173e-07
## [6,] 8.434833e+00 -3.678398e-01 -6.830323e-08 -2.019510e-08 8.108250e-08
## [7,] -3.678398e-01 1.193362e+01 1.362458e-07 3.586428e-08 -3.054613e-08
## [8,] -6.830323e-08 1.362458e-07 1.300760e+01 -2.998959e-01 -3.091102e-01
## [9,] -2.019510e-08 3.586428e-08 -2.998959e-01 1.099669e+01 -5.484860e-01
## [10,] 8.108250e-08 -3.054613e-08 -3.091102e-01 -5.484860e-01 1.212450e+01
## [11,] -3.082308e-02 -6.105782e-02 -3.424316e-07 -8.235563e-08 1.374512e-07
## [12,] -3.189537e-02 -2.063602e-02 5.045476e-10 -4.580776e-09 1.894908e-08
## [13,] -1.794652e-02 -5.552779e-02 6.073195e-08 1.203298e-08 -1.820422e-08
## [14,] -4.027986e-08 7.321332e-08 -2.475766e-02 -1.637855e-02 -3.362429e-02
## [15,] -1.585963e-08 8.337106e-09 -2.366714e-02 -2.242471e-02 -1.526583e-02
## [16,] 1.412255e-08 -1.289654e-10 -1.068923e-02 -5.424499e-03 -3.289267e-02
## [,11] [,12] [,13] [,14] [,15]
## [1,] 1.460777e-01 8.416323e-03 -2.978886e-02 -2.561218e-06 -4.419886e-07
## [2,] -1.339621e-05 -5.912161e-07 2.181263e-06 4.459439e-02 6.805948e-03
## [3,] 5.475124e-05 8.215272e-06 -8.809437e-06 6.223590e-06 4.275629e-06
## [4,] 2.903803e-05 7.494480e-06 -5.027699e-06 2.416267e-08 -1.616029e-08
## [5,] -4.583996e-02 -2.994154e-02 -3.515644e-02 -1.886760e-07 -1.943994e-08
## [6,] -3.082308e-02 -3.189537e-02 -1.794652e-02 -4.027986e-08 -1.585963e-08
## [7,] -6.105782e-02 -2.063602e-02 -5.552779e-02 7.321332e-08 8.337106e-09
## [8,] -3.424316e-07 5.045476e-10 6.073195e-08 -2.475766e-02 -2.366714e-02
## [9,] -8.235563e-08 -4.580776e-09 1.203298e-08 -1.637855e-02 -2.242471e-02
## [10,] 1.374512e-07 1.894908e-08 -1.820422e-08 -3.362429e-02 -1.526583e-02
## [11,] 4.297629e+00 -2.462764e-02 -1.542217e-01 -1.774305e-07 -3.417733e-08
## [12,] -2.462764e-02 2.060116e+00 -1.788043e-02 -1.789512e-09 -4.008611e-09
## [13,] -1.542217e-01 -1.788043e-02 3.164284e+00 3.205089e-08 4.303616e-09
## [14,] -1.774305e-07 -1.789512e-09 3.205089e-08 2.392398e+00 -7.592085e-02
## [15,] -3.417733e-08 -4.008611e-09 4.303616e-09 -7.592085e-02 1.481738e+00
## [16,] 8.484911e-09 4.086340e-09 -4.811789e-10 -1.819315e-01 8.522086e-02
## [,16]
## [1,] 5.099160e-08
## [2,] -1.223311e-02
## [3,] -4.769035e-06
## [4,] 2.563815e-08
## [5,] 8.134387e-09
## [6,] 1.412255e-08
## [7,] -1.289655e-10
## [8,] -1.068923e-02
## [9,] -5.424499e-03
## [10,] -3.289267e-02
## [11,] 8.484911e-09
## [12,] 4.086340e-09
## [13,] -4.811789e-10
## [14,] -1.819315e-01
## [15,] 8.522086e-02
## [16,] 1.864432e+00
##
## Model rank = 1346 / 1346
##
## 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):vowelContrastFLEECE 11.00 6.65 1 0.52
## s(Time):vowelContrastKIT 11.00 7.84 1 0.53
## s(Time):isAustralian 12.00 2.02 1 0.51
## s(Time):isAustralianFleece 10.00 2.01 1 0.54
## s(Time,pp):vowelContrastFLEECE 450.00 113.58 1 0.55
## s(Time,pp):vowelContrastKIT 450.00 122.15 1 0.53
## s(Time,item):vowelContrastFLEECE 200.00 32.98 1 0.49
## s(Time,item):vowelContrastKIT 200.00 20.55 1 0.53We now examine this model to determine whether there is an overall effect of contrast.
summary(interaction.target_kitFleece.gamm1)##
## Family: binomial
## Link function: logit
##
## Formula:
## position ~ vowelContrast + s(Time, by = vowelContrast, k = 12) +
## s(Time, by = isAustralian, k = 12) + s(Time, by = isAustralianFleece) +
## 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.19097 0.09737 1.961 0.0499 *
## vowelContrastKIT 0.10656 0.12749 0.836 0.4032
## ---
## 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):vowelContrastFLEECE 6.655 8.298 282.646 <2e-16 ***
## s(Time):vowelContrastKIT 7.837 9.474 375.653 <2e-16 ***
## s(Time):isAustralian 2.017 2.033 1.341 0.441
## s(Time):isAustralianFleece 2.010 2.020 5.197 0.099 .
## s(Time,pp):vowelContrastFLEECE 113.582 448.000 374.373 <2e-16 ***
## s(Time,pp):vowelContrastKIT 122.151 448.000 429.997 <2e-16 ***
## s(Time,item):vowelContrastFLEECE 32.983 98.000 139.294 <2e-16 ***
## s(Time,item):vowelContrastKIT 20.545 98.000 79.259 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.25 Deviance explained = 20%
## fREML = 1.0578e+05 Scale est. = 1 n = 320222The model shows nonlinear time effects and substantial participant and item variability, but no significant main effect of vowel contrast and no reliable accent-related effects, indicating that fixation trajectories change over time but do not significantly differ between KIT and FLEECE or across accents.
smooth <- plot_smooth(interaction.target_kitFleece.gamm1, view = "Time",
plot_all = "vowelContrast",
rug = FALSE, rm.ranef = TRUE, lwd = 4)## Summary:
## * vowelContrast : factor; set to the value(s): FLEECE, KIT.
## * Time : numeric predictor; with 30 values ranging from 10.000000 to 890.000000.
## * isAustralian : numeric predictor; set to the value(s): 0.
## * isAustralianFleece : 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): beetle. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(Time,pp):vowelContrastFLEECE,s(Time,pp):vowelContrastKIT,s(Time,item):vowelContrastFLEECE,s(Time,item):vowelContrastKIT
## 
We now proceed to model comparison by successively removing effects to evaluate their contribution to model fit.
reducedModelFile = "target_kitFleece_reducedModel.Rdata"
if (file.exists(reducedModelFile))
{
load(reducedModelFile)
}else{
mainEffect.target_kitFleece.gamm <- bam( position ~ vowelContrast+
s(Time, by = vowelContrast) +
#s(Time, by = isAustralian) +
#s(Time, by = isAustralianFleece) +
s(Time, pp, by = vowelContrast, bs="fs") +
s(Time, item, by = vowelContrast, bs="fs"), family = "binomial",
data=gamData, rho = myRhoVal,
AR.start = gamData$start.event, nThreads = 3, discrete = T)
save(list= c("mainEffect.target_kitFleece.gamm"), file = reducedModelFile)
}summary(mainEffect.target_kitFleece.gamm)##
## Family: binomial
## Link function: logit
##
## Formula:
## position ~ vowelContrast + s(Time, by = vowelContrast) + 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.17262 0.09188 1.879 0.0603 .
## vowelContrastKIT 0.14288 0.11890 1.202 0.2295
## ---
## 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):vowelContrastFLEECE 6.120 7.343 363.1 <2e-16 ***
## s(Time):vowelContrastKIT 6.958 8.121 415.2 <2e-16 ***
## s(Time,pp):vowelContrastFLEECE 112.706 448.000 690.1 <2e-16 ***
## s(Time,pp):vowelContrastKIT 121.372 448.000 708.4 <2e-16 ***
## s(Time,item):vowelContrastFLEECE 32.751 99.000 197.1 <2e-16 ***
## s(Time,item):vowelContrastKIT 20.477 98.000 250.0 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.252 Deviance explained = 20.2%
## fREML = 1.0551e+05 Scale est. = 1 n = 320222There is no significant overall effect of vowel contrast (KIT vs. FLEECE), as the parametric term for vowelContrastKIT is not significant (p = .23). This indicates that, averaged across the time window, fixation proportions do not differ reliably between KIT and FLEECE items.
smooth <- plot_smooth(mainEffect.target_kitFleece.gamm, view = "Time",
plot_all = "vowelContrast",
rug = FALSE, rm.ranef = TRUE, lwd = 4)## Summary:
## * vowelContrast : factor; set to the value(s): FLEECE, KIT.
## * Time : numeric predictor; with 30 values ranging from 10.000000 to 890.000000.
## * pp : factor; set to the value(s): CAT01. (Might be canceled as random effect, check below.)
## * item : factor; set to the value(s): beetle. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(Time,pp):vowelContrastFLEECE,s(Time,pp):vowelContrastKIT,s(Time,item):vowelContrastFLEECE,s(Time,item):vowelContrastKIT
## 
3.2.2.2. Competitor Fixations
A. TIME-WINDOW ANALYSIS
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("KIT|FLEECE", 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
## 1799 1799trialDataE$contrastC = -0.5 + as.numeric(trialDataE$contrast == "FLEECE")
table(trialDataE$contrastC)##
## -0.5 0.5
## 1799 1799useData = cbind(trialDataE,gazeDataE)We first need to reorganize the data.
windowData = useData[ ,c(1,2,4:8,10:134)] %>%
pivot_longer(., cols = starts_with("T_"),
names_to = "Time",
values_to = "position") %>% as.data.frame
windowData$position = as.numeric(windowData$position)
windowData$Time = gsub("T_","", windowData$Time) %>% as.numeric()
windowData = windowData[ order(windowData$pp,windowData$trial,windowData$Time), ]
windowData$pp = as.factor(windowData$pp)
windowData$item = as.factor(windowData$item)
windowData$isAustralian = ifelse(windowData$accent == "AusE",1,0)
windowData$isAustralianFleece = windowData$isAustralian *
as.numeric(windowData$contrast == "FLEECE")
windowData$vowelContrast = windowData$contrast %>% as.factor()
#restricting the data to 0-900 after vowel onset
windowData = windowData %>% filter(Time > 0 & Time < 900)
windowData <- start_event(windowData, event=c("pp","trial"))
windowDataPlot = windowData %>% group_by(vowelContrast, Time, accent) %>%
summarize(meanFix = mean(position))
Fig.compet_kitFleece <- ggplot(windowDataPlot, aes(x = Time, y = meanFix, col = vowelContrast)) + geom_line() + facet_wrap(~accent)
print(Fig.compet_kitFleece)
In both accents (AusE and SEng), fixation proportions to the competitor rise to a peak around 250–300 ms and then decline, with FLEECE showing consistently slightly higher competitor fixation proportions than KIT across most of the time course, and only minimal differences between accents.
lmerData = windowData %>% filter(Time > 150 & Time < 600) %>%
group_by(pp, item, trial, isAustralian, contrastC) %>% summarize(meanLooks = mean(position))
lmerData$accentC = lmerData$isAustralian - 0.5
outcomes = sort(unique(lmerData$meanLooks))
if (outcomes[1] == 0){
correction = outcomes[2] * 0.5
}
lmerData$logitLooks = ifelse(lmerData$meanLooks == 0, qlogis(correction),
ifelse(lmerData$meanLooks == 1, qlogis(1- correction),
qlogis(lmerData$meanLooks))
)
plot(lmerData$meanLooks, lmerData$logitLooks)
lmerData$trialScaled = scale(lmerData$trial)
compet_fleeceKit.lmer = lmer(logitLooks ~ contrastC*accentC*trialScaled +
(1+contrastC+contrastC||pp) +
(1+accentC*trialScaled||item), data = lmerData,
control=lmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))
summary(compet_fleeceKit.lmer)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logitLooks ~ contrastC * accentC * trialScaled + (1 + contrastC +
## contrastC || pp) + (1 + accentC * trialScaled || item)
## Data: lmerData
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## REML criterion at convergence: 16295.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4813 -0.9655 0.1864 0.8378 3.1328
##
## Random effects:
## Groups Name Variance Std.Dev.
## pp (Intercept) 0.05718 0.2391
## pp.1 contrastC 0.05394 0.2322
## item (Intercept) 0.06489 0.2547
## item.1 accentC 0.12455 0.3529
## item.2 trialScaled 0.02393 0.1547
## item.3 accentC:trialScaled 0.03659 0.1913
## Residual 5.27238 2.2962
## Number of obs: 3598, groups: pp, 45; item, 20
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.07524 0.07737 25.47510 -26.823 <2e-16 ***
## contrastC 0.39072 0.14163 19.25931 2.759 0.0124 *
## accentC -0.03409 0.11003 18.01479 -0.310 0.7602
## trialScaled -0.05497 0.05181 17.98592 -1.061 0.3028
## contrastC:accentC 0.16490 0.22007 18.01478 0.749 0.4633
## contrastC:trialScaled 0.03384 0.10367 18.01836 0.326 0.7478
## accentC:trialScaled 0.07335 0.08833 16.74850 0.830 0.4180
## contrastC:accentC:trialScaled -0.14668 0.17666 16.75335 -0.830 0.4180
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cntrsC accntC trlScl cntC:C cntC:S accC:S
## contrastC 0.000
## accentC 0.000 0.000
## trialScaled 0.000 0.001 -0.014
## cntrstC:ccC 0.000 0.000 0.000 -0.013
## cntrstC:trS 0.001 0.000 -0.013 0.002 -0.014
## accntC:trlS -0.011 -0.012 0.001 0.000 0.002 -0.017
## cntrstC:C:S -0.011 -0.012 0.002 -0.017 0.001 0.001 0.003There is a significant main effect of contrast (p = .0137), indicating that logit-transformed looks differ reliably between the two vowel contrasts, but there are no significant effects of accent, trial, or any interactions, meaning the contrast effect does not vary across accents or over trials.
B. Generalized Additive Mixed Model
gamData <- windowDatafullModelFile = "compet_kitFleece_fullModels.Rdata"
if (file.exists(fullModelFile))
{
load(fullModelFile)
}else{
interaction.compet_kitFleece.gamm0 <- bam( position ~ vowelContrast+
s(Time, by = vowelContrast) +
s(Time, by = isAustralian) +
s(Time, by = isAustralianFleece) +
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.compet_kitFleece.gamm0)
myRhoVal = myAutoCorr[2]
interaction.compet_kitFleece.gamm1 <- bam( position ~ vowelContrast+
s(Time, by = vowelContrast, k=12) +
s(Time, by = isAustralian, k= 12) +
s(Time, by = isAustralianFleece) +
s(Time, pp, by = vowelContrast, bs="fs") +
s(Time, item, by = vowelContrast, bs="fs"), family = "binomial",
data=gamData, rho = myRhoVal,
AR.start = gamData$start.event, nThreads = 3, discrete = T)
save(list= c("interaction.compet_kitFleece.gamm0",
"interaction.compet_kitFleece.gamm1"),
file = fullModelFile)
}Model criticism: How are we doing regarding auto-correlation?
myAutoCorr =acf_resid(interaction.compet_kitFleece.gamm1)
Is the default for the smooth (k = 10) sufficient?
gam.check(interaction.compet_kitFleece.gamm1)
##
## Method: fREML Optimizer: perf chol
## $grad
## [1] 3.488698e-10 -5.758505e-12 -2.080648e-03 -1.490578e-03 5.505285e-11
## [6] 1.487330e-08 -3.012283e-07 -1.595239e-09 5.075282e-10 -9.554002e-11
## [11] 8.704504e-11 5.817733e-10 4.767715e-10 6.443379e-11 -3.904450e-09
## [16] -3.568913e-04
##
## $hess
## [,1] [,2] [,3] [,4] [,5]
## [1,] 3.488878e+00 -7.118330e-05 1.883164e-04 2.086634e-04 2.167088e-02
## [2,] -7.118330e-05 3.575813e+00 1.686063e-04 -2.082987e-07 -2.111883e-07
## [3,] 1.883164e-04 1.686063e-04 2.075738e-03 -1.576270e-06 -4.557843e-07
## [4,] 2.086634e-04 -2.082987e-07 -1.576270e-06 1.488017e-03 -1.190370e-06
## [5,] 2.167088e-02 -2.111883e-07 -4.557843e-07 -1.190370e-06 4.348542e+00
## [6,] -2.969726e-02 8.170725e-08 -1.072585e-06 -1.411217e-06 -3.540216e-01
## [7,] -1.232164e-02 2.251190e-07 -6.386259e-07 -7.342654e-07 2.626423e-02
## [8,] -3.876821e-07 4.655849e-02 3.288944e-07 1.321195e-09 -2.374359e-08
## [9,] 1.586598e-06 -1.017152e-02 -4.960540e-06 1.224518e-08 -1.777662e-08
## [10,] 7.072524e-08 -3.039729e-03 -2.175744e-07 4.662594e-10 -9.964996e-10
## [11,] 9.883974e-02 -2.100073e-06 4.650470e-06 4.748376e-06 -2.874770e-02
## [12,] -9.672676e-03 4.173984e-08 -8.184191e-07 -1.253865e-06 -1.359073e-02
## [13,] -1.559061e-02 3.535761e-07 -9.443235e-07 -1.106510e-06 -2.066704e-02
## [14,] -7.443209e-07 4.232443e-02 1.555389e-06 -2.105582e-09 -6.523833e-09
## [15,] 1.731081e-07 -3.814853e-03 -5.098902e-07 1.004885e-09 -1.211248e-09
## [16,] -2.161490e-09 5.364334e-05 4.601878e-09 -1.030844e-11 -4.210237e-11
## [,6] [,7] [,8] [,9] [,10]
## [1,] -2.969726e-02 -1.232164e-02 -3.876821e-07 1.586598e-06 7.072524e-08
## [2,] 8.170725e-08 2.251190e-07 4.655849e-02 -1.017152e-02 -3.039729e-03
## [3,] -1.072585e-06 -6.386259e-07 3.288944e-07 -4.960540e-06 -2.175744e-07
## [4,] -1.411217e-06 -7.342654e-07 1.321195e-09 1.224518e-08 4.662594e-10
## [5,] -3.540216e-01 2.626423e-02 -2.374359e-08 -1.777662e-08 -9.964995e-10
## [6,] 7.666502e+00 -4.441255e-01 -7.843717e-09 -2.837885e-08 -9.067094e-10
## [7,] -4.441255e-01 2.361734e+00 1.572049e-09 -5.117327e-09 -2.365266e-10
## [8,] -7.843717e-09 1.572049e-09 4.904163e+00 -5.875864e-01 1.540887e-02
## [9,] -2.837885e-08 -5.117327e-09 -5.875864e-01 1.047881e+01 -4.271804e-01
## [10,] -9.067095e-10 -2.365266e-10 1.540887e-02 -4.271804e-01 5.181825e-01
## [11,] -1.716472e-02 -1.542121e-02 -3.435465e-08 3.982536e-08 1.254875e-09
## [12,] -2.324477e-02 -7.385966e-03 -6.764359e-09 -2.193674e-08 -7.124700e-10
## [13,] -1.570228e-02 -1.097573e-02 2.274705e-09 -6.669225e-09 -3.537795e-10
## [14,] 2.169655e-10 2.590221e-09 -1.140304e-02 -1.622682e-02 -6.996366e-04
## [15,] -1.961956e-09 -5.337287e-10 -6.512315e-03 -1.355881e-02 6.580466e-05
## [16,] 1.468009e-11 9.562038e-12 -6.481737e-05 -4.328781e-05 -2.642882e-05
## [,11] [,12] [,13] [,14] [,15]
## [1,] 9.883974e-02 -9.672676e-03 -1.559061e-02 -7.443209e-07 1.731081e-07
## [2,] -2.100073e-06 4.173984e-08 3.535761e-07 4.232443e-02 -3.814853e-03
## [3,] 4.650470e-06 -8.184191e-07 -9.443235e-07 1.555389e-06 -5.098902e-07
## [4,] 4.748376e-06 -1.253865e-06 -1.106510e-06 -2.105582e-09 1.004885e-09
## [5,] -2.874770e-02 -1.359073e-02 -2.066704e-02 -6.523833e-09 -1.211248e-09
## [6,] -1.716472e-02 -2.324477e-02 -1.570228e-02 2.169655e-10 -1.961956e-09
## [7,] -1.542121e-02 -7.385966e-03 -1.097573e-02 2.590221e-09 -5.337287e-10
## [8,] -3.435465e-08 -6.764359e-09 2.274705e-09 -1.140304e-02 -6.512315e-03
## [9,] 3.982536e-08 -2.193674e-08 -6.669225e-09 -1.622682e-02 -1.355881e-02
## [10,] 1.254875e-09 -7.124701e-10 -3.537795e-10 -6.996366e-04 6.580466e-05
## [11,] 2.962036e+00 7.366857e-02 9.789900e-02 -3.024778e-08 4.672609e-09
## [12,] 7.366857e-02 2.947388e+00 -2.420712e-01 -4.748003e-10 -1.515864e-09
## [13,] 9.789900e-02 -2.420712e-01 3.108273e+00 4.003006e-09 -7.624866e-10
## [14,] -3.024778e-08 -4.748003e-10 4.003006e-09 1.034150e+00 -3.205520e-02
## [15,] 4.672609e-09 -1.515864e-09 -7.624866e-10 -3.205520e-02 9.420318e-01
## [16,] -1.547068e-10 5.904136e-12 1.255166e-11 6.604327e-04 -2.008086e-03
## [,16]
## [1,] -2.161490e-09
## [2,] 5.364334e-05
## [3,] 4.601878e-09
## [4,] -1.030844e-11
## [5,] -4.210237e-11
## [6,] 1.468009e-11
## [7,] 9.562038e-12
## [8,] -6.481737e-05
## [9,] -4.328781e-05
## [10,] -2.642882e-05
## [11,] -1.547068e-10
## [12,] 5.904135e-12
## [13,] 1.255166e-11
## [14,] 6.604327e-04
## [15,] -2.008086e-03
## [16,] 1.875207e-03
##
## Model rank = 1346 / 1346
##
## 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):vowelContrastFLEECE 11.0 7.9 0.96 0.28
## s(Time):vowelContrastKIT 11.0 7.8 0.96 0.29
## s(Time):isAustralian 12.0 2.0 0.96 0.38
## s(Time):isAustralianFleece 10.0 2.0 0.96 0.35
## s(Time,pp):vowelContrastFLEECE 450.0 76.0 0.96 0.30
## s(Time,pp):vowelContrastKIT 450.0 70.5 0.96 0.33
## s(Time,item):vowelContrastFLEECE 200.0 30.9 0.96 0.34
## s(Time,item):vowelContrastKIT 200.0 13.8 0.96 0.33We now examine this model to determine whether there is an overall effect of contrast.
summary(interaction.compet_kitFleece.gamm1)##
## Family: binomial
## Link function: logit
##
## Formula:
## position ~ vowelContrast + s(Time, by = vowelContrast, k = 12) +
## s(Time, by = isAustralian, k = 12) + s(Time, by = isAustralianFleece) +
## 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.44524 0.09834 -14.70 <2e-16 ***
## vowelContrastKIT -0.13714 0.11237 -1.22 0.222
## ---
## 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):vowelContrastFLEECE 7.901 9.481 119.945 <2e-16 ***
## s(Time):vowelContrastKIT 7.799 9.435 163.014 <2e-16 ***
## s(Time):isAustralian 2.005 2.010 1.451 0.488
## s(Time):isAustralianFleece 2.004 2.007 0.930 0.582
## s(Time,pp):vowelContrastFLEECE 75.973 448.000 346.049 <2e-16 ***
## s(Time,pp):vowelContrastKIT 70.463 449.000 331.189 <2e-16 ***
## s(Time,item):vowelContrastFLEECE 30.876 98.000 141.688 <2e-16 ***
## s(Time,item):vowelContrastKIT 13.791 98.000 84.218 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0875 Deviance explained = 9.41%
## fREML = 85505 Scale est. = 1 n = 320222The model shows significant nonlinear time effects and substantial participant and item variability, but no significant main effect of vowel contrast and no significant accent-related smooths, indicating that fixation trajectories change over time but do not reliably differ between KIT and FLEECE or across accents.
smooth <- plot_smooth(interaction.compet_kitFleece.gamm1, view = "Time",
plot_all = "vowelContrast",
rug = FALSE, rm.ranef = TRUE, lwd = 4)## Summary:
## * vowelContrast : factor; set to the value(s): FLEECE, KIT.
## * Time : numeric predictor; with 30 values ranging from 10.000000 to 890.000000.
## * isAustralian : numeric predictor; set to the value(s): 0.
## * isAustralianFleece : 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): beetle. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(Time,pp):vowelContrastFLEECE,s(Time,pp):vowelContrastKIT,s(Time,item):vowelContrastFLEECE,s(Time,item):vowelContrastKIT
## 
We now proceed to model comparison by successively removing effects to evaluate their contribution to model fit.
reducedModelFile = "compet_kitFleece_reducedModel.Rdata"
if (file.exists(reducedModelFile))
{
load(reducedModelFile)
}else{
mainEffect.compet_kitFleece.gamm <- bam( position ~ vowelContrast+
s(Time, by = vowelContrast) +
#s(Time, by = isAustralian) +
#s(Time, by = isAustralianFleece) +
s(Time, pp, by = vowelContrast, bs="fs") +
s(Time, item, by = vowelContrast, bs="fs"), family = "binomial",
data=gamData, rho = myRhoVal,
AR.start = gamData$start.event, nThreads = 3, discrete = T)
save(list= c("mainEffect.compet_kitFleece.gamm"), file = reducedModelFile)
}summary(mainEffect.compet_kitFleece.gamm)##
## Family: binomial
## Link function: logit
##
## Formula:
## position ~ vowelContrast + s(Time, by = vowelContrast) + 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.43386 0.09112 -15.735 <2e-16 ***
## vowelContrastKIT -0.17853 0.09986 -1.788 0.0738 .
## ---
## 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):vowelContrastFLEECE 7.172 8.234 125.0 <2e-16 ***
## s(Time):vowelContrastKIT 7.107 8.241 165.5 <2e-16 ***
## s(Time,pp):vowelContrastFLEECE 75.625 448.000 578.8 <2e-16 ***
## s(Time,pp):vowelContrastKIT 70.291 448.000 529.8 <2e-16 ***
## s(Time,item):vowelContrastFLEECE 30.732 99.000 163.3 <2e-16 ***
## s(Time,item):vowelContrastKIT 13.667 98.000 158.2 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0899 Deviance explained = 9.61%
## fREML = 85259 Scale est. = 1 n = 320222This model shows strong nonlinear time effects and substantial participant and item variability, with only a marginal main effect of vowel contrast (p = .074), suggesting weak evidence that KIT and FLEECE differ overall, but no indication of additional effects beyond contrast-specific time trajectories.
smooth <- plot_smooth(mainEffect.compet_kitFleece.gamm, view = "Time",
plot_all = "vowelContrast",
rug = FALSE, rm.ranef = TRUE, lwd = 4)## Summary:
## * vowelContrast : factor; set to the value(s): FLEECE, KIT.
## * Time : numeric predictor; with 30 values ranging from 10.000000 to 890.000000.
## * pp : factor; set to the value(s): CAT01. (Might be canceled as random effect, check below.)
## * item : factor; set to the value(s): beetle. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(Time,pp):vowelContrastFLEECE,s(Time,pp):vowelContrastKIT,s(Time,item):vowelContrastFLEECE,s(Time,item):vowelContrastKIT
## 
3.2.3 DRESS-TRAP contrast
3.2.3.1. Target Fixations
A. TIME-WINDOW ANALYSIS
rm(list = ls())
library(lmerTest)
library(stringr)
library(tidyverse)
library(mgcv)
library(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("DRESS|TRAP", 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
## 1797 1799trialDataE$contrastC = -0.5 + as.numeric(trialDataE$contrast == "TRAP")
table(trialDataE$contrastC)##
## -0.5 0.5
## 1798 1798useData = cbind(trialDataE,gazeDataE)
We first need to reorganize the data.
windowData = useData[ ,c(1,2,4:8,10:134)] %>%
pivot_longer(., cols = starts_with("T_"),
names_to = "Time",
values_to = "position") %>% as.data.frame
windowData$position = as.numeric(windowData$position)
windowData$Time = gsub("T_","", windowData$Time) %>% as.numeric()
windowData = windowData[ order(windowData$pp,windowData$trial,windowData$Time), ]
windowData$pp = as.factor(windowData$pp)
windowData$item = as.factor(windowData$item)
windowData$isAustralian = ifelse(windowData$accent == "AusE",1,0)
windowData$isAustralianTrap = windowData$isAustralian *
as.numeric(windowData$contrast == "TRAP")
windowData$vowelContrast = windowData$contrast %>% as.factor()
#restricting the data to 0-900 after vowel onset
windowData = windowData %>% filter(Time > 0 & Time < 900)
windowData <- start_event(windowData, event=c("pp","trial"))
windowDataPlot = windowData %>% group_by(vowelContrast, Time, accent) %>%
summarize(meanFix = mean(position))
Fig.target_dressTrap <- ggplot(windowDataPlot, aes(x = Time, y = meanFix, col = vowelContrast)) + geom_line() + facet_wrap(~accent)
print(Fig.target_dressTrap)
In both accents (AusE and SEng), fixation proportions increase steadily over time toward the target, reaching high levels by around 700–800 ms, with only minimal differences between accents.
aggregatedData = windowData %>% filter(Time > 200 & Time < 600) %>%
group_by(pp, item, contrast) %>% summarize(meanLooks = mean(position)) %>%
ungroup()
meanLooksValues = sort(unique(aggregatedData$meanLooks))
correction = meanLooksValues[2]/2
aggregatedData$logitLooks = ifelse( aggregatedData$meanLooks == 0,
qlogis(correction),
ifelse(aggregatedData$meanLooks == 1,
qlogis(1 - correction),
qlogis(aggregatedData$meanLooks)
)
)
plot(aggregatedData$meanLooks, aggregatedData$logitLooks)
library(lmerTest)
target.dressTrap.lmer = lmer(meanLooks ~ contrast +
(1+contrast|pp)+
(1|item), data = aggregatedData)summary(target.dressTrap.lmer)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanLooks ~ contrast + (1 + contrast | pp) + (1 | item)
## Data: aggregatedData
##
## REML criterion at convergence: -450.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.58405 -0.61820 0.01342 0.70206 2.59147
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## pp (Intercept) 0.0101361 0.100678
## contrastTRAP 0.0000113 0.003362 1.00
## item (Intercept) 0.0022677 0.047620
## Residual 0.0306898 0.175185
## Number of obs: 900, groups: pp, 45; item, 20
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.49136 0.02281 40.08042 21.543 <2e-16 ***
## contrastTRAP 0.03712 0.02429 18.01392 1.528 0.144
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## contrstTRAP -0.519
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')This model shows no significant effect of contrast (p = .144), indicating that mean target looks do not reliably differ between DRESS and TRAP.
B. Generalized Additive Mixed Model
gamData <- windowDatafullModelFile = "target_dressTrap_fullModels.Rdata"
if (file.exists(fullModelFile))
{
load(fullModelFile)
}else{
interaction.target_dressTrap.gamm0 <- bam( position ~ vowelContrast+
s(Time, by = vowelContrast) +
s(Time, by = isAustralian) +
s(Time, by = isAustralianTrap) +
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.target_dressTrap.gamm0)
myRhoVal = myAutoCorr[2]
interaction.target_dressTrap.gamm1 <- bam( position ~ vowelContrast+
s(Time, by = vowelContrast, k=12) +
s(Time, by = isAustralian, k= 12) +
s(Time, by = isAustralianTrap) +
s(Time, pp, by = vowelContrast, bs="fs") +
s(Time, item, by = vowelContrast, bs="fs"), family = "binomial",
data=gamData, rho = myRhoVal,
AR.start = gamData$start.event, nThreads = 3, discrete = T)
save(list= c("interaction.target_dressTrap.gamm0","interaction.target_dressTrap.gamm1"), file = fullModelFile)
}Model criticism: How are we doing regarding auto-correlation?
myAutoCorr =acf_resid(interaction.target_dressTrap.gamm1)
Is the default for the smooth (k = 10) sufficient?
gam.check(interaction.target_dressTrap.gamm1)
##
## Method: fREML Optimizer: perf chol
## $grad
## [1] 1.124132e-10 8.562928e-12 2.434613e-10 -2.620573e-09 -2.240462e-09
## [6] 9.517720e-11 7.313439e-10 9.979573e-11 -2.366566e-09 6.884875e-10
## [11] -3.523208e-10 2.073630e-11 -1.455225e-03 -2.800586e-10 -4.784861e-10
## [16] 1.989444e-10
##
## $hess
## [,1] [,2] [,3] [,4] [,5]
## [1,] 2.574389e+00 -5.372136e-02 2.671644e-01 -9.291128e-03 1.885889e-01
## [2,] -5.372136e-02 2.441941e+00 1.102698e-01 9.249520e-03 -4.604970e-03
## [3,] 2.671644e-01 1.102698e-01 9.345165e-01 -7.267178e-02 1.018538e-02
## [4,] -9.291128e-03 9.249520e-03 -7.267178e-02 2.803835e-01 -3.588512e-03
## [5,] 1.885889e-01 -4.604970e-03 1.018538e-02 -3.588512e-03 1.271000e+01
## [6,] 5.117274e-02 3.542331e-04 -1.391567e-05 5.980944e-04 1.622817e-01
## [7,] -6.315414e-02 9.081398e-04 -1.149494e-02 -3.302645e-03 -4.124627e-01
## [8,] -3.812244e-03 1.883458e-01 8.225007e-03 2.363819e-03 -3.871683e-04
## [9,] -9.275696e-05 9.007293e-03 -5.712868e-04 -2.505065e-03 4.953290e-05
## [10,] 1.585711e-03 -4.470818e-02 -5.825819e-03 -7.007061e-04 1.268153e-05
## [11,] 6.780741e-02 -1.582060e-03 4.075920e-03 -8.099319e-04 -5.025720e-02
## [12,] 3.728936e-03 2.057554e-04 -6.253283e-04 1.334829e-04 -4.135785e-02
## [13,] -4.580549e-06 2.165997e-07 2.317746e-07 1.419592e-07 -6.433368e-06
## [14,] -1.664102e-03 5.877961e-02 3.591199e-03 1.195273e-03 -1.588594e-04
## [15,] -5.970585e-05 2.553971e-03 4.344918e-05 -4.339907e-04 7.463628e-06
## [16,] 3.553481e-04 -1.120619e-02 -1.134819e-03 7.170709e-05 -7.017299e-06
## [,6] [,7] [,8] [,9] [,10]
## [1,] 5.117274e-02 -6.315414e-02 -3.812244e-03 -9.275696e-05 1.585711e-03
## [2,] 3.542331e-04 9.081398e-04 1.883458e-01 9.007293e-03 -4.470818e-02
## [3,] -1.391567e-05 -1.149494e-02 8.225007e-03 -5.712868e-04 -5.825819e-03
## [4,] 5.980944e-04 -3.302645e-03 2.363819e-03 -2.505065e-03 -7.007061e-04
## [5,] 1.622817e-01 -4.124627e-01 -3.871683e-04 4.953290e-05 1.268153e-05
## [6,] 8.096357e+00 -7.784636e-01 3.263589e-05 -1.236014e-05 1.147544e-05
## [7,] -7.784636e-01 1.308235e+01 2.719411e-05 3.590109e-05 -7.706644e-05
## [8,] 3.263589e-05 2.719411e-05 1.769747e+01 1.001442e-01 -5.841955e-01
## [9,] -1.236014e-05 3.590109e-05 1.001442e-01 6.024813e+00 -6.316207e-01
## [10,] 1.147544e-05 -7.706644e-05 -5.841955e-01 -6.316207e-01 1.193955e+01
## [11,] -2.834672e-02 -4.203819e-02 -9.287350e-05 2.253274e-06 2.848508e-05
## [12,] -2.767136e-02 -2.278868e-02 1.758904e-05 -2.543764e-06 -3.229051e-07
## [13,] -5.125626e-06 -2.333097e-06 1.679765e-08 1.115155e-09 -3.037668e-09
## [14,] 1.305264e-05 1.173333e-05 -7.453526e-02 -2.417521e-02 -5.676286e-02
## [15,] -2.784932e-06 8.171771e-06 -4.288636e-02 -2.131897e-02 -1.581359e-02
## [16,] 4.256593e-06 -2.222244e-05 -6.513843e-02 -1.542573e-02 -6.140306e-02
## [,11] [,12] [,13] [,14] [,15]
## [1,] 6.780741e-02 3.728936e-03 -4.580549e-06 -1.664102e-03 -5.970585e-05
## [2,] -1.582060e-03 2.057554e-04 2.165997e-07 5.877961e-02 2.553971e-03
## [3,] 4.075920e-03 -6.253283e-04 2.317746e-07 3.591199e-03 4.344918e-05
## [4,] -8.099319e-04 1.334829e-04 1.419592e-07 1.195273e-03 -4.339907e-04
## [5,] -5.025720e-02 -4.135785e-02 -6.433368e-06 -1.588594e-04 7.463628e-06
## [6,] -2.834672e-02 -2.767136e-02 -5.125626e-06 1.305264e-05 -2.784932e-06
## [7,] -4.203819e-02 -2.278868e-02 -2.333097e-06 1.173333e-05 8.171771e-06
## [8,] -9.287350e-05 1.758904e-05 1.679765e-08 -7.453526e-02 -4.288636e-02
## [9,] 2.253274e-06 -2.543764e-06 1.115155e-09 -2.417521e-02 -2.131897e-02
## [10,] 2.848508e-05 -3.229051e-07 -3.037668e-09 -5.676286e-02 -1.581359e-02
## [11,] 2.707940e+00 -8.205933e-02 -6.407937e-05 -4.174864e-05 9.380925e-08
## [12,] -8.205933e-02 2.014569e+00 -6.392426e-05 7.003773e-06 -3.864735e-07
## [13,] -6.407937e-05 -6.392426e-05 1.455041e-03 6.825663e-09 4.827229e-10
## [14,] -4.174864e-05 7.003773e-06 6.825663e-09 3.343737e+00 -6.860880e-02
## [15,] 9.380925e-08 -3.864735e-07 4.827229e-10 -6.860880e-02 1.507825e+00
## [16,] 6.428046e-06 3.674116e-07 -6.235596e-10 -1.157565e-01 -3.277858e-02
## [,16]
## [1,] 3.553481e-04
## [2,] -1.120619e-02
## [3,] -1.134819e-03
## [4,] 7.170709e-05
## [5,] -7.017299e-06
## [6,] 4.256593e-06
## [7,] -2.222244e-05
## [8,] -6.513843e-02
## [9,] -1.542573e-02
## [10,] -6.140306e-02
## [11,] 6.428046e-06
## [12,] 3.674116e-07
## [13,] -6.235596e-10
## [14,] -1.157565e-01
## [15,] -3.277858e-02
## [16,] 3.151849e+00
##
## Model rank = 1346 / 1346
##
## 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):vowelContrastDRESS 11.00 7.36 0.98 0.100 .
## s(Time):vowelContrastTRAP 11.00 8.16 0.98 0.120
## s(Time):isAustralian 12.00 5.68 0.98 0.115
## s(Time):isAustralianTrap 10.00 3.28 0.98 0.090 .
## s(Time,pp):vowelContrastDRESS 450.00 116.62 0.98 0.115
## s(Time,pp):vowelContrastTRAP 450.00 124.31 0.98 0.080 .
## s(Time,item):vowelContrastDRESS 200.00 17.26 0.98 0.085 .
## s(Time,item):vowelContrastTRAP 200.00 26.29 0.98 0.125
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We now examine this model to determine whether there is an overall effect of contrast.
summary(interaction.target_dressTrap.gamm1)##
## Family: binomial
## Link function: logit
##
## Formula:
## position ~ vowelContrast + s(Time, by = vowelContrast, k = 12) +
## s(Time, by = isAustralian, k = 12) + s(Time, by = isAustralianTrap) +
## 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.16876 0.07426 2.273 0.0231 *
## vowelContrastTRAP 0.12578 0.12375 1.016 0.3094
## ---
## 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):vowelContrastDRESS 7.357 8.984 342.773 <2e-16 ***
## s(Time):vowelContrastTRAP 8.164 9.634 417.975 <2e-16 ***
## s(Time):isAustralian 5.676 7.068 13.509 0.0728 .
## s(Time):isAustralianTrap 3.277 3.946 3.601 0.4964
## s(Time,pp):vowelContrastDRESS 116.619 448.000 670.491 <2e-16 ***
## s(Time,pp):vowelContrastTRAP 124.308 448.000 629.734 <2e-16 ***
## s(Time,item):vowelContrastDRESS 17.260 98.000 344.925 <2e-16 ***
## s(Time,item):vowelContrastTRAP 26.295 98.000 180.866 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.271 Deviance explained = 21.8%
## fREML = 1.0623e+05 Scale est. = 1 n = 320044The model shows strong nonlinear time effects and substantial participant and item variability, but no significant main effect of vowel contrast* and no reliable accent-related effects** (aside from a marginal overall time-by-accent smooth), indicating that while fixation trajectories change robustly over time, they do not significantly differ between DRESS and TRAP or across accents.
smooth <- plot_smooth(interaction.target_dressTrap.gamm1, view = "Time",
plot_all = "vowelContrast",
rug = FALSE, rm.ranef = TRUE, lwd = 4)## Summary:
## * vowelContrast : factor; set to the value(s): DRESS, TRAP.
## * Time : numeric predictor; with 30 values ranging from 10.000000 to 890.000000.
## * isAustralian : numeric predictor; set to the value(s): 0.
## * isAustralianTrap : 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):vowelContrastDRESS,s(Time,pp):vowelContrastTRAP,s(Time,item):vowelContrastDRESS,s(Time,item):vowelContrastTRAP
## 
We now proceed to model comparison by successively removing effects to evaluate their contribution to model fit.
reducedModelFile = "target_dressTrap_reducedModel.Rdata"
if (file.exists(reducedModelFile))
{
load(reducedModelFile)
}else{
mainEffect.target_dressTrap.gamm <- bam( position ~ vowelContrast+
s(Time, by = vowelContrast) +
#s(Time, by = isAustralian) +
#s(Time, by = isAustralianTrap) +
s(Time, pp, by = vowelContrast, bs="fs") +
s(Time, item, by = vowelContrast, bs="fs"), family = "binomial",
data=gamData, rho = myRhoVal,
AR.start = gamData$start.event, nThreads = 3, discrete = T)
save(list= c("mainEffect.target_dressTrap.gamm"), file = reducedModelFile)
}summary(mainEffect.target_dressTrap.gamm)##
## Family: binomial
## Link function: logit
##
## Formula:
## position ~ vowelContrast + s(Time, by = vowelContrast) + 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.20684 0.06785 3.048 0.0023 **
## vowelContrastTRAP 0.11507 0.11593 0.992 0.3210
## ---
## 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):vowelContrastDRESS 7.284 8.341 444.6 <2e-16 ***
## s(Time):vowelContrastTRAP 7.264 8.288 556.0 <2e-16 ***
## s(Time,pp):vowelContrastDRESS 116.393 448.000 671.3 <2e-16 ***
## s(Time,pp):vowelContrastTRAP 124.500 448.000 630.6 <2e-16 ***
## s(Time,item):vowelContrastDRESS 17.155 99.000 345.8 <2e-16 ***
## s(Time,item):vowelContrastTRAP 26.364 99.000 183.3 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.271 Deviance explained = 21.8%
## fREML = 1.0633e+05 Scale est. = 1 n = 320044This model shows strong nonlinear time effects and substantial participant and item variability, but no significant main effect of vowel contrast, indicating that although fixation proportions change robustly over time, they do not reliably differ between DRESS and TRAP overall.
smooth <- plot_smooth(mainEffect.target_dressTrap.gamm, view = "Time",
plot_all = "vowelContrast",
rug = FALSE, rm.ranef = TRUE, lwd = 4)## Summary:
## * vowelContrast : factor; set to the value(s): DRESS, TRAP.
## * Time : numeric predictor; with 30 values ranging from 10.000000 to 890.000000.
## * 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):vowelContrastDRESS,s(Time,pp):vowelContrastTRAP,s(Time,item):vowelContrastDRESS,s(Time,item):vowelContrastTRAP
## 
3.2.3.2. Competitor Fixations
A. TIME-WINDOW ANALYSIS
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("DRESS|TRAP", 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
## 1797 1799trialDataE$contrastC = -0.5 + as.numeric(trialDataE$contrast == "TRAP")
table(trialDataE$contrastC)##
## -0.5 0.5
## 1798 1798useData = cbind(trialDataE,gazeDataE)We first need to reorganize the data.
windowData = useData[ ,c(1,2,4:8,10:134)] %>%
pivot_longer(., cols = starts_with("T_"),
names_to = "Time",
values_to = "position") %>% as.data.frame
windowData$position = as.numeric(windowData$position)
windowData$Time = gsub("T_","", windowData$Time) %>% as.numeric()
windowData = windowData[ order(windowData$pp,windowData$trial,windowData$Time), ]
windowData$pp = as.factor(windowData$pp)
windowData$item = as.factor(windowData$item)
windowData$isAustralian = ifelse(windowData$accent == "AusE",1,0)
windowData$isAustralianTrap = windowData$isAustralian *
as.numeric(windowData$contrast == "TRAP")
windowData$vowelContrast = windowData$contrast %>% as.factor()
#restricting the data to 0-900 after vowel onset
windowData = windowData %>% filter(Time > 0 & Time < 900)
windowData <- start_event(windowData, event=c("pp","trial"))
windowDataPlot = windowData %>% group_by(vowelContrast, Time, accent) %>%
summarize(meanFix = mean(position))
Fig.compet_dressTrap <- ggplot(windowDataPlot, aes(x = Time, y = meanFix, col = vowelContrast)) + geom_line() + facet_wrap(~accent)
print(Fig.compet_dressTrap)
In both accents (AusE and SEng), competitor fixations rise to a peak around 250–300 ms and then decline over time, with DRESS showing slightly higher fixation proportions than TRAP (particularly in SEng), while overall accent differences appear small.
lmerData = windowData %>% filter(Time > 200 & Time < 600) %>%
group_by(pp, item, trial, isAustralian, contrastC) %>% summarize(meanLooks = mean(position))
lmerData$accentC = lmerData$isAustralian - 0.5
outcomes = sort(unique(lmerData$meanLooks))
if (outcomes[1] == 0){
correction = outcomes[2] * 0.5
}
lmerData$logitLooks = ifelse(lmerData$meanLooks == 0, qlogis(correction),
ifelse(lmerData$meanLooks == 1, qlogis(1- correction),
qlogis(lmerData$meanLooks))
)
plot(lmerData$meanLooks, lmerData$logitLooks)
library(lmerTest)
lmerData$trialScaled = scale(lmerData$trial)
compet_dressTrap.lmer = lmer(logitLooks ~ contrastC*accentC*trialScaled +
(1+contrastC+contrastC||pp) +
(1+accentC||item), data = lmerData,
control=lmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))
summary(compet_dressTrap.lmer)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logitLooks ~ contrastC * accentC * trialScaled + (1 + contrastC +
## contrastC || pp) + (1 + accentC || item)
## Data: lmerData
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## REML criterion at convergence: 16048.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.1663 -0.8789 -0.6926 0.9179 3.1333
##
## Random effects:
## Groups Name Variance Std.Dev.
## pp (Intercept) 0.061053 0.24709
## pp.1 contrastC 0.002189 0.04678
## item (Intercept) 0.011749 0.10839
## item.1 accentC 0.011225 0.10595
## Residual 4.998325 2.23569
## Number of obs: 3596, groups: pp, 45; item, 20
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.34776 0.05776 28.51522 -40.650 <2e-16
## contrastC -0.20083 0.08924 15.04838 -2.250 0.0398
## accentC -0.09765 0.07828 18.00636 -1.247 0.2282
## trialScaled -0.00319 0.03747 3564.41139 -0.085 0.9321
## contrastC:accentC 0.29730 0.15655 18.00601 1.899 0.0737
## contrastC:trialScaled 0.02034 0.07495 3564.77276 0.271 0.7861
## accentC:trialScaled 0.05534 0.07495 3569.75753 0.738 0.4603
## contrastC:accentC:trialScaled 0.15509 0.15003 3573.31192 1.034 0.3013
##
## (Intercept) ***
## contrastC *
## accentC
## trialScaled
## contrastC:accentC .
## contrastC:trialScaled
## accentC:trialScaled
## contrastC:accentC:trialScaled
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cntrsC accntC trlScl cntC:C cntC:S accC:S
## contrastC 0.000
## accentC 0.001 0.000
## trialScaled 0.001 -0.015 -0.025
## cntrstC:ccC 0.000 0.001 0.000 -0.007
## cntrstC:trS -0.012 0.001 -0.008 -0.028 -0.025
## accntC:trlS -0.017 -0.007 0.001 -0.005 -0.017 -0.026
## cntrstC:C:S -0.005 -0.022 -0.017 -0.027 0.001 -0.003 -0.029This model shows a significant main effect of contrast (p = .040), indicating that competitor looks differ reliably between DRESS and TRAP, but no significant effects of accent, trial, or higher-order interactions, suggesting that the contrast effect does not vary across accents or over trials (with only a marginal contrast-by-accent interaction, p = .074).
B. Generalized Additive Mixed Model
gamData <- windowDatafullModelFile = "compet_dressTrap_fullModels.Rdata"
if (file.exists(fullModelFile))
{
load(fullModelFile)
}else{
interaction.compet_dressTrap.gamm0 <- bam( position ~ vowelContrast+
s(Time, by = vowelContrast) +
s(Time, by = isAustralian) +
s(Time, by = isAustralianTrap) +
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.compet_dressTrap.gamm0)
myRhoVal = myAutoCorr[2]
interaction.compet_dressTrap.gamm1 <- bam( position ~ vowelContrast+
s(Time, by = vowelContrast, k=12) +
s(Time, by = isAustralian, k= 12) +
s(Time, by = isAustralianTrap) +
s(Time, pp, by = vowelContrast, bs="fs") +
s(Time, item, by = vowelContrast, bs="fs"), family = "binomial",
data=gamData, rho = myRhoVal,
AR.start = gamData$start.event, nThreads = 3, discrete = T)
save(list= c("interaction.compet_dressTrap.gamm0",
"interaction.compet_dressTrap.gamm1"),
file = fullModelFile)
}Model criticism: How are we doing regarding auto-correlation?
myAutoCorr =acf_resid(interaction.compet_dressTrap.gamm1)
Is the default for the smooth (k = 10) sufficient?
gam.check(interaction.compet_dressTrap.gamm1)
##
## Method: fREML Optimizer: perf chol
## $grad
## [1] -1.178316e-05 -6.967773e-06 -7.489831e-04 -9.684682e-06 -3.303210e-05
## [6] -3.812109e-06 -5.411167e-04 3.319510e-05 5.845322e-06 -2.178664e-05
## [11] -1.613847e-06 -2.299812e-06 -1.823013e-06 -7.301697e-04 4.639945e-06
## [16] 1.196212e-06
##
## $hess
## [,1] [,2] [,3] [,4] [,5]
## [1,] 3.133878e+00 -2.809268e-05 5.098853e-05 -9.433452e-05 8.272094e-02
## [2,] -2.809268e-05 1.888164e+00 7.637095e-05 2.485599e-01 -8.285359e-07
## [3,] 5.098853e-05 7.637095e-05 7.487630e-04 -7.925943e-05 -1.565917e-06
## [4,] -9.433452e-05 2.485599e-01 -7.925943e-05 6.384780e-01 -2.688339e-06
## [5,] 8.272094e-02 -8.285359e-07 -1.565917e-06 -2.688339e-06 7.497056e+00
## [6,] -3.875749e-02 2.940057e-07 1.119619e-05 2.379910e-06 -9.649535e-01
## [7,] -5.072738e-06 4.830259e-11 3.767417e-09 6.183713e-10 1.800560e-05
## [8,] -1.573085e-06 7.978232e-02 5.603121e-06 3.459462e-02 -7.652765e-08
## [9,] 6.017420e-08 -3.021403e-02 2.434089e-06 1.931896e-03 5.325669e-09
## [10,] -1.725266e-10 3.810211e-05 -2.691839e-09 4.540213e-06 -3.677829e-11
## [11,] 2.861680e-02 -3.139979e-07 3.895279e-07 -1.132166e-06 -6.470501e-03
## [12,] -2.382350e-03 3.152343e-08 4.398105e-06 4.486090e-07 -1.505899e-02
## [13,] -2.709489e-04 2.918004e-09 3.040644e-07 3.745140e-08 8.073588e-04
## [14,] -1.935033e-09 1.511147e-04 7.514155e-09 2.418647e-05 -6.347924e-11
## [15,] -4.227225e-08 -1.536106e-02 1.291319e-06 8.372008e-04 1.079039e-09
## [16,] 9.620642e-08 -1.012532e-02 3.537853e-07 -1.987040e-03 4.008188e-09
## [,6] [,7] [,8] [,9] [,10]
## [1,] -3.875749e-02 -5.072738e-06 -1.573085e-06 6.017420e-08 -1.725266e-10
## [2,] 2.940057e-07 4.830259e-11 7.978232e-02 -3.021403e-02 3.810211e-05
## [3,] 1.119619e-05 3.767417e-09 5.603121e-06 2.434089e-06 -2.691839e-09
## [4,] 2.379910e-06 6.183713e-10 3.459462e-02 1.931896e-03 4.540213e-06
## [5,] -9.649535e-01 1.800560e-05 -7.652765e-08 5.325669e-09 -3.677829e-11
## [6,] 8.988545e+00 -1.924923e-03 1.623265e-08 -1.574602e-08 2.140269e-11
## [7,] -1.924923e-03 5.542436e-04 3.234351e-12 -4.524284e-12 5.256780e-15
## [8,] 1.623265e-08 3.234351e-12 8.143346e+00 -2.290968e-01 9.276075e-04
## [9,] -1.574602e-08 -4.524284e-12 -2.290968e-01 7.048228e+00 -8.555279e-03
## [10,] 2.140269e-11 5.256780e-15 9.276075e-04 -8.555279e-03 4.446642e-04
## [11,] -5.726829e-03 -1.848140e-06 -2.031621e-08 5.100882e-10 -2.588485e-12
## [12,] -2.662642e-02 5.906348e-06 -5.716248e-09 -6.220002e-09 5.187373e-12
## [13,] -7.293254e-04 1.635124e-06 -2.078163e-10 -4.077588e-10 3.705165e-13
## [14,] 1.174257e-11 2.735617e-15 2.937778e-05 1.768113e-05 -6.508111e-08
## [15,] -7.008085e-09 -2.425713e-12 -1.533530e-02 -1.796480e-02 1.241710e-06
## [16,] -5.609351e-09 -1.730943e-12 -1.320097e-02 -1.106668e-02 -3.321188e-05
## [,11] [,12] [,13] [,14] [,15]
## [1,] 2.861680e-02 -2.382350e-03 -2.709489e-04 -1.935033e-09 -4.227225e-08
## [2,] -3.139979e-07 3.152343e-08 2.918004e-09 1.511147e-04 -1.536106e-02
## [3,] 3.895279e-07 4.398105e-06 3.040644e-07 7.514155e-09 1.291319e-06
## [4,] -1.132166e-06 4.486090e-07 3.745140e-08 2.418647e-05 8.372008e-04
## [5,] -6.470501e-03 -1.505899e-02 8.073588e-04 -6.347924e-11 1.079039e-09
## [6,] -5.726829e-03 -2.662642e-02 -7.293254e-04 1.174257e-11 -7.008085e-09
## [7,] -1.848140e-06 5.906348e-06 1.635124e-06 2.735617e-15 -2.425713e-12
## [8,] -2.031621e-08 -5.716248e-09 -2.078163e-10 2.937778e-05 -1.533530e-02
## [9,] 5.100882e-10 -6.220002e-09 -4.077588e-10 1.768113e-05 -1.796480e-02
## [10,] -2.588485e-12 5.187373e-12 3.705165e-13 -6.508111e-08 1.241710e-06
## [11,] 8.147358e-01 -1.597045e-01 1.704864e-02 -2.399692e-11 -1.814001e-10
## [12,] -1.597045e-01 2.414046e+00 -3.720758e-02 -2.922371e-12 -2.739914e-09
## [13,] 1.704864e-02 -3.720758e-02 4.352086e-02 -5.370356e-14 -1.892795e-10
## [14,] -2.399692e-11 -2.922372e-12 -5.370356e-14 7.566872e-04 -1.590628e-03
## [15,] -1.814001e-10 -2.739914e-09 -1.892795e-10 -1.590628e-03 2.552177e+00
## [16,] 1.299943e-09 -1.452803e-09 -1.101203e-10 1.376592e-04 -3.953442e-01
## [,16]
## [1,] 9.620642e-08
## [2,] -1.012532e-02
## [3,] 3.537853e-07
## [4,] -1.987040e-03
## [5,] 4.008188e-09
## [6,] -5.609351e-09
## [7,] -1.730943e-12
## [8,] -1.320097e-02
## [9,] -1.106668e-02
## [10,] -3.321188e-05
## [11,] 1.299943e-09
## [12,] -1.452803e-09
## [13,] -1.101203e-10
## [14,] 1.376592e-04
## [15,] -3.953442e-01
## [16,] 2.461612e+00
##
## Model rank = 1346 / 1346
##
## 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):vowelContrastDRESS 11.00 7.70 0.96 0.65
## s(Time):vowelContrastTRAP 11.00 5.55 0.96 0.59
## s(Time):isAustralian 12.00 2.00 0.96 0.62
## s(Time):isAustralianTrap 10.00 5.19 0.96 0.63
## s(Time,pp):vowelContrastDRESS 450.00 88.73 0.96 0.68
## s(Time,pp):vowelContrastTRAP 450.00 69.98 0.96 0.60
## s(Time,item):vowelContrastDRESS 200.00 13.40 0.96 0.60
## s(Time,item):vowelContrastTRAP 200.00 13.45 0.96 0.68We now examine this model to determine whether there is an overall effect of contrast.
summary(interaction.compet_dressTrap.gamm1)##
## Family: binomial
## Link function: logit
##
## Formula:
## position ~ vowelContrast + s(Time, by = vowelContrast, k = 12) +
## s(Time, by = isAustralian, k = 12) + s(Time, by = isAustralianTrap) +
## 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.54766 0.05171 -29.927 <2e-16 ***
## vowelContrastTRAP -0.12903 0.09696 -1.331 0.183
## ---
## 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):vowelContrastDRESS 7.698 9.346 119.976 < 2e-16 ***
## s(Time):vowelContrastTRAP 5.545 7.127 70.882 < 2e-16 ***
## s(Time):isAustralian 2.004 2.007 2.691 0.261
## s(Time):isAustralianTrap 5.191 6.375 9.032 0.206
## s(Time,pp):vowelContrastDRESS 88.729 448.000 188.679 < 2e-16 ***
## s(Time,pp):vowelContrastTRAP 69.978 449.000 147.293 < 2e-16 ***
## s(Time,item):vowelContrastDRESS 13.396 98.000 34.496 4.92e-05 ***
## s(Time,item):vowelContrastTRAP 13.447 99.000 44.940 3.11e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0602 Deviance explained = 6.92%
## fREML = 82243 Scale est. = 1 n = 320044This model shows significant nonlinear time effects and substantial participant and item variability, but no significant main effect of vowel contrast (p = .183) and no significant accent-related smooths, indicating that competitor fixation trajectories change over time but do not reliably differ overall between DRESS and TRAP or across accent.
smooth <- plot_smooth(interaction.compet_dressTrap.gamm1, view = "Time",
plot_all = "vowelContrast",
rug = FALSE, rm.ranef = TRUE, lwd = 4)## Summary:
## * vowelContrast : factor; set to the value(s): DRESS, TRAP.
## * Time : numeric predictor; with 30 values ranging from 10.000000 to 890.000000.
## * isAustralian : numeric predictor; set to the value(s): 0.
## * isAustralianTrap : 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):vowelContrastDRESS,s(Time,pp):vowelContrastTRAP,s(Time,item):vowelContrastDRESS,s(Time,item):vowelContrastTRAP
## 
We now proceed to model comparison by successively removing effects to evaluate their contribution to model fit.
reducedModelFile = "compet_dressTrap_reducedModel.Rdata"
if (file.exists(reducedModelFile))
{
load(reducedModelFile)
}else{
mainEffect.compet_dressTrap.gamm <- bam( position ~ vowelContrast+
s(Time, by = vowelContrast) +
#s(Time, by = isAustralian) +
#s(Time, by = isAustralianTrap) +
s(Time, pp, by = vowelContrast, bs="fs") +
s(Time, item, by = vowelContrast, bs="fs"), family = "binomial",
data=gamData, rho = myRhoVal,
AR.start = gamData$start.event, nThreads = 3, discrete = T)
save(list= c("mainEffect.compet_dressTrap.gamm"), file = reducedModelFile)
}summary(mainEffect.compet_dressTrap.gamm)##
## Family: binomial
## Link function: logit
##
## Formula:
## position ~ vowelContrast + s(Time, by = vowelContrast) + 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.59969 0.03859 -41.45 <2e-16 ***
## vowelContrastTRAP -0.10061 0.08315 -1.21 0.226
## ---
## 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):vowelContrastDRESS 6.991 8.148 125.66 < 2e-16 ***
## s(Time):vowelContrastTRAP 6.100 7.455 116.86 < 2e-16 ***
## s(Time,pp):vowelContrastDRESS 88.634 448.000 188.34 < 2e-16 ***
## s(Time,pp):vowelContrastTRAP 70.198 448.000 148.01 < 2e-16 ***
## s(Time,item):vowelContrastDRESS 13.369 98.000 34.38 5.11e-05 ***
## s(Time,item):vowelContrastTRAP 13.486 98.000 45.05 3.06e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0596 Deviance explained = 6.86%
## fREML = 82316 Scale est. = 1 n = 320044This reduced model shows strong nonlinear time effects and substantial participant and item variability, but no significant main effect of vowel contrast (p = .226), indicating that competitor fixation trajectories change over time but do not reliably differ overall between DRESS and TRAP.
smooth <- plot_smooth(mainEffect.compet_dressTrap.gamm, view = "Time",
plot_all = "vowelContrast",
rug = FALSE, rm.ranef = TRUE, lwd = 4)## Summary:
## * vowelContrast : factor; set to the value(s): DRESS, TRAP.
## * Time : numeric predictor; with 30 values ranging from 10.000000 to 890.000000.
## * 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):vowelContrastDRESS,s(Time,pp):vowelContrastTRAP,s(Time,item):vowelContrastDRESS,s(Time,item):vowelContrastTRAP
## 