Disclaimer: This short overview of a GAMM analysis is somewhat outdated, because some newer features of the packages mcgv and itsadug are not yet included. Note, however, that the basics (covered here) are still the same. I'm planning to update this summary, but in the meantime you may want to look at more recent analysis scripts here.

Data and R packages

To give an overview of a GAMM analysis, we use the simulation data simdat of the package itsadug. First load the packages mgcv and itsadug.

How to install the R packages

Both packages could be installed from CRAN using the commandline in R:

# installing packages from CRAN:
install.packages('mgcv')
install.packages('itsadug')

To install the newest version of itsadug directly from GitHub:

# installing the package devtools from CRAN:
install.packages('devtools')
# load the package devtools:
library(devtools)
# install the newest version from itsadug:
devtools::install_github("vr-vr/itsadug", build_vignettes=TRUE)

Load packages

# itsadug automatically loads mgcv:
library(itsadug)
## Loading required package: mgcv
## Loading required package: nlme
## This is mgcv 1.8-3. For overview type 'help("mgcv-package")'.
## Loaded package itsadug 0.92 (see 'help("itsadug")' ).

As the packages mgcv and itsadug are often updated, please check whether you have the latest version:

packageVersion('mgcv')
## [1] '1.8.3'
packageVersion('itsadug')
## [1] '0.92'

Data

In this tutorial we use the simulated data simdat from itsadug:

data(simdat)
# inspect data:
head(simdat)
str(simdat)
# see for more info:
help(simdat)

However, you could also load your own data in R. If you have your data in Microsoft Excel, the easiest way is to export the data as .csv file:

File \(\rightarrow\) Save As, then choose ‘Comma Separated Values (.csv)’.

Then load your data in R as follows:

dat <- read.table(file="mydata.csv", sep=",", header=TRUE)

If you did not use column names set the argument header to FALSE. The same function could be used to read a text file. The argument sep specifies the character being used to separate the columns. In text files this is often a tab, which could be specified with sep="\t".

Analysis

For the current example, our question is whether adults and children reaction differently in their ‘Y’ response on the stimulus onset asynchrony, captured by the column Condition. Condition is manipulated between participants. The data set does not contain information about the different items, so for simplicity we assume that the experiment contained just one single item.

gam() or bam()

There are two functions for implementing a GAMM model: gam() and bam(). There are largely similar. The most important difference is that bam() is optimized for big data sets (e.g., more than 10.000 data points)

nrow(simdat)
## [1] 75600

Smooths terms

To model a potentially nonlinear smooth or surface, three different smooth functions are available:
  • s() :

    • for modeling a 1-dimensional smooth, or
    • for modeling isotropic interactions (variables are measured in same units and on same scale)
  • te(): for modeling 2- or n-dimensional interaction surfaces of variables that are not isotropic (but see info about d parameter below). Includes ‘main’ effects.

  • ti(): for modeling 2- or n-dimensional interaction surfaces that do not include the ‘main effects’.

Example models:

# 1-dimensional smooth of Time:
m1 <- bam(Y ~ s(Time), data=simdat)
# interaction with non-isotropicterms:
m2 <- bam(Y ~ te(Time, Trial), data=simdat)
# decompose the same interaction in main effects + interaction:
# (note: include main effects)
m3 <- bam(Y ~ s(Time) + s(Trial) + ti(Time, Trial), data=simdat)

Parameters of smooth functions

The smooth functions have several parameters that could be set to change their behavior. The most often-used parameters are listed here:

  • k : number of ‘knots’. This parameter determines the upper bound of the number of underlying base functions being used to build up the curve. Thus, this parameter constraints the wigglyness of a smooth, or - as a metaphor - the number of bowpoints of a curve.
    • Note that the model will base the number of base functions (reflected in the edf of the summary) on the data with the setting for k as upper bound. By default, the value of k for s() is around 9, and for te() and ti() 5 per dimension. Importantly, the value of k should be at most one less than the number of unique data points, otherwise it will fit the density of that predictor.
  • d : for specifiying that predictors in the interaction are on the same scale or dimension (only used in te() and ti()). For example, in te(Time, width, height, d=c(1,2)), with width and height reflecting the picture size measured in pixels, we specify that Time is on a different dimension than the next two variables. By default, the value would be d=c(1,1,1) in this case.

  • bs: specifies the type of underlying base functios. For s() this defaults to "tp" (thin plate regression spline) and for te() and ti() this defaults to "cr" (cubic regression spline). For random intercepts and linear random slopes use bs="re", but for random smooths use bs="fs" (see below).

