Vowel contrast: GAMM analysis

Note that this Markdown requires the file “mouse_data.txt”, which is generated by running the preprocessing markdown as well as the mousePreProcess.R script.

TARGET FIXATIONS

rm(list = ls())
library(lmerTest)

## Loading required package: lme4

## Loading required package: Matrix

## 
## Attaching package: 'lmerTest'

## The following object is masked from 'package:lme4':
## 
##     lmer

## The following object is masked from 'package:stats':
## 
##     step

library(stringr)
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.0     ✔ purrr     1.2.1
## ✔ forcats   1.0.1     ✔ readr     2.2.0
## ✔ ggplot2   4.0.2     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ tidyr::pack()   masks Matrix::pack()
## ✖ tidyr::unpack() masks Matrix::unpack()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(mgcv)

## Loading required package: nlme
## 
## Attaching package: 'nlme'
## 
## The following object is masked from 'package:dplyr':
## 
##     collapse
## 
## The following object is masked from 'package:lme4':
## 
##     lmList
## 
## This is mgcv 1.9-3. For overview type 'help("mgcv-package")'.

library(itsadug)

## Loading required package: plotfunctions
## 
## Attaching package: 'plotfunctions'
## 
## The following object is masked from 'package:ggplot2':
## 
##     alpha
## 
## Loaded package itsadug 2.5 (see 'help("itsadug")' ).

options(warn = 0)
gazeData = read.table("dataCatalan/gazeTarget.txt", T, sep = "\t", quote = "")
trialData = read.table("dataCatalan/trialInfo.txt", T, sep = "\t", quote = "\"")
expFilter = grepl("_", trialData$condition)

gazeDataE = gazeData %>% filter(expFilter)
trialDataE = trialData %>% filter(expFilter)
timeline = -200 + (0:120) * 10
colnames(gazeDataE) = paste0("T_",timeline)

#contrast coding
temp = str_split_fixed(trialDataE$condition,"_",2) 
trialDataE$accent = temp[, 1]
trialDataE$contrast = temp[, 2]

trialDataE$accentC = -0.5 + as.numeric(trialDataE$accent == "SEng") 
table(trialDataE$accentC)

## 
## -0.5  0.5 
## 3596 3598

trialDataE$contrastC = -0.5 + as.numeric(trialDataE$contrast == "DRESS" |trialDataE$contrast == "TRAP" )
table(trialDataE$contrastC)

## 
## -0.5  0.5 
## 3598 3596

useData = cbind(trialDataE,gazeDataE)

Analysis of the time course using GAMMs

We first need to reorganize the data.

gamData = useData[ ,c(1,2,4,11:14,15:135)] %>% 
               pivot_longer(., cols = starts_with("T_"),
                            names_to = "Time",
                            values_to = "position") %>% as.data.frame

gamData$position = as.numeric(gamData$position)
gamData$Time = gsub("T_","", gamData$Time) %>% as.numeric() 
gamData = gamData[ order(gamData$pp,gamData$trial,gamData$Time), ]


gamData$pp = as.factor(gamData$pp)
gamData$item = as.factor(gamData$item)

gamData$isAustralian = ifelse(gamData$accent == "AusE",1,0)
gamData$isAustralianDress = gamData$isAustralian *
                          as.numeric(gamData$contrast == "DRESS" |gamData$contrast == "TRAP" ) 

gamData$vowelContrast = ifelse(gamData$contrast == "DRESS" |gamData$contrast == "TRAP",
                               "midFront","tensenessI") %>% as.factor()

#restricting the data to 0-900 after vowel onset
gamData = gamData %>% filter(Time > 0 & Time < 900)
gamData <- start_event(gamData, event=c("pp","trial"))

gamDataPlot = gamData %>% group_by(vowelContrast, Time, accent)  %>%
                summarize(meanFix = mean(position))

## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by vowelContrast, Time, and accent.
## ℹ Output is grouped by vowelContrast and Time.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(vowelContrast, Time, accent))` for per-operation
##   grouping (`?dplyr::dplyr_by`) instead.

ggplot(gamDataPlot, aes(x = Time, y = meanFix,  col = vowelContrast)) + 
  geom_line() + facet_wrap(~accent)

This is the descriptive data

Building an initial GAMM

fullModelFile = "target_overall_fullModels.Rdata" 

if (file.exists(fullModelFile))
{
  load(fullModelFile)
}else{
  interaction.gamm0 <- bam( position ~  vowelContrast+
                       s(Time, by =  vowelContrast) + 
                       s(Time, by = isAustralian) +
                       s(Time, by = isAustralianDress) +
                       s(Time, pp, by =  vowelContrast, bs="fs") +
                       s(Time, item, by =  vowelContrast, bs="fs"), 
                       family = "binomial",
                       data=gamData, discrete = T, nThreads = 3)
  #dealing with autocorrelation
  myAutoCorr  =acf_resid(interaction.gamm0)
  myRhoVal = myAutoCorr[2]
  interaction.gamm1 <-  bam( position ~  vowelContrast+
                       s(Time, by =  vowelContrast, k=12) + 
                       s(Time, by = isAustralian, k= 12) +
                       s(Time, by = isAustralianDress) +
                       s(Time, pp, by =  vowelContrast, bs="fs") +
                       s(Time, item, by =  vowelContrast, bs="fs"), 
                       family = "binomial",
                       data=gamData, rho = myRhoVal, 
                       AR.start = gamData$start.event, nThreads = 3, discrete = T)
 save(list= c("interaction.gamm0","interaction.gamm1"), file = fullModelFile)
}

Model criticism

How are we doing regarding auto-correlation?

myAutoCorr  =acf_resid(interaction.gamm1)

Is the default for the smooth (k = 10) sufficient?

