Packages & data

For creating animations you need to install the package animation from CRAN:

# when being asked for a mirrir to choose, 
# you could choose the first one:
install.packages("animation")

For this example you will need the newest version of the package itsadug, at least version 2.2 (on CRAN) or newer versions, see my website:

install.packages("itsadug", repos="http://cran.us.r-project.org")

Load the packages:

library(animation)
library(itsadug)
library(mgcv)
# set info messages to true:
infoMessages("on")

For the EEG example I used the data of one participant of the Kryuchkova et al. (2012) study. Note that this example is only for illustration purposes and does not present the best analysis. You can download a simplified version of the data here or run the following code:

download.file("http://jacolienvanrij.com/Tutorials/Files/EEG_example.rds", "EEG_example.rds", mode="wb")

Load the data:

dat <- readRDS("EEG_example.rds")

head(dat,3)
##   Electrode Item Trial      Time Animate      Amp   X  Y
## 1       AF3   w1   116 -192.1788     yes 9.312411 -16 31
## 2       AF3   w1   116 -184.3575     yes 6.989169 -16 31
## 3       AF3   w1   116 -176.5363     yes 7.051056 -16 31

Example: A GAMM model to investigate the pattern over different electrodes

Generalized Additive Mixed Modeling (GAMM; see Wood (2006) for more info) is a nonlinear regression technique that allows to investigate the changes in voltages measured by different electrodes over time.

We run a relatively simple nonlinear regression model that only includes Time and an interaction between X (approximate position of the electrode on the scalp, determining on which hemisphere the electrode is located) and Y (approximate position of the electrode on the scalp, determining how frontal / backward the electrode is located).

# may take some time to run:
m1<- bam(Amp ~ te(Time, X, Y, k=c(20,5,5), d=c(1,2)),
         data=dat, discrete=TRUE)
We can inspect the model by printing the summary and plotting the modelterms using the following R code. Note: In R markdown add ```{r, results='asis'} on top of the R block to show the output as HTML table. In addition, some HTML code was included in the markdown file to increase the padding of the table cells: <style>td { padding:2pt 10pt;}</style> 
gamtabs(m1, type='HTML')
A. parametric coefficients Estimate Std. Error t-value p-value
(Intercept) -0.2208 0.0305 -7.2395 < 0.0001
B. smooth terms edf Ref.df F-value p-value
te(Time,X,Y) 100.1025 113.7483 28.6043 < 0.0001

In the package mgcv that implements GAMMs interaction surfaces are by default plotted using a contour plot generated by the function plot. However, a three-way interaction is not possible (i.e., plot(m1, select=1) will generate an error). Instead, we could plot the partial effect of the interaction with pvisgam(), which uses colors to improve the readability (left plot below).

pvisgam(m1, view=c("X","Y"), zlim=c(-3,2))

Alternatively, we could use the function fvisgam(). In contrast with pvisgam, the function fvisgam visualizes the fitted effects. In other words, all modelterms (Intercept, s(Time) and s(X,Y)) are summed together and are plotted (right plot below). As the intercept is close to 0 here, the functions do not result in very different plots.

Note that in both plots the predictor Time is automatically set to the value 50, but we can change this value with the argument cond.

# Note the scale in z-scale:
fvisgam(m1, view=c("X","Y"), cond=list(Time=150), zlim=c(-3,2))

With the function plot_topo we could visualize the electrode positions on the scalp. This function requires the electrode positions as a list with ‘X’ and ‘Y’ coordinates, which we could retrieve from the data:

electrodes <- dat[,c("X","Y", "Electrode")]
electrodes <- as.list(electrodes[!duplicated(electrodes),])
str(electrodes)
## List of 3
##  $ X        : int [1:32] -16 16 -21 21 -10 10 -30 30 0 -18 ...
##  $ Y        : int [1:32] 31 31 0 0 -11 -11 -12 -12 0 22 ...
##  $ Electrode: chr [1:32] "AF3" "AF4" "C3" "C4" ...

The function topo_plot calls the functions pvisgam, fvisgam or plot_diff2, so their arguments need to be provided. Make sure to make the size of the plot regions squared.

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

plot_topo(m1, view=c("X","Y"), fun='pvisgam', zlim=c(-3,2),
          el.pos=electrodes, cond=list(Time=150), 
          main='pvisgam')
# automatically the function fvisgam is being selected:
plot_topo(m1, view=c("X","Y"), zlim=c(-3,2),
          el.pos=electrodes, cond=list(Time=150),
          main='fvisgam')

(In R Studio the plot window might give problems. If you do not see a plot appearing, type dev.off() and see if R STudio saved the plot as pdf in your working directory (“Rplots.pdf”). )

Instead of just selecting one value for Time, we could loop over different values for Time in afor-loop. The following code implements a for-loop to inspect the fitted (summed) effects of X and Y over time. A movie basically combines such a series of plots, and presents them consequently in quick pace.

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

for(i in seq(0,300, length=6)){
    # don't forget to fix the zlimits. Why is this important?
    plot_topo(m1, view=c("X","Y"), zlim=c(-6,4),
          el.pos=electrodes, cond=list(Time=i),
          main=paste("Time:",i))
}

Making a movie

The functions saveGIF and saveHTML from the package animation allow for storing the sequence of plots as a movie. Some of these functions require additional software. Please inspect the help pages of these functions, e.g. help(saveGIF).

library(animation)
saveHTML({
  for(i in seq(-50,300, length=50)){
      # don't forget to fix the zlimits. Why is this important?
      plot_topo(m1, view=c("X","Y"), zlim=c(-6,4),
            el.pos=electrodes, cond=list(Time=i),
            main=paste("Time:",round(i)))
     # add labels:
      text(electrodes[['X']], electrodes[['Y']],
        labels=electrodes[['Electrode']], cex=.75,
        pos=3, xpd=TRUE)     
  }
},ani.width = 300, ani.height=300)

Including the movie in a R Markdown file

If you are using R Markdown and would like to include a movie in your markdown file, the saveGIF function is easiest:

saveGIF({
  for(i in seq(-200,300, length=10)){
      # don't forget to fix the zlimits. Why is this important?
      plot_topo(m1, view=c("X","Y"), zlim=c(-6,4),
            el.pos=electrodes, cond=list(Time=i),
            main=paste("Time:",round(i)))
      
  }
})
  • And include the GIF in your markdown: <img src="animation.gif"/>

Alternative

See https://github.com/yihui/knitr/issues/883 for another solution:

  • add <link rel="stylesheet" href="http://vis.supstat.com/assets/themes/dinky/css/scianimator.css"> in your R markdown file (outside a R block).

  • start the following R block with ```{r,fig.width=4, fig.height=4, dpi=50, fig.show='animate', fig.path="img//"}

  • include the following code (so don’t use the saveGIF function):

for(i in seq(-50,300, length=50)){
    # don't forget to fix the zlimits. Why is this important?
    plot_topo(m1, view=c("X","Y"), zlim=c(-6,4),
          el.pos=electrodes, cond=list(Time=i),
          main=paste("Time:",round(i)))
}

References:

Kryuchkova, T., B. Tucker, Lee Wurm, and R.H. Baayen. 2012. “Danger and Usefulness in Auditory Lexical Processing: Evidence from Electroencephalography.” Brain and Language 122: 81–91.

Wood, Simon N. 2006. Generalized Additive Models: An Introduction with R. Chapman; Hall/CRC.