Different types of interactions

For modeling potentially different trends over time for different age groups (CHildren, Adults), the by argument in the smooth functions could be used.

How to set up the interaction depends on the type of grouping predictor:

  • with factor include intercept difference: Group + s(Time, by=Group).

  • with ordered factor include intercept difference and reference smooth: Group + s(Time) + s(Time, by=Group).

  • with binary predictor include reference smooth: s(Time) + s(Time, by=IsGroupChildren).

An intercept term needs to be included when using factors and ordered factors, because smooths with factors are centered. Including the intercept terms provides more flexibility for the model to fit overall intercept differences between factor levels, and avoids artifacts due to the centering constraints. Binary predictors are modeled as separate smooths, and hence are not centered around zero. As a result, binary predictors may yield different estimates that factors and ordered factors.

Example factor:

# make predictor as factor:
simdat$Group <- as.factor(simdat$Group)
# inspect contrasts:
contrasts(simdat$Group)
##          Adults
## Children      0
## Adults        1
# change reference level: Adults is reference
simdat$Group <- relevel(simdat$Group, ref="Adults")
# inspect contrasts:
contrasts(simdat$Group)
##          Children
## Adults          0
## Children        1
# run model:
m.factor <- bam(Y ~ Group + s(Time, by=Group), data=simdat)
# inspect summary:
summary(m.factor)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ Group + s(Time, by = Group)
## 
## Parametric coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    5.26277    0.03241  162.39   <2e-16 ***
## GroupChildren -3.18626    0.04583  -69.52   <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):GroupAdults   8.586  8.951 5367  <2e-16 ***
## s(Time):GroupChildren 8.124  8.793 2924  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.509   Deviance explained =   51%
## fREML = 2.4646e+05  Scale est. = 39.69     n = 75600

The summary of m.factor shows:

  • Intercept adjustment for Children in comparison with Adults in parametric summary.

  • The smooth terms summary shows two smooths that are significantly different from 0. We cannot conclude from the summary that the two curves are different from each other.

Example ordered factor:

simdat$OFGroup <- as.factor(simdat$Group)
# change factor to ordered factor:
simdat$OFGroup <- as.ordered(simdat$OFGroup)
# change contrast to treatment coding (difference curves)
contrasts(simdat$OFGroup) <- 'contr.treatment'
# Inspect contrasts:
contrasts(simdat$OFGroup)
##          Children
## Adults          0
## Children        1
# run model:
m.orderedfactor <- bam(Y ~ OFGroup + s(Time) + s(Time, by=OFGroup), data=simdat)
# inspect summary:
summary(m.orderedfactor)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ OFGroup + s(Time) + s(Time, by = OFGroup)
## 
## Parametric coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      5.26276    0.03241  162.39   <2e-16 ***
## OFGroupChildren -3.18621    0.04583  -69.52   <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)                 8.662  8.926 5382  <2e-16 ***
## s(Time):OFGroupChildren 7.408  8.347  595  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.509   Deviance explained =   51%
## fREML = 2.4646e+05  Scale est. = 39.689    n = 75600

The summary of m.orderedfactor shows:

  • Intercept adjustment for Children in comparison with Adults in parametric summary.

  • The smooth terms summary shows two smooths that are significantly different from 0: the reference smooth (for the Adults) and the difference smooth (Children - Adults). When the difference smooth is significantly different from 0, we can conclude that the two groups follow a different trend over Time.

Example binary:

# create binary predictor:
simdat$IsChildren <- with(simdat, ifelse(Group=="Children", 1, 0))
# run model:
m.binary <- bam(Y ~ s(Time) + s(Time, by=IsChildren), data=simdat)
# inspect summary:
summary(m.binary)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ s(Time) + s(Time, by = IsChildren)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.26264    0.03241   162.4   <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)            8.663  8.926 5382  <2e-16 ***
## s(Time):IsChildren 8.408  9.347 1040  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.509   Deviance explained =   51%
## fREML = 2.4646e+05  Scale est. = 39.689    n = 75600