gam.check(interaction.gamm1)

## 
## Method: fREML   Optimizer: perf chol
## $grad
##  [1]  0.0069406337  0.0024180552 -0.0022415194 -0.0019898415  0.0314626836
##  [6]  0.0153926254  0.0354470408  0.0445616939  0.0195573102  0.0377280812
## [11]  0.0048887509  0.0066139624  0.0002983201  0.0079791663  0.0053302942
## [16]  0.0053596874
## 
## $hess
##                [,1]          [,2]          [,3]          [,4]          [,5]
##  [1,]  3.340892e+00 -1.798569e-04  9.577181e-04  2.132693e-01  2.147290e-01
##  [2,] -1.798569e-04  3.053790e+00  2.571085e-03 -7.880004e-04 -4.850833e-06
##  [3,]  9.577181e-04  2.571085e-03  2.311638e-03  2.669913e-03  3.149169e-05
##  [4,]  2.132693e-01 -7.880004e-04  2.669913e-03  7.238158e-01  7.384829e-03
##  [5,]  2.147290e-01 -4.850833e-06  3.149169e-05  7.384829e-03  2.474965e+01
##  [6,]  1.514245e-02 -4.698830e-07  4.643919e-06 -4.979736e-04  2.700870e-01
##  [7,] -3.966766e-02  4.387085e-06 -2.615465e-05 -8.140549e-03 -4.895610e-01
##  [8,] -1.141793e-05  1.428172e-01  8.473813e-05 -2.807682e-05 -6.192438e-07
##  [9,]  7.968583e-07  4.237776e-02  2.217294e-05 -5.403145e-06  1.060016e-07
## [10,]  1.838134e-06 -4.716512e-02 -5.144949e-05  1.067204e-05  3.869345e-08
## [11,]  7.242385e-02 -2.295805e-06  1.297510e-05  1.937084e-03 -7.943797e-02
## [12,]  1.582360e-03  1.666258e-07  3.421469e-07 -4.328364e-04 -6.932169e-02
## [13,] -1.222761e-02  1.182772e-06 -7.226080e-06 -2.399865e-03 -4.764104e-02
## [14,] -5.168744e-06  6.239933e-02  5.907710e-05 -1.778927e-05 -2.573258e-07
## [15,]  3.484781e-07  8.836137e-03  4.127427e-07 -1.037451e-06  3.380736e-08
## [16,]  7.547161e-07 -3.197547e-02 -3.306720e-05  7.285799e-06 -7.595624e-09
##                [,6]          [,7]          [,8]          [,9]         [,10]
##  [1,]  1.514245e-02 -3.966766e-02 -1.141793e-05  7.968583e-07  1.838134e-06
##  [2,] -4.698830e-07  4.387085e-06  1.428172e-01  4.237776e-02 -4.716512e-02
##  [3,]  4.643919e-06 -2.615465e-05  8.473813e-05  2.217294e-05 -5.144949e-05
##  [4,] -4.979736e-04 -8.140549e-03 -2.807682e-05 -5.403145e-06  1.067204e-05
##  [5,]  2.700870e-01 -4.895610e-01 -6.192438e-07  1.060016e-07  3.869345e-08
##  [6,]  1.108683e+01 -7.695496e-01  5.950836e-08 -2.945867e-08  1.745073e-08
##  [7,] -7.695496e-01  1.662183e+01  1.237807e-07  4.088049e-08 -6.623594e-08
##  [8,]  5.950836e-08  1.237807e-07  2.069404e+01 -3.315983e-01 -6.395533e-01
##  [9,] -2.945867e-08  4.088049e-08 -3.315983e-01  1.399071e+01 -5.609156e-01
## [10,]  1.745073e-08 -6.623594e-08 -6.395533e-01 -5.609156e-01  1.570437e+01
## [11,] -3.558707e-02 -5.819340e-02 -2.064639e-07  2.566554e-08  2.269722e-08
## [12,] -3.876018e-02 -2.676317e-02  3.358399e-08 -8.862064e-09  1.955237e-09
## [13,] -1.524750e-02 -5.407243e-02  2.361057e-08  1.411396e-08 -1.895829e-08
## [14,]  1.108600e-08  8.871286e-08 -8.433821e-02 -4.216662e-02 -8.025604e-02
## [15,] -7.207956e-09  7.932692e-09 -5.520802e-02 -4.532412e-02 -2.846319e-02
## [16,]  1.572271e-08 -4.623023e-08 -4.613904e-02 -2.238309e-02 -6.695430e-02
##               [,11]         [,12]         [,13]         [,14]         [,15]
##  [1,]  7.242385e-02  1.582360e-03 -1.222761e-02 -5.168744e-06  3.484781e-07
##  [2,] -2.295805e-06  1.666258e-07  1.182772e-06  6.239933e-02  8.836137e-03
##  [3,]  1.297510e-05  3.421469e-07 -7.226080e-06  5.907710e-05  4.127427e-07
##  [4,]  1.937084e-03 -4.328364e-04 -2.399865e-03 -1.778927e-05 -1.037451e-06
##  [5,] -7.943797e-02 -6.932169e-02 -4.764104e-02 -2.573258e-07  3.380736e-08
##  [6,] -3.558707e-02 -3.876018e-02 -1.524750e-02  1.108600e-08 -7.207956e-09
##  [7,] -5.819340e-02 -2.676317e-02 -5.407243e-02  8.871286e-08  7.932692e-09
##  [8,] -2.064639e-07  3.358399e-08  2.361057e-08 -8.433821e-02 -5.520802e-02
##  [9,]  2.566554e-08 -8.862064e-09  1.411396e-08 -4.216662e-02 -4.532412e-02
## [10,]  2.269722e-08  1.955237e-09 -1.895829e-08 -8.025604e-02 -2.846319e-02
## [11,]  6.437702e+00 -1.742088e-01 -2.092983e-01 -8.719682e-08  8.840444e-09
## [12,] -1.742088e-01  3.545384e+00 -6.236342e-02  1.148506e-08 -2.224622e-09
## [13,] -2.092983e-01 -6.236342e-02  4.774143e+00  2.164207e-08  2.697182e-09
## [14,] -8.719682e-08  1.148506e-08  2.164207e-08  8.121517e+00 -1.602243e-01
## [15,]  8.840444e-09 -2.224622e-09  2.697182e-09 -1.602243e-01  3.910068e+00
## [16,]  4.075270e-09  3.163292e-09 -1.405559e-08 -3.110829e-01  4.323348e-02
##               [,16]
##  [1,]  7.547161e-07
##  [2,] -3.197547e-02
##  [3,] -3.306720e-05
##  [4,]  7.285799e-06
##  [5,] -7.595624e-09
##  [6,]  1.572271e-08
##  [7,] -4.623023e-08
##  [8,] -4.613904e-02
##  [9,] -2.238309e-02
## [10,] -6.695430e-02
## [11,]  4.075270e-09
## [12,]  3.163292e-09
## [13,] -1.405559e-08
## [14,] -3.110829e-01
## [15,]  4.323348e-02
## [16,]  6.009899e+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):vowelContrastmidFront         11.00   8.67    0.97   0.040 *
## s(Time):vowelContrasttensenessI       11.00   8.15    0.97   0.040 *
## s(Time):isAustralian                  12.00   2.03    0.97   0.035 *
## s(Time):isAustralianDress             10.00   5.16    0.97   0.030 *
## s(Time,pp):vowelContrastmidFront     450.00 152.01    0.97   0.045 *
## s(Time,pp):vowelContrasttensenessI   450.00 157.48    0.97   0.045 *
## s(Time,item):vowelContrastmidFront   400.00  51.94    0.97   0.040 *
## s(Time,item):vowelContrasttensenessI 400.00  61.49    0.97   0.025 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Let’s have a look at this model no whether there is an overall effect of contrast

