Fit multilevel model to evoked responses using lme4

Denis A. Engemann

INRIA Saclay, Parietal
Source: vignettes/plot_evoked_multilevel_model.Rmd

Here we will use a multilevel model to assess the difference between two conditions over time. For now we will stick with the lme4 optimization-based implementation and use the parametric bootstrap for uncertainty estimation. Bayesian inference will be covered in a future example.

library(tidyverse)
library(lme4)
library(merTools)
library(mne)
## Importing MNE version=0.18.dev0, path='/Users/dengeman/github/mne-python/mne'

Let’s read in the raw data.

We can now go ahead and compute evokeds.

Let’s get the epochs into a dataframe.

Now we can set up a simple varying effects model:

\[ y_i \sim N(\alpha + \alpha[j] + T\beta[j]_i, \sigma^2) \] Where \(y_i\) is the MEG observation at a given epoch and time point, \(\alpha\) is the global intercept, \(\alpha[j]\) the intercept of each tome point, \(T\) is our treatment information, and \(\beta[j]\) is the estimated treatment effect over all \(j\) time points of the epoch. The multilevel model will learn a population-shrinkage that depewnds on the variance over the time.

One of the assumptions of this model is that all \(J\) coefficients come from a common distribution:

\[ \beta_{j} \sim N(\mu_{\beta[j]}, \sigma_{j}^2) \] This is what makes this model one model instead of \(j\) separate models. In this model, all \(j\) coefficients are shrunk to \(\mu_j\) to an extent that is proportional to the variance of the group effects.

mod1 <- lmer(observation ~ 1 + condition + (1 + condition | time),
             data = epochs_df)
mod1 %>% summary() %>% print()
## Linear mixed model fit by REML ['lmerMod']
## Formula: observation ~ 1 + condition + (1 + condition | time)
##    Data: epochs_df
## 
## REML criterion at convergence: 177903.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.9168 -0.6324 -0.0012  0.6395  4.2404 
## 
## Random effects:
##  Groups   Name           Variance Std.Dev. Corr 
##  time     (Intercept)    2003     44.76         
##           conditionvis/l 2506     50.06    -0.94
##  Residual                6032     77.67         
## Number of obs: 15370, groups:  time, 106
## 
## Fixed effects:
##                Estimate Std. Error t value
## (Intercept)      10.050      4.437   2.265
## conditionvis/l  -21.756      5.021  -4.333
## 
## Correlation of Fixed Effects:
##             (Intr)
## conditnvs/l -0.927

As often, looking at the coefficients and text summaries of such models is not tremendously helpful. We will instead generate predictions.

We can see that the model predictions are recapitulating the temporal structure of the data. But how much uncertainty do we have about those predicrions?

In the absence of a full Bayesian analysis we can do an informed parametric bootstrap that takes into account all sources of uncertainty of our multilevel model, which yields more conservative estimates of compatibility intervals. We will use the merTools package for this purpose.

We’ll put a few custom options which should speak for themselves.

Let’s first look at the compatibility intervals.

Now let’s plot the single simulations on top of it. We first have to extract the simulations from the merTools output.

Now we can plot the prediction from each bootstrap replica. In order to not clutter things we will make generous use of transparency.

While we see a clear condition-relative modulation of the signal, at the same time we have to appreciate considerable uncertainty.