The summary of m.binary shows:

  • The smooth terms summary shows two smooths that are significantly different from 0: the reference smooth (for the Adults) and the difference smooth (Children - Adults). When the difference smooth is significantly different from 0, we can conclude that the two groups follow a different trend over Time. However, we do not know whether the trend over time is different, or whether the groups are overall showing a difference (intercept difference).

Random effects

Three different types of random effects are distinghuished when using GAMMs:

  • random intercepts adjust the height of other modelterms with a constant value: s(Subject, bs="re").

  • random slopes adjust the slope ofthe trend of a numeric predictor: s(Subject, Time, bs="re").

  • random smooths adjust the trend of a numeric predictor in a nonlinear way: s(Time, Subject, bs="fs", m=1).

Notes:

  • Random intercepts and random slopes could be combined, but the random smooths already include random intercepts and random slope effects.

  • The argument m=1 sets a heavier penalty for the smooth moving away from 0, causing shrinkage to the mean.

Example of random intercept:

# first make sure Subject is a factor:
simdat$Subject <- as.factor(simdat$Subject)
# run model with random intercept:
m1 <- bam(Y ~ Group + s(Time, by=Group)
      + s(Subject, bs="re"), data=simdat)
summary(m1)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ Group + s(Time, by = Group) + s(Subject, bs = "re")
## 
## Parametric coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     5.2614     0.6903   7.621 2.54e-14 ***
## GroupChildren  -3.1895     0.9763  -3.267  0.00109 ** 
## ---
## 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):GroupAdults    8.66  8.966 6730  <2e-16 ***
## s(Time):GroupChildren  8.26  8.850 3649  <2e-16 ***
## s(Subject)            33.94 34.000  569  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.609   Deviance explained =   61%
## fREML = 2.3795e+05  Scale est. = 31.601    n = 75600

Plot of distribution of random intercepts:

plot(m1, select=3)

Plot the summed effects for the adults without random effects (Left), and for two specific subjects (Right):

# specify two plot panels:
par(mfrow=c(1,2), cex=1.1)
# summed effects without random effects:
# This is actually: intercept + s(Time):GroupAdults
plot_smooth(m1, view="Time", cond=list(Group="Adults"), 
      rm.ranef=TRUE, main="intercept+s(Time):GroupAdults", rug=FALSE)

# Plot summed effects for Subject "a01":
plot_smooth(m1, view="Time", cond=list(Group="Adults", Subject="a01"), 
      main="... + s(Subject)", col='red',
      ylim=c(-15,20), rug=FALSE)
# Add Subject "a18":
plot_smooth(m1, view="Time", cond=list(Group="Adults", Subject="a18"), 
      add=TRUE, col='blue')

Example of random slope:

# first make sure Subject is a factor:
simdat$Subject <- as.factor(simdat$Subject)
# run model with random intercept:
m2 <- bam(Y ~ Group + s(Time, by=Group)
      + s(Subject, Time, bs="re"), data=simdat)