summary(interaction.gamm1)

## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## position ~ vowelContrast + s(Time, by = vowelContrast, k = 12) + 
##     s(Time, by = isAustralian, k = 12) + s(Time, by = isAustralianDress) + 
##     s(Time, pp, by = vowelContrast, bs = "fs") + s(Time, item, 
##     by = vowelContrast, bs = "fs")
## 
## Parametric coefficients:
##                         Estimate Std. Error z value Pr(>|z|)   
## (Intercept)              0.23484    0.07650   3.070  0.00214 **
## vowelContrasttensenessI  0.01743    0.10737   0.162  0.87103   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                          edf  Ref.df  Chi.sq p-value    
## s(Time):vowelContrastmidFront          8.672  10.003 628.280  <2e-16 ***
## s(Time):vowelContrasttensenessI        8.146   9.644 472.640  <2e-16 ***
## s(Time):isAustralian                   2.032   2.059   0.389   0.714    
## s(Time):isAustralianDress              5.159   6.370   5.516   0.511    
## s(Time,pp):vowelContrastmidFront     152.015 449.000 714.637  <2e-16 ***
## s(Time,pp):vowelContrasttensenessI   157.481 449.000 714.361  <2e-16 ***
## s(Time,item):vowelContrastmidFront    51.936 199.000 210.473  <2e-16 ***
## s(Time,item):vowelContrasttensenessI  61.485 199.000 269.413  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.254   Deviance explained = 20.3%
## fREML = 2.1037e+05  Scale est. = 1         n = 640266

There is a main effect of competitor, meaning that, overall, the mouse position is closer to target when the competitor is dissimilar to the target. Note that the significant effects of the tensors only indicate changes over time, they do not constitute evidence for a significant difference between conditions.

smooth <-   plot_smooth(interaction.gamm1, view = "Time", 
                        plot_all = "vowelContrast", 
                        rug = FALSE, rm.ranef = TRUE, lwd = 4)

## Summary:
##  * vowelContrast : factor; set to the value(s): midFront, tensenessI. 
##  * Time : numeric predictor; with 30 values ranging from 10.000000 to 890.000000. 
##  * isAustralian : numeric predictor; set to the value(s): 0. 
##  * isAustralianDress : numeric predictor; set to the value(s): 0. 
##  * pp : factor; set to the value(s): CAT01. (Might be canceled as random effect, check below.) 
##  * item : factor; set to the value(s): baggage. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(Time,pp):vowelContrastmidFront,s(Time,pp):vowelContrasttensenessI,s(Time,item):vowelContrastmidFront,s(Time,item):vowelContrasttensenessI
## 

Now model comparison

We succesively take effects out to prepare model comparisons.

reducedModelFile = "target_overall_reducedModels.Rdata" 

