During the lecture we have looked at the effect of Correctness and the effect of Age of Arrival. In this lab session we will look at another additional factor: Structure. Research question: Is there an effect of the distance between determiner and noun in the violation? We compare the DN (determiner-noun, e.g. der Garten/*das Garten) structure with the DAN (determiner-adjective-noun, e.g. das frische Gras/*der frische Gras) structure.
Load the required packages and download the required files.
# test if the packages (mgcv, car and itsadug) are installed
# if not, install them
packages <- c("mgcv","car","itsadug")
if (length(setdiff(packages, rownames(installed.packages()))) > 0) {
install.packages(setdiff(packages, rownames(installed.packages())))
}
# load the packages
library(mgcv)
library(itsadug)
infoMessages("on") # set info messages of itsadug on
library(car)
# display version information
R.version.string## [1] "R version 3.2.1 (2015-06-18)"packageVersion('mgcv')## [1] '1.8.6'packageVersion('car')## [1] '2.0.25'packageVersion('itsadug')## [1] '1.0.1'# download and load the following files containing the dataset and some functions:
download.file('www.sfs.uni-tuebingen.de/~jvanrij/LSA2015/Data/dat.rda', 'dat.rda')
load('dat.rda') # This dataset is a subset of 15 subjectsLook at the data and sort it (necessary for correcting for autocorrelation).
head(dat)
str(dat)
summary(dat)
dim(dat)
dat = start_event(dat, event=c("Subject","TrialNr","Roi"))
head(dat)We would like to know whether the P600 effect is dependent on the structure of the sentence (Determiner-Noun: DN vs. Determiner-Adjective-Noun: DAN). Remember that the size of the P600 is determined by the difference between incorrect and correct. Consequently, we are interested in assessing if that difference varies between the levels of Structure, and we will make a model where we look at the interaction between Correctness and Structure. Note that to prevent lengthy computations during this lab session, we will not take into account factor smooths per word, but only per subject.
dat$CorStruct = interaction(dat$Correctness,dat$Structure) # create a new variable with 4 levels
levels(dat$CorStruct) # show levels
# also create a new predictor to model variability per subject
dat$SubjectCorStruct = interaction(dat$SubjectCor,dat$Structure)
# determine rho value:
m1_forrho <- bam(uV ~ s(Time,by=CorStruct) + CorStruct + s(Time,SubjectCorStruct,bs='fs',m=1),
data=dat)
rhoval = acf(resid(m1_forrho))$acf[2]rhoval
# rerun model with rhoval, which we keep constant for the rest
# of the lecture for simplicity.
m1 <- bam(uV ~ s(Time,by=CorStruct) + CorStruct + s(Time,SubjectCorStruct,bs='fs',m=1),
data=dat, rho = rhoval, AR.start=start.event)
summary(m1)# plot the separate smooths
par(mfrow=c(1,3))
plot_smooth(m1,"Time",plot_all="CorStruct",rm.ranef=T,eegAxis=T,rug=F)
# plot the differences between correct and incorrect for the two cases (DN and DAN)
plot_diff(m1, "Time", comp=list(CorStruct=c("incor.DAN","cor.DAN")), rm.ranef=T, eegAxis=T)
plot_diff(m1, "Time", comp=list(CorStruct=c("incor.DN","cor.DN")), rm.ranef=T, eegAxis=T)For this exercise, please assess if the distinction between DAN and DN (stored in the predictor Structure) is necessary. Put your code in the following code block, generate the output html file (see below for instructions), and put your conclusion below.
# Your code here.What do you conclude?
In the following, we will investigate how the DN vs. DAN patterns vary depending on the AoArr of the participants. Below we will fit the model and show the commands to visualize these difference surfaces. Note that we only use a subset of the data here (N = 15), so there are not so many different values of AoArr.
m2 = bam(uV ~ te(Time,AoArr,by=CorStruct) + CorStruct + s(Time,SubjectCorStruct,bs='fs',m=1), data=dat,
rho=rhoval, AR.start=start.event)
summary(m2)
# plot the differences between correct and incorrect for the two cases (DN and DAN)
par(mfrow=c(1,2))
plot_diff2(m2, view=c("Time","AoArr"), comp=list(CorStruct=c("incor.DAN","cor.DAN")), rm.ranef=T)
fadeRug(dat$Time,dat$AoArr) # make areas lighter where there is no data
plot_diff2(m2, view=c("Time","AoArr"), comp=list(CorStruct=c("incor.DN","cor.DN")), rm.ranef=T)
fadeRug(dat$Time,dat$AoArr) # make areas lighter where there is no dataFor this exercise, please assess if a simpler model would be a better choice than sticking with m2. First assess if a model with a smooth over Time and a smooth over AoArr separated by Corstruct would be better. Then proceed with assessing if a model with only one smooth would be a better option (test both smooths separately). Note that each model will take up to 5 minutes to run.
# Your code here.What do your conclude?
Model criticism using model m2 shows the familiar pattern:
par(mfrow=c(1,2))
qqp(resid_gam(m2)) # from library(car), combines qqnorm and qqline
hist(resid_gam(m2))check_resid(m2) # from library(itsadug)The code to run the scaled-t analysis is shown below. However, as this analysis is very time consuming (taking many hours), do not run it here. We therefore only show the commands for reference.
download.file('http://www.let.rug.nl/wieling/LSA/bam.art.fit.R', 'bam.art.fit.R')
source('bam.art.fit.R')
g0 <- gam(uV ~ 1, data=dat, method='REML', family=scat(rho=rhoval, AR.start=dat$SeqStart))
g1 <- gam(uV ~ te(Time,AoArr,by=CorStruct) + CorStruct, data=dat, method='REML',
family=scat(rho=rhoval, AR.start=dat$SeqStart))
theta0 = g0$family$getTheta(TRUE)
theta1 = g1$family$getTheta(TRUE)
theta0 - theta1 # check if difference is small
m3 = bam(uV ~ te(Time,Aoarr,by=CorStruct) + CorStruct + s(Time,SubjectCorStruct,bs='fs',m=1), data=dat
method='fREML', family=art(theta=theta1,rho=rhoval), AR.start=SeqStart)There are two other types of predictors which we did not cover. These two types allow us to investigate difference curves directly (rather than using model comparison). The two predictors are binary predictors and ordered factors and model the difference between different levels of a factorial predictor directly as a non-linear patter. You can find more information about these types in another online lecture of Martijn (using the same data covered here) available at http://www.let.rug.nl/~wieling/statscourse/lecture4/presentation.pdf, and an online tutorial of Jacolien available at http://www.sfs.uni-tuebingen.de/~jvanrij/Tutorial/GAMM.html. Some example code is shown below, which you can run yourself.
# Our standard approach, including intercept difference
m <- bam(uV ~ s(Time,by=Correctness) + Correctness + s(Time,Subject,bs='fs',m=1),
data=dat, rho = rhoval, AR.start=start.event)
summary(m)
# A binary curve is only active when the binary predictor equals 1:
# therefore it models the difference between the reference smooth: s(Time)
# and the other level. Note that a binary curve is not centered,
# so an intercept difference is not necessary
dat$IsIncorrect = (dat$Correctness=='incor')*1 # 0 for correct, 1 for incorrect
m_binary <- bam(uV ~ s(Time) + s(Time,by=IsIncorrect) + s(Time,Subject,bs='fs',m=1),
data=dat, rho = rhoval, AR.start=start.event)
summary(m_binary)
# An ordered factor is similar as a binary curve, but it can also contrast
# factors with more than two levels directly and needs a separate intercept
# difference, as these smooths are centered.
dat$CorrectnessO = as.ordered(dat$Correctness)
contrasts(dat$CorrectnessO) = 'contr.treatment' # contrast treatment
m_ordfac <- bam(uV ~ s(Time) + s(Time,by=CorrectnessO) + CorrectnessO +
s(Time,Subject,bs='fs',m=1), data=dat, rho = rhoval, AR.start=start.event)
summary(m_ordfac)
par(mfrow=c(1,3))
plot_diff(m,view='Time',comp=list(Correctness=c('incor','cor')),rm.ranef=T,
eegAxis=T,ylim=c(-6,6), main='m: calculated difference')
plot(m_binary,select=2, rug=F, shade=T,ylim=c(6,-6), main='difference curve (binary: not centered)')
abline(h=0)
plot(m_ordfac,select=2, rug=F, shade=T,ylim=c(6,-6), main='difference curve (ordered factor: centered)')
abline(h=0)Given that generalized additive models have only relatively recently been applied to analyzing time series data in linguistics, the method will be unfamiliar to many potential reviewers. Consequently, when using GAMs, you need to put some more effort in describing the method than if you would use a more traditional (but less adequate!) approach. Fortunately, more and more publications are coming out which have used GAMs in their analysis (and provide a detailed explanation of the method) to which you can refer. Some examples are given below:
To generate the output of the analysis presented above, first install Pandoc. Then copy the following lines to the most recent version of R. Note that you need to remove the results="hide" and the fig.show="hide" parts from the r blocks to see the output.
# install rmarkdown package if not installed
if(!"rmarkdown" %in% rownames(installed.packages())) {
install.packages("rmarkdown")
}
library(rmarkdown) # load rmarkdown package
# download original file if not already exists (to prevent overwriting)
if (!file.exists('analysisEEG.Rmd')) {
download.file('http://www.let.rug.nl/wieling/LSA/analysisEEG.Rmd', 'analysisEEG.Rmd')
}
# generate output
render('analysisEEG.Rmd') # generates html file with results
# view output in browser
browseURL(paste('file://', file.path(getwd(),'analysisEEG.html'), sep='')) # shows result