summary(m2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ Group + s(Time, by = Group) + s(Subject, Time, bs = "re")
## 
## Parametric coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     5.2613     0.2345  22.436   <2e-16 ***
## GroupChildren  -3.1895     0.3316  -9.617   <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):GroupAdults    8.597  8.953 3338.56  <2e-16 ***
## s(Time):GroupChildren  8.144  8.802  729.72  <2e-16 ***
## s(Subject,Time)       33.527 34.000   70.92  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.525   Deviance explained = 52.5%
## fREML = 2.4534e+05  Scale est. = 38.463    n = 75600

Plot of distribution of random intercepts:

plot(m2, select=3)

Plot the summed effects for the adults without random effects (Left), and for two specific subjects (Right):

# specify two plot panels:
par(mfrow=c(1,2), cex=1.1)
# summed effects without random effects:
# This is actually: intercept + s(Time):GroupAdults
plot_smooth(m2, view="Time", cond=list(Group="Adults"), 
      rm.ranef=TRUE, main="intercept+s(Time):GroupAdults", rug=FALSE)

# Plot summed effects for Subject "a01":
plot_smooth(m2, view="Time", cond=list(Group="Adults", Subject="a01"), 
      main="... + s(Subject)", col='red',
      ylim=c(-15,20), rug=FALSE)
# Add Subject "a04":
plot_smooth(m2, view="Time", cond=list(Group="Adults", Subject="a18"), 
      add=TRUE, col='blue')

Example of random intercepts and slopes:

# first make sure Subject is a factor:
simdat$Subject <- as.factor(simdat$Subject)
# run model with random intercept:
m3 <- bam(Y ~ Group + s(Time, by=Group)
      + s(Subject, bs="re") + s(Subject, Time, bs="re"), data=simdat)
summary(m3)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ Group + s(Time, by = Group) + s(Subject, bs = "re") + s(Subject, 
##     Time, bs = "re")
## 
## Parametric coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)      5.261      2.951   1.783   0.0746 .
## GroupChildren   -3.189      4.174  -0.764   0.4447  
## ---
## 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):GroupAdults    8.837  8.992      8798  <2e-16 ***
## s(Time):GroupChildren  8.620  8.958      1630  <2e-16 ***
## s(Subject)            33.991 34.000 125335343  <2e-16 ***
## s(Subject,Time)       33.988 34.000 125318455  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.827   Deviance explained = 82.7%
## fREML = 2.0739e+05  Scale est. = 14.006    n = 75600

Plot of distribution of random intercepts:

# specify two plot panels:
par(mfrow=c(1,2), cex=1.1)
# random intercepts:
plot(m3, select=3)
# random slopes:
plot(m3, select=4)

Plot the summed effects for the adults without random effects (Left), and for two specific subjects (Right):

# specify two plot panels:
par(mfrow=c(1,2), cex=1.1)
# summed effects without random effects:
# This is actually: intercept + s(Time):GroupAdults
plot_smooth(m3, view="Time", cond=list(Group="Adults"), 
      rm.ranef=TRUE, main="intercept+s(Time):GroupAdults", rug=FALSE)

# Plot summed effects for Subject "a01":
plot_smooth(m3, view="Time", cond=list(Group="Adults", Subject="a01"), 
      main="... + s(Subject) + s(Subject, Time)", col='red',
      ylim=c(-15,20), rug=FALSE)
# Add Subject "a04":
plot_smooth(m3, view="Time", cond=list(Group="Adults", Subject="a18"), 
      add=TRUE, col='blue', xpd=TRUE)

Example of random smooths:

# first make sure Subject is a factor:
simdat$Subject <- as.factor(simdat$Subject)
# run model with random intercept:
m4 <- bam(Y ~ Group + s(Time, by=Group)
      + s(Time, Subject, bs="fs", m=1), data=simdat)
summary(m4)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ Group + s(Time, by = Group) + s(Time, Subject, bs = "fs", 
##     m = 1)
## 
## Parametric coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     4.9320     0.7286   6.770  1.3e-11 ***
## GroupChildren  -2.9868     1.0231  -2.919  0.00351 ** 
## ---
## 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):GroupAdults     7.710   7.837   69.87  <2e-16 ***
## s(Time):GroupChildren   6.804   7.003   15.83  <2e-16 ***
## s(Time,Subject)       309.911 322.000 1753.11  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.942   Deviance explained = 94.2%
## fREML = 1.6661e+05  Scale est. = 4.6865    n = 75600

Plot of distribution of random intercepts:

# random smooths:
plot(m4, select=3)
# add horizontal line:
abline(h=0)

Plot the summed effects for the adults without random effects (Left), and for two specific subjects (Right):

# specify two plot panels:
par(mfrow=c(1,2), cex=1.1)
# summed effects without random effects:
# This is actually: intercept + s(Time):GroupAdults
plot_smooth(m4, view="Time", cond=list(Group="Adults"), 
      rm.ranef=TRUE, main="intercept+s(Time):GroupAdults", rug=FALSE)

# Plot summed effects for Subject "a01":
plot_smooth(m4, view="Time", cond=list(Group="Adults", Subject="a01"), 
      main="... + s(Time, Subject)", col='red',
      ylim=c(-15,20), rug=FALSE)
# Add Subject "a04":
plot_smooth(m4, view="Time", cond=list(Group="Adults", Subject="a18"), 
      add=TRUE, col='blue', xpd=TRUE)

More about random effects

In time series data it is important to include random effects for Subjects and Items, because these cluster the observations. The model does not automatically know that a series of observations is produces by the same participant in a particular trial.

As a result, the estimates for the fixed effects may become less confident. However, the residuals of the model improve: the autocorrelation is reduced.

par(mfrow=c(1,3), cex=1.1)

acf_resid(m1, split_pred="Subject", main="ACF resid(m1)")
acf_resid(m3, split_pred="Subject", main="ACF resid(m3)")
acf_resid(m4, split_pred="Subject", main="ACF resid(m4)")

Setting up a GAMM model

Before the different components of a GAMM model were summarized. In this section we focus on finding the model that best fits the data.

We start with the full model. Note that we did not include interactions with the predictor Trial, because it is not the focus of our simulated experiment, and it takes much longer to run. Therefore, we only included Trial as random effect. We did not include Item effects, because these were not specified in this data set.

# convert to factors:
simdat$Subject <- as.factor(simdat$Subject)
simdat$Group <- as.factor(simdat$Group)

# This model takes a long time to run, therefore it is saved for later use:
if(F){
  m1 <- bam(Y ~ Group 
        + s(Time, by=Group) 
        + s(Condition, by=Group, k=5) 
        + ti(Time, Condition, by=Group)
        + s(Time, Subject, bs='fs', m=1)
        + s(Trial, Subject, bs='fs', m=1),
        data=simdat)
  save(m1, file="Models/m1.rda", compress='xz')
}
# load the model:
load('Models//m1.rda')
# summary:
summary(m1)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ Group + s(Time, by = Group) + s(Condition, by = Group, k = 5) + 
##     ti(Time, Condition, by = Group) + s(Time, Subject, bs = "fs", 
##     m = 1) + s(Trial, Subject, bs = "fs", m = 1)
## 
## Parametric coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    5.23738    0.04578  114.41   <2e-16 ***
## GroupChildren -3.16723    0.06474  -48.93   <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):GroupAdults                8.943   8.991 9602.908  <2e-16 ***
## s(Time):GroupChildren              8.860   8.976 3109.113  <2e-16 ***
## s(Condition):GroupAdults           3.908   3.936 1354.915  <2e-16 ***
## s(Condition):GroupChildren         3.855   3.898  649.766  <2e-16 ***
## ti(Time,Condition):GroupAdults    15.776  15.836 2232.226  <2e-16 ***
## ti(Time,Condition):GroupChildren  15.816  15.865 2135.490  <2e-16 ***
## s(Time,Subject)                  108.351 320.000    2.243  <2e-16 ***
## s(Trial,Subject)                 259.367 320.000   35.889  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =   0.95   Deviance explained =   95%
## fREML = 1.6089e+05  Scale est. = 4.0555    n = 75600

Then we check the distribution and autocorrelation structure of the residuals:

check_resid(m1, split_by=c("Subject", "Trial"))

The distribution is not completely normal, but rather approaching a t-distribution. The autocorrelation of the residuals is too high.

Including AR1 model

We try to account for the autocorrelation by including an AR1 model, specified by the parameter rho. When including rho in time series data, the model automatically treats the residuals as a single time series. However, this is clearly not the case:

Instead of treating the data as a single trial, the model should treat each individual trial as a unique time series, as indicated by the blue axis. To do that, one needs to specify in the data what is the first data point of each time series. In the plot the first data point of each time series is marked with a blue cross. We use the function start_event to determine the start of each time series.

simdat <- start_event(simdat, column="Time", event=c("Subject", "Trial"))
head(simdat)
##    Group      Time Trial Condition Subject         Y OFGroup IsChildren
## 1 Adults   0.00000   -10        -1     a01 0.7554469  Adults          0
## 2 Adults  20.20202   -10        -1     a01 2.7834759  Adults          0
## 3 Adults  40.40404   -10        -1     a01 1.9696963  Adults          0
## 4 Adults  60.60606   -10        -1     a01 0.6814298  Adults          0
## 5 Adults  80.80808   -10        -1     a01 1.6939195  Adults          0
## 6 Adults 101.01010   -10        -1     a01 2.3651969  Adults          0
##   start.event
## 1        TRUE
## 2       FALSE
## 3       FALSE
## 4       FALSE
## 5       FALSE
## 6       FALSE

The second number we need is the value of the autocorrelation parameter rho. We will take as starting value the ‘lag1’ value. However, it is better to use model comparisons to determine the correct value (not shown here).

(valRho <- acf(resid(m1), plot=FALSE)$acf[2])
## [1] 0.7375302
# This model takes a long time to run, therefore it is saved for later use:
if(F){
  m1.rho <- bam(Y ~ Group 
        + s(Time, by=Group) 
        + s(Condition, by=Group, k=5) 
        + ti(Time, Condition, by=Group)
        + s(Time, Subject, bs='fs', m=1)
        + s(Trial, Subject, bs='fs', m=1),
        data=simdat,
        AR.start=simdat$start.event, rho=valRho)
  save(m1.rho, file="Models/m1rho.rda", compress='xz')
}
# load the model:
load('Models//m1rho.rda')
# print summary:
summary(m1.rho)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ Group + s(Time, by = Group) + s(Condition, by = Group, k = 5) + 
##     ti(Time, Condition, by = Group) + s(Time, Subject, bs = "fs", 
##     m = 1) + s(Trial, Subject, bs = "fs", m = 1)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.06267    0.05621   36.70   <2e-16 ***
## GroupAdults  3.19516    0.07950   40.19   <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):GroupChildren              8.614   8.953 2621.146  < 2e-16 ***
## s(Time):GroupAdults                8.828   8.988 6354.578  < 2e-16 ***
## s(Condition):GroupChildren         3.761   3.877  419.494  < 2e-16 ***
## s(Condition):GroupAdults           3.853   3.928  894.686  < 2e-16 ***
## ti(Time,Condition):GroupChildren  15.702  15.919 1792.453  < 2e-16 ***
## ti(Time,Condition):GroupAdults    15.646  15.900 1470.099  < 2e-16 ***
## s(Time,Subject)                   32.460 320.000    0.195 1.04e-07 ***
## s(Trial,Subject)                 145.514 320.000    4.216  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.949   Deviance explained = 94.9%
## fREML = 1.2974e+05  Scale est. = 3.9107    n = 75600

Visualize the corrected residuals:

check_resid(m1.rho, split_by=c("Subject", "Trial"))

The red lines in the ACF plot are the corrected residuals after including an AR1 model. As is clear from this plot, the value of rho is probably too high, because the lag1 value is too low. To keep this tutorial short, we continue with this rho - although this is not an optimal value.

Testing for significance

The question is whether Children and Adults differ in their trend of Y over Time by Condition.

Three important methods to determine the significance of predictors:

  1. Model-comparison procedure (using function compareML).

  2. Test statistics of the smooth terms as presented in the summary.

  3. Visualization of the smooth terms (e.g., difference plots and difference surfaces).

Model-comparison procedure

# This model takes a long time to run, therefore it is saved for later use.
if(F){
  # Note: here the same rho is used, 
  # but it's better to set an optimal rho for this model.
  # Although here it probably does not differ so much, 
  # because we do not change the random effect structure.
  m2a.rho <- bam(Y ~ Group 
        + s(Time, by=Group) 
        + s(Condition, by=Group, k=5) 
        + ti(Time, Condition) # check: no difference between groups?
        + s(Time, Subject, bs='fs', m=1)
        + s(Trial, Subject, bs='fs', m=1),
        data=simdat,
        AR.start=simdat$start.event, rho=valRho)
  save(m2a.rho, file="Models/m2arho.rda", compress='xz')
}
load("Models/m2arho.rda")
# The interaction is different between groups:
compareML(m1.rho, m2a.rho)
## m1.rho: Y ~ Group + s(Time, by = Group) + s(Condition, by = Group, k = 5) + 
##      ti(Time, Condition, by = Group) + s(Time, Subject, bs = "fs", 
##      m = 1) + s(Trial, Subject, bs = "fs", m = 1)
## 
## m2a.rho: Y ~ Group + s(Time, by = Group) + s(Condition, by = Group, k = 5) + 
##      ti(Time, Condition) + s(Time, Subject, bs = "fs", m = 1) + 
##      s(Trial, Subject, bs = "fs", m = 1)
## 
## Chi-square test of fREML scores
## -----
##     Model    Score Edf   Chisq    Df  p.value Sig.
## 1 m2a.rho 129953.4  17                            
## 2  m1.rho 129744.3  20 209.132 3.000  < 2e-16  ***
## Warning in compareML(m1.rho, m2a.rho): AIC is not reliable, because an AR1
## model is included (rho1 = 0.737530, rho2 = 0.737530).
if(F){
  m2b.rho <- bam(Y ~ Group 
        + s(Time, by=Group) 
        + s(Condition, by=Group, k=5) 
        + s(Time, Subject, bs='fs', m=1)
        + s(Trial, Subject, bs='fs', m=1),
        data=simdat,
        AR.start=simdat$start.event, rho=valRho)
  save(m2b.rho, file="Models/m2brho.rda", compress='xz')  
}
load("Models/m2brho.rda")
# The interaction needs to be included:
compareML(m1.rho, m2b.rho)
## m1.rho: Y ~ Group + s(Time, by = Group) + s(Condition, by = Group, k = 5) + 
##      ti(Time, Condition, by = Group) + s(Time, Subject, bs = "fs", 
##      m = 1) + s(Trial, Subject, bs = "fs", m = 1)
## 
## m2b.rho: Y ~ Group + s(Time, by = Group) + s(Condition, by = Group, k = 5) + 
##      s(Time, Subject, bs = "fs", m = 1) + s(Trial, Subject, bs = "fs", 
##      m = 1)
## 
## Chi-square test of fREML scores
## -----
##     Model    Score Edf   Chisq    Df  p.value Sig.
## 1 m2b.rho 130343.5  14                            
## 2  m1.rho 129744.3  20 599.233 6.000  < 2e-16  ***
## Warning in compareML(m1.rho, m2b.rho): AIC is not reliable, because an AR1
## model is included (rho1 = 0.737530, rho2 = 0.737530).

Summary

With factors such as Group we cannot use the summary to test differences between groups. However, to investigate differences between groups, we could use difference smooths with ordered factors or binary predictors.

simdat$OFGroup <- as.factor(simdat$Group)
# change factor to ordered factor:
simdat$OFGroup <- as.ordered(simdat$OFGroup)
# change contrast to treatment coding (difference curves)
contrasts(simdat$OFGroup) <- 'contr.treatment'
# Inspect contrasts:
contrasts(simdat$OFGroup)
##          Children
## Adults          0
## Children        1
# This model takes a long time to run, therefore it is saved for later use:
if(F){
  # Note: here the same rho is used, 
  # but it's better to set an optimal rho for this model.
  # Although here it probably does not differ so much, 
  # because we do not change the random effect structure.
  m3.rho <- bam(Y ~ OFGroup 
        # reference curves:
        + s(Time) + s(Condition, k=5) 
        # difference curves:
        + s(Time, by=OFGroup) + s(Condition, by=OFGroup, k=5) 
        # reference surface:
        + ti(Time, Condition)
        # difference surface:
        + ti(Time, Condition, by=OFGroup)
        + s(Time, Subject, bs='fs', m=1)
        + s(Trial, Subject, bs='fs', m=1),
        data=simdat,
        AR.start=simdat$start.event, rho=valRho)
  save(m3.rho, file="Models/m3rho.rda", compress='xz')
}
# load the model:
load('Models//m3rho.rda')
# print summary:
summary(m3.rho)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ OFGroup + s(Time) + s(Condition, k = 5) + s(Time, by = OFGroup) + 
##     s(Condition, by = OFGroup, k = 5) + ti(Time, Condition) + 
##     ti(Time, Condition, by = OFGroup) + s(Time, Subject, bs = "fs", 
##     m = 1) + s(Trial, Subject, bs = "fs", m = 1)
## 
## Parametric coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    2.06316    0.05609   36.78   <2e-16 ***
## OFGroupAdults  3.19444    0.07933   40.27   <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)                            8.652   8.932 2657.910  < 2e-16 ***
## s(Condition)                       3.821   3.899  424.367  < 2e-16 ***
## s(Time):OFGroupAdults              8.116   8.762  850.445  < 2e-16 ***
## s(Condition):OFGroupAdults         2.952   3.228   55.770  < 2e-16 ***
## ti(Time,Condition)                15.750  15.915 1815.941  < 2e-16 ***
## ti(Time,Condition):OFGroupAdults  13.906  15.051  189.400  < 2e-16 ***
## s(Time,Subject)                   32.426 320.000    0.193 1.05e-07 ***
## s(Trial,Subject)                 145.714 320.000    4.217  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.949   Deviance explained = 94.9%
## fREML = 1.2973e+05  Scale est. = 3.9106    n = 75600

The summary now provides the significance of the difference terms: Children and adults show a significantly different trend over Time (F=850.445, edf=8.116, p<.001), a significantly different trend over Condition (F=55.770, edf=2.952, p<.001) and a significantly different interaction between Condition and Time (F=.193, edf=32.426, p<.001).

Plots

To find in which time window and on which condition children and adults differ, we can plot the model predictions using the functions fvisgam and plot_diff2.

# set two plot panels:
par(mfrow=c(1,2), cex=1.1)

# Plot difference surface:
plot_diff2(m1.rho, view=c("Time", "Condition"), 
    comp=list(Group=c("Adults", "Children")), 
    plotCI=TRUE, rm.ranef=TRUE,
    main="Adults-Children",
    zlim=c(-13,13))

# plot differences per condition:
for(i in -1:4){
  add=TRUE
  if(i==-1){
    add = FALSE
  }
  p <- plot_diff(m1.rho, view="Time", 
    comp=list(Group=c("Adults", "Children")),
    cond=list(Condition=i),
    ylim=c(-13,13),
    main="Adults-Children",
    col=rainbow(6)[i+2], add=add,
    rm.ranef=TRUE)
  text(2000, p$est[length(p$est)], pos=4, labels=i, col=rainbow(6)[i+2],
       font=2, xpd=TRUE)
}

The plots indicate that the children differ from adults in the conditions 2, 3, and 4 (longer stimulus onset asynchrony), but that the differences on the other conditions are not significant.

Visualizing the model predictions

There are many options for plotting the model’s predictions. Below are different examples.

Model terms / partial effects

mgcv’s plot.gam function allows for plotting the model terms. Additionally, pvisgam from itsadug provides more colorful plots, and also allows for plotting n-dimensional surfaces (e.g., ti(Time, Trial, Condition)).

par(mfrow=c(1,3), cex=1.1)

plot(m1.rho, select=2, scale=0, 
     shade=TRUE, rug=FALSE, bty='n')
# add horizontal line:
abline(h=0)
# add title:
title("s(Time):GroupAdults")

plot(m1.rho, select=6, 
    rug=FALSE)

pvisgam(m1.rho, view=c("Time", "Condition"), select=6,
  zlim=c(-20,10), main="pvisgam plot")

Summed effects with / without random effects

The function fvisgam plots the summed effects, without random effects (rm.ranef=TRUE) are plotted below. The random effects are small in comparison with the fixed effects, so the plots do not differ much with random effect included.

par(mfrow=c(1,2), cex=1.1)

fvisgam(m1.rho, view=c("Time", "Condition"),
        cond=list(Group="Adults"),
        main="Adults excl. random",
        rm.ranef=TRUE,
        zlim=c(-30,30))
fvisgam(m1.rho, view=c("Time", "Condition"),
        cond=list(Group="Children"),
        main="Children  excl. random",
        rm.ranef=TRUE,
        zlim=c(-30,30))

Random effects

The function plot.gam plots the random effects in a single plot. Sometimes it might be interesting to compare the different participants.

The function get_modelterm extracts the predicted values for a specific model term. These could be plotted with your preferred plot functions, for example using the package lattice.

# extract random effects for Trial:
pp <- get_random(m1.rho, cond=list(Time=1000))
# rename column:
names(pp)[2] <- 'fit'
require(lattice)
## Loading required package: lattice
lattice::xyplot(fit~Trial|Subject,
    data=pp, type="l",
    xlab="Time", ylab="Partial effect")

Or comparing children and adults:

# extract random effects for Trial:
pp <- get_random(m1.rho, cond=list(Trial=0))

emptyPlot(c(0,2000), range(pp$'fit.s(Time,Subject)'),
    h0=0, main="s(Time, Subject)")
for(i in sort(unique(pp$Subject))){
  tmp <- pp[pp$Subject==i,]
  lines(tmp$Time,tmp$'fit.s(Time,Subject)',
      col=ifelse(substr(i,1,1)=="a", alpha('blue'), alpha('red'))) 
}

# add legend:
legend(getFigCoords()[1], getFigCoords()[4],
    legend=c("Children", "Adults"),
    fill=c(alpha('blue'), alpha('red')),
    border=0, bty='n', xpd=TRUE)

Conclusion

On the basis of this (somewhat incomplete) analysis, we could argue that adults and children differ in their Y response over time on the larger values of Condition: children have a more ‘shallow’ or flat response than adults.