if (file.exists(reducedModelFile))
{
  load(reducedModelFile)
}else{

  minusInteraction.gamm <-  bam( position ~  vowelContrast+
                       s(Time, by =  vowelContrast) + 
                       s(Time, by = isAustralian) +
                       #s(Time, by = isAustralianDress) +
                       s(Time, pp, by =  vowelContrast, bs="fs") +
                       s(Time, item, by =  vowelContrast, bs="fs"), 
                       family = "binomial",
                       data=gamData, rho = myRhoVal,
                       AR.start = gamData$start.event, nThreads = 3, discrete = T)
  
  minusAccent.gamm <-  bam( position ~  vowelContrast+
                       s(Time, by =  vowelContrast) + 
                       #s(Time, by = isAustralian) +
                       #s(Time, by = isAustralianDress) +
                       s(Time, pp, by =  vowelContrast, bs="fs") +
                       s(Time, item, by =  vowelContrast, bs="fs"), 
                       family = "binomial",
                       data=gamData, rho = myRhoVal, 
                       AR.start = gamData$start.event, nThreads = 3, discrete = T)
  
  timeOnly.gamm  <-  bam( position ~  1 +
                       s(Time) + 
                       #s(Time, by = isAustralian) +
                       #s(Time, by = isAustralianDress) +
                       s(Time, pp, by =  vowelContrast, bs="fs") +
                       s(Time, item, by =  vowelContrast, bs="fs"), family = "binomial",
                       data=gamData, rho = myRhoVal, 
                       AR.start = gamData$start.event, nThreads = 3, discrete = T)
  
 save(list= c("minusInteraction.gamm","minusAccent.gamm","timeOnly.gamm"), file = reducedModelFile)
}

Model comparisons

Full model versus no interaction

compareML(interaction.gamm1, minusInteraction.gamm)

## interaction.gamm1: position ~ vowelContrast + s(Time, by = vowelContrast, k = 12) + 
##     s(Time, by = isAustralian, k = 12) + s(Time, by = isAustralianDress) + 
##     s(Time, pp, by = vowelContrast, bs = "fs") + s(Time, item, 
##     by = vowelContrast, bs = "fs")
## 
## minusInteraction.gamm: position ~ vowelContrast + s(Time, by = vowelContrast) + s(Time, 
##     by = isAustralian) + s(Time, pp, by = vowelContrast, bs = "fs") + 
##     s(Time, item, by = vowelContrast, bs = "fs")
## 
## Model minusInteraction.gamm preferred: lower fREML score (145.423), and lower df (1.000).
## -----
##                   Model    Score Edf Difference    Df
## 1     interaction.gamm1 210367.6  16                 
## 2 minusInteraction.gamm 210222.2  15   -145.423 1.000
## 
## AIC difference: 492.28, model minusInteraction.gamm has lower AIC.

Both main effects versus minus accent

compareML(minusInteraction.gamm, minusAccent.gamm)

## minusInteraction.gamm: position ~ vowelContrast + s(Time, by = vowelContrast) + s(Time, 
##     by = isAustralian) + s(Time, pp, by = vowelContrast, bs = "fs") + 
##     s(Time, item, by = vowelContrast, bs = "fs")
## 
## minusAccent.gamm: position ~ vowelContrast + s(Time, by = vowelContrast) + s(Time, 
##     pp, by = vowelContrast, bs = "fs") + s(Time, item, by = vowelContrast, 
##     bs = "fs")
## 
## Model minusAccent.gamm preferred: lower fREML score (219.395), and lower df (1.000).
## -----
##                   Model    Score Edf Difference    Df
## 1 minusInteraction.gamm 210222.2  15                 
## 2      minusAccent.gamm 210002.8  14   -219.395 1.000
## 
## AIC difference: 972.60, model minusAccent.gamm has lower AIC.

Main effect of Contrast versus time-only

compareML(minusAccent.gamm, timeOnly.gamm)

## minusAccent.gamm: position ~ vowelContrast + s(Time, by = vowelContrast) + s(Time, 
##     pp, by = vowelContrast, bs = "fs") + s(Time, item, by = vowelContrast, 
##     bs = "fs")
## 
## timeOnly.gamm: position ~ 1 + s(Time) + s(Time, pp, by = vowelContrast, bs = "fs") + 
##     s(Time, item, by = vowelContrast, bs = "fs")
## 
## Model timeOnly.gamm preferred: lower fREML score (137.133), and lower df (1.000).
## -----
##              Model    Score Edf Difference    Df
## 1 minusAccent.gamm 210002.8  14                 
## 2    timeOnly.gamm 209865.7  13   -137.133 1.000
## 
## AIC difference: 613.35, model timeOnly.gamm has lower AIC.

COMPETITOR FIXATIONS

Setting up the environment and then Getting all the mouse positions as a matrix

rm(list = ls())
library(lmerTest)
library(stringr)
library(tidyverse)
library(mgcv)
library(itsadug)
options(warn = 0)
gazeData = read.table("dataCatalan/gazeCompet.txt", T, sep = "\t", quote = "")
trialData = read.table("dataCatalan/trialInfo.txt", T, sep = "\t", quote = "\"")
expFilter = grepl("_", trialData$condition)

gazeDataE = gazeData %>% filter(expFilter)
trialDataE = trialData %>% filter(expFilter)
timeline = -200 + (0:120) * 10
colnames(gazeDataE) = paste0("T_",timeline)

#contrast coding
temp = str_split_fixed(trialDataE$condition,"_",2) 
trialDataE$accent = temp[, 1]
trialDataE$contrast = temp[, 2]

trialDataE$accentC = -0.5 + as.numeric(trialDataE$accent == "SEng") 
table(trialDataE$accentC)

## 
## -0.5  0.5 
## 3596 3598

trialDataE$contrastC = -0.5 + as.numeric(trialDataE$contrast == "DRESS" |trialDataE$contrast == "TRAP" )
table(trialDataE$contrastC)

## 
## -0.5  0.5 
## 3598 3596

useData = cbind(trialDataE,gazeDataE)

Analysis of the time course using GAMMs

We first need to reorganize the data.

gamData = useData[ ,c(1,2,4,11:14,15:135)] %>% 
               pivot_longer(., cols = starts_with("T_"),
                            names_to = "Time",
                            values_to = "position") %>% as.data.frame

gamData$position = as.numeric(gamData$position)
gamData$Time = gsub("T_","", gamData$Time) %>% as.numeric() 
gamData = gamData[ order(gamData$pp,gamData$trial,gamData$Time), ]


