Multilevel Model Example

simulated rat data

Take a look at simulated data

rbind(head(d), tail(d))

    seizure neuron region rat      rate
1    before      1      1   1 12.236493
2    during      1      1   1 15.072837
3    before      2      1   1 12.092462
4    during      2      1   1 16.449610
5    before      3      1   1  8.312029
6    during      3      1   1 14.339275
595  before      8      3  10 12.346769
596  during      8      3  10 18.994619
597  before      9      3  10 11.043983
598  during      9      3  10 14.348025
599  before     10      3  10  9.935970
600  during     10      3  10 14.255610

visualize data

library(ggplot2)
ggplot(d) +
  aes(x = seizure, y = rate) +
  geom_jitter(width = 0.1, height = 0) +
  facet_wrap(~ rat, labeller = "label_both")

Model data

Model rate as a function of seizure level (before/during) controlling for rat, regions within rat, and neurons within regions within rats.

library(lme4)
m <- lmer(rate ~ seizure + (1 | rat/region/neuron), data = d)
summary(m, corr = FALSE)

Linear mixed model fit by REML ['lmerMod']
Formula: rate ~ seizure + (1 | rat/region/neuron)
   Data: d

REML criterion at convergence: 2226

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.3149 -0.6325 -0.0110  0.6253  3.3229 

Random effects:
 Groups              Name        Variance Std.Dev.
 neuron:(region:rat) (Intercept) 0.4871   0.6979  
 region:rat          (Intercept) 1.1862   1.0891  
 rat                 (Intercept) 0.6731   0.8204  
 Residual                        1.6643   1.2901  
Number of obs: 600, groups:  neuron:(region:rat), 300; region:rat, 30; rat, 10

Fixed effects:
              Estimate Std. Error t value
(Intercept)    10.1589     0.3377   30.09
seizureduring   5.0490     0.1053   47.93

Expected rate goes up by about 5 during seizure. Notice there are four sources of variation:

1. between neurons, which are within regions, which are within rats
2. between regions, which are within rats
3. betweem rats
4. within rats (Residual)

Calculate 95% confidence interval on the seizure coefficient.

confint(m, parm = "seizureduring")

Computing profile confidence intervals ...

                 2.5 %   97.5 %
seizureduring 4.842273 5.255804

Rate increases about 4.8 to 5.3 during seizure.

visualize model

library(ggeffects)
plot(ggpredict(m, terms = "seizure"), connect.lines = TRUE)

basic model diagnostics

Check for constant variance of residuals. Looks good

plot(m)

Check for constant variance of residuals within rat. Looks good.

plot(m, rat ~ resid(., scaled=TRUE))

Check normality of residuals. Looks good.

lattice::qqmath(m)

how data was simulated

This may be of interest. Notice how closely the model coefficients (under fixed effects) and standard deviations (under random effects) come close to the “true” values used to simulate the data below. For example, below the “true” between rat standard deviation is set to 0.8. Above the model estimated that standard deviation between rats to be 0.8204. Likewise, the “true” increase in rate during seizure is set to 5 below. The model output above estimates it to be 5.0490. Of course the model performs well because I fit the “correct” model (since I knew how I generated the data).

# generate levels that define unique observations
rat <- 1:10
region <- 1:3
neuron <- 1:10
seizure <- c("before","during")

# combine levels into a data frame
d <- expand.grid(seizure = seizure, neuron = neuron, 
                 region = region, rat = rat)

# generate random effects;
set.seed(1)
e_rat <- rnorm(n = 10, mean = 0, sd = 0.8) # rat
e_region <- rnorm(n = 10*3, mean = 0, sd = 1.3) # region:rat
e_neuron <- rnorm(n = 10*3*10, mean = 0, sd = 0.7) # neuron:(region:rat) 
e <- rnorm(nrow(d), mean = 0, sd = 1.2) # Residual

# generate nesting identifiers
d$region_rat <- as.numeric(interaction(d$region, d$rat))
d$region_rat_neuron <- as.numeric(interaction(d$region, d$rat, d$neuron))

# generate rate:
# 10 is true "before" seizure rate (ie, intercept)
# 5 is true change in rate "during" seizure
d$rate <- (10 + e_rat[d$rat] + 
             e_region[d$region_rat] + 
             e_neuron[d$region_rat_neuron]) + 
  5*(d$seizure == "during") + e

# tidy up and declare rat, region, and neuron as factors
d$region_rat <- NULL
d$region_rat_neuron <- NULL
d$rat <- factor(d$rat)
d$region <- factor(d$region)
d$neuron <- factor(d$neuron)