Data, documentation and analysis scripts for Analyzing dynamic phonetic data using generalized additive mixed modeling: a tutorial focusing on articulatory differences between L1 and L2 speakers of English
Martijn Wieling (University of Groningen, Haskins Laboratories, http://www.martijnwieling.nl)
Generated on: Sat Jun 5 20:03:02 2021
1 Abstract
In phonetics, many datasets are encountered which deal with dynamic data collected over time. Examples include diphthongal formant trajectories and articulator trajectories observed using electromagnetic articulography. Traditional approaches for analyzing this type of data generally aggregate data over a certain timespan, or only include measurements at a fixed time point (e.g., formant measurements at the midpoint of a vowel). In this paper, I discuss generalized additive modeling, a non-linear regression method which does not require aggregation or the pre-selection of a fixed time point. Instead, the method is able to identify general patterns over dynamically varying data, while simultaneously accounting for subject and item-related variability. An advantage of this approach is that patterns may be discovered which are hidden when data is aggregated or when a single time point is selected. A corresponding disadvantage is that these analyses are generally more time consuming and complex. This tutorial aims to overcome this disadvantage by providing a hands-on introduction to generalized additive modeling using articulatory trajectories from L1 and L2 speakers of English within the freely available R environment. All data and R code is made available to reproduce the analysis presented in this paper.
Journal: Journal of Phonetics (https://www.sciencedirect.com/science/article/abs/pii/S0095447017301377)
Preprint: http://www.martijnwieling.nl/files/GAM-tutorial-Wieling.pdf
Keywords: Generalized additive modeling, Tutorial, Articulography, Second language acquisition
2 Important note
A recent publication (Sóskuthy, 2021) has shown that the random-effects specification in the models below (from model m6 onwards) results in overly conservative \(p\)-values (i.e. \(p\)-values which are higher than they should be) when determining the fixed-effect differences between two levels (such as Word in model m6, or Sound in model ffmc0s). The reason for this is that the factor smooth does not take the dependency between the subjects with respect to the two levels of the fixed-effect factor into account, thereby moving some of the differences between the two levels to the random effects (and hence an estimated uncertainty about the difference in the fixed-effects which is too high, i.e. confidence bands which are too high). To obtain a more appropriate estimated significance of the difference curve, the random-effects specification (for m6) s(Time, Subject, by=Word, bs="fs", m=1) should be replaced by (with WordO as an orderdered factor, see Section 6.4.3) s(Time, Subject, bs="fs", m=1) + s(Time, Subject, by=WordO, bs="fs", m=1). Note, however, that this model is suitable to estimate the difference between the two words (Tenth and Tent), but should not be used to make inferences about the word-specific non-linear trajectories, as the confidence bands for the non-reference level will be overly large (overly conservative). Instead, model m6 should be used for that purpose. An example of this approach (for all words ending in "th") is shown in Section 7.4. Note that the other models in this tutorial which focus on assessing whether differences are significant are (still) overly conservative, and would need to be refit using the adjusted random-effects specification. This was not done here (except for the additional model shown in Section 7.4) to not deviate from the results discussed in the paper (which this tutorial accompanies).
3 Libraries and functions
The following commands load the necessary functions and libraries and show the version information.
# install packages if not yet installed
packages <- c("mgcv","itsadug","lme4","sp")
if (length(setdiff(packages, rownames(installed.packages()))) > 0) {
install.packages(setdiff(packages, rownames(installed.packages())))
}
# load required packages
library(mgcv)
library(itsadug)
library(sp) # for colors which also print well in grayscale
library(lme4)
# version information
R.version.string# [1] "R version 4.1.0 (2021-05-18)"cat(paste("mgcv version:",packageVersion("mgcv")))# mgcv version: 1.8.36cat(paste("itsadug version:",packageVersion("itsadug")))# itsadug version: 2.44 Dataset
The following shows the columns of the full dataset and their explanation.
if (!file.exists("full.rda")) {
download.file("http://www.let.rug.nl/wieling/Tutorial/full.rda","full.rda") # 2 MB
}
load("full.rda")4.1 Column names
The dataset consists of 126177 rows and 9 columns with the following column names:
str(full)# 'data.frame': 126177 obs. of 9 variables:
# $ Speaker: Factor w/ 42 levels "VENI_EN_1","VENI_EN_10",..: 1 1 1 1 1 1 1 1 1 1 ...
# $ Lang : Factor w/ 2 levels "EN","NL": 1 1 1 1 1 1 1 1 1 1 ...
# $ Word : Factor w/ 20 levels "faith","fate",..: 18 18 18 18 18 18 18 18 18 18 ...
# $ Sound : Factor w/ 2 levels "T","TH": 1 1 1 1 1 1 1 1 1 1 ...
# $ Loc : Factor w/ 2 levels "Final","Init": 2 2 2 2 2 2 2 2 2 2 ...
# $ Trial : int 1 1 1 1 1 1 1 1 1 1 ...
# $ Time : num 0 0.0161 0.0323 0.0484 0.0645 ...
# $ Pos : num -0.392 -0.44 -0.44 -0.503 -0.513 ...
# $ Pos01 : num 0.434 0.425 0.425 0.412 0.41 ...4.2 Data description
Speaker- ID of the speakerLang- Native language of the speaker ("NL"for Dutch, or"EN"for * English)Word- The label of the wordSound- The sound contrast ("TH"for words with the dental fricative,"T"for words with the stop)Loc- The location where in the word the sound contrasts occurs ("Init"when it occurs at the beginning of the word or"Final"when it occurs at the back of the word)Trial- The trial number of the word for each speakerTime- The normalized (between 0: beginning of the word, to 1: end of the word)Pos- The standardized (mean 0, standard deviation 1) position for each speaker of the T1 sensor in the anterior-posterior direction (higher values, more anterior)Pos01- The normalized (1: most anterior, 0: most posterior) position for each speaker of the T1 sensor
4.3 Example data
head(full)# Speaker Lang Word Sound Loc Trial Time Pos Pos01
# 1 VENI_EN_1 EN tick T Init 1 0.0000 -0.392 0.434
# 2 VENI_EN_1 EN tick T Init 1 0.0161 -0.440 0.425
# 3 VENI_EN_1 EN tick T Init 1 0.0323 -0.440 0.425
# 4 VENI_EN_1 EN tick T Init 1 0.0484 -0.503 0.412
# 5 VENI_EN_1 EN tick T Init 1 0.0645 -0.513 0.410
# 6 VENI_EN_1 EN tick T Init 1 0.0806 -0.677 0.3785 Visualizing basis functions
# Summary:
# * Word : factor; set to the value(s): tenth.
# * Time : numeric predictor; with 100 values ranging from 0.000000 to 1.000000.
6 Comparing "tenth" and "tent"
In the following, several models are fitted for the two words "tenth" and "tent".
words = c("tent","tenth")
dat = droplevels(full[full$Word %in% words,])
dat$Word = relevel(dat$Word,"tent") # tent is set as the reference level
# track first observation per trial (to correct for autocorrelation later on)
# and sort per trajectory
dat$start.event = (dat$Time == 0)
dat = dat[order(dat$Speaker,dat$Trial,dat$Time),]
# dat <- start_event(dat,event=c("Speaker","Trial")) # bug in itsadug 2.46.1 Fitting a linear model
m1 <- bam(Pos ~ Word, data=dat, method="fREML")
(smry1 <- summary(m1))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ Word
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 0.0654 0.0117 5.57 2.5e-08 ***
# Wordtenth 0.6642 0.0164 40.41 < 2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
#
# R-sq.(adj) = 0.113 Deviance explained = 11.3%
# -REML = 17307 Scale est. = 0.86694 n = 128396.2 Fitting a non-linear model
m2 <- bam(Pos ~ Word + s(Time, by=Word, bs="tp", k=10), data=dat)
(smry2 <- summary(m2))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ Word + s(Time, by = Word, bs = "tp", k = 10)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 0.0651 0.0107 6.11 1e-09 ***
# Wordtenth 0.6600 0.0149 44.18 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):Wordtent 7.52 8.46 28.4 <2e-16 ***
# s(Time):Wordtenth 8.55 8.94 276.2 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.267 Deviance explained = 26.8%
# fREML = 16112 Scale est. = 0.71584 n = 12839gam.check(m2)#
# Method: fREML Optimizer: perf newton
# full convergence after 9 iterations.
# Gradient range [-4.63e-07,3.87e-07]
# (score 16112 & scale 0.716).
# Hessian positive definite, eigenvalue range [2.95,6418].
# Model rank = 20 / 20
#
# 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):Wordtent 9.00 7.52 1 0.58
# s(Time):Wordtenth 9.00 8.55 1 0.606.3 Visualizing the results
par(mfrow=c(1,2))
plot(m2, ylim=c(-1,2), select=1)
abline(h=0)
plot(m2, ylim=c(-1,2), select=2)
abline(h=0)
two.colors = bpy.colors(n = 2, cutoff.tails = 0.7, alpha = 1.0)
plot_smooth(m2, view="Time", plot_all="Word", main="", rug=FALSE, ylim=c(-1,2), col=two.colors) # Summary:
# * Word : factor; set to the value(s): tent, tenth.
# * Time : numeric predictor; with 30 values ranging from 0.000000 to 1.000000.
# * NOTE : No random effects in the model to cancel.
# 
plot_diff(m2, view="Time", comp=list(Word=c("tenth","tent")), ylim=c(-1,2))# Summary:
# * Time : numeric predictor; with 100 values ranging from 0.000000 to 1.000000.
# * NOTE : No random effects in the model to cancel.
# 
#
# Time window(s) of significant difference(s):
# 0.070707 - 1.0000006.4 Complexity necessary?
6.4.1 Model comparison
m2a.ml <- bam(Pos ~ Word + s(Time), data=dat, method="ML")
m2b.ml <- bam(Pos ~ Word + s(Time, by=Word), data=dat, method="ML")
compareML(m2a.ml, m2b.ml)# m2a.ml: Pos ~ Word + s(Time)
#
# m2b.ml: Pos ~ Word + s(Time, by = Word)# Warning in sprintf("***", h1): one argument not used by format '***'#
# Chi-square test of ML scores
# -----
# Model Score Edf Difference Df p.value Sig.
# 1 m2a.ml 16505 4
# 2 m2b.ml 16103 6 401.810 2.000 < 2e-16 ***
#
# AIC difference: 823.83, model m2b.ml has lower AIC.6.4.2 Binary difference smooth
dat$IsTenth <- (dat$Word == "tenth")*1
m2.bin <- bam(Pos ~ s(Time) + s(Time,by=IsTenth), data=dat)
(smry2bin <- summary(m2.bin))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ s(Time) + s(Time, by = IsTenth)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 0.0650 0.0107 6.09 1.1e-09 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time) 7.69 8.49 28.8 <2e-16 ***
# s(Time):IsTenth 9.02 9.66 293.9 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.267 Deviance explained = 26.8%
# fREML = 16111 Scale est. = 0.71583 n = 12839plot(m2.bin, select=2, shade=TRUE, ylim=c(-1,2))
abline(h=0)
6.4.2.1 Binary predictors: use only once!
summary(bam(Pos ~ s(Time) + s(Time,by=IsTenth) + s(Trial) + s(Trial,by=IsTenth), data=dat))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ s(Time) + s(Time, by = IsTenth) + s(Trial) + s(Trial, by = IsTenth)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.0253 0.0261 -0.97 0.33
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time) 7.76 8.53 31.4 <2e-16 ***
# s(Time):IsTenth 8.08 8.70 108.3 <2e-16 ***
# s(Trial) 8.97 9.00 75.7 <2e-16 ***
# s(Trial):IsTenth 9.97 10.00 101.4 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Rank: 38/39
# R-sq.(adj) = 0.329 Deviance explained = 33.1%
# fREML = 15622 Scale est. = 0.65564 n = 12839summary(bam(Pos ~ s(Trial) + s(Trial,by=IsTenth) + s(Time) + s(Time,by=IsTenth), data=dat))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ s(Trial) + s(Trial, by = IsTenth) + s(Time) + s(Time, by = IsTenth)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.0253 0.0261 -0.97 0.33
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Trial) 8.97 9.00 75.7 <2e-16 ***
# s(Trial):IsTenth 8.97 9.00 109.7 <2e-16 ***
# s(Time) 7.76 8.53 31.4 <2e-16 ***
# s(Time):IsTenth 9.08 9.70 156.8 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Rank: 38/39
# R-sq.(adj) = 0.329 Deviance explained = 33.1%
# fREML = 15622 Scale est. = 0.65564 n = 12839Note that the edf value for the binary smooths will differ by exactly 1 comparing these two models. The reason for this problematic order effect is that the model is not able to determine which of the binary difference smooths will include the constant difference between the two words.
6.4.3 Ordered factor difference
dat$WordO <- as.ordered(dat$Word)
contrasts(dat$WordO) <- "contr.treatment"
m2.ord <- bam(Pos ~ s(Time) + s(Time,by=WordO) + WordO, data=dat)
(smrym2.ord <- summary(m2.ord))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ s(Time) + s(Time, by = WordO) + WordO
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 0.0651 0.0107 6.1 1.1e-09 ***
# WordOtenth 0.6601 0.0149 44.2 < 2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time) 7.69 8.48 28.8 <2e-16 ***
# s(Time):WordOtenth 8.02 8.66 99.8 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.267 Deviance explained = 26.8%
# fREML = 16111 Scale est. = 0.71583 n = 12839plot(m2.ord, select=2, shade=TRUE, ylim=c(-1,2))
abline(h=0)
6.5 Model criticism
Model criticism shows the residuals are fine.
gam.check(m2)
#
# Method: fREML Optimizer: perf newton
# full convergence after 9 iterations.
# Gradient range [-4.63e-07,3.87e-07]
# (score 16112 & scale 0.716).
# Hessian positive definite, eigenvalue range [2.95,6418].
# Model rank = 20 / 20
#
# 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):Wordtent 9.00 7.52 1.02 0.91
# s(Time):Wordtenth 9.00 8.55 1.02 0.926.6 Mixed-effects
6.6.1 Random intercepts
m3 <- bam(Pos ~ Word + s(Time, by=Word) + s(Speaker,bs="re"), data=dat)
(smrym3 <- summary(m3))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ Word + s(Time, by = Word) + s(Speaker, bs = "re")
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 0.0914 0.0680 1.35 0.18
# Wordtenth 0.6782 0.0134 50.79 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):Wordtent 7.76 8.61 36.3 <2e-16 ***
# s(Time):Wordtenth 8.64 8.96 352.8 <2e-16 ***
# s(Speaker) 40.58 41.00 86.9 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.427 Deviance explained = 42.9%
# fREML = 14634 Scale est. = 0.56012 n = 12839par(mfrow=c(1,2))
plot_smooth(m3, view="Time", plot_all="Word", main="m3", rug=FALSE, rm.ranef=T, ylim=c(-1,2),
col=two.colors)# Summary:
# * Word : factor; set to the value(s): tent, tenth.
# * Time : numeric predictor; with 30 values ranging from 0.000000 to 1.000000.
# * Speaker : factor; set to the value(s): VENI_NL_5. (Might be canceled as random effect, check below.)
# * NOTE : The following random effects columns are canceled: s(Speaker)
# plot_diff(m3, view="Time", comp=list(Word=c("tenth","tent")), rm.ranef=T, ylim=c(-1,2))# Summary:
# * Time : numeric predictor; with 100 values ranging from 0.000000 to 1.000000.
# * Speaker : factor; set to the value(s): VENI_NL_5. (Might be canceled as random effect, check below.)
# * NOTE : The following random effects columns are canceled: s(Speaker)
# 
#
# Time window(s) of significant difference(s):
# 0.050505 - 1.0000006.6.2 Random slopes
m4 <- bam(Pos ~ Word + s(Time, by=Word) + s(Speaker,bs="re") + s(Speaker,Word,bs="re"), data=dat)
(smrym4 <- summary(m4))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ Word + s(Time, by = Word) + s(Speaker, bs = "re") + s(Speaker,
# Word, bs = "re")
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 0.1087 0.0828 1.31 0.19
# Wordtenth 0.6179 0.1032 5.99 2.2e-09 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):Wordtent 7.95 8.71 44.6 < 2e-16 ***
# s(Time):Wordtenth 8.70 8.97 433.0 < 2e-16 ***
# s(Speaker) 15.48 41.00 1080.0 0.11
# s(Speaker,Word) 64.59 81.00 960.4 3.3e-05 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.534 Deviance explained = 53.7%
# fREML = 13397 Scale est. = 0.45546 n = 12839compareML(m3,m4)# m3: Pos ~ Word + s(Time, by = Word) + s(Speaker, bs = "re")
#
# m4: Pos ~ Word + s(Time, by = Word) + s(Speaker, bs = "re") + s(Speaker,
# Word, bs = "re")# Warning in sprintf("***", h1): one argument not used by format '***'#
# Chi-square test of fREML scores
# -----
# Model Score Edf Difference Df p.value Sig.
# 1 m3 14634 7
# 2 m4 13397 8 1237.194 1.000 < 2e-16 ***
#
# AIC difference: 2616.08, model m4 has lower AIC.par(mfrow=c(1,2))
plot_smooth(m4, view="Time", plot_all="Word", main="m4", rug=FALSE, rm.ranef=T, ylim=c(-1,2),
col=two.colors) # Summary:
# * Word : factor; set to the value(s): tent, tenth. (Might be canceled as random effect, check below.)
# * Time : numeric predictor; with 30 values ranging from 0.000000 to 1.000000.
# * Speaker : factor; set to the value(s): VENI_NL_5. (Might be canceled as random effect, check below.)
# * NOTE : The following random effects columns are canceled: s(Speaker),s(Speaker,Word)
# plot_diff(m4, view="Time", comp=list(Word=c("tenth","tent")), rm.ranef=T, ylim=c(-1,2))# Summary:
# * Time : numeric predictor; with 100 values ranging from 0.000000 to 1.000000.
# * Speaker : factor; set to the value(s): VENI_NL_5. (Might be canceled as random effect, check below.)
# * NOTE : The following random effects columns are canceled: s(Speaker),s(Speaker,Word)
# 
#
# Time window(s) of significant difference(s):
# 0.151515 - 0.282828
# 0.393939 - 1.0000006.6.3 Factor smooths
m5 <- bam(Pos ~ Word + s(Time, by=Word) + s(Speaker,Word,bs="re") + s(Time,Speaker,bs="fs",m=1),
data=dat)
(smrym5 <- summary(m5))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ Word + s(Time, by = Word) + s(Speaker, Word, bs = "re") +
# s(Time, Speaker, bs = "fs", m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 0.1090 0.0829 1.31 0.19
# Wordtenth 0.6176 0.1032 5.99 2.2e-09 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):Wordtent 7.41 7.97 9.47 <2e-16 ***
# s(Time):Wordtenth 8.57 8.80 43.72 <2e-16 ***
# s(Speaker,Word) 64.79 81.00 52.32 <2e-16 ***
# s(Time,Speaker) 296.36 377.00 249.79 0.075 .
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.67 Deviance explained = 68%
# fREML = 11596 Scale est. = 0.32251 n = 12839plot(m5, select=4)
compareML(m4,m5)# m4: Pos ~ Word + s(Time, by = Word) + s(Speaker, bs = "re") + s(Speaker,
# Word, bs = "re")
#
# m5: Pos ~ Word + s(Time, by = Word) + s(Speaker, Word, bs = "re") +
# s(Time, Speaker, bs = "fs", m = 1)# Warning in sprintf("***", h1): one argument not used by format '***'#
# Chi-square test of fREML scores
# -----
# Model Score Edf Difference Df p.value Sig.
# 1 m4 13397 8
# 2 m5 11596 9 1801.628 1.000 < 2e-16 ***
#
# AIC difference: 4156.35, model m5 has lower AIC.par(mfrow=c(1,2))
plot_smooth(m5, view="Time", plot_all="Word", main="m5", rug=FALSE, rm.ranef=T, ylim=c(-1,2),
col=two.colors) # Summary:
# * Word : factor; set to the value(s): tent, tenth. (Might be canceled as random effect, check below.)
# * Time : numeric predictor; with 30 values ranging from 0.000000 to 1.000000.
# * Speaker : factor; set to the value(s): VENI_NL_5. (Might be canceled as random effect, check below.)
# * NOTE : The following random effects columns are canceled: s(Speaker,Word),s(Time,Speaker)
# plot_diff(m5, view="Time", comp=list(Word=c("tenth","tent")), rm.ranef=T, ylim=c(-1,2))# Summary:
# * Time : numeric predictor; with 100 values ranging from 0.000000 to 1.000000.
# * Speaker : factor; set to the value(s): VENI_NL_5. (Might be canceled as random effect, check below.)
# * NOTE : The following random effects columns are canceled: s(Speaker,Word),s(Time,Speaker)
# 
#
# Time window(s) of significant difference(s):
# 0.141414 - 0.292929
# 0.393939 - 1.000000m6 <- bam(Pos ~ Word + s(Time, by=Word) + s(Time,Speaker,by=Word,bs="fs",m=1), data=dat)# Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated 1-d smooths of same
# variable.(smrym6 <- summary(m6))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ Word + s(Time, by = Word) + s(Time, Speaker, by = Word,
# bs = "fs", m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 0.109 0.094 1.16 0.25
# Wordtenth 0.612 0.117 5.22 1.9e-07 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):Wordtent 7.54 7.95 9.31 <2e-16 ***
# s(Time):Wordtenth 8.38 8.54 22.13 <2e-16 ***
# s(Time,Speaker):Wordtent 316.83 377.00 38.13 <2e-16 ***
# s(Time,Speaker):Wordtenth 328.72 368.00 43.27 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.782 Deviance explained = 79.3%
# fREML = 9386.6 Scale est. = 0.21284 n = 12839compareML(m5,m6)# m5: Pos ~ Word + s(Time, by = Word) + s(Speaker, Word, bs = "re") +
# s(Time, Speaker, bs = "fs", m = 1)
#
# m6: Pos ~ Word + s(Time, by = Word) + s(Time, Speaker, by = Word,
# bs = "fs", m = 1)# Warning in sprintf("***", h1): one argument not used by format '***'#
# Chi-square test of fREML scores
# -----
# Model Score Edf Difference Df p.value Sig.
# 1 m5 11596 9
# 2 m6 9387 10 2209.002 1.000 < 2e-16 ***
#
# AIC difference: 5063.34, model m6 has lower AIC.par(mfrow=c(1,2))
plot_smooth(m6, view="Time", plot_all="Word", main="m6", rug=FALSE, rm.ranef=T, ylim=c(-1,2),
col=two.colors) # Summary:
# * Word : factor; set to the value(s): tent, tenth.
# * Time : numeric predictor; with 30 values ranging from 0.000000 to 1.000000.
# * Speaker : factor; set to the value(s): VENI_NL_5. (Might be canceled as random effect, check below.)
# * NOTE : The following random effects columns are canceled: s(Time,Speaker):Wordtent,s(Time,Speaker):Wordtenth
# plot_diff(m6, view="Time", comp=list(Word=c("tenth","tent")), rm.ranef=T, ylim=c(-1,2))# Summary:
# * Time : numeric predictor; with 100 values ranging from 0.000000 to 1.000000.
# * Speaker : factor; set to the value(s): VENI_NL_5. (Might be canceled as random effect, check below.)
# * NOTE : The following random effects columns are canceled: s(Time,Speaker):Wordtent,s(Time,Speaker):Wordtenth
# 
#
# Time window(s) of significant difference(s):
# 0.404040 - 1.0000006.7 Autocorrelation
m6acf <- acf_resid(m6)
(rhoval <- m6acf[2]) # autocorrelation at lag 1# 1
# 0.912m7 <- bam(Pos ~ Word + s(Time, by=Word) + s(Time,Speaker,by=Word,bs="fs",m=1), data=dat,
rho=rhoval, AR.start=dat$start.event) # Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated 1-d smooths of same
# variable.acf_resid(m7)
summary(m7)#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ Word + s(Time, by = Word) + s(Time, Speaker, by = Word,
# bs = "fs", m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 0.1093 0.0896 1.22 0.22
# Wordtenth 0.6137 0.1126 5.45 5.1e-08 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):Wordtent 7.4 8.07 8.88 <2e-16 ***
# s(Time):Wordtenth 8.3 8.58 21.29 <2e-16 ***
# s(Time,Speaker):Wordtent 231.2 377.00 2.86 <2e-16 ***
# s(Time,Speaker):Wordtenth 271.1 368.00 3.86 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.762 Deviance explained = 77.2%
# fREML = -3293.3 Scale est. = 0.18777 n = 12839par(mfrow=c(1,2))
plot_smooth(m7, view="Time", plot_all="Word", main="m7", rug=FALSE, rm.ranef=T, ylim=c(-1,2),
col=two.colors) # Summary:
# * Word : factor; set to the value(s): tent, tenth.
# * Time : numeric predictor; with 30 values ranging from 0.000000 to 1.000000.
# * Speaker : factor; set to the value(s): VENI_NL_5. (Might be canceled as random effect, check below.)
# * NOTE : The following random effects columns are canceled: s(Time,Speaker):Wordtent,s(Time,Speaker):Wordtenth
# plot_diff(m7, view="Time", comp=list(Word=c("tenth","tent")), rm.ranef=T, ylim=c(-1,2))# Summary:
# * Time : numeric predictor; with 100 values ranging from 0.000000 to 1.000000.
# * Speaker : factor; set to the value(s): VENI_NL_5. (Might be canceled as random effect, check below.)
# * NOTE : The following random effects columns are canceled: s(Time,Speaker):Wordtent,s(Time,Speaker):Wordtenth
# 
#
# Time window(s) of significant difference(s):
# 0.404040 - 1.000000m7.alt <- bam(Pos ~ Word + s(Time, by=Word) + s(Time,Speaker,by=Word,bs="fs",m=1), data=dat,
rho=rhoval - 0.1, AR.start=dat$start.event) # Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated 1-d smooths of same
# variable.acf_resid(m7.alt)
m7.alt2 <- bam(Pos ~ Word + s(Time, by=Word) + s(Time,Speaker,by=Word,bs="fs",m=1), data=dat,
rho=0.99, AR.start=dat$start.event)
acf_resid(m7.alt2)par(mfrow=c(1,2))
plot_smooth(m7.alt2, view="Time", plot_all="Word", main="m7.alt2", rug=FALSE, rm.ranef=T, ylim=c(-1,2),
col=two.colors)
plot_diff(m7.alt2, view="Time", comp=list(Word=c("tenth","tent")), rm.ranef=T, ylim=c(-1,2))6.8 Tensor product interaction
m8 <- bam(Pos ~ Word + te(Time, Trial, by=Word) + s(Time,Speaker,by=Word,bs="fs",m=1), data=dat,
rho=rhoval, AR.start=dat$start.event)
(smrym8 <- summary(m8))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ Word + te(Time, Trial, by = Word) + s(Time, Speaker, by = Word,
# bs = "fs", m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 0.0895 0.0937 0.96 0.34
# Wordtenth 0.6636 0.1184 5.60 2.2e-08 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# te(Time,Trial):Wordtent 9.09 9.92 8.07 <2e-16 ***
# te(Time,Trial):Wordtenth 8.55 8.77 15.61 <2e-16 ***
# s(Time,Speaker):Wordtent 232.93 377.00 2.88 <2e-16 ***
# s(Time,Speaker):Wordtenth 282.68 368.00 4.15 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.776 Deviance explained = 78.5%
# fREML = -3290 Scale est. = 0.18737 n = 12839par(mfrow=c(2,2))
fvisgam(m8, view=c("Time","Trial"), cond=list(Word=c("tent")), main="m8: tent", rm.ranef=T,
zlim=c(-0.9,1.6), ylim=c(0,600), color="gray") # Summary:
# * Word : factor; set to the value(s): tent.
# * Time : numeric predictor; with 30 values ranging from 0.000000 to 1.000000.
# * Trial : numeric predictor; with 30 values ranging from 0.000000 to 600.000000.
# * Speaker : factor; set to the value(s): VENI_NL_5. (Might be canceled as random effect, check below.)
# * NOTE : The following random effects columns are canceled: s(Time,Speaker):Wordtent,s(Time,Speaker):Wordtenth
# # Warning in gradientLegend(c(min.z, max.z), n.seg = 3, pos = 0.875, color = pal, : Increase right
# margin to fit labels or decrease the number of decimals, see help(gradientLegend).fvisgam(m8, view=c("Time","Trial"), cond=list(Word=c("tenth")), main="m8: tenth", rm.ranef=T,
zlim=c(-0.9,1.6), ylim=c(0,600), color="gray") # Summary:
# * Word : factor; set to the value(s): tenth.
# * Time : numeric predictor; with 30 values ranging from 0.000000 to 1.000000.
# * Trial : numeric predictor; with 30 values ranging from 0.000000 to 600.000000.
# * Speaker : factor; set to the value(s): VENI_NL_5. (Might be canceled as random effect, check below.)
# * NOTE : The following random effects columns are canceled: s(Time,Speaker):Wordtent,s(Time,Speaker):Wordtenth
# # Warning in gradientLegend(c(min.z, max.z), n.seg = 3, pos = 0.875, color = pal, : Increase right
# margin to fit labels or decrease the number of decimals, see help(gradientLegend).plot_diff2(m8, view=c("Time","Trial"), comp=list(Word=c("tenth","tent")), rm.ranef=T, se=0,
main="Difference tenth - tent", zlim=c(-0.1,1.2), ylim=c(0,600), color="gray")# Summary:
# * Time : numeric predictor; with 30 values ranging from 0.000000 to 1.000000.
# * Trial : numeric predictor; with 30 values ranging from 0.000000 to 600.000000.
# * Speaker : factor; set to the value(s): VENI_NL_5. (Might be canceled as random effect, check below.)
# * NOTE : The following random effects columns are canceled: s(Time,Speaker):Wordtent,s(Time,Speaker):Wordtenth
# # Warning in gradientLegend(zlim, n.seg = 3, pos = 0.875, dec = dec, color = color, : Increase right
# margin to fit labels or decrease the number of decimals, see help(gradientLegend).
par(mfcol=c(3,2))
plot_diff2(m8, view=c("Time","Trial"), comp=list(Word=c("tenth","tent")), rm.ranef=T, se=0,
main="Difference tenth - tent", zlim=c(-0.1,1.2), ylim=c(0,600), color="gray",
cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5) # Summary:
# * Time : numeric predictor; with 30 values ranging from 0.000000 to 1.000000.
# * Trial : numeric predictor; with 30 values ranging from 0.000000 to 600.000000.
# * Speaker : factor; set to the value(s): VENI_NL_5. (Might be canceled as random effect, check below.)
# * NOTE : The following random effects columns are canceled: s(Time,Speaker):Wordtent,s(Time,Speaker):Wordtenth
# # Warning in gradientLegend(zlim, n.seg = 3, pos = 0.875, dec = dec, color = color, : Increase right
# margin to fit labels or decrease the number of decimals, see help(gradientLegend).abline(h=500,lty=2,lwd=2,col="white")
plot_diff2(m8, view=c("Time","Trial"), comp=list(Word=c("tenth","tent")), rm.ranef=T, se=0,
main="Difference tenth - tent", zlim=c(-0.1,1.2), ylim=c(0,600), color="gray",
cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5) # Summary:
# * Time : numeric predictor; with 30 values ranging from 0.000000 to 1.000000.
# * Trial : numeric predictor; with 30 values ranging from 0.000000 to 600.000000.
# * Speaker : factor; set to the value(s): VENI_NL_5. (Might be canceled as random effect, check below.)
# * NOTE : The following random effects columns are canceled: s(Time,Speaker):Wordtent,s(Time,Speaker):Wordtenth
# # Warning in gradientLegend(zlim, n.seg = 3, pos = 0.875, dec = dec, color = color, : Increase right
# margin to fit labels or decrease the number of decimals, see help(gradientLegend).abline(h=300,lty=2,lwd=2,col="white")
plot_diff2(m8, view=c("Time","Trial"), comp=list(Word=c("tenth","tent")), rm.ranef=T, se=0,
main="Difference tenth - tent", zlim=c(-0.1,1.2), ylim=c(0,600), color="gray",
cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5) # Summary:
# * Time : numeric predictor; with 30 values ranging from 0.000000 to 1.000000.
# * Trial : numeric predictor; with 30 values ranging from 0.000000 to 600.000000.
# * Speaker : factor; set to the value(s): VENI_NL_5. (Might be canceled as random effect, check below.)
# * NOTE : The following random effects columns are canceled: s(Time,Speaker):Wordtent,s(Time,Speaker):Wordtenth
# # Warning in gradientLegend(zlim, n.seg = 3, pos = 0.875, dec = dec, color = color, : Increase right
# margin to fit labels or decrease the number of decimals, see help(gradientLegend).abline(h=100,lty=2,lwd=2,col="white")
plot_diff(m8, view="Time", comp=list(Word=c("tenth","tent")), cond=list(Trial=500), rm.ranef=T,
main="Difference for trial 500", ylim=c(-1,2), cex.lab=1.5, cex.axis=1.5,
cex.main=1.5, cex.sub=1.5) # Summary:
# * Time : numeric predictor; with 100 values ranging from 0.000000 to 1.000000.
# * Trial : numeric predictor; set to the value(s): 500.
# * Speaker : factor; set to the value(s): VENI_NL_5. (Might be canceled as random effect, check below.)
# * NOTE : The following random effects columns are canceled: s(Time,Speaker):Wordtent,s(Time,Speaker):Wordtenth
# #
# Time window(s) of significant difference(s):
# 0.333333 - 1.000000abline(h=1.035,lty=3)
abline(v=0.745,lty=3)
plot_diff(m8, view="Time", comp=list(Word=c("tenth","tent")), cond=list(Trial=300), rm.ranef=T,
main="Difference for trial 300", ylim=c(-1,2), cex.lab=1.5, cex.axis=1.5,
cex.main=1.5, cex.sub=1.5)# Summary:
# * Time : numeric predictor; with 100 values ranging from 0.000000 to 1.000000.
# * Trial : numeric predictor; set to the value(s): 300.
# * Speaker : factor; set to the value(s): VENI_NL_5. (Might be canceled as random effect, check below.)
# * NOTE : The following random effects columns are canceled: s(Time,Speaker):Wordtent,s(Time,Speaker):Wordtenth
# #
# Time window(s) of significant difference(s):
# 0.272727 - 1.000000abline(h=1.035,lty=3)
abline(v=0.745,lty=3)
plot_diff(m8, view="Time", comp=list(Word=c("tenth","tent")), cond=list(Trial=100), rm.ranef=T,
main="Difference for trial 100", ylim=c(-1,2), cex.lab=1.5, cex.axis=1.5,
cex.main=1.5, cex.sub=1.5)# Summary:
# * Time : numeric predictor; with 100 values ranging from 0.000000 to 1.000000.
# * Trial : numeric predictor; set to the value(s): 100.
# * Speaker : factor; set to the value(s): VENI_NL_5. (Might be canceled as random effect, check below.)
# * NOTE : The following random effects columns are canceled: s(Time,Speaker):Wordtent,s(Time,Speaker):Wordtenth
# #
# Time window(s) of significant difference(s):
# 0.282828 - 1.000000abline(h=1.035,lty=3)
abline(v=0.745,lty=3)
6.9 Decomposition of tensor
m8.dc <- bam(Pos ~ Word + s(Time, by=Word) + s(Trial, by=Word) + ti(Time, Trial, by=Word) +
s(Time,Speaker,by=Word,bs="fs",m=1), data=dat, rho=rhoval,
AR.start=dat$start.event)
(smrym8.dc <- summary(m8.dc))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ Word + s(Time, by = Word) + s(Trial, by = Word) + ti(Time,
# Trial, by = Word) + s(Time, Speaker, by = Word, bs = "fs",
# m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 0.0877 0.0937 0.94 0.35
# Wordtenth 0.6636 0.1181 5.62 2e-08 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):Wordtent 7.40 8.06 9.02 < 2e-16 ***
# s(Time):Wordtenth 8.29 8.57 21.47 < 2e-16 ***
# s(Trial):Wordtent 1.00 1.00 14.15 0.00017 ***
# s(Trial):Wordtenth 1.68 1.97 7.27 0.00048 ***
# ti(Time,Trial):Wordtent 1.00 1.00 3.50 0.06142 .
# ti(Time,Trial):Wordtenth 1.13 1.23 3.54 0.06351 .
# s(Time,Speaker):Wordtent 230.98 377.00 2.88 < 2e-16 ***
# s(Time,Speaker):Wordtenth 272.15 368.00 3.94 < 2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.777 Deviance explained = 78.6%
# fREML = -3307 Scale est. = 0.18707 n = 12839par(mfrow=c(2,2))
plot(m8.dc, select=1, shade=T, rug=F, cex.lab=1.35, cex.axis=1.35, cex.main=1.35, cex.sub=1.35,
ylim=c(-1.5,1.5))
abline(h=0)
plot(m8.dc, select=2, shade=T, rug=F, cex.lab=1.35, cex.axis=1.35, cex.main=1.35, cex.sub=1.35,
ylim=c(-1.5,1.5))
abline(h=0)
plot(m8.dc, select=3, shade=T, rug=F, cex.lab=1.35, cex.axis=1.35, cex.main=1.35, cex.sub=1.35,
ylim=c(-1.5,1.5))
abline(h=0)
plot(m8.dc, select=4, shade=T, rug=F, cex.lab=1.35, cex.axis=1.35, cex.main=1.35, cex.sub=1.35,
ylim=c(-1.5,1.5))
abline(h=0)
#plot(m8.dc, select=5, scheme=3, cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5,
# contour.col='black', zlim=c(-0.2,0.2)) # non-significant: excluded
#plot(m8.dc, select=6, scheme=3, cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5,
# contour.col='black', zlim=c(-0.2,0.2)) # non-significant: excluded6.10 Including language difference
dat$WordLang <- interaction(dat$Word, dat$Lang)
m9 <- bam(Pos ~ WordLang + s(Time, by=WordLang) + s(Time,Speaker,by=Word,bs="fs",m=1), data=dat,
rho=rhoval, AR.start=dat$start.event)
(smrym9 <- summary(m9))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ WordLang + s(Time, by = WordLang) + s(Time, Speaker, by = Word,
# bs = "fs", m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.0686 0.1197 -0.57 0.566
# WordLangtenth.EN 0.7315 0.1515 4.83 1.4e-06 ***
# WordLangtent.NL 0.3755 0.1739 2.16 0.031 *
# WordLangtenth.NL 0.8636 0.1566 5.51 3.6e-08 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):WordLangtent.EN 3.86 4.56 2.47 0.042 *
# s(Time):WordLangtenth.EN 7.93 8.34 15.48 <2e-16 ***
# s(Time):WordLangtent.NL 7.45 8.12 10.74 <2e-16 ***
# s(Time):WordLangtenth.NL 7.69 8.16 11.56 <2e-16 ***
# s(Time,Speaker):Wordtent 219.72 376.00 2.65 <2e-16 ***
# s(Time,Speaker):Wordtenth 258.63 367.00 3.51 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.762 Deviance explained = 77.1%
# fREML = -3305.3 Scale est. = 0.1878 n = 12839par(mfrow=c(2,2))
plot_smooth(m9, view="Time", cond=list(WordLang=c("tent.EN")), main="m9: English speakers",
rug=FALSE, rm.ranef=T, ylim=c(-1,2), col=two.colors[1]) # Summary:
# * WordLang : factor; set to the value(s): tent.EN.
# * Time : numeric predictor; with 30 values ranging from 0.000000 to 1.000000.
# * Speaker : factor; set to the value(s): VENI_NL_5. (Might be canceled as random effect, check below.)
# * Word : factor; set to the value(s): tenth.
# * NOTE : The following random effects columns are canceled: s(Time,Speaker):Wordtent,s(Time,Speaker):Wordtenth
# plot_smooth(m9, view="Time", cond=list(WordLang=c("tenth.EN")), rm.ranef=T, add=T,
col=two.colors[2], rug=FALSE) # Summary:
# * WordLang : factor; set to the value(s): tenth.EN.
# * Time : numeric predictor; with 30 values ranging from 0.000000 to 1.000000.
# * Speaker : factor; set to the value(s): VENI_NL_5. (Might be canceled as random effect, check below.)
# * Word : factor; set to the value(s): tenth.
# * NOTE : The following random effects columns are canceled: s(Time,Speaker):Wordtent,s(Time,Speaker):Wordtenth
# # add legend (from code of plot_smooth with parameter plot_all)
gfc <- getFigCoords()
legend(gfc[2], gfc[4], legend = c("tent","tenth"), text.col = two.colors, text.font = 2,
xjust = 1, yjust = 1, bty = "n", xpd = TRUE)
plot_diff(m9, view="Time", comp=list(WordLang=c("tenth.EN","tent.EN")), rm.ranef=T,
ylim=c(-1,2), main="Difference tenth - tent")# Summary:
# * Time : numeric predictor; with 100 values ranging from 0.000000 to 1.000000.
# * Speaker : factor; set to the value(s): VENI_NL_5. (Might be canceled as random effect, check below.)
# * Word : factor; set to the value(s): tenth.
# * NOTE : The following random effects columns are canceled: s(Time,Speaker):Wordtent,s(Time,Speaker):Wordtenth
# #
# Time window(s) of significant difference(s):
# 0.424242 - 1.000000plot_smooth(m9, view="Time", cond=list(WordLang=c("tent.NL")), main="m9: Dutch speakers",
rug=FALSE, rm.ranef=T, ylim=c(-1,2), col=two.colors[1]) # Summary:
# * WordLang : factor; set to the value(s): tent.NL.
# * Time : numeric predictor; with 30 values ranging from 0.000000 to 1.000000.
# * Speaker : factor; set to the value(s): VENI_NL_5. (Might be canceled as random effect, check below.)
# * Word : factor; set to the value(s): tenth.
# * NOTE : The following random effects columns are canceled: s(Time,Speaker):Wordtent,s(Time,Speaker):Wordtenth
# plot_smooth(m9, view="Time", cond=list(WordLang=c("tenth.NL")), rm.ranef=T, add=T,
col=two.colors[2], rug=FALSE) # Summary:
# * WordLang : factor; set to the value(s): tenth.NL.
# * Time : numeric predictor; with 30 values ranging from 0.000000 to 1.000000.
# * Speaker : factor; set to the value(s): VENI_NL_5. (Might be canceled as random effect, check below.)
# * Word : factor; set to the value(s): tenth.
# * NOTE : The following random effects columns are canceled: s(Time,Speaker):Wordtent,s(Time,Speaker):Wordtenth
# # add legend
gfc <- getFigCoords()
legend(gfc[2], gfc[4], legend = c("tent","tenth"), text.col = two.colors, text.font = 2,
xjust = 1, yjust = 1, bty = "n", xpd = TRUE)
plot_diff(m9, view="Time", comp=list(WordLang=c("tenth.NL","tent.NL")), rm.ranef=T,
ylim=c(-1,2), main="Difference tenth - tent")# Summary:
# * Time : numeric predictor; with 100 values ranging from 0.000000 to 1.000000.
# * Speaker : factor; set to the value(s): VENI_NL_5. (Might be canceled as random effect, check below.)
# * Word : factor; set to the value(s): tenth.
# * NOTE : The following random effects columns are canceled: s(Time,Speaker):Wordtent,s(Time,Speaker):Wordtenth
# 
#
# Time window(s) of significant difference(s):
# 0.222222 - 0.262626
# 0.434343 - 0.8686876.10.1 Model comparison
m7.ml <- bam(Pos ~ Word + s(Time, by=Word) + s(Time,Speaker,by=Word,bs="fs",m=1), data=dat,
rho=rhoval, AR.start=dat$start.event, method='ML') # Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated 1-d smooths of same
# variable.m9.ml <- bam(Pos ~ WordLang + s(Time, by=WordLang) + s(Time,Speaker,by=Word,bs="fs",m=1), data=dat,
rho=rhoval, AR.start=dat$start.event, method='ML')
compareML(m7.ml, m9.ml)# m7.ml: Pos ~ Word + s(Time, by = Word) + s(Time, Speaker, by = Word,
# bs = "fs", m = 1)
#
# m9.ml: Pos ~ WordLang + s(Time, by = WordLang) + s(Time, Speaker, by = Word,
# bs = "fs", m = 1)# Warning in sprintf("***", h1): one argument not used by format '***'#
# Chi-square test of ML scores
# -----
# Model Score Edf Difference Df p.value Sig.
# 1 m7.ml -3299 10
# 2 m9.ml -3314 16 15.217 6.000 3.251e-05 ***
#
# AIC difference: 4.51, model m9.ml has lower AIC.6.10.2 Ordered factor
dat$ENTenthO <- as.ordered(dat$Lang == "EN" & dat$Word == "tenth")
contrasts(dat$ENTenthO) <- "contr.treatment"
dat$NLTenthO <- as.ordered(dat$Lang == "NL" & dat$Word == "tenth")
contrasts(dat$NLTenthO) <- "contr.treatment"
m9.ord <- bam(Pos ~ Lang + ENTenthO + NLTenthO + s(Time,by=Lang) + s(Time,by=ENTenthO) +
s(Time,by=NLTenthO) + s(Time,Speaker,by=Word,bs="fs",m=1), data=dat,
rho=rhoval, AR.start=dat$start.event)
(smrym9.ord <- summary(m9.ord))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ Lang + ENTenthO + NLTenthO + s(Time, by = Lang) + s(Time,
# by = ENTenthO) + s(Time, by = NLTenthO) + s(Time, Speaker,
# by = Word, bs = "fs", m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.0681 0.1197 -0.57 0.5692
# LangNL 0.3774 0.1738 2.17 0.0300 *
# ENTenthOTRUE 0.7312 0.1514 4.83 1.4e-06 ***
# NLTenthOTRUE 0.4840 0.1615 3.00 0.0027 **
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):LangEN 4.29 5.01 2.46 0.03 *
# s(Time):LangNL 8.04 8.47 15.96 <2e-16 ***
# s(Time):ENTenthOTRUE 7.63 8.08 8.61 <2e-16 ***
# s(Time):NLTenthOTRUE 2.11 2.42 0.79 0.55
# s(Time,Speaker):Wordtent 220.88 376.00 2.66 <2e-16 ***
# s(Time,Speaker):Wordtenth 261.90 367.00 3.55 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.762 Deviance explained = 77.1%
# fREML = -3312.1 Scale est. = 0.1877 n = 12839par(mfrow=c(1,2))
plot(m9.ord,select=3,shade=T,rug=F, ylim=c(-1.5,1.5), main="m9.ord: English difference")
abline(h=0)
plot(m9.ord,select=4,shade=T,rug=F, ylim=c(-1.5,1.5), main="m9.ord: Dutch difference")
abline(h=0)
6.10.2.1 Increasing k for NLTenthO smooth
m9.ord20 <- bam(Pos ~ Lang + ENTenthO + NLTenthO + s(Time,by=Lang) + s(Time,by=ENTenthO) +
s(Time,by=NLTenthO,k=20) + s(Time,Speaker,by=Word,bs="fs",m=1), data=dat,
rho=rhoval, AR.start=dat$start.event)
(smrym9.ord20 <- summary(m9.ord20))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ Lang + ENTenthO + NLTenthO + s(Time, by = Lang) + s(Time,
# by = ENTenthO) + s(Time, by = NLTenthO, k = 20) + s(Time,
# Speaker, by = Word, bs = "fs", m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.0681 0.1196 -0.57 0.5692
# LangNL 0.3754 0.1738 2.16 0.0308 *
# ENTenthOTRUE 0.7312 0.1514 4.83 1.4e-06 ***
# NLTenthOTRUE 0.4873 0.1616 3.02 0.0026 **
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):LangEN 4.29 5.02 2.46 0.03 *
# s(Time):LangNL 7.81 8.29 12.26 <2e-16 ***
# s(Time):ENTenthOTRUE 7.65 8.09 8.73 <2e-16 ***
# s(Time):NLTenthOTRUE 7.72 10.11 1.09 0.43
# s(Time,Speaker):Wordtent 220.24 376.00 2.65 <2e-16 ***
# s(Time,Speaker):Wordtenth 259.11 367.00 3.51 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.762 Deviance explained = 77.1%
# fREML = -3312.6 Scale est. = 0.18764 n = 12839m9.ord.ml <- bam(Pos ~ Lang + ENTenthO + NLTenthO + s(Time,by=Lang) + s(Time,by=ENTenthO) +
s(Time,by=NLTenthO) + s(Time,Speaker,by=Word,bs="fs",m=1), data=dat,
rho=rhoval, AR.start=dat$start.event, method='ML')
m9.ord20.ml <- bam(Pos ~ Lang + ENTenthO + NLTenthO + s(Time,by=Lang) + s(Time,by=ENTenthO) +
s(Time,by=NLTenthO,k=20) + s(Time,Speaker,by=Word,bs="fs",m=1), data=dat,
rho=rhoval, AR.start=dat$start.event, method='ML')
compareML(m9.ord.ml,m9.ord20.ml)# m9.ord.ml: Pos ~ Lang + ENTenthO + NLTenthO + s(Time, by = Lang) + s(Time,
# by = ENTenthO) + s(Time, by = NLTenthO) + s(Time, Speaker,
# by = Word, bs = "fs", m = 1)
#
# m9.ord20.ml: Pos ~ Lang + ENTenthO + NLTenthO + s(Time, by = Lang) + s(Time,
# by = ENTenthO) + s(Time, by = NLTenthO, k = 20) + s(Time,
# Speaker, by = Word, bs = "fs", m = 1)
#
# Model m9.ord.ml preferred: lower ML score (0.055), and equal df (0.000).
# -----
# Model Score Edf Difference Df
# 1 m9.ord20.ml -3321 16
# 2 m9.ord.ml -3321 16 0.055 0.000
#
# AIC difference: -5.14, model m9.ord.ml has lower AIC.# Warning in compareML(m9.ord.ml, m9.ord20.ml): Only small difference in ML...par(mfrow=c(1,2))
plot(m9.ord20,select=3,shade=T,rug=F, ylim=c(-1.5,1.5), main="m9.ord: English difference")
abline(h=0)
plot(m9.ord20,select=4,shade=T,rug=F, ylim=c(-1.5,1.5), main="m9.ord: Dutch difference")
abline(h=0)
6.10.2.2 Excluding the NLTenthO smooth
m9.ord2 <- bam(Pos ~ Lang + ENTenthO + NLTenthO + s(Time,by=Lang) + s(Time,by=ENTenthO) +
s(Time,Speaker,by=Word,bs="fs",m=1), data=dat,
rho=rhoval, AR.start=dat$start.event)
(smrym9.ord2 <- summary(m9.ord2))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ Lang + ENTenthO + NLTenthO + s(Time, by = Lang) + s(Time,
# by = ENTenthO) + s(Time, Speaker, by = Word, bs = "fs", m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.0681 0.1197 -0.57 0.5694
# LangNL 0.3809 0.1739 2.19 0.0285 *
# ENTenthOTRUE 0.7312 0.1515 4.83 1.4e-06 ***
# NLTenthOTRUE 0.4768 0.1614 2.95 0.0031 **
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):LangEN 4.29 5.02 2.46 0.03 *
# s(Time):LangNL 8.08 8.50 20.45 <2e-16 ***
# s(Time):ENTenthOTRUE 7.62 8.07 8.56 <2e-16 ***
# s(Time,Speaker):Wordtent 221.02 376.00 2.66 <2e-16 ***
# s(Time,Speaker):Wordtenth 263.60 367.00 3.58 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.762 Deviance explained = 77.1%
# fREML = -3312.9 Scale est. = 0.18768 n = 128396.10.3 Binary smooth
dat$IsENTenth <- (dat$Lang == "EN" & dat$Word == "tenth")*1
dat$IsNLTenth <- (dat$Lang == "NL" & dat$Word == "tenth")*1
m9.bin <- bam(Pos ~ Lang + s(Time,by=Lang) + s(Time,by=IsENTenth) + s(Time,by=IsNLTenth,k=20) +
s(Time,Speaker,by=Word,bs="fs",m=1), data=dat, rho=rhoval,
AR.start=dat$start.event)
(smrym9.bin <- summary(m9.bin))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ Lang + s(Time, by = Lang) + s(Time, by = IsENTenth) + s(Time,
# by = IsNLTenth, k = 20) + s(Time, Speaker, by = Word, bs = "fs",
# m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.0681 0.1196 -0.57 0.569
# LangNL 0.3751 0.1738 2.16 0.031 *
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):LangEN 4.29 5.02 2.46 0.030 *
# s(Time):LangNL 7.81 8.29 12.26 <2e-16 ***
# s(Time):IsENTenth 8.65 9.09 10.15 <2e-16 ***
# s(Time):IsNLTenth 8.72 11.11 1.78 0.065 .
# s(Time,Speaker):Wordtent 220.25 376.00 2.65 <2e-16 ***
# s(Time,Speaker):Wordtenth 259.11 367.00 3.51 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.762 Deviance explained = 77.1%
# fREML = -3312.6 Scale est. = 0.18764 n = 12839par(mfrow=c(2,2))
plot(m9.ord20,select=3,shade=T,rug=F, ylim=c(-1.5,1.5), main="m9.ord: English difference")
abline(h=0)
plot(m9.ord20,select=4,shade=T,rug=F, ylim=c(-1.5,1.5), main="m9.ord: Dutch difference")
abline(h=0)
plot(m9.bin,select=3,shade=T,rug=F, ylim=c(-1,1.999), main="m9.bin: English difference")
abline(h=0)
plot(m9.bin,select=4,shade=T,rug=F, ylim=c(-1,1.999), main="m9.bin: Dutch difference")
abline(h=0)
6.11 Differences different?
6.11.1 Binary smooth
dat$IsTenth <- (dat$Word == "tenth")*1
m9b.bin <- bam(Pos ~ Lang + s(Time,by=Lang) + s(Time,by=IsTenth) + s(Time,by=IsNLTenth) +
s(Time,Speaker,by=Word,bs="fs",m=1), data=dat, rho=rhoval,
AR.start=dat$start.event)
(smrym9b.bin <- summary(m9b.bin))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ Lang + s(Time, by = Lang) + s(Time, by = IsTenth) + s(Time,
# by = IsNLTenth) + s(Time, Speaker, by = Word, bs = "fs",
# m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.0671 0.1196 -0.56 0.575
# LangNL 0.3723 0.1737 2.14 0.032 *
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):LangEN 4.75 5.52 3.01 0.0075 **
# s(Time):LangNL 7.68 8.23 11.89 <2e-16 ***
# s(Time):IsTenth 8.91 9.27 10.32 <2e-16 ***
# s(Time):IsNLTenth 4.22 4.69 1.77 0.1386
# s(Time,Speaker):Wordtent 219.97 376.00 2.65 <2e-16 ***
# s(Time,Speaker):Wordtenth 260.48 367.00 3.53 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.762 Deviance explained = 77.1%
# fREML = -3313.5 Scale est. = 0.18767 n = 12839par(mfrow=c(1,2))
plot(m9b.bin,select=3,shade=T,rug=F, ylim=c(-1,2), main="m9b.bin: English difference")
abline(h=0)
plot(m9b.bin,select=4,shade=T,rug=F, ylim=c(-1.5,1.5), main="m9b.bin: Dutch vs. English difference")
abline(h=0)
6.11.1.1 Model criticism
gam.check(m9b.bin)
#
# Method: fREML Optimizer: perf newton
# full convergence after 14 iterations.
# Gradient range [-6.07e-06,2.19e-06]
# (score -3313 & scale 0.188).
# Hessian positive definite, eigenvalue range [0.568,6419].
# Model rank = 796 / 796
#
# 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):LangEN 9.00 4.75 1.02 0.94
# s(Time):LangNL 9.00 7.68 1.02 0.88
# s(Time):IsTenth 10.00 8.91 1.02 0.87
# s(Time):IsNLTenth 10.00 4.22 1.02 0.83
# s(Time,Speaker):Wordtent 378.00 219.97 1.02 0.89
# s(Time,Speaker):Wordtenth 378.00 260.48 1.02 0.926.11.2 Ordered factor
dat$IsTenthO <- as.ordered(dat$Word == "tenth")
contrasts(dat$IsTenthO) <- "contr.treatment"
dat$IsNLTenthO <- as.ordered(dat$Lang == "NL" & dat$Word == "tenth")
contrasts(dat$IsNLTenthO) <- "contr.treatment"
m9b.ord <- bam(Pos ~ Lang + s(Time,by=Lang) + s(Time,by=IsTenthO) + IsTenthO +
s(Time,by=IsNLTenthO) + IsNLTenthO +
s(Time,Speaker,by=Word,bs="fs",m=1), data=dat, rho=rhoval,
AR.start=dat$start.event)
(smrym9b.ord <- summary(m9b.ord))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ Lang + s(Time, by = Lang) + s(Time, by = IsTenthO) + IsTenthO +
# s(Time, by = IsNLTenthO) + IsNLTenthO + s(Time, Speaker,
# by = Word, bs = "fs", m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.0671 0.1196 -0.56 0.575
# LangNL 0.3727 0.1738 2.14 0.032 *
# IsTenthOTRUE 0.7291 0.1514 4.82 1.5e-06 ***
# IsNLTenthOTRUE -0.2377 0.2213 -1.07 0.283
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):LangEN 4.75 5.52 3.01 0.0075 **
# s(Time):LangNL 7.68 8.23 11.89 <2e-16 ***
# s(Time):IsTenthOTRUE 7.91 8.27 8.95 <2e-16 ***
# s(Time):IsNLTenthOTRUE 3.22 3.69 1.97 0.1184
# s(Time,Speaker):Wordtent 219.97 376.00 2.65 <2e-16 ***
# s(Time,Speaker):Wordtenth 260.49 367.00 3.53 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.762 Deviance explained = 77.1%
# fREML = -3313.4 Scale est. = 0.18767 n = 12839m9b.ord20 <- bam(Pos ~ Lang + s(Time,by=Lang) + s(Time,by=IsTenthO) + IsTenthO +
s(Time,by=IsNLTenthO,k=20) + IsNLTenthO +
s(Time,Speaker,by=Word,bs="fs",m=1), data=dat, rho=rhoval,
AR.start=dat$start.event)
(smrym9b.ord20 <- summary(m9b.ord20))par(mfrow=c(1,2))
plot(m9b.ord,select=3,shade=T,rug=F, ylim=c(-1.5,1.5), main="m9b.ord: English difference")
abline(h=0)
plot(m9b.ord,select=4,shade=T,rug=F, ylim=c(-1.5,1.5), main="m9b.ord: Dutch vs. English difference")
abline(h=0)
6.11.3 Model comparison
m9b.bin.ml <- bam(Pos ~ Lang + s(Time,by=Lang) + s(Time,by=IsTenth) + s(Time,by=IsNLTenth) +
s(Time,Speaker,by=Word,bs="fs",m=1), data=dat,
rho=rhoval, AR.start=dat$start.event, method='ML')
m9a.bin.ml <- bam(Pos ~ Lang + s(Time,by=Lang) + s(Time,by=IsTenth) +
s(Time,Speaker,by=Word,bs="fs",m=1), data=dat,
rho=rhoval, AR.start=dat$start.event, method='ML')
compareML(m9a.bin.ml, m9b.bin.ml)# m9a.bin.ml: Pos ~ Lang + s(Time, by = Lang) + s(Time, by = IsTenth) + s(Time,
# Speaker, by = Word, bs = "fs", m = 1)
#
# m9b.bin.ml: Pos ~ Lang + s(Time, by = Lang) + s(Time, by = IsTenth) + s(Time,
# by = IsNLTenth) + s(Time, Speaker, by = Word, bs = "fs",
# m = 1)# Warning in sprintf(" ", h1): one argument not used by format ' '#
# Chi-square test of ML scores
# -----
# Model Score Edf Difference Df p.value Sig.
# 1 m9a.bin.ml -3317 13
# 2 m9b.bin.ml -3320 16 3.162 3.000 0.097
#
# AIC difference: -35.61, model m9a.bin.ml has lower AIC.# Warning in compareML(m9a.bin.ml, m9b.bin.ml): Only small difference in ML...6.11.4 Full binary model
The following model shows that the difference between the Dutch and the English group for the word 'tent' is significant (s(Time,by=IsNL)).
dat$IsNL <- (dat$Lang == "NL")*1
m9c.bin <- bam(Pos ~ Lang + s(Time) + s(Time,by=IsNL) + s(Time,by=IsTenth) + s(Time,by=IsNLTenth) +
s(Time,Speaker,by=Word,bs="fs",m=1), data=dat, rho=rhoval,
AR.start=dat$start.event)
(smrym9c.bin <- summary(m9c.bin))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ Lang + s(Time) + s(Time, by = IsNL) + s(Time, by = IsTenth) +
# s(Time, by = IsNLTenth) + s(Time, Speaker, by = Word, bs = "fs",
# m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.0656 0.1196 -0.55 0.58
# LangNL 1.3837 0.3341 4.14 3.5e-05 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time) 5.93 6.69 3.38 0.0012 **
# s(Time):IsNL 6.82 7.47 5.98 <2e-16 ***
# s(Time):IsTenth 8.79 9.19 9.77 <2e-16 ***
# s(Time):IsNLTenth 4.24 4.70 1.70 0.1447
# s(Time,Speaker):Wordtent 218.91 376.00 2.63 <2e-16 ***
# s(Time,Speaker):Wordtenth 260.66 367.00 3.50 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Rank: 796/797
# R-sq.(adj) = 0.762 Deviance explained = 77.1%
# fREML = -3314.8 Scale est. = 0.18767 n = 128396.12 Differences different (at end)?
6.12.1 Binary smooth
dat50 <- dat[dat$Time >= 0.5,] # only use time points in second half of word
m9d.bin <- bam(Pos ~ Lang + s(Time,by=Lang) + s(Time,by=IsTenth) + s(Time,by=IsNLTenth) +
s(Time,Speaker,by=Word,bs="fs",m=1), data=dat50, rho=rhoval,
AR.start=dat50$start.event)
(smrym9d.bin <- summary(m9d.bin))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ Lang + s(Time, by = Lang) + s(Time, by = IsTenth) + s(Time,
# by = IsNLTenth) + s(Time, Speaker, by = Word, bs = "fs",
# m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.0534 0.1256 -0.43 0.671
# LangNL 0.3829 0.1824 2.10 0.036 *
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):LangEN 4.28 5.17 4.37 0.00053 ***
# s(Time):LangNL 5.08 6.03 9.24 < 2e-16 ***
# s(Time):IsTenth 6.61 7.48 10.37 < 2e-16 ***
# s(Time):IsNLTenth 3.98 4.53 2.42 0.03927 *
# s(Time,Speaker):Wordtent 174.73 376.00 3.30 < 2e-16 ***
# s(Time,Speaker):Wordtenth 192.73 367.00 4.75 < 2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.763 Deviance explained = 77.7%
# fREML = -2161.4 Scale est. = 0.15694 n = 6463par(mfrow=c(1,2))
plot(m9d.bin,select=3,shade=T,rug=F, ylim=c(-1,2), main="m9d.bin: English difference")
abline(h=0)
plot(m9d.bin,select=4,shade=T,rug=F, ylim=c(-1.5,1.5), main="m9d.bin: Dutch vs. English difference")
abline(h=0)
6.12.1.1 Model criticism
gam.check(m9d.bin)
#
# Method: fREML Optimizer: perf newton
# full convergence after 6 iterations.
# Gradient range [-2.68e-08,2.22e-08]
# (score -2161 & scale 0.157).
# Hessian positive definite, eigenvalue range [0.848,3231].
# Model rank = 796 / 796
#
# 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):LangEN 9.00 4.28 0.96 <2e-16 ***
# s(Time):LangNL 9.00 5.08 0.96 <2e-16 ***
# s(Time):IsTenth 10.00 6.61 0.96 <2e-16 ***
# s(Time):IsNLTenth 10.00 3.98 0.96 <2e-16 ***
# s(Time,Speaker):Wordtent 378.00 174.73 0.96 <2e-16 ***
# s(Time,Speaker):Wordtenth 378.00 192.73 0.96 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 16.12.2 Ordered factor
m9d.ord <- bam(Pos ~ Lang + s(Time,by=Lang) + s(Time,by=IsTenthO) + IsTenthO +
s(Time,by=IsNLTenthO) + IsNLTenthO +
s(Time,Speaker,by=Word,bs="fs",m=1), data=dat50, rho=rhoval,
AR.start=dat50$start.event)
(smrym9d.ord <- summary(m9d.ord))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ Lang + s(Time, by = Lang) + s(Time, by = IsTenthO) + IsTenthO +
# s(Time, by = IsNLTenthO) + IsNLTenthO + s(Time, Speaker,
# by = Word, bs = "fs", m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.0534 0.1256 -0.43 0.671
# LangNL 0.3829 0.1824 2.10 0.036 *
# IsTenthOTRUE 1.2421 0.1879 6.61 4.2e-11 ***
# IsNLTenthOTRUE -0.5674 0.2745 -2.07 0.039 *
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):LangEN 4.28 5.17 4.37 0.00053 ***
# s(Time):LangNL 5.08 6.03 9.24 < 2e-16 ***
# s(Time):IsTenthOTRUE 5.61 6.48 6.57 < 2e-16 ***
# s(Time):IsNLTenthOTRUE 2.98 3.53 2.25 0.06942 .
# s(Time,Speaker):Wordtent 174.73 376.00 3.30 < 2e-16 ***
# s(Time,Speaker):Wordtenth 192.73 367.00 4.75 < 2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.763 Deviance explained = 77.7%
# fREML = -2161.4 Scale est. = 0.15694 n = 6463m9d.ord20 <- bam(Pos ~ Lang + s(Time,by=Lang) + s(Time,by=IsTenthO) + IsTenthO +
s(Time,by=IsNLTenthO,k=20) + IsNLTenthO +
s(Time,Speaker,by=Word,bs="fs",m=1), data=dat50, rho=rhoval,
AR.start=da60t$start.event)
(smrym9d.ord20 <- summary(m9d.ord20))par(mfrow=c(1,2))
plot(m9d.ord,select=3,shade=T,rug=F, ylim=c(-1.5,1.5), main="m9d.ord: English difference")
abline(h=0)
plot(m9d.ord,select=4,shade=T,rug=F, ylim=c(-1.5,1.5), main="m9d.ord: Dutch vs. English difference")
abline(h=0)
6.12.3 Model comparison
m9d.bin.ml <- bam(Pos ~ Lang + s(Time,by=Lang) + s(Time,by=IsTenth) + s(Time,by=IsNLTenth) +
s(Time,Speaker,by=Word,bs="fs",m=1), data=dat50,
rho=rhoval, AR.start=dat50$start.event, method='ML')
m9d9.bin.ml <- bam(Pos ~ Lang + s(Time,by=Lang) + s(Time,by=IsTenth) +
s(Time,Speaker,by=Word,bs="fs",m=1), data=dat50,
rho=rhoval, AR.start=dat50$start.event, method='ML')
compareML(m9d.bin.ml, m9d9.bin.ml)# m9d.bin.ml: Pos ~ Lang + s(Time, by = Lang) + s(Time, by = IsTenth) + s(Time,
# by = IsNLTenth) + s(Time, Speaker, by = Word, bs = "fs",
# m = 1)
#
# m9d9.bin.ml: Pos ~ Lang + s(Time, by = Lang) + s(Time, by = IsTenth) + s(Time,
# Speaker, by = Word, bs = "fs", m = 1)# Warning in sprintf(" ", h1): one argument not used by format ' '#
# Chi-square test of ML scores
# -----
# Model Score Edf Difference Df p.value Sig.
# 1 m9d9.bin.ml -2168 13
# 2 m9d.bin.ml -2171 16 3.310 3.000 0.085
#
# AIC difference: -3.93, model m9d.bin.ml has lower AIC.# Warning in compareML(m9d.bin.ml, m9d9.bin.ml): Only small difference in ML...6.12.4 Full binary model
The following model shows that the difference between the Dutch and the English group for the word 'tent' is significant (s(Time,by=IsNL)).
m9e.bin <- bam(Pos ~ Lang + s(Time) + s(Time,by=IsNL) + s(Time,by=IsTenth) + s(Time,by=IsNLTenth) +
s(Time,Speaker,by=Word,bs="fs",m=1), data=dat50, rho=rhoval,
AR.start=dat50$start.event)
(smrym9e.bin <- summary(m9e.bin))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ Lang + s(Time) + s(Time, by = IsNL) + s(Time, by = IsTenth) +
# s(Time, by = IsNLTenth) + s(Time, Speaker, by = Word, bs = "fs",
# m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.047 0.126 -0.37 0.71
# LangNL 0.000 0.000 NaN NaN
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time) 4.93 5.83 5.91 8.7e-06 ***
# s(Time):IsNL 3.71 4.28 1.94 0.096 .
# s(Time):IsTenth 6.40 7.27 10.13 < 2e-16 ***
# s(Time):IsNLTenth 3.81 4.34 2.17 0.065 .
# s(Time,Speaker):Wordtent 175.08 376.00 3.30 < 2e-16 ***
# s(Time,Speaker):Wordtenth 192.86 367.00 4.75 < 2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Rank: 796/797
# R-sq.(adj) = 0.763 Deviance explained = 77.7%
# fREML = -2163.3 Scale est. = 0.15697 n = 64636.13 Faster computation
Note that the benefit of using more threads is greater for models which take longer to fit.
system.time(
m9b.bin <- bam(Pos ~ Lang + s(Time,by=Lang) + s(Time,by=IsTenth) + s(Time,by=IsNLTenth) +
s(Time,Speaker,by=Word,bs="fs",m=1), data=dat, rho=rhoval,
AR.start=dat$start.event)
)# user system elapsed
# 71.84 95.09 8.71system.time(
m9b.bin.discrete <- bam(Pos ~ Lang + s(Time,by=Lang) + s(Time,by=IsTenth) +
s(Time,by=IsNLTenth) + s(Time,Speaker,by=Word,bs="fs",m=1),
data=dat, rho=rhoval, AR.start=dat$start.event, discrete=TRUE)
)# user system elapsed
# 15.91 25.06 3.66system.time(
m9b.bin.discrete2 <- bam(Pos ~ Lang + s(Time,by=Lang) + s(Time,by=IsTenth) +
s(Time,by=IsNLTenth) + s(Time,Speaker,by=Word,bs="fs",m=1),
data=dat, rho=rhoval, AR.start=dat$start.event, discrete=TRUE,
nthreads=2)
)# user system elapsed
# 43.52 89.98 5.23(smrym9b.bin.discrete <- summary(m9b.bin.discrete))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ Lang + s(Time, by = Lang) + s(Time, by = IsTenth) + s(Time,
# by = IsNLTenth) + s(Time, Speaker, by = Word, bs = "fs",
# m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.0633 0.1201 -0.53 0.598
# LangNL 0.3820 0.1745 2.19 0.029 *
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):LangEN 4.76 5.54 2.98 0.0079 **
# s(Time):LangNL 7.71 8.26 11.96 <2e-16 ***
# s(Time):IsTenth 8.92 9.29 10.30 <2e-16 ***
# s(Time):IsNLTenth 4.29 4.77 1.82 0.1252
# s(Time,Speaker):Wordtent 217.88 377.00 2.61 <2e-16 ***
# s(Time,Speaker):Wordtenth 258.08 367.00 3.46 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.761 Deviance explained = 77%
# fREML = -3302.5 Scale est. = 0.18819 n = 12839par(mfrow=c(1,2))
plot(m9b.bin.discrete,select=3,shade=T,rug=F, ylim=c(-1,2),
main="m9b.bin.discrete: English difference", cex.main=0.9)
abline(h=0)
plot(m9b.bin.discrete,select=4,shade=T,rug=F, ylim=c(-1.5,1.5),
main="m9b.bin.discrete: Dutch vs. English difference", cex.main=0.9)
abline(h=0)
7 Full model: word-final contrast
fullfinal <- droplevels(full[full$Loc == "Final" & full$Time >= 0.5,])
fullfinal$IsTH = (fullfinal$Sound == "TH")*1
fullfinal$IsNLTH = (fullfinal$Sound == "TH" & fullfinal$Lang == "NL")*1
# track first observation per trial (to correct for autocorrelation later on)
# and sort per trajectory
fullfinal$start.event = (fullfinal$Time == 0.5)
fullfinal = fullfinal[order(fullfinal$Speaker,fullfinal$Trial,fullfinal$Time),]
#fullfinal <- start_event(fullfinal,event=c("Speaker","Trial")) # bug in itsadug 2.4
system.time(
ffmc0 <- bam(Pos ~ Lang + s(Time, by=Lang) + s(Time, by=IsTH) + s(Time, by=IsNLTH) +
s(Time,Speaker,by=Sound,bs="fs",m=1) + s(Time,Word,by=Lang,bs="fs",m=1),
data=fullfinal, discrete=TRUE)
)# Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated 1-d smooths of same
# variable.# user system elapsed
# 23.7 39.5 4.4ffmc0acf <- acf_resid(ffmc0)
(rhoval <- ffmc0acf[2]) # autocorrelation at lag 1# 1
# 0.951system.time(
ffmc1 <- bam(Pos ~ Lang + s(Time, by=Lang) + s(Time, by=IsTH) + s(Time, by=IsNLTH) +
s(Time,Speaker,by=Sound,bs="fs",m=1) + s(Time,Word,by=Lang,bs="fs",m=1),
data=fullfinal, discrete=TRUE, rho=rhoval, AR.start=fullfinal$start.event)
)# Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated 1-d smooths of same
# variable.# user system elapsed
# 39.5 61.8 10.9acf_resid(ffmc1) 
(smryffmc1 <- summary(ffmc1))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ Lang + s(Time, by = Lang) + s(Time, by = IsTH) + s(Time,
# by = IsNLTH) + s(Time, Speaker, by = Sound, bs = "fs", m = 1) +
# s(Time, Word, by = Lang, bs = "fs", m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.137 0.249 -0.55 0.58
# LangNL 0.175 0.452 0.39 0.70
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):LangEN 3.34 3.75 3.44 0.0117 *
# s(Time):LangNL 5.01 5.43 2.53 0.0224 *
# s(Time):IsTH 5.28 5.71 4.54 0.0002 ***
# s(Time):IsNLTH 2.00 2.00 0.92 0.3988
# s(Time,Speaker):SoundT 224.41 377.00 5.20 <2e-16 ***
# s(Time,Speaker):SoundTH 245.82 376.00 7.88 <2e-16 ***
# s(Time,Word):LangEN 71.81 88.00 57.64 <2e-16 ***
# s(Time,Word):LangNL 75.34 90.00 85.55 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.626 Deviance explained = 63.3%
# fREML = -4928.1 Scale est. = 0.41178 n = 31225par(mfrow=c(1,2))
plot(ffmc1,select=3,shade=T,rug=F, ylim=c(-0.8,2.0), main="ffmc1: English difference")
abline(h=0)
plot(ffmc1,select=4,shade=T,rug=F, ylim=c(-1.8,1.0), main="ffmc1: Dutch vs. English difference")
abline(h=0)
7.1 Model criticism
gam.check(ffmc1)
#
# Method: fREML Optimizer: perf chol
# $grad
# [1] -5.06e-08 -7.20e-07 8.03e-08 -1.93e-05 9.80e-08 -1.26e-09 2.51e-08 1.45e-10 -3.23e-09
# [10] -4.71e-12 1.11e-06 -5.42e-10 -9.51e-07
#
# $hess
# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
# 9.98e-01 -0.065977 0.191355 1.15e-04 -2.10e-02 -2.12e-05 3.11e-02 1.00e-03 4.30e-02
# -6.60e-02 1.475521 0.149282 -3.90e-05 7.85e-02 -3.64e-03 -6.05e-02 -9.40e-04 -1.09e-01
# 1.91e-01 0.149282 1.204532 -1.75e-04 -2.12e-02 1.09e-03 5.79e-02 -2.62e-03 3.05e-01
# 1.15e-04 -0.000039 -0.000175 1.95e-05 9.43e-06 -3.05e-06 1.73e-06 3.18e-06 1.25e-04
# -2.10e-02 0.078525 -0.021200 9.43e-06 7.00e+01 -1.55e+00 -8.62e-02 -2.08e-01 -1.24e-01
# -2.12e-05 -0.003644 0.001094 -3.05e-06 -1.55e+00 1.59e+01 -1.38e-01 -7.60e-03 -1.01e-02
# 3.11e-02 -0.060452 0.057880 1.73e-06 -8.62e-02 -1.38e-01 8.01e+01 -7.41e-01 -1.49e-01
# 1.00e-03 -0.000940 -0.002616 3.18e-06 -2.08e-01 -7.60e-03 -7.41e-01 1.82e+01 -1.38e-02
# 4.30e-02 -0.108533 0.304655 1.25e-04 -1.24e-01 -1.01e-02 -1.49e-01 -1.38e-02 3.06e+01
# 1.13e-03 -0.000059 0.000471 2.47e-07 -6.09e-03 -4.20e-03 1.72e-02 -1.83e-04 1.15e-01
# 5.33e-02 0.354251 -0.237855 -4.80e-05 -7.57e-02 9.45e-03 -1.74e-01 -3.21e-02 4.58e-02
# -3.89e-04 -0.000325 0.000662 -9.11e-07 1.90e-02 -4.61e-03 6.36e-03 8.72e-03 -5.24e-04
# d -1.17e+00 -2.002814 -1.640884 -3.77e-04 -9.41e+01 -1.81e+01 -1.04e+02 -1.91e+01 -3.19e+01
# [,10] [,11] [,12] [,13]
# 1.13e-03 5.33e-02 -3.89e-04 -1.17e+00
# -5.90e-05 3.54e-01 -3.25e-04 -2.00e+00
# 4.71e-04 -2.38e-01 6.62e-04 -1.64e+00
# 2.47e-07 -4.80e-05 -9.11e-07 -3.77e-04
# -6.09e-03 -7.57e-02 1.90e-02 -9.41e+01
# -4.20e-03 9.45e-03 -4.61e-03 -1.81e+01
# 1.72e-02 -1.74e-01 6.36e-03 -1.04e+02
# -1.83e-04 -3.21e-02 8.72e-03 -1.91e+01
# 1.15e-01 4.58e-02 -5.24e-04 -3.19e+01
# 3.93e+00 4.84e-05 -7.83e-07 -3.96e+00
# 4.84e-05 3.29e+01 -4.69e-02 -3.37e+01
# -7.83e-07 -4.69e-02 3.97e+00 -3.98e+00
# d -3.96e+00 -3.37e+01 -3.98e+00 1.56e+04
#
# Model rank = 976 / 976
#
# 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):LangEN 9.00 3.34 0.99 0.26
# s(Time):LangNL 9.00 5.01 0.99 0.19
# s(Time):IsTH 10.00 5.28 0.99 0.23
# s(Time):IsNLTH 10.00 2.00 0.99 0.30
# s(Time,Speaker):SoundT 378.00 224.41 0.99 0.22
# s(Time,Speaker):SoundTH 378.00 245.82 0.99 0.24
# s(Time,Word):LangEN 90.00 71.81 0.99 0.24
# s(Time,Word):LangNL 90.00 75.34 0.99 0.28Model criticism shows that the distribution of the residuals is not normal. As it is similar to a scaled-\(t\) model, we refit the model with setting family="scat". Note that this is substantially slower (taking about 25-30 times as long).
7.2 Ordered factor model
fullfinal$IsTHO <- as.ordered(fullfinal$Sound == "TH")
contrasts(fullfinal$IsTHO) <- "contr.treatment"
fullfinal$IsNLTHO <- as.ordered(fullfinal$Sound == "TH" & fullfinal$Lang == "NL")
contrasts(fullfinal$IsNLTHO) <- "contr.treatment"
system.time(
ffmo0 <- bam(Pos ~ Lang + s(Time,by=Lang) + s(Time,by=IsTHO) + IsTHO +
s(Time,by=IsNLTHO) + IsNLTHO +
s(Time,Speaker,by=Sound,bs="fs",m=1) + s(Time,Word,by=Lang,bs="fs",m=1),
data=fullfinal, discrete=TRUE)
)# Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated 1-d smooths of same
# variable.# user system elapsed
# 22.85 40.78 4.41ffmo0acf <- acf_resid(ffmo0)
(rhoval <- ffmo0acf[2]) # autocorrelation at lag 1# 1
# 0.951system.time(
ffmo1 <- bam(Pos ~ Lang + s(Time,by=Lang) + s(Time,by=IsTHO) + IsTHO +
s(Time,by=IsNLTHO) + IsNLTHO +
s(Time,Speaker,by=Sound,bs="fs",m=1) + s(Time,Word,by=Lang,bs="fs",m=1),
data=fullfinal, discrete=TRUE, rho=rhoval, AR.start=fullfinal$start.event)
)# Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated 1-d smooths of same
# variable.# user system elapsed
# 35.5 51.9 10.4acf_resid(ffmo1)
(smryffmo1 <- summary(ffmo1))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ Lang + s(Time, by = Lang) + s(Time, by = IsTHO) + IsTHO +
# s(Time, by = IsNLTHO) + IsNLTHO + s(Time, Speaker, by = Sound,
# bs = "fs", m = 1) + s(Time, Word, by = Lang, bs = "fs", m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.137 0.249 -0.55 0.5831
# LangNL 0.175 0.452 0.39 0.6979
# IsTHOTRUE 1.029 0.364 2.83 0.0047 **
# IsNLTHOTRUE -0.420 0.653 -0.64 0.5194
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):LangEN 3.34 3.75 3.44 0.0117 *
# s(Time):LangNL 5.01 5.43 2.53 0.0225 *
# s(Time):IsTHOTRUE 4.28 4.71 3.91 0.0021 **
# s(Time):IsNLTHOTRUE 1.00 1.01 1.40 0.2349
# s(Time,Speaker):SoundT 224.41 376.00 4.50 <2e-16 ***
# s(Time,Speaker):SoundTH 245.82 376.00 7.11 <2e-16 ***
# s(Time,Word):LangEN 71.81 88.00 57.52 <2e-16 ***
# s(Time,Word):LangNL 75.34 89.00 83.80 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.626 Deviance explained = 63.3%
# fREML = -4928.1 Scale est. = 0.41178 n = 31225par(mfrow=c(1,2))
plot(ffmo1,select=3,shade=T,rug=F, ylim=c(-1.3,1.5), main="ffmo1: English difference")
abline(h=0)
plot(ffmo1,select=4,shade=T,rug=F, ylim=c(-1.3,1.5), main="ffmo1: Dutch vs. English difference")
abline(h=0)
7.3 Scaled-\(t\) model
system.time(
ffmc0s <- bam(Pos ~ Lang + s(Time, by=Lang) + s(Time, by=IsTH) + s(Time, by=IsNLTH) +
s(Time,Speaker,by=Sound,bs="fs",m=1) + s(Time,Word,by=Lang,bs="fs",m=1),
data=fullfinal, discrete=TRUE, family="scat")
)# Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated 1-d smooths of same
# variable.# user system elapsed
# 59.0 56.2 34.2ffmc0sacf <- acf_resid(ffmc0s)
(rhoval <- ffmc0sacf[2]) # autocorrelation at lag 1# 1
# 0.946system.time(
ffmc1s <- bam(Pos ~ Lang + s(Time, by=Lang) + s(Time, by=IsTH) + s(Time, by=IsNLTH) +
s(Time,Speaker,by=Sound,bs="fs",m=1) + s(Time,Word,by=Lang,bs="fs",m=1),
data=fullfinal, discrete=TRUE, rho=rhoval, AR.start=fullfinal$start.event,
family="scat")
)# Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated 1-d smooths of same
# variable.# user system elapsed
# 183 87 143acf_resid(ffmc1s)
(smryffmc1s <- summary(ffmc1s))#
# Family: Scaled t(3.263,0.518)
# Link function: identity
#
# Formula:
# Pos ~ Lang + s(Time, by = Lang) + s(Time, by = IsTH) + s(Time,
# by = IsNLTH) + s(Time, Speaker, by = Sound, bs = "fs", m = 1) +
# s(Time, Word, by = Lang, bs = "fs", m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.165 0.241 -0.68 0.49
# LangNL 0.213 0.432 0.49 0.62
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):LangEN 3.00 3.37 2.59 0.045 *
# s(Time):LangNL 5.15 5.54 2.90 0.010 *
# s(Time):IsTH 5.73 6.13 4.79 5.8e-05 ***
# s(Time):IsNLTH 2.93 3.18 1.04 0.388
# s(Time,Speaker):SoundT 229.22 376.00 4.38 < 2e-16 ***
# s(Time,Speaker):SoundTH 238.07 376.00 11.86 < 2e-16 ***
# s(Time,Word):LangEN 74.14 88.00 58.81 < 2e-16 ***
# s(Time,Word):LangNL 77.04 89.00 82.00 < 2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.614 Deviance explained = 50.9%
# fREML = 18505 Scale est. = 1 n = 31225par(mfrow=c(1,2))
plot(ffmc1s,select=3,shade=T,rug=F, ylim=c(-0.8,2.0), main="ffmc1s: English difference")
abline(h=0)
plot(ffmc1s,select=4,shade=T,rug=F, ylim=c(-1.8,1.0), main="ffmc1s: Dutch vs. English difference")
abline(h=0)
7.3.1 Model criticism
gam.check(ffmc1s)
#
# Method: fREML Optimizer: perf chol
# $grad
# [1] -1.15e-12 -8.04e-10 -3.33e-11 -3.09e-11 3.17e-11 7.12e-11 -8.28e-11 1.02e-11 7.89e-11
# [10] 2.67e-12 2.41e-09 1.42e-09
#
# $hess
# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
# [1,] 6.15e-01 -7.11e-02 0.240260 0.074660 0.018701 5.46e-05 0.015027 0.000104 0.23595
# [2,] -7.11e-02 1.42e+00 0.222243 -0.009241 0.099547 -3.25e-03 -0.038215 -0.000443 -0.12104
# [3,] 2.40e-01 2.22e-01 1.407399 -0.088046 -0.034967 6.48e-04 0.048968 -0.001024 0.26812
# [4,] 7.47e-02 -9.24e-03 -0.088046 0.261340 -0.005872 -2.06e-03 -0.008372 0.001742 0.04258
# [5,] 1.87e-02 9.95e-02 -0.034967 -0.005872 69.536816 -1.76e+00 0.122283 -0.375058 -0.09976
# [6,] 5.46e-05 -3.25e-03 0.000648 -0.002062 -1.763025 1.67e+01 -0.189981 -0.004590 -0.00186
# [7,] 1.50e-02 -3.82e-02 0.048968 -0.008372 0.122283 -1.90e-01 75.413663 -0.530734 -0.06231
# [8,] 1.04e-04 -4.43e-04 -0.001024 0.001742 -0.375058 -4.59e-03 -0.530734 18.411703 -0.00551
# [9,] 2.36e-01 -1.21e-01 0.268120 0.042582 -0.099756 -1.86e-03 -0.062312 -0.005506 30.96991
# [10,] 7.85e-04 -4.86e-05 0.000457 0.000210 -0.011468 -4.83e-03 0.017514 -0.000713 0.14256
# [11,] 3.08e-02 5.50e-01 -0.256232 -0.167618 0.033335 7.32e-03 -0.039243 -0.026909 0.02706
# [12,] -2.09e-04 -4.08e-04 0.000366 -0.000736 -0.000105 -2.02e-03 -0.000939 0.010330 -0.00021
# [,10] [,11] [,12]
# [1,] 7.85e-04 3.08e-02 -2.09e-04
# [2,] -4.86e-05 5.50e-01 -4.08e-04
# [3,] 4.57e-04 -2.56e-01 3.66e-04
# [4,] 2.10e-04 -1.68e-01 -7.36e-04
# [5,] -1.15e-02 3.33e-02 -1.05e-04
# [6,] -4.83e-03 7.32e-03 -2.02e-03
# [7,] 1.75e-02 -3.92e-02 -9.39e-04
# [8,] -7.13e-04 -2.69e-02 1.03e-02
# [9,] 1.43e-01 2.71e-02 -2.10e-04
# [10,] 3.92e+00 3.11e-06 -4.90e-07
# [11,] 3.11e-06 3.29e+01 -3.00e-02
# [12,] -4.90e-07 -3.00e-02 3.96e+00
#
# Model rank = 976 / 976
#
# 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):LangEN 9.00 3.00 1.03 0.95
# s(Time):LangNL 9.00 5.15 1.03 0.96
# s(Time):IsTH 10.00 5.73 1.03 0.96
# s(Time):IsNLTH 10.00 2.93 1.03 0.98
# s(Time,Speaker):SoundT 378.00 229.22 1.03 0.96
# s(Time,Speaker):SoundTH 378.00 238.07 1.03 0.98
# s(Time,Word):LangEN 90.00 74.14 1.03 0.96
# s(Time,Word):LangNL 90.00 77.04 1.03 0.967.4 NEW random effects
In this section, the model is fitted with a random-effects specification, which does not result in overly conservative (i.e. too high \(p\)-values) for the difference between the two sounds and languages. This approach follows Sóskuthy (2021).
# create ordered factors for appropriate random-effects specification (not overly conservative; Soskuthy, 2021)
fullfinal$SoundO <- as.ordered(fullfinal$Sound)
contrasts(fullfinal$SoundO) <- "contr.treatment"
fullfinal$LangO <- as.ordered(fullfinal$Lang)
contrasts(fullfinal$LangO) <- "contr.treatment"
system.time(
ffmc0s2 <- bam(Pos ~ Lang + s(Time, by=Lang) + s(Time, by=IsTH) + s(Time, by=IsNLTH) +
s(Time,Speaker,bs="fs",m=1) + s(Time,Speaker,by=SoundO,bs="fs",m=1) +
s(Time,Word,bs="fs",m=1) + s(Time,Word,by=LangO,bs="fs",m=1),
data=fullfinal, discrete=TRUE, family="scat")
)# Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated 1-d smooths of same
# variable.# user system elapsed
# 53.1 53.0 29.8ffmc0sacf2 <- acf_resid(ffmc0s2)
(rhoval <- ffmc0sacf2[2]) # autocorrelation at lag 1# 1
# 0.946system.time(
ffmc1s2 <- bam(Pos ~ Lang + s(Time, by=Lang) + s(Time, by=IsTH) + s(Time, by=IsNLTH) +
s(Time,Speaker,bs="fs",m=1) + s(Time,Speaker,by=SoundO,bs="fs",m=1) +
s(Time,Word,bs="fs",m=1) + s(Time,Word,by=LangO,bs="fs",m=1),
data=fullfinal, discrete=TRUE, rho=rhoval, AR.start=fullfinal$start.event,
family="scat")
)# Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated 1-d smooths of same
# variable.# user system elapsed
# 85.3 41.0 66.5acf_resid(ffmc1s2)
(smryffmc1s2 <- summary(ffmc1s2))#
# Family: Scaled t(3.797,0.543)
# Link function: identity
#
# Formula:
# Pos ~ Lang + s(Time, by = Lang) + s(Time, by = IsTH) + s(Time,
# by = IsNLTH) + s(Time, Speaker, bs = "fs", m = 1) + s(Time,
# Speaker, by = SoundO, bs = "fs", m = 1) + s(Time, Word, bs = "fs",
# m = 1) + s(Time, Word, by = LangO, bs = "fs", m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.133 0.250 -0.53 0.59
# LangNL 0.184 0.175 1.05 0.29
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):LangEN 2.29 2.62 2.82 0.02957 *
# s(Time):LangNL 5.38 5.92 3.23 0.00377 **
# s(Time):IsTH 5.34 5.74 4.43 0.00025 ***
# s(Time):IsNLTH 3.08 3.37 2.72 0.03741 *
# s(Time,Speaker) 215.97 377.00 4.81 < 2e-16 ***
# s(Time,Speaker):SoundOTH 218.67 376.00 8.11 < 2e-16 ***
# s(Time,Word) 75.20 89.00 55.48 < 2e-16 ***
# s(Time,Word):LangONL 51.22 89.00 4.76 < 2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.627 Deviance explained = 52.1%
# fREML = 17310 Scale est. = 1 n = 31225par(mfrow=c(1,2))
plot(ffmc1s2,select=3,shade=T,rug=F, ylim=c(-0.8,2.0), main="ffmc1s2: English difference")
abline(h=0)
plot(ffmc1s2,select=4,shade=T,rug=F, ylim=c(-1.8,1.0), main="ffmc1s2: Dutch vs. English difference")
abline(h=0)
7.4.1 Model criticism
gam.check(ffmc1s)
#
# Method: fREML Optimizer: perf chol
# $grad
# [1] -1.15e-12 -8.04e-10 -3.33e-11 -3.09e-11 3.17e-11 7.12e-11 -8.28e-11 1.02e-11 7.89e-11
# [10] 2.67e-12 2.41e-09 1.42e-09
#
# $hess
# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
# [1,] 6.15e-01 -7.11e-02 0.240260 0.074660 0.018701 5.46e-05 0.015027 0.000104 0.23595
# [2,] -7.11e-02 1.42e+00 0.222243 -0.009241 0.099547 -3.25e-03 -0.038215 -0.000443 -0.12104
# [3,] 2.40e-01 2.22e-01 1.407399 -0.088046 -0.034967 6.48e-04 0.048968 -0.001024 0.26812
# [4,] 7.47e-02 -9.24e-03 -0.088046 0.261340 -0.005872 -2.06e-03 -0.008372 0.001742 0.04258
# [5,] 1.87e-02 9.95e-02 -0.034967 -0.005872 69.536816 -1.76e+00 0.122283 -0.375058 -0.09976
# [6,] 5.46e-05 -3.25e-03 0.000648 -0.002062 -1.763025 1.67e+01 -0.189981 -0.004590 -0.00186
# [7,] 1.50e-02 -3.82e-02 0.048968 -0.008372 0.122283 -1.90e-01 75.413663 -0.530734 -0.06231
# [8,] 1.04e-04 -4.43e-04 -0.001024 0.001742 -0.375058 -4.59e-03 -0.530734 18.411703 -0.00551
# [9,] 2.36e-01 -1.21e-01 0.268120 0.042582 -0.099756 -1.86e-03 -0.062312 -0.005506 30.96991
# [10,] 7.85e-04 -4.86e-05 0.000457 0.000210 -0.011468 -4.83e-03 0.017514 -0.000713 0.14256
# [11,] 3.08e-02 5.50e-01 -0.256232 -0.167618 0.033335 7.32e-03 -0.039243 -0.026909 0.02706
# [12,] -2.09e-04 -4.08e-04 0.000366 -0.000736 -0.000105 -2.02e-03 -0.000939 0.010330 -0.00021
# [,10] [,11] [,12]
# [1,] 7.85e-04 3.08e-02 -2.09e-04
# [2,] -4.86e-05 5.50e-01 -4.08e-04
# [3,] 4.57e-04 -2.56e-01 3.66e-04
# [4,] 2.10e-04 -1.68e-01 -7.36e-04
# [5,] -1.15e-02 3.33e-02 -1.05e-04
# [6,] -4.83e-03 7.32e-03 -2.02e-03
# [7,] 1.75e-02 -3.92e-02 -9.39e-04
# [8,] -7.13e-04 -2.69e-02 1.03e-02
# [9,] 1.43e-01 2.71e-02 -2.10e-04
# [10,] 3.92e+00 3.11e-06 -4.90e-07
# [11,] 3.11e-06 3.29e+01 -3.00e-02
# [12,] -4.90e-07 -3.00e-02 3.96e+00
#
# Model rank = 976 / 976
#
# 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):LangEN 9.00 3.00 1 0.33
# s(Time):LangNL 9.00 5.15 1 0.34
# s(Time):IsTH 10.00 5.73 1 0.36
# s(Time):IsNLTH 10.00 2.93 1 0.39
# s(Time,Speaker):SoundT 378.00 229.22 1 0.36
# s(Time,Speaker):SoundTH 378.00 238.07 1 0.35
# s(Time,Word):LangEN 90.00 74.14 1 0.37
# s(Time,Word):LangNL 90.00 77.04 1 0.418 Full model
In this section, we will assess if Dutch speakers contrast /θ/-words from /t/-words less strongly than English speakers for both word groups. Since we do not wish to explicitly contrast word-initial from word-final patterns, we create a model with four reference levels (the four possible combinations of language and location, combined into the variable LangLoc, below) and four binary difference smooths (IsTHInit, IsTHFinal, IsNLTHInit, and IsNLTHFinal, created below). The first two binary difference smooths will inform us about the difference between the /θ/-sound and /t/-sound for English speakers, separately for words in which the two sounds are word-initial and word-final. The final two difference smooths will indicate if these differences are substantially smaller for the Dutch speakers. As there might be individual speaker variability in how the sounds /θ/ and /t/ are pronounced depending on their location (word-initial or word-final), this will need to be included in the by-speaker factor smooth (using the variable SoundLoc, below).
full$LangLoc <- interaction(full$Lang, full$Loc)
full$IsTHInit <- (full$Sound == "TH" & full$Loc == "Init")*1
full$IsTHFinal <- (full$Sound == "TH" & full$Loc == "Final")*1
full$IsNLTHInit <- (full$Lang == "NL" & full$Sound == "TH" & full$Loc == "Init")*1
full$IsNLTHFinal <- (full$Lang == "NL" & full$Sound == "TH" & full$Loc == "Final")*1
full$SoundLoc <- interaction(full$Sound, full$Loc)
# track first observation per trial (to correct for autocorrelation later on)
# and sort per trajectory
full$start.event = (full$Time == 0)
full = full[order(full$Speaker,full$Trial,full$Time),]
#full <- start_event(full,event=c("Speaker","Trial")) # bug in itsadug 2.4
system.time(
fmc0 <- bam(Pos ~ LangLoc + s(Time,by=LangLoc) + s(Time,by=IsTHInit) + s(Time,by=IsTHFinal) +
s(Time,by=IsNLTHInit) + s(Time,by=IsNLTHFinal) +
s(Time,Speaker,by=SoundLoc,bs="fs",m=1) + s(Time,Word,by=Lang,bs="fs",m=1),
data=full, discrete=TRUE)
)# user system elapsed
# 128.5 155.8 52.6fmc0acf <- acf_resid(fmc0)
(rhoval <- fmc0acf[2]) # autocorrelation at lag 1# 1
# 0.963system.time(
fmc1 <- bam(Pos ~ LangLoc + s(Time,by=LangLoc) + s(Time,by=IsTHInit) + s(Time,by=IsTHFinal) +
s(Time,by=IsNLTHInit) + s(Time,by=IsNLTHFinal) +
s(Time,Speaker,by=SoundLoc,bs="fs",m=1) + s(Time,Word,by=Lang,bs="fs",m=1),
data=full, discrete=TRUE, rho=rhoval, AR.start=full$start.event)
)# user system elapsed
# 183 137 116acf_resid(fmc1)
(smryfmc1 <- summary(fmc1))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ LangLoc + s(Time, by = LangLoc) + s(Time, by = IsTHInit) +
# s(Time, by = IsTHFinal) + s(Time, by = IsNLTHInit) + s(Time,
# by = IsNLTHFinal) + s(Time, Speaker, by = SoundLoc, bs = "fs",
# m = 1) + s(Time, Word, by = Lang, bs = "fs", m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.0669 0.1667 -0.40 0.69
# LangLocNL.Final 0.1167 0.2661 0.44 0.66
# LangLocEN.Init 0.1324 0.2356 0.56 0.57
# LangLocNL.Init 0.3470 0.2661 1.30 0.19
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):LangLocEN.Final 3.34 3.70 2.89 0.04564 *
# s(Time):LangLocNL.Final 5.91 6.21 2.43 0.02526 *
# s(Time):LangLocEN.Init 7.02 7.24 3.98 0.00011 ***
# s(Time):LangLocNL.Init 6.89 7.12 3.56 0.00049 ***
# s(Time):IsTHInit 6.66 6.97 5.92 5e-06 ***
# s(Time):IsTHFinal 6.79 7.10 3.40 0.00129 **
# s(Time):IsNLTHInit 2.00 2.00 0.87 0.42090
# s(Time):IsNLTHFinal 2.00 2.00 0.55 0.57466
# s(Time,Speaker):SoundLocT.Final 266.94 376.00 3.07 < 2e-16 ***
# s(Time,Speaker):SoundLocTH.Final 286.99 376.00 4.62 < 2e-16 ***
# s(Time,Speaker):SoundLocT.Init 307.41 376.00 5.82 < 2e-16 ***
# s(Time,Speaker):SoundLocTH.Init 302.64 376.00 5.28 < 2e-16 ***
# s(Time,Word):LangEN 153.34 176.00 16.20 < 2e-16 ***
# s(Time,Word):LangNL 154.76 176.00 19.96 < 2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.593 Deviance explained = 59.8%
# fREML = -49820 Scale est. = 0.34505 n = 126177par(mfrow=c(2,2))
plot(fmc1,select=5,shade=T,rug=F, ylim=c(-0.8,2.0), main="fmc1: English difference initial")
abline(h=0)
plot(fmc1,select=6,shade=T,rug=F, ylim=c(-0.8,2.0), main="fmc1: English difference final")
abline(h=0)
plot(fmc1,select=7,shade=T,rug=F, ylim=c(-1.8,1.0), main="fmc1: Dutch vs. English difference initial")
abline(h=0)
plot(fmc1,select=8,shade=T,rug=F, ylim=c(-1.8,1.0), main="fmc1: Dutch vs. English difference final")
abline(h=0)
8.1 Model criticism
gam.check(fmc1)
#
# Method: fREML Optimizer: perf chol
# $grad
# [1] 3.81e-02 2.74e-01 -2.22e-01 4.77e-01 4.25e-01 -2.56e-01 -8.93e-04 5.04e-05 3.62e-01
# [10] 3.33e-02 4.20e-01 5.13e-01 2.60e+00 -1.93e-02 7.22e-01 2.25e-01 4.11e+00 7.47e-01
# [19] 4.26e+00 2.60e+00 -1.07e+02
#
# $hess
# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
# 7.65e-01 -7.12e-02 1.74e-06 -9.55e-07 -5.08e-07 1.43e-01 1.13e-08 6.77e-05 2.55e-02
# -7.12e-02 1.07e+00 -3.75e-07 -1.25e-06 1.47e-08 2.45e-01 -4.06e-09 -1.06e-05 1.78e-01
# 1.74e-06 -3.75e-07 1.91e+00 -1.89e-01 2.76e-01 2.32e-06 -3.84e-05 3.30e-09 9.07e-05
# -9.55e-07 -1.25e-06 -1.89e-01 1.99e+00 3.25e-01 -1.25e-06 7.13e-06 -1.46e-09 4.57e-05
# -5.08e-07 1.47e-08 2.76e-01 3.25e-01 2.38e+00 -6.78e-07 -2.33e-04 -9.41e-10 -2.04e-05
# 1.43e-01 2.45e-01 2.32e-06 -1.25e-06 -6.78e-07 1.57e+00 1.51e-08 -1.76e-05 -2.54e-02
# 1.13e-08 -4.06e-09 -3.84e-05 7.13e-06 -2.33e-04 1.51e-08 8.92e-04 2.19e-11 6.94e-07
# 6.77e-05 -1.06e-05 3.30e-09 -1.46e-09 -9.41e-10 -1.76e-05 2.19e-11 -5.03e-05 4.18e-06
# 2.55e-02 1.78e-01 9.07e-05 4.57e-05 -2.04e-05 -2.54e-02 6.94e-07 4.18e-06 8.68e+01
# -4.47e-04 -2.52e-04 -1.76e-05 -8.90e-06 3.87e-06 1.16e-04 -1.38e-07 -3.69e-06 -1.54e-01
# 1.83e-02 9.27e-03 -9.37e-05 8.34e-05 2.94e-05 1.01e-01 -5.75e-07 -3.13e-06 -5.36e-03
# 7.87e-04 -1.45e-04 1.72e-05 -1.53e-05 -5.42e-06 -1.85e-04 1.05e-07 1.48e-06 9.52e-05
# 5.13e-05 -1.36e-04 2.10e-01 2.82e-01 -1.33e-01 7.07e-05 6.53e-05 1.28e-07 1.07e-02
# 2.25e-05 -6.11e-05 -2.16e-04 -6.28e-04 9.07e-04 3.10e-05 -3.46e-06 5.63e-08 4.59e-03
# -4.65e-04 -2.71e-04 4.76e-02 -4.21e-02 6.96e-02 -6.16e-04 4.04e-05 -7.93e-07 1.81e-05
# -3.43e-04 -2.00e-04 -1.99e-04 -5.52e-04 6.60e-04 -4.55e-04 2.59e-06 -5.85e-07 -9.02e-06
# 2.97e-01 -1.42e-01 8.22e-01 -2.11e-02 1.85e-01 4.66e-01 2.29e-04 3.79e-05 1.74e-02
# -1.42e-03 1.88e-05 7.95e-06 3.79e-06 -8.07e-05 1.20e-04 -1.07e-07 7.36e-08 -6.32e-03
# -1.27e-02 1.58e+00 -3.97e-02 8.56e-01 -3.69e-01 -3.57e-02 1.94e-04 -3.78e-05 2.13e-01
# -8.49e-05 2.48e-04 -3.37e-05 -3.00e-05 -3.44e-05 -8.68e-05 -3.55e-07 -2.43e-07 -1.71e-02
# d -1.21e+00 -2.73e+00 -2.79e+00 -3.42e+00 -2.75e+00 -2.14e+00 -4.06e-04 -2.41e-04 -1.18e+02
# [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18]
# -4.47e-04 1.83e-02 7.87e-04 5.13e-05 2.25e-05 -4.65e-04 -3.43e-04 2.97e-01 -1.42e-03
# -2.52e-04 9.27e-03 -1.45e-04 -1.36e-04 -6.11e-05 -2.71e-04 -2.00e-04 -1.42e-01 1.88e-05
# -1.76e-05 -9.37e-05 1.72e-05 2.10e-01 -2.16e-04 4.76e-02 -1.99e-04 8.22e-01 7.95e-06
# -8.90e-06 8.34e-05 -1.53e-05 2.82e-01 -6.28e-04 -4.21e-02 -5.52e-04 -2.11e-02 3.79e-06
# 3.87e-06 2.94e-05 -5.42e-06 -1.33e-01 9.07e-04 6.96e-02 6.60e-04 1.85e-01 -8.07e-05
# 1.16e-04 1.01e-01 -1.85e-04 7.07e-05 3.10e-05 -6.16e-04 -4.55e-04 4.66e-01 1.20e-04
# -1.38e-07 -5.75e-07 1.05e-07 6.53e-05 -3.46e-06 4.04e-05 2.59e-06 2.29e-04 -1.07e-07
# -3.69e-06 -3.13e-06 1.48e-06 1.28e-07 5.63e-08 -7.93e-07 -5.85e-07 3.79e-05 7.36e-08
# -1.54e-01 -5.36e-03 9.52e-05 1.07e-02 4.59e-03 1.81e-05 -9.02e-06 1.74e-02 -6.32e-03
# 1.22e+01 -4.94e-05 -5.22e-05 -2.34e-03 -1.27e-03 1.35e-05 1.43e-05 -3.58e-03 -3.17e-03
# -4.94e-05 9.77e+01 -2.09e-02 -4.55e-05 4.89e-06 3.30e-02 2.44e-02 1.19e-01 -5.76e-04
# -5.22e-05 -2.09e-02 1.66e+01 9.30e-06 -4.71e-07 -6.16e-03 -4.54e-03 -2.52e-03 -8.91e-03
# -2.34e-03 -4.55e-05 9.30e-06 1.17e+02 -8.48e-01 1.27e-02 8.90e-04 1.06e-01 4.57e-04
# -1.27e-03 4.89e-06 -4.71e-07 -8.48e-01 1.27e+01 -5.61e-04 -2.97e-05 -1.63e-03 6.93e-04
# 1.35e-05 3.30e-02 -6.16e-03 1.27e-02 -5.61e-04 1.09e+02 -2.88e-01 2.60e-01 7.99e-04
# 1.43e-05 2.44e-02 -4.54e-03 8.90e-04 -2.97e-05 -2.88e-01 1.57e+01 -6.73e-04 -8.32e-04
# -3.58e-03 1.19e-01 -2.52e-03 1.06e-01 -1.63e-03 2.60e-01 -6.73e-04 6.54e+01 -2.28e-01
# -3.17e-03 -5.76e-04 -8.91e-03 4.57e-04 6.93e-04 7.99e-04 -8.32e-04 -2.28e-01 8.26e+00
# -6.20e-03 -9.48e-02 -2.69e-03 9.94e-01 -1.55e-02 2.92e-01 3.74e-03 2.73e-01 -8.90e-05
# -3.20e-02 -6.82e-04 2.68e-02 -5.92e-03 -1.01e-02 -2.18e-03 6.05e-03 -1.31e-04 -2.19e-07
# d -1.62e+01 -1.25e+02 -1.90e+01 -1.40e+02 -1.61e+01 -1.34e+02 -1.83e+01 -7.30e+01 -8.52e+00
# [,19] [,20] [,21]
# -1.27e-02 -8.49e-05 -1.21e+00
# 1.58e+00 2.48e-04 -2.73e+00
# -3.97e-02 -3.37e-05 -2.79e+00
# 8.56e-01 -3.00e-05 -3.42e+00
# -3.69e-01 -3.44e-05 -2.75e+00
# -3.57e-02 -8.68e-05 -2.14e+00
# 1.94e-04 -3.55e-07 -4.06e-04
# -3.78e-05 -2.43e-07 -2.41e-04
# 2.13e-01 -1.71e-02 -1.18e+02
# -6.20e-03 -3.20e-02 -1.62e+01
# -9.48e-02 -6.82e-04 -1.25e+02
# -2.69e-03 2.68e-02 -1.90e+01
# 9.94e-01 -5.92e-03 -1.40e+02
# -1.55e-02 -1.01e-02 -1.61e+01
# 2.92e-01 -2.18e-03 -1.34e+02
# 3.74e-03 6.05e-03 -1.83e+01
# 2.73e-01 -1.31e-04 -7.30e+01
# -8.90e-05 -2.19e-07 -8.52e+00
# 6.49e+01 -2.62e-01 -7.38e+01
# -2.62e-01 1.01e+01 -1.04e+01
# d -7.38e+01 -1.04e+01 6.32e+04
#
# Model rank = 1952 / 1952
#
# 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):LangLocEN.Final 9.00 3.34 1.01 0.71
# s(Time):LangLocNL.Final 9.00 5.91 1.01 0.75
# s(Time):LangLocEN.Init 9.00 7.02 1.01 0.72
# s(Time):LangLocNL.Init 9.00 6.89 1.01 0.70
# s(Time):IsTHInit 10.00 6.66 1.01 0.70
# s(Time):IsTHFinal 10.00 6.79 1.01 0.70
# s(Time):IsNLTHInit 10.00 2.00 1.01 0.70
# s(Time):IsNLTHFinal 10.00 2.00 1.01 0.70
# s(Time,Speaker):SoundLocT.Final 378.00 266.94 1.01 0.80
# s(Time,Speaker):SoundLocTH.Final 378.00 286.99 1.01 0.69
# s(Time,Speaker):SoundLocT.Init 378.00 307.41 1.01 0.69
# s(Time,Speaker):SoundLocTH.Init 378.00 302.64 1.01 0.71
# s(Time,Word):LangEN 180.00 153.34 1.01 0.74
# s(Time,Word):LangNL 180.00 154.76 1.01 0.77Model criticism shows that the distribution of the residuals is not normal. As it is similar to a scaled-\(t\) model, we refit the model with setting family="scat". Note that this is substantially slower (taking about 25-30 times as long).
8.2 Ordered factor model
full$IsTHInitO <- as.ordered(full$Sound == "TH" & full$Loc == "Init")
contrasts(full$IsTHInitO) <- "contr.treatment"
full$IsTHFinalO <- as.ordered(full$Sound == "TH" & full$Loc == "Final")
contrasts(full$IsTHFinalO) <- "contr.treatment"
full$IsNLTHInitO <- as.ordered(full$Sound == "TH" & full$Loc == "Init" & full$Lang == "NL")
contrasts(full$IsNLTHInitO) <- "contr.treatment"
full$IsNLTHFinalO <- as.ordered(full$Sound == "TH" & full$Loc == "Final" & full$Lang == "NL")
contrasts(full$IsNLTHFinalO) <- "contr.treatment"
system.time(
fmo0 <- bam(Pos ~ LangLoc + s(Time,by=LangLoc) + s(Time,by=IsTHInitO) + IsTHInitO +
s(Time,by=IsTHFinalO) + IsTHFinalO + s(Time,by=IsNLTHInitO) + IsNLTHInitO +
s(Time,by=IsNLTHFinalO) + IsNLTHFinalO +
s(Time,Speaker,by=SoundLoc,bs="fs",m=1) + s(Time,Word,by=Lang,bs="fs",m=1),
data=full, discrete=TRUE)
)# user system elapsed
# 150.6 195.8 54.7fmo0acf <- acf_resid(fmo0)
(rhoval <- fmo0acf[2]) # autocorrelation at lag 1# 1
# 0.963system.time(
fmo1 <- bam(Pos ~ LangLoc + s(Time,by=LangLoc) + s(Time,by=IsTHInitO) + IsTHInitO +
s(Time,by=IsTHFinalO) + IsTHFinalO + s(Time,by=IsNLTHInitO) + IsNLTHInitO +
s(Time,by=IsNLTHFinalO) + IsNLTHFinalO +
s(Time,Speaker,by=SoundLoc,bs="fs",m=1) + s(Time,Word,by=Lang,bs="fs",m=1),
data=full, discrete=TRUE, rho=rhoval, AR.start=full$start.event)
)# user system elapsed
# 218 201 121acf_resid(fmo1)
(smryfmo1 <- summary(fmo1))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ LangLoc + s(Time, by = LangLoc) + s(Time, by = IsTHInitO) +
# IsTHInitO + s(Time, by = IsTHFinalO) + IsTHFinalO + s(Time,
# by = IsNLTHInitO) + IsNLTHInitO + s(Time, by = IsNLTHFinalO) +
# IsNLTHFinalO + s(Time, Speaker, by = SoundLoc, bs = "fs",
# m = 1) + s(Time, Word, by = Lang, bs = "fs", m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.0667 0.1738 -0.38 0.7010
# LangLocNL.Final 0.1169 0.2939 0.40 0.6908
# LangLocEN.Init 0.1320 0.2456 0.54 0.5910
# LangLocNL.Init 0.3468 0.2939 1.18 0.2380
# IsTHInitOTRUE 0.8050 0.2534 3.18 0.0015 **
# IsTHFinalOTRUE 0.5833 0.2583 2.26 0.0239 *
# IsNLTHInitOTRUE -0.4892 0.4252 -1.15 0.2500
# IsNLTHFinalOTRUE -0.2351 0.4314 -0.54 0.5858
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):LangLocEN.Final 3.41 3.77 2.89 0.04351 *
# s(Time):LangLocNL.Final 6.18 6.46 2.41 0.02551 *
# s(Time):LangLocEN.Init 6.73 6.97 3.67 0.00064 ***
# s(Time):LangLocNL.Init 7.04 7.25 3.48 0.00052 ***
# s(Time):IsTHInitOTRUE 5.88 6.18 4.74 4.3e-05 ***
# s(Time):IsTHFinalOTRUE 5.37 5.70 2.93 0.00968 **
# s(Time):IsNLTHInitOTRUE 1.06 1.10 0.14 0.78273
# s(Time):IsNLTHFinalOTRUE 1.00 1.00 0.75 0.38757
# s(Time,Speaker):SoundLocT.Final 266.98 377.00 3.06 < 2e-16 ***
# s(Time,Speaker):SoundLocTH.Final 287.17 376.00 4.62 < 2e-16 ***
# s(Time,Speaker):SoundLocT.Init 308.14 376.00 5.82 < 2e-16 ***
# s(Time,Speaker):SoundLocTH.Init 302.81 376.00 5.28 < 2e-16 ***
# s(Time,Word):LangEN 154.26 177.00 16.15 < 2e-16 ***
# s(Time,Word):LangNL 155.16 178.00 19.77 < 2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.593 Deviance explained = 59.8%
# fREML = -49821 Scale est. = 0.34553 n = 126177par(mfrow=c(2,2))
plot(fmo1,select=5,shade=T,rug=F, ylim=c(-1.3,1.5), main="fmo1: English difference initial")
abline(h=0)
plot(fmo1,select=6,shade=T,rug=F, ylim=c(-1.3,1.5), main="fmo1: English difference final")
abline(h=0)
plot(fmo1,select=7,shade=T,rug=F, ylim=c(-1.3,1.5), main="fmo1: Dutch vs. English difference initial")
abline(h=0)
plot(fmo1,select=8,shade=T,rug=F, ylim=c(-1.3,1.5), main="fmo1: Dutch vs. English difference final")
abline(h=0)
8.3 Model comparison
system.time(
fmc1.ml <- bam(Pos ~ LangLoc + s(Time,by=LangLoc) + s(Time,by=IsTHInit) + s(Time,by=IsTHFinal) +
s(Time,by=IsNLTHInit) + s(Time,by=IsNLTHFinal) +
s(Time,Speaker,by=SoundLoc,bs="fs",m=1) + s(Time,Word,by=Lang,bs="fs",m=1),
data=full, method="ML", rho=rhoval, AR.start=full$start.event)
)# user system elapsed
# 21701 28737 4524system.time(
fmc1a.ml <- bam(Pos ~ LangLoc + s(Time,by=LangLoc) + s(Time,by=IsTHInit) + s(Time,by=IsTHFinal) +
s(Time,by=IsNLTHInit) +
s(Time,Speaker,by=SoundLoc,bs="fs",m=1) + s(Time,Word,by=Lang,bs="fs",m=1),
data=full, method="ML", rho=rhoval, AR.start=full$start.event)
)# user system elapsed
# 18096 23479 3569system.time(
fmc1b.ml <- bam(Pos ~ LangLoc + s(Time,by=LangLoc) + s(Time,by=IsTHInit) + s(Time,by=IsTHFinal) +
s(Time,by=IsNLTHFinal) +
s(Time,Speaker,by=SoundLoc,bs="fs",m=1) + s(Time,Word,by=Lang,bs="fs",m=1),
data=full, method="ML", rho=rhoval, AR.start=full$start.event)
)# user system elapsed
# 16720 21454 3165system.time(
fmc1c.ml <- bam(Pos ~ LangLoc + s(Time,by=LangLoc) + s(Time,by=IsTHInit) + s(Time,by=IsTHFinal) +
s(Time,Speaker,by=SoundLoc,bs="fs",m=1) + s(Time,Word,by=Lang,bs="fs",m=1),
data=full, method="ML", rho=rhoval, AR.start=full$start.event)
)# user system elapsed
# 8640 10058 2035compareML(fmc1a.ml,fmc1.ml)# fmc1a.ml: Pos ~ LangLoc + s(Time, by = LangLoc) + s(Time, by = IsTHInit) +
# s(Time, by = IsTHFinal) + s(Time, by = IsNLTHInit) + s(Time,
# Speaker, by = SoundLoc, bs = "fs", m = 1) + s(Time, Word,
# by = Lang, bs = "fs", m = 1)
#
# fmc1.ml: Pos ~ LangLoc + s(Time, by = LangLoc) + s(Time, by = IsTHInit) +
# s(Time, by = IsTHFinal) + s(Time, by = IsNLTHInit) + s(Time,
# by = IsNLTHFinal) + s(Time, Speaker, by = SoundLoc, bs = "fs",
# m = 1) + s(Time, Word, by = Lang, bs = "fs", m = 1)# Warning in sprintf(" ", h1): one argument not used by format ' '#
# Chi-square test of ML scores
# -----
# Model Score Edf Difference Df p.value Sig.
# 1 fmc1a.ml -49854 33
# 2 fmc1.ml -49855 36 0.582 3.000 0.762
#
# AIC difference: -1.19, model fmc1a.ml has lower AIC.# Warning in compareML(fmc1a.ml, fmc1.ml): Only small difference in ML...compareML(fmc1b.ml,fmc1.ml)# fmc1b.ml: Pos ~ LangLoc + s(Time, by = LangLoc) + s(Time, by = IsTHInit) +
# s(Time, by = IsTHFinal) + s(Time, by = IsNLTHFinal) + s(Time,
# Speaker, by = SoundLoc, bs = "fs", m = 1) + s(Time, Word,
# by = Lang, bs = "fs", m = 1)
#
# fmc1.ml: Pos ~ LangLoc + s(Time, by = LangLoc) + s(Time, by = IsTHInit) +
# s(Time, by = IsTHFinal) + s(Time, by = IsNLTHInit) + s(Time,
# by = IsNLTHFinal) + s(Time, Speaker, by = SoundLoc, bs = "fs",
# m = 1) + s(Time, Word, by = Lang, bs = "fs", m = 1)# Warning in sprintf(" ", h1): one argument not used by format ' '#
# Chi-square test of ML scores
# -----
# Model Score Edf Difference Df p.value Sig.
# 1 fmc1b.ml -49854 33
# 2 fmc1.ml -49855 36 0.805 3.000 0.657
#
# AIC difference: -2.47, model fmc1b.ml has lower AIC.# Warning in compareML(fmc1b.ml, fmc1.ml): Only small difference in ML...compareML(fmc1c.ml,fmc1.ml)# fmc1c.ml: Pos ~ LangLoc + s(Time, by = LangLoc) + s(Time, by = IsTHInit) +
# s(Time, by = IsTHFinal) + s(Time, Speaker, by = SoundLoc,
# bs = "fs", m = 1) + s(Time, Word, by = Lang, bs = "fs", m = 1)
#
# fmc1.ml: Pos ~ LangLoc + s(Time, by = LangLoc) + s(Time, by = IsTHInit) +
# s(Time, by = IsTHFinal) + s(Time, by = IsNLTHInit) + s(Time,
# by = IsNLTHFinal) + s(Time, Speaker, by = SoundLoc, bs = "fs",
# m = 1) + s(Time, Word, by = Lang, bs = "fs", m = 1)# Warning in sprintf(" ", h1): one argument not used by format ' '#
# Chi-square test of ML scores
# -----
# Model Score Edf Difference Df p.value Sig.
# 1 fmc1c.ml -49854 30
# 2 fmc1.ml -49855 36 1.395 6.000 0.835
#
# AIC difference: -3.15, model fmc1c.ml has lower AIC.# Warning in compareML(fmc1c.ml, fmc1.ml): Only small difference in ML...8.4 k set to 20
system.time(
fmk1 <- bam(Pos ~ LangLoc + s(Time,by=LangLoc) + s(Time,by=IsTHInit) + s(Time,by=IsTHFinal) +
s(Time, by=IsNLTHInit, k=20) + s(Time, by=IsNLTHFinal, k=20) +
s(Time,Speaker,by=SoundLoc,bs="fs",m=1) + s(Time,Word,by=Lang,bs="fs",m=1),
data=full, discrete=TRUE, rho=0.963, AR.start=full$start.event)
)# user system elapsed
# 202 168 118acf_resid(fmk1)
(smryfmk1 <- summary(fmk1))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ LangLoc + s(Time, by = LangLoc) + s(Time, by = IsTHInit) +
# s(Time, by = IsTHFinal) + s(Time, by = IsNLTHInit, k = 20) +
# s(Time, by = IsNLTHFinal, k = 20) + s(Time, Speaker, by = SoundLoc,
# bs = "fs", m = 1) + s(Time, Word, by = Lang, bs = "fs", m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.0674 0.1738 -0.39 0.70
# LangLocNL.Final 0.1184 0.2940 0.40 0.69
# LangLocEN.Init 0.1325 0.2457 0.54 0.59
# LangLocNL.Init 0.3478 0.2939 1.18 0.24
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):LangLocEN.Final 3.30 3.65 2.57 0.07692 .
# s(Time):LangLocNL.Final 6.17 6.45 2.46 0.02314 *
# s(Time):LangLocEN.Init 6.73 6.97 3.67 0.00065 ***
# s(Time):LangLocNL.Init 7.01 7.21 3.35 0.00081 ***
# s(Time):IsTHInit 6.88 7.16 5.27 6.7e-06 ***
# s(Time):IsTHFinal 6.39 6.71 3.30 0.00229 **
# s(Time):IsNLTHInit 4.25 5.52 0.33 0.91196
# s(Time):IsNLTHFinal 2.76 3.24 0.40 0.75891
# s(Time,Speaker):SoundLocT.Final 267.13 377.00 3.07 < 2e-16 ***
# s(Time,Speaker):SoundLocTH.Final 287.27 377.00 4.63 < 2e-16 ***
# s(Time,Speaker):SoundLocT.Init 308.21 378.00 5.80 < 2e-16 ***
# s(Time,Speaker):SoundLocTH.Init 302.84 377.00 5.28 < 2e-16 ***
# s(Time,Word):LangEN 154.14 179.00 16.01 < 2e-16 ***
# s(Time,Word):LangNL 154.79 180.00 19.61 < 2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.593 Deviance explained = 59.8%
# fREML = -49814 Scale est. = 0.34349 n = 126177par(mfrow=c(2,2))
plot(fmk1,select=5,shade=T,rug=F, ylim=c(-0.8,2.0), main="fmk1: English difference initial")
abline(h=0)
plot(fmk1,select=6,shade=T,rug=F, ylim=c(-0.8,2.0), main="fmk1: English difference final")
abline(h=0)
plot(fmk1,select=7,shade=T,rug=F, ylim=c(-1.8,1.0), main="fmk1: Dutch vs. English difference initial")
abline(h=0)
plot(fmk1,select=8,shade=T,rug=F, ylim=c(-1.8,1.0), main="fmk1: Dutch vs. English difference final")
abline(h=0)
8.5 Per language
full$IsENTHInit <- (full$Lang == "EN" & full$Sound == "TH" & full$Loc == "Init")*1
full$IsENTHFinal <- (full$Lang == "EN" & full$Sound == "TH" & full$Loc == "Final")*1
system.time(
fm0 <- bam(Pos ~ LangLoc + s(Time,by=LangLoc) + s(Time,by=IsENTHInit) + s(Time,by=IsENTHFinal) +
s(Time,by=IsNLTHInit) + s(Time,by=IsNLTHFinal) +
s(Time,Speaker,by=SoundLoc,bs="fs",m=1) + s(Time,Word,by=Lang,bs="fs",m=1),
data=full, discrete=TRUE)
)# user system elapsed
# 170.5 232.1 55.4fm0acf <- acf_resid(fm0)
(rhoval <- fm0acf[2]) # autocorrelation at lag 1# 1
# 0.963system.time(
fm1 <- bam(Pos ~ LangLoc + s(Time,by=LangLoc) + s(Time,by=IsENTHInit) + s(Time,by=IsENTHFinal) +
s(Time,by=IsNLTHInit) + s(Time,by=IsNLTHFinal) +
s(Time,Speaker,by=SoundLoc,bs="fs",m=1) + s(Time,Word,by=Lang,bs="fs",m=1),
data=full, discrete=TRUE, rho=rhoval, AR.start=full$start.event)
)# user system elapsed
# 206 181 118acf_resid(fm1)
(smryfm1 <- summary(fm1))#
# Family: gaussian
# Link function: identity
#
# Formula:
# Pos ~ LangLoc + s(Time, by = LangLoc) + s(Time, by = IsENTHInit) +
# s(Time, by = IsENTHFinal) + s(Time, by = IsNLTHInit) + s(Time,
# by = IsNLTHFinal) + s(Time, Speaker, by = SoundLoc, bs = "fs",
# m = 1) + s(Time, Word, by = Lang, bs = "fs", m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.0694 0.1738 -0.40 0.69
# LangLocNL.Final 0.1227 0.2939 0.42 0.68
# LangLocEN.Init 0.1342 0.2456 0.55 0.58
# LangLocNL.Init 0.3509 0.2939 1.19 0.23
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):LangLocEN.Final 2.86 3.20 1.63 0.24651
# s(Time):LangLocNL.Final 6.73 6.97 3.34 0.00175 **
# s(Time):LangLocEN.Init 6.82 7.05 3.61 0.00044 ***
# s(Time):LangLocNL.Init 7.29 7.47 4.47 2.8e-05 ***
# s(Time):IsENTHInit 6.43 6.75 4.92 1.9e-05 ***
# s(Time):IsENTHFinal 5.95 6.30 3.15 0.00326 **
# s(Time):IsNLTHInit 5.19 5.54 2.51 0.08728 .
# s(Time):IsNLTHFinal 3.36 3.68 0.55 0.69526
# s(Time,Speaker):SoundLocT.Final 266.97 378.00 3.05 < 2e-16 ***
# s(Time,Speaker):SoundLocTH.Final 287.19 376.00 4.62 < 2e-16 ***
# s(Time,Speaker):SoundLocT.Init 307.85 378.00 5.79 < 2e-16 ***
# s(Time,Speaker):SoundLocTH.Init 302.62 376.00 5.28 < 2e-16 ***
# s(Time,Word):LangEN 152.99 178.00 16.02 < 2e-16 ***
# s(Time,Word):LangNL 154.62 178.00 19.76 < 2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.593 Deviance explained = 59.8%
# fREML = -49817 Scale est. = 0.34556 n = 126177par(mfrow=c(2,2))
plot(fm1,select=5,shade=T,rug=F, ylim=c(-0.8,2.0), main="fm1: English difference initial")
abline(h=0)
plot(fm1,select=6,shade=T,rug=F, ylim=c(-0.8,2.0), main="fm1: English difference final")
abline(h=0)
plot(fm1,select=7,shade=T,rug=F, ylim=c(-0.8,2.0), main="fm1: Dutch difference initial")
abline(h=0)
plot(fm1,select=8,shade=T,rug=F, ylim=c(-0.8,2.0), main="fm1: Dutch difference final")
abline(h=0)
9 Full model: scaled-\(t\)
system.time(
fmc0s <- bam(Pos ~ LangLoc + s(Time,by=LangLoc) + s(Time,by=IsTHInit) +
s(Time,by=IsTHFinal) + s(Time,by=IsNLTHInit) + s(Time,by=IsNLTHFinal) +
s(Time,Speaker,by=SoundLoc,bs="fs",m=1) + s(Time,Word,by=Lang,bs="fs",m=1),
data=full, discrete=TRUE, family="scat")
)# user system elapsed
# 761 221 664fmc0sacf <- acf_resid(fmc0s)
(rhoval <- fmc0sacf[2]) # autocorrelation at lag 1# 1
# 0.963system.time(
fmc1s <- bam(Pos ~ LangLoc + s(Time,by=LangLoc) + s(Time,by=IsTHInit) +
s(Time,by=IsTHFinal) + s(Time,by=IsNLTHInit) + s(Time,by=IsNLTHFinal) +
s(Time,Speaker,by=SoundLoc,bs="fs",m=1) + s(Time,Word,by=Lang,bs="fs",m=1),
data=full, discrete=TRUE, rho=rhoval, AR.start=full$start.event, family="scat")
)# user system elapsed
# 4266 483 4064acf_resid(fmc1s)
(smryfmc1s <- summary(fmc1s))#
# Family: Scaled t(4.163,0.506)
# Link function: identity
#
# Formula:
# Pos ~ LangLoc + s(Time, by = LangLoc) + s(Time, by = IsTHInit) +
# s(Time, by = IsTHFinal) + s(Time, by = IsNLTHInit) + s(Time,
# by = IsNLTHFinal) + s(Time, Speaker, by = SoundLoc, bs = "fs",
# m = 1) + s(Time, Word, by = Lang, bs = "fs", m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.0734 0.1743 -0.42 0.67
# LangLocNL.Final 0.1445 0.2792 0.52 0.60
# LangLocEN.Init 0.1441 0.2445 0.59 0.56
# LangLocNL.Init 0.3613 0.2773 1.30 0.19
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):LangLocEN.Final 2.90 3.24 1.45 0.27781
# s(Time):LangLocNL.Final 6.73 6.96 3.04 0.00400 **
# s(Time):LangLocEN.Init 7.16 7.35 4.53 2.3e-05 ***
# s(Time):LangLocNL.Init 7.58 7.73 4.60 1.4e-05 ***
# s(Time):IsTHInit 6.53 6.84 4.73 2.6e-05 ***
# s(Time):IsTHFinal 6.69 6.99 3.71 0.00054 ***
# s(Time):IsNLTHInit 2.00 2.00 0.79 0.45527
# s(Time):IsNLTHFinal 3.15 3.44 0.64 0.64260
# s(Time,Speaker):SoundLocT.Final 270.98 377.00 3.29 < 2e-16 ***
# s(Time,Speaker):SoundLocTH.Final 275.88 376.00 4.12 < 2e-16 ***
# s(Time,Speaker):SoundLocT.Init 306.76 376.00 5.71 < 2e-16 ***
# s(Time,Speaker):SoundLocTH.Init 304.95 376.00 5.51 < 2e-16 ***
# s(Time,Word):LangEN 155.64 177.00 18.07 < 2e-16 ***
# s(Time,Word):LangNL 155.71 176.00 21.79 < 2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.58 Deviance explained = 49.3%
# fREML = 42155 Scale est. = 1 n = 126177par(mfrow=c(2,2))
plot(fmc1s,select=5,shade=T,rug=F, ylim=c(-0.8,2.0), main="fm1s: English difference initial")
abline(h=0)
plot(fmc1s,select=6,shade=T,rug=F, ylim=c(-0.8,2.0), main="fm1s: English difference final")
abline(h=0)
plot(fmc1s,select=7,shade=T,rug=F, ylim=c(-1.8,1.0), main="fm1s: Dutch vs. English difference initial")
abline(h=0)
plot(fmc1s,select=8,shade=T,rug=F, ylim=c(-1.8,1.0), main="fm1s: Dutch vs. English difference final")
abline(h=0)
9.1 Model criticism
gam.check(fmc1s)
#
# Method: fREML Optimizer: perf chol
# $grad
# [1] -7.39e-12 4.91e-12 4.70e-13 1.82e-11 -1.36e-12 5.06e-12 -9.41e-04 6.30e-12 1.43e-10
# [10] 1.41e-10 7.91e-11 1.85e-11 3.32e-09 -2.09e-11 1.36e-10 2.03e-11 9.67e-11 7.34e-10
# [19] 1.36e-10 -1.41e-11
#
# $hess
# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
# [1,] 5.99e-01 -3.86e-02 5.21e-07 -4.26e-07 9.92e-09 1.42e-01 1.38e-09 1.01e-01 1.94e-02
# [2,] -3.86e-02 1.44e+00 1.25e-06 -1.05e-06 1.46e-08 2.62e-01 3.29e-09 -8.42e-03 1.29e-01
# [3,] 5.21e-07 1.25e-06 2.03e+00 -1.49e-01 3.01e-01 1.37e-06 -2.63e-05 3.19e-06 6.29e-05
# [4,] -4.26e-07 -1.05e-06 -1.49e-01 2.43e+00 2.60e-01 -6.89e-07 -1.21e-06 -1.82e-06 1.21e-06
# [5,] 9.92e-09 1.46e-08 3.01e-01 2.60e-01 1.90e+00 1.81e-07 -1.50e-04 3.42e-07 2.00e-05
# [6,] 1.42e-01 2.62e-01 1.37e-06 -6.89e-07 1.81e-07 1.49e+00 4.02e-09 -6.77e-02 -2.06e-02
# [7,] 1.38e-09 3.29e-09 -2.63e-05 -1.21e-06 -1.50e-04 4.02e-09 9.39e-04 9.13e-09 2.10e-07
# [8,] 1.01e-01 -8.42e-03 3.19e-06 -1.82e-06 3.42e-07 -6.77e-02 9.13e-09 3.90e-01 -1.08e-03
# [9,] 1.94e-02 1.29e-01 6.29e-05 1.21e-06 2.00e-05 -2.06e-02 2.10e-07 -1.08e-03 8.90e+01
# [10,] 8.05e-04 -7.25e-04 -1.63e-05 -3.31e-07 -5.23e-06 -1.21e-04 -5.75e-08 -4.85e-03 -4.06e-02
# [11,] 8.85e-03 2.10e-02 -6.92e-05 6.23e-05 6.91e-07 1.07e-01 -1.80e-07 -1.14e-02 -3.78e-03
# [12,] 2.66e-04 -6.28e-04 1.31e-05 -1.18e-05 -1.53e-07 -2.00e-04 3.33e-08 1.71e-03 1.04e-05
# [13,] 4.43e-06 7.93e-06 1.57e-01 1.17e-01 -1.01e-01 5.44e-05 1.73e-05 1.05e-04 5.56e-03
# [14,] 6.21e-07 8.15e-07 -2.63e-04 -5.50e-04 9.44e-04 8.57e-06 -1.75e-06 1.63e-05 6.72e-04
# [15,] -1.35e-04 -3.31e-04 6.81e-02 -4.88e-02 7.46e-02 -2.21e-04 8.78e-06 -5.82e-04 6.91e-05
# [16,] -1.51e-04 -3.71e-04 -1.71e-04 -4.96e-04 2.36e-04 -2.48e-04 1.15e-06 -6.53e-04 3.57e-05
# [17,] 2.47e-01 -1.87e-01 9.07e-01 -3.73e-02 2.72e-01 7.23e-01 1.07e-04 1.55e-02 2.45e-02
# [18,] -1.45e-03 -2.04e-05 3.00e-05 -6.59e-06 -2.32e-05 1.43e-04 -2.28e-08 5.66e-05 -6.42e-03
# [19,] -1.91e-02 1.44e+00 -3.95e-02 2.49e-01 -3.62e-01 1.09e-01 3.38e-05 -1.27e-01 9.04e-02
# [20,] -1.44e-05 2.06e-04 -1.85e-05 1.13e-05 -3.32e-05 -9.81e-05 -9.79e-08 -3.09e-04 -1.13e-02
# [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18]
# [1,] 8.05e-04 8.85e-03 2.66e-04 4.43e-06 6.21e-07 -1.35e-04 -1.51e-04 2.47e-01 -1.45e-03
# [2,] -7.25e-04 2.10e-02 -6.28e-04 7.93e-06 8.15e-07 -3.31e-04 -3.71e-04 -1.87e-01 -2.04e-05
# [3,] -1.63e-05 -6.92e-05 1.31e-05 1.57e-01 -2.63e-04 6.81e-02 -1.71e-04 9.07e-01 3.00e-05
# [4,] -3.31e-07 6.23e-05 -1.18e-05 1.17e-01 -5.50e-04 -4.88e-02 -4.96e-04 -3.73e-02 -6.59e-06
# [5,] -5.23e-06 6.91e-07 -1.53e-07 -1.01e-01 9.44e-04 7.46e-02 2.36e-04 2.72e-01 -2.32e-05
# [6,] -1.21e-04 1.07e-01 -2.00e-04 5.44e-05 8.57e-06 -2.21e-04 -2.48e-04 7.23e-01 1.43e-04
# [7,] -5.75e-08 -1.80e-07 3.33e-08 1.73e-05 -1.75e-06 8.78e-06 1.15e-06 1.07e-04 -2.28e-08
# [8,] -4.85e-03 -1.14e-02 1.71e-03 1.05e-04 1.63e-05 -5.82e-04 -6.53e-04 1.55e-02 5.66e-05
# [9,] -4.06e-02 -3.78e-03 1.04e-05 5.56e-03 6.72e-04 6.91e-05 3.57e-05 2.45e-02 -6.42e-03
# [10,] 1.31e+01 1.05e-04 -2.01e-05 -1.70e-03 -5.60e-04 -6.81e-06 3.21e-06 -1.11e-03 -2.65e-03
# [11,] 1.05e-04 8.88e+01 4.96e-02 -2.73e-05 7.84e-06 1.96e-02 2.20e-02 2.74e-02 -1.32e-03
# [12,] -2.01e-05 4.96e-02 1.65e+01 4.46e-06 -1.60e-06 -3.80e-03 -4.23e-03 -2.34e-03 -7.75e-03
# [13,] -1.70e-03 -2.73e-05 4.46e-06 1.16e+02 -7.78e-01 4.52e-04 5.72e-04 4.71e-02 6.62e-04
# [14,] -5.60e-04 7.84e-06 -1.60e-06 -7.78e-01 1.21e+01 -4.79e-04 -2.71e-05 -2.19e-03 7.43e-04
# [15,] -6.81e-06 1.96e-02 -3.80e-03 4.52e-04 -4.79e-04 1.10e+02 -3.89e-01 4.66e-01 1.04e-03
# [16,] 3.21e-06 2.20e-02 -4.23e-03 5.72e-04 -2.71e-05 -3.89e-01 1.60e+01 -1.10e-03 -1.13e-03
# [17,] -1.11e-03 2.74e-02 -2.34e-03 4.71e-02 -2.19e-03 4.66e-01 -1.10e-03 6.32e+01 -1.57e-01
# [18,] -2.65e-03 -1.32e-03 -7.75e-03 6.62e-04 7.43e-04 1.04e-03 -1.13e-03 -1.57e-01 7.58e+00
# [19,] -4.22e-03 -8.84e-03 -5.02e-03 8.86e-01 -1.21e-02 4.57e-01 8.05e-04 1.85e-01 -1.32e-04
# [20,] -2.93e-02 -3.13e-03 1.58e-02 -2.83e-03 -6.55e-03 -2.80e-03 4.34e-03 1.27e-05 -1.10e-07
# [,19] [,20]
# [1,] -1.91e-02 -1.44e-05
# [2,] 1.44e+00 2.06e-04
# [3,] -3.95e-02 -1.85e-05
# [4,] 2.49e-01 1.13e-05
# [5,] -3.62e-01 -3.32e-05
# [6,] 1.09e-01 -9.81e-05
# [7,] 3.38e-05 -9.79e-08
# [8,] -1.27e-01 -3.09e-04
# [9,] 9.04e-02 -1.13e-02
# [10,] -4.22e-03 -2.93e-02
# [11,] -8.84e-03 -3.13e-03
# [12,] -5.02e-03 1.58e-02
# [13,] 8.86e-01 -2.83e-03
# [14,] -1.21e-02 -6.55e-03
# [15,] 4.57e-01 -2.80e-03
# [16,] 8.05e-04 4.34e-03
# [17,] 1.85e-01 1.27e-05
# [18,] -1.32e-04 -1.10e-07
# [19,] 6.22e+01 -1.61e-01
# [20,] -1.61e-01 7.68e+00
#
# Model rank = 1952 / 1952
#
# 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):LangLocEN.Final 9.00 2.90 0.99 0.15
# s(Time):LangLocNL.Final 9.00 6.73 0.99 0.14
# s(Time):LangLocEN.Init 9.00 7.16 0.99 0.17
# s(Time):LangLocNL.Init 9.00 7.58 0.99 0.20
# s(Time):IsTHInit 10.00 6.53 0.99 0.19
# s(Time):IsTHFinal 10.00 6.69 0.99 0.18
# s(Time):IsNLTHInit 10.00 2.00 0.99 0.17
# s(Time):IsNLTHFinal 10.00 3.15 0.99 0.20
# s(Time,Speaker):SoundLocT.Final 378.00 270.98 0.99 0.23
# s(Time,Speaker):SoundLocTH.Final 378.00 275.88 0.99 0.17
# s(Time,Speaker):SoundLocT.Init 378.00 306.76 0.99 0.14
# s(Time,Speaker):SoundLocTH.Init 378.00 304.95 0.99 0.16
# s(Time,Word):LangEN 180.00 155.64 0.99 0.20
# s(Time,Word):LangNL 180.00 155.71 0.99 0.179.2 Per language
system.time(
fm0s <- bam(Pos ~ LangLoc + s(Time,by=LangLoc) + s(Time,by=IsTHInit) +
s(Time,by=IsTHFinal) + s(Time,by=IsNLTHInit) + s(Time,by=IsNLTHFinal) +
s(Time,Speaker,by=SoundLoc,bs="fs",m=1) + s(Time,Word,by=Lang,bs="fs",m=1),
data=full, discrete=TRUE, family="scat")
)# user system elapsed
# 739 189 654fm0sacf <- acf_resid(fm0s)
(rhoval <- fm0sacf[2]) # autocorrelation at lag 1# 1
# 0.963system.time(
fm1s <- bam(Pos ~ LangLoc + s(Time,by=LangLoc) + s(Time,by=IsENTHInit) +
s(Time,by=IsENTHFinal) + s(Time,by=IsNLTHInit) + s(Time,by=IsNLTHFinal) +
s(Time,Speaker,by=SoundLoc,bs="fs",m=1) + s(Time,Word,by=Lang,bs="fs",m=1),
data=full, discrete=TRUE, rho=rhoval, AR.start=full$start.event, family="scat")
)# user system elapsed
# 4831 581 4581acf_resid(fm1s)
(smryfm1s <- summary(fm1s)) #
# Family: Scaled t(4.165,0.506)
# Link function: identity
#
# Formula:
# Pos ~ LangLoc + s(Time, by = LangLoc) + s(Time, by = IsENTHInit) +
# s(Time, by = IsENTHFinal) + s(Time, by = IsNLTHInit) + s(Time,
# by = IsNLTHFinal) + s(Time, Speaker, by = SoundLoc, bs = "fs",
# m = 1) + s(Time, Word, by = Lang, bs = "fs", m = 1)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.0749 0.1743 -0.43 0.67
# LangLocNL.Final 0.1476 0.2793 0.53 0.60
# LangLocEN.Init 0.1453 0.2445 0.59 0.55
# LangLocNL.Init 0.3647 0.2774 1.31 0.19
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(Time):LangLocEN.Final 2.62 2.95 0.98 0.4335
# s(Time):LangLocNL.Final 7.34 7.52 4.66 1.5e-05 ***
# s(Time):LangLocEN.Init 7.16 7.35 4.38 3.4e-05 ***
# s(Time):LangLocNL.Init 7.74 7.87 5.38 2.8e-06 ***
# s(Time):IsENTHInit 6.18 6.52 4.48 7.4e-05 ***
# s(Time):IsENTHFinal 6.16 6.50 3.36 0.0016 **
# s(Time):IsNLTHInit 4.73 5.08 2.12 0.1008
# s(Time):IsNLTHFinal 2.00 2.00 0.75 0.4713
# s(Time,Speaker):SoundLocT.Final 271.03 377.00 3.29 < 2e-16 ***
# s(Time,Speaker):SoundLocTH.Final 276.08 376.00 4.13 < 2e-16 ***
# s(Time,Speaker):SoundLocT.Init 306.56 376.00 5.70 < 2e-16 ***
# s(Time,Speaker):SoundLocTH.Init 304.83 376.00 5.51 < 2e-16 ***
# s(Time,Word):LangEN 154.55 177.00 18.05 < 2e-16 ***
# s(Time,Word):LangNL 156.30 177.00 21.76 < 2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# R-sq.(adj) = 0.58 Deviance explained = 49.3%
# fREML = 42149 Scale est. = 1 n = 126177par(mfrow=c(2,2))
plot(fm1s,select=5,shade=T,rug=F, ylim=c(-0.8,2.0), main="fm1s: English difference initial")
abline(h=0)
plot(fm1s,select=6,shade=T,rug=F, ylim=c(-0.8,2.0), main="fm1s: English difference final")
abline(h=0)
plot(fm1s,select=7,shade=T,rug=F, ylim=c(-0.8,2.0), main="fm1s: Dutch difference initial")
abline(h=0)
plot(fm1s,select=8,shade=T,rug=F, ylim=c(-0.8,2.0), main="fm1s: Dutch difference final")
abline(h=0)
gam.check(fm1s)10 Speedup with nthreads
# Using 2 processors
system.time(
tmp <- bam(Pos ~ LangLoc + s(Time,by=LangLoc) + s(Time,by=IsTHInit) +
s(Time,by=IsTHFinal) + s(Time,by=IsNLTHInit) + s(Time,by=IsNLTHFinal) +
s(Time,Speaker,by=SoundLoc,bs="fs",m=1) + s(Time,Word,by=Lang,bs="fs",m=1),
data=full, discrete=TRUE, rho=rhoval, AR.start=full$start.event, family="scat",
nthreads=2)
)# user system elapsed
# 4678 1060 2153# Using 4 processors
system.time(
tmp <- bam(Pos ~ LangLoc + s(Time,by=LangLoc) + s(Time,by=IsTHInit) +
s(Time,by=IsTHFinal) + s(Time,by=IsNLTHInit) + s(Time,by=IsNLTHFinal) +
s(Time,Speaker,by=SoundLoc,bs="fs",m=1) + s(Time,Word,by=Lang,bs="fs",m=1),
data=full, discrete=TRUE, rho=rhoval, AR.start=full$start.event, family="scat",
nthreads=4)
)# user system elapsed
# 5493 1837 1273# Using 8 processors
system.time(
tmp <- bam(Pos ~ LangLoc + s(Time,by=LangLoc) + s(Time,by=IsTHInit) +
s(Time,by=IsTHFinal) + s(Time,by=IsNLTHInit) + s(Time,by=IsNLTHFinal) +
s(Time,Speaker,by=SoundLoc,bs="fs",m=1) + s(Time,Word,by=Lang,bs="fs",m=1),
data=full, discrete=TRUE, rho=rhoval, AR.start=full$start.event, family="scat",
nthreads=8)
)# user system elapsed
# 7051 4040 869# Using 16 processors
system.time(
tmp <- bam(Pos ~ LangLoc + s(Time,by=LangLoc) + s(Time,by=IsTHInit) +
s(Time,by=IsTHFinal) + s(Time,by=IsNLTHInit) + s(Time,by=IsNLTHFinal) +
s(Time,Speaker,by=SoundLoc,bs="fs",m=1) + s(Time,Word,by=Lang,bs="fs",m=1),
data=full, discrete=TRUE, rho=rhoval, AR.start=full$start.event, family="scat",
nthreads=16)
)# user system elapsed
# 6553 1532 468# Using 32 processors
system.time(
tmp <- bam(Pos ~ LangLoc + s(Time,by=LangLoc) + s(Time,by=IsTHInit) +
s(Time,by=IsTHFinal) + s(Time,by=IsNLTHInit) + s(Time,by=IsNLTHFinal) +
s(Time,Speaker,by=SoundLoc,bs="fs",m=1) + s(Time,Word,by=Lang,bs="fs",m=1),
data=full, discrete=TRUE, rho=rhoval, AR.start=full$start.event, family="scat",
nthreads=32)
)# user system elapsed
# 7309 2042 333# Linear mixed model fit by REML ['lmerMod']
# Formula: Pos ~ Lang * Sound + (1 + Lang | Word) + (1 + Sound | Speaker)
# Data: agg
#
# REML criterion at convergence: 2148
#
# Scaled residuals:
# Min 1Q Median 3Q Max
# -4.735 -0.560 0.023 0.575 3.677
#
# Random effects:
# Groups Name Variance Std.Dev. Corr
# Speaker (Intercept) 0.0658 0.256
# SoundTH 0.2150 0.464 -0.39
# Word (Intercept) 0.1375 0.371
# LangNL 0.0315 0.178 0.64
# Residual 0.1814 0.426
# Number of obs: 1618, groups: Speaker, 42; Word, 20
#
# Fixed effects:
# Estimate Std. Error t value
# (Intercept) -0.00289 0.13097 -0.02
# LangNL 0.17649 0.10182 1.73
# SoundTH 0.69663 0.19520 3.57
# LangNL:SoundTH -0.37260 0.16939 -2.20
#
# Correlation of Fixed Effects:
# (Intr) LangNL SondTH
# LangNL 0.061
# SoundTH -0.637 -0.085
# LngNL:SndTH -0.076 -0.496 -0.065
11 Replication
To replicate the analysis presented above, you can just copy the following lines to the most recent version of R. Please note that you first need to install Pandoc.
download.file("http://www.let.rug.nl/wieling/Tutorial/analysis.Rmd", "analysis.Rmd")
if (length(setdiff("rmarkdown", rownames(installed.packages()))) > 0) {
install.packages("rmarkdown")
}
library(rmarkdown)
render("analysis.Rmd") # generates html file with results
browseURL(paste("file://", file.path(getwd(),"analysis.html"), sep="")) # shows result