gamData$pp = as.factor(gamData$pp)
gamData$item = as.factor(gamData$item)

gamData$isAustralian = ifelse(gamData$accent == "AusE",1,0)
gamData$isAustralianDress = gamData$isAustralian *
                          as.numeric(gamData$contrast == "DRESS" |gamData$contrast == "TRAP" ) 

gamData$vowelContrast = ifelse(gamData$contrast == "DRESS" |gamData$contrast == "TRAP",
                               "midFront","tensenessI") %>% as.factor()

#restricting the data to 0-900 after vowel onset
gamData = gamData %>% filter(Time > 0 & Time < 900)
gamData <- start_event(gamData, event=c("pp","trial"))

gamDataPlot = gamData %>% group_by(vowelContrast, Time, accent)  %>%
                summarize(meanFix = mean(position))

## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by vowelContrast, Time, and accent.
## ℹ Output is grouped by vowelContrast and Time.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(vowelContrast, Time, accent))` for per-operation
##   grouping (`?dplyr::dplyr_by`) instead.

ggplot(gamDataPlot, aes(x = Time, y = meanFix,  col = vowelContrast)) + 
  geom_line() + facet_wrap(~accent)

This is the descriptive data

Building an initial GAMM

fullModelFile = "competitor_overall_fullModels.Rdata" 

if (file.exists(fullModelFile))
{
  load(fullModelFile)
}else{
  interaction.Comp.gamm0 <- bam( position ~  vowelContrast+
                       s(Time, by =  vowelContrast) + 
                       s(Time, by = isAustralian) +
                       s(Time, by = isAustralianDress) +
                       s(Time, pp, by =  vowelContrast, bs="fs") +
                       s(Time, item, by =  vowelContrast, bs="fs"), 
                       family = "binomial",
                       data=gamData, discrete = T, nThreads = 3)
  #dealing with autocorrelation
  myAutoCorr  =acf_resid(interaction.Comp.gamm0)
  myRhoVal = myAutoCorr[2]
  interaction.Comp.gamm1 <-  bam( position ~  vowelContrast+
                       s(Time, by =  vowelContrast, k=12) + 
                       s(Time, by = isAustralian, k= 12) +
                       s(Time, by = isAustralianDress) +
                       s(Time, pp, by =  vowelContrast, bs="fs") +
                       s(Time, item, by =  vowelContrast, bs="fs"), 
                       family = "binomial",
                       data=gamData, rho = myRhoVal,
                       AR.start = gamData$start.event, nThreads = 3, discrete = T)
 save(list= c("interaction.Comp.gamm0","interaction.Comp.gamm1"), file = fullModelFile)
}

Model criticism

How are we doing regarding auto-correlation?

myAutoCorr  =acf_resid(interaction.Comp.gamm1)

Is the default for the smooth (k = 10) sufficient?

gam.check(interaction.Comp.gamm1)

## 
## Method: fREML   Optimizer: perf chol
## $grad
##  [1]  4.344936e-09 -6.668310e-11 -3.617088e-03 -2.361695e-03 -1.437751e-08
##  [6]  4.397094e-09 -1.403622e-07 -3.270612e-08  1.320329e-08 -1.223585e-07
## [11] -3.114113e-10 -9.537457e-09 -3.080776e-08 -5.261800e-09 -1.131491e-08
## [16] -3.280221e-09
## 
## $hess
##                [,1]          [,2]          [,3]          [,4]          [,5]
##  [1,]  3.635382e+00 -1.376290e-04  8.619980e-04  8.586773e-04  1.138355e-01
##  [2,] -1.376290e-04  3.880525e+00  7.466937e-04 -1.558598e-06 -4.275009e-06
##  [3,]  8.619980e-04  7.466937e-04  3.597892e-03 -4.456656e-06  3.025956e-05
##  [4,]  8.586773e-04 -1.558598e-06 -4.456656e-06  2.349919e-03  3.030709e-05
##  [5,]  1.138355e-01 -4.275009e-06  3.025956e-05  3.030709e-05  1.292494e+01
##  [6,] -3.700456e-02  2.276860e-06 -1.621932e-05 -1.369423e-05 -9.831543e-01
##  [7,] -7.906943e-03  2.875409e-07 -2.106202e-06 -2.068759e-06  1.255385e-01
##  [8,] -1.602161e-06  3.773216e-02  6.891158e-06 -8.127057e-09 -7.072587e-08
##  [9,]  1.045072e-06 -1.702599e-02 -6.433961e-06  3.504716e-08  6.856727e-08
## [10,]  4.089598e-07 -1.273209e-02 -2.499181e-06  7.974350e-09  2.095914e-08
## [11,]  2.629652e-02 -5.265767e-07  3.552015e-06  4.319390e-06 -8.710586e-03
## [12,] -5.603244e-03  1.024203e-06 -7.201901e-06 -5.408272e-06 -1.884878e-02
## [13,] -7.166771e-03  2.895458e-07 -1.989687e-06 -1.841388e-06 -1.822608e-02
## [14,] -2.529303e-06  5.614114e-02  1.268922e-05 -2.578752e-08 -1.389110e-07
## [15,]  6.679473e-07 -7.128001e-03 -4.334080e-06  1.994009e-08  3.573136e-08
## [16,]  2.588168e-07 -1.940266e-02 -1.579194e-06 -1.970647e-09  1.683026e-09
##                [,6]          [,7]          [,8]          [,9]         [,10]
##  [1,] -3.700456e-02 -7.906943e-03 -1.602161e-06  1.045072e-06  4.089598e-07
##  [2,]  2.276860e-06  2.875409e-07  3.773216e-02 -1.702599e-02 -1.273209e-02
##  [3,] -1.621932e-05 -2.106202e-06  6.891158e-06 -6.433961e-06 -2.499181e-06
##  [4,] -1.369423e-05 -2.068759e-06 -8.127057e-09  3.504716e-08  7.974350e-09
##  [5,] -9.831543e-01  1.255385e-01 -7.072587e-08  6.856727e-08  2.095914e-08
##  [6,]  1.261535e+01 -9.316762e-01  1.032658e-08 -7.923694e-08 -1.505821e-08
##  [7,] -9.316762e-01  2.657327e+00  2.262021e-09 -6.058032e-09 -1.609988e-09
##  [8,]  1.032658e-08  2.262021e-09  1.015215e+01 -1.013432e+00 -8.810777e-02
##  [9,] -7.923694e-08 -6.058032e-09 -1.013432e+00  1.368019e+01 -1.156318e+00
## [10,] -1.505821e-08 -1.609988e-09 -8.810777e-02 -1.156318e+00  5.072175e+00
## [11,] -9.425580e-03 -2.471591e-03 -1.264574e-08 -4.784365e-09  9.058073e-10
## [12,] -3.369710e-02 -4.802659e-03  6.349869e-11 -4.345179e-08 -7.465993e-09
## [13,] -1.632764e-02 -6.561908e-03  1.299653e-09 -5.977144e-09 -1.389967e-09
## [14,]  4.858758e-08  6.524317e-09 -4.271093e-02 -2.416663e-02 -1.883049e-02
## [15,] -4.188182e-08 -3.449844e-09 -1.691265e-02 -3.940526e-02 -1.346172e-02
## [16,]  9.881644e-09  1.251340e-10 -4.067646e-02 -2.646676e-02 -1.665774e-02
##               [,11]         [,12]         [,13]         [,14]         [,15]
##  [1,]  2.629652e-02 -5.603244e-03 -7.166771e-03 -2.529303e-06  6.679473e-07
##  [2,] -5.265767e-07  1.024203e-06  2.895458e-07  5.614114e-02 -7.128001e-03
##  [3,]  3.552015e-06 -7.201901e-06 -1.989687e-06  1.268922e-05 -4.334080e-06
##  [4,]  4.319390e-06 -5.408272e-06 -1.841388e-06 -2.578752e-08  1.994009e-08
##  [5,] -8.710586e-03 -1.884878e-02 -1.822608e-02 -1.389110e-07  3.573136e-08
##  [6,] -9.425580e-03 -3.369710e-02 -1.632764e-02  4.858758e-08 -4.188182e-08
##  [7,] -2.471591e-03 -4.802659e-03 -6.561908e-03  6.524317e-09 -3.449844e-09
##  [8,] -1.264574e-08  6.349865e-11  1.299653e-09 -4.271093e-02 -1.691265e-02
##  [9,] -4.784365e-09 -4.345179e-08 -5.977144e-09 -2.416663e-02 -3.940526e-02
## [10,]  9.058073e-10 -7.465993e-09 -1.389967e-09 -1.883049e-02 -1.346172e-02
## [11,]  1.097995e+00 -2.705718e-01  5.232447e-02 -1.882327e-08 -2.334746e-09
## [12,] -2.705718e-01  5.151456e+00 -6.577427e-01  1.434919e-08 -2.315007e-08
## [13,]  5.232447e-02 -6.577427e-01  3.283239e+00  3.786390e-09 -3.514842e-09
## [14,] -1.882327e-08  1.434919e-08  3.786390e-09  6.572665e+00 -3.769508e-03
## [15,] -2.334746e-09 -2.315007e-08 -3.514842e-09 -3.769508e-03  4.759903e+00
## [16,]  3.295020e-09  6.427254e-09  2.429660e-10  1.597832e-01 -4.868616e-01
##               [,16]
##  [1,]  2.588168e-07
##  [2,] -1.940266e-02
##  [3,] -1.579194e-06
##  [4,] -1.970647e-09
##  [5,]  1.683026e-09
##  [6,]  9.881644e-09
##  [7,]  1.251340e-10
##  [8,] -4.067646e-02
##  [9,] -2.646676e-02
## [10,] -1.665774e-02
## [11,]  3.295020e-09
## [12,]  6.427254e-09
## [13,]  2.429660e-10
## [14,]  1.597832e-01
## [15,] -4.868616e-01
## [16,]  5.361220e+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):vowelContrastmidFront         11.00   8.11    0.96    0.17
## s(Time):vowelContrasttensenessI       11.00   8.79    0.96    0.23
## s(Time):isAustralian                  12.00   2.01    0.96    0.20
## s(Time):isAustralianDress             10.00   2.01    0.96    0.20
## s(Time,pp):vowelContrastmidFront     450.00 121.50    0.96    0.20
## s(Time,pp):vowelContrasttensenessI   450.00 120.12    0.96    0.17
## s(Time,item):vowelContrastmidFront   400.00  35.44    0.96    0.20
## s(Time,item):vowelContrasttensenessI 400.00  61.45    0.96    0.23

Let’s have a look at this model no whether there is an overall effect of contrast

summary(interaction.Comp.gamm1)

## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## position ~ vowelContrast + s(Time, by = vowelContrast, k = 12) + 
##     s(Time, by = isAustralian, k = 12) + s(Time, by = isAustralianDress) + 
##     s(Time, pp, by = vowelContrast, bs = "fs") + s(Time, item, 
##     by = vowelContrast, bs = "fs")
## 
## Parametric coefficients:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             -1.61502    0.05197 -31.076   <2e-16 ***
## vowelContrasttensenessI  0.06879    0.08114   0.848    0.397    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                          edf  Ref.df  Chi.sq p-value    
## s(Time):vowelContrastmidFront          8.114   9.673 168.256  <2e-16 ***
## s(Time):vowelContrasttensenessI        8.790  10.120 229.128  <2e-16 ***
## s(Time):isAustralian                   2.013   2.025   0.704   0.772    
## s(Time):isAustralianDress              2.009   2.018   0.327   0.880    
## s(Time,pp):vowelContrastmidFront     121.499 449.000 352.710  <2e-16 ***
## s(Time,pp):vowelContrasttensenessI   120.124 448.000 354.939  <2e-16 ***
## s(Time,item):vowelContrastmidFront    35.440 199.000 146.981  <2e-16 ***
## s(Time,item):vowelContrasttensenessI  61.451 198.000 209.584  <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.6334e+05  Scale est. = 1         n = 640266

There is a main effect of competitor, meaning that, overall, the mouse position is closer to target when the competitor is dissimilar to the target. Note that the significant effects of the tensors only indicate changes over time, they do not constitute evidence for a significant difference between conditions.

smooth <-   plot_smooth(interaction.Comp.gamm1, view = "Time", 
                        plot_all = "vowelContrast", 
                        rug = FALSE, rm.ranef = TRUE, lwd = 4)

## Summary:
##  * vowelContrast : factor; set to the value(s): midFront, tensenessI. 
##  * Time : numeric predictor; with 30 values ranging from 10.000000 to 890.000000. 
##  * isAustralian : numeric predictor; set to the value(s): 0. 
##  * isAustralianDress : numeric predictor; set to the value(s): 0. 
##  * pp : factor; set to the value(s): CAT01. (Might be canceled as random effect, check below.) 
##  * item : factor; set to the value(s): baggage. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(Time,pp):vowelContrastmidFront,s(Time,pp):vowelContrasttensenessI,s(Time,item):vowelContrastmidFront,s(Time,item):vowelContrasttensenessI
## 

Now model comparison

We succesively take effects out to prepare model comparisons.

reducedModelFile = "competitor_overall_reducedModels.Rdata" 

if (file.exists(reducedModelFile))
{
  load(reducedModelFile)
}else{

  minusInteraction.Comp.gamm <-  bam( position ~  vowelContrast+
                       s(Time, by =  vowelContrast) + 
                       s(Time, by = isAustralian) +
                       #s(Time, by = isAustralianDress) +
                       s(Time, pp, by =  vowelContrast, bs="fs") +
                       s(Time, item, by =  vowelContrast, bs="fs"), 
                       family = "binomial",
                       data=gamData, rho = myRhoVal, 
                       AR.start = gamData$start.event, nThreads = 3, discrete = T)
  
  minusAccent.Comp.gamm <-  bam( position ~  vowelContrast+
                       s(Time, by =  vowelContrast) + 
                       #s(Time, by = isAustralian) +
                       #s(Time, by = isAustralianDress) +
                       s(Time, pp, by =  vowelContrast, bs="fs") +
                       s(Time, item, by =  vowelContrast, bs="fs"), 
                       family = "binomial",
                       data=gamData, rho = myRhoVal, 
                       AR.start = gamData$start.event, nThreads = 3, discrete = T)
  
  timeOnly.Comp.gamm  <-  bam( position ~  1 +
                       s(Time) + 
                       #s(Time, by = isAustralian) +
                       #s(Time, by = isAustralianDress) +
                       s(Time, pp, by =  vowelContrast, bs="fs") +
                       s(Time, item, by =  vowelContrast, bs="fs"), 
                       family = "binomial",
                       data=gamData, rho = myRhoVal,
                       AR.start = gamData$start.event, nThreads = 3, discrete = T)
  
 save(list= c("minusInteraction.Comp.gamm","minusAccent.Comp.gamm","timeOnly.Comp.gamm"), file = reducedModelFile)
}

Model comparisons

Full model versus no interaction

compareML(interaction.Comp.gamm1, minusInteraction.Comp.gamm)

## interaction.Comp.gamm1: position ~ vowelContrast + s(Time, by = vowelContrast, k = 12) + 
##     s(Time, by = isAustralian, k = 12) + s(Time, by = isAustralianDress) + 
##     s(Time, pp, by = vowelContrast, bs = "fs") + s(Time, item, 
##     by = vowelContrast, bs = "fs")
## 
## minusInteraction.Comp.gamm: position ~ vowelContrast + s(Time, by = vowelContrast) + s(Time, 
##     by = isAustralian) + s(Time, pp, by = vowelContrast, bs = "fs") + 
##     s(Time, item, by = vowelContrast, bs = "fs")
## 
## Model minusInteraction.Comp.gamm preferred: lower fREML score (111.479), and lower df (1.000).
## -----
##                        Model    Score Edf Difference    Df
## 1     interaction.Comp.gamm1 163337.1  16                 
## 2 minusInteraction.Comp.gamm 163225.6  15   -111.479 1.000
## 
## AIC difference: 221.24, model minusInteraction.Comp.gamm has lower AIC.

Both main effects versus minus accent

compareML(minusInteraction.Comp.gamm, minusAccent.Comp.gamm)

## minusInteraction.Comp.gamm: position ~ vowelContrast + s(Time, by = vowelContrast) + s(Time, 
##     by = isAustralian) + s(Time, pp, by = vowelContrast, bs = "fs") + 
##     s(Time, item, by = vowelContrast, bs = "fs")
## 
## minusAccent.Comp.gamm: position ~ vowelContrast + s(Time, by = vowelContrast) + s(Time, 
##     pp, by = vowelContrast, bs = "fs") + s(Time, item, by = vowelContrast, 
##     bs = "fs")
## 
## Model minusAccent.Comp.gamm preferred: lower fREML score (14.380), and lower df (1.000).
## -----
##                        Model    Score Edf Difference    Df
## 1 minusInteraction.Comp.gamm 163225.6  15                 
## 2      minusAccent.Comp.gamm 163211.3  14    -14.380 1.000
## 
## AIC difference: 374.14, model minusAccent.Comp.gamm has lower AIC.

Main effect of Contrast versus time-only

compareML(minusAccent.Comp.gamm, timeOnly.Comp.gamm)

## minusAccent.Comp.gamm: position ~ vowelContrast + s(Time, by = vowelContrast) + s(Time, 
##     pp, by = vowelContrast, bs = "fs") + s(Time, item, by = vowelContrast, 
##     bs = "fs")
## 
## timeOnly.Comp.gamm: position ~ 1 + s(Time) + s(Time, pp, by = vowelContrast, bs = "fs") + 
##     s(Time, item, by = vowelContrast, bs = "fs")
## 
## Model timeOnly.Comp.gamm preferred: lower fREML score (436.489), and lower df (1.000).
## -----
##                   Model    Score Edf Difference    Df
## 1 minusAccent.Comp.gamm 163211.3  14                 
## 2    timeOnly.Comp.gamm 162774.8  13   -436.489 1.000
## 
## AIC difference: 1703.05, model timeOnly.Comp.gamm has lower AIC.

FILLER ITEMS

Target fixations

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
useData$condition <- ifelse(
  useData$condition == "filler" & grepl("oo", useData$item, ignore.case = TRUE),
  "GOOSE",
  ifelse(useData$condition == "filler", "LOT", useData$condition)
)

useData$condition <- sub(".*_", "", useData$condition)

gamData = useData[ ,c(1,2,4:131)] %>% 
               pivot_longer(., cols = starts_with("T_"),
                            names_to = "Time",
                            values_to = "position") %>% as.data.frame

gamData$position = as.numeric(gamData$position)
gamData$Time = gsub("T_","", gamData$Time) %>% as.numeric() 

# Sort (like original)
gamData <- gamData[order(gamData$pp, gamData$trial, gamData$Time), ]

# Factors (like original)
gamData$pp   <- as.factor(gamData$pp)
gamData$item <- as.factor(gamData$item)

# Time window
gamData <- gamData %>% filter(Time > 0 & Time < 900)
gamData <- start_event(gamData, event = c("pp","trial"))

gamDataPlot = gamData %>% group_by(condition, Time)  %>%
                summarize(meanFix = mean(position))

## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by condition and Time.
## ℹ Output is grouped by condition.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(condition, Time))` for per-operation grouping
##   (`?dplyr::dplyr_by`) instead.

ggplot(gamDataPlot, aes(x = Time, y = meanFix,  col = condition)) + 
  geom_line() 

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.

Competitor fixations

library(tidyverse)
library(mgcv)
library(itsadug)
library(stringr)

options(warn = 0)

# Load data
gazeCompet <- read.table("dataCatalan/gazeCompet.txt",
                         header = TRUE, sep = "\t", quote = "")
trialData  <- read.table("dataCatalan/trialInfo.txt",
                         header = TRUE, sep = "\t", quote = "\"")

# Timeline
timeline <- -200 + (0:120) * 10
colnames(gazeCompet) <- paste0("T_", timeline)


# Combine datasets (same style as original)
useData <- cbind(trialData, gazeCompet) %>% as.data.frame()

# Remove practice trials
useData <- useData %>%
  filter(useData$condition != "practice") %>%
  droplevels()

# Rename levels of condition
useData$condition <- ifelse(
  useData$condition == "filler" & grepl("oo", useData$item, ignore.case = TRUE),
  "GOOSE",
  ifelse(useData$condition == "filler", "LOT", useData$condition)
)

useData$condition <- sub(".*_", "", useData$condition)

gamData = useData[ ,c(1,2,4:131)] %>% 
               pivot_longer(., cols = starts_with("T_"),
                            names_to = "Time",
                            values_to = "position") %>% as.data.frame

gamData$position = as.numeric(gamData$position)
gamData$Time = gsub("T_","", gamData$Time) %>% as.numeric() 

# Sort (like original)
gamData <- gamData[order(gamData$pp, gamData$trial, gamData$Time), ]

# Factors (like original)
gamData$pp   <- as.factor(gamData$pp)
gamData$item <- as.factor(gamData$item)

# Time window
gamData <- gamData %>% filter(Time > 0 & Time < 900)
gamData <- start_event(gamData, event = c("pp","trial"))

gamDataPlot = gamData %>% group_by(condition, Time)  %>%
                summarize(meanFix = mean(position))

## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by condition and Time.
## ℹ Output is grouped by condition.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(condition, Time))` for per-operation grouping
##   (`?dplyr::dplyr_by`) instead.

ggplot(gamDataPlot, aes(x = Time, y = meanFix,  col = condition)) + 
  geom_line() 

# STIMULUS MATERIALS



``` r
library(dplyr)
library(ggplot2)
library(ggrepel)

Stimuli_Acoustics <- read.csv("WP3_Stimuli_Acoustics.csv", row.names = NULL, stringsAsFactors = TRUE)

stimuli_plot <- Stimuli_Acoustics %>%
  mutate(
    Accent = factor(Accent, levels = c("S.Eng", "AusE")),
    Vowel = factor(toupper(Vowel))
  )

stimuli_plot <- stimuli_plot %>%
  filter(stimuli_plot$Contrast != "GOOSE-LOT") %>%
  droplevels()

means_va <- stimuli_plot %>%
  group_by(Accent, Vowel) %>%
  summarise(
    F1_mean = mean(F1, na.rm = TRUE),
    F2_mean = mean(F2, na.rm = TRUE),
    .groups = "drop"
  )

Fig.Stimuli <- ggplot(stimuli_plot, 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)