Introduction

This script will simulate a very simple experiment with two conditions, show how to model the data using brms, and plotting the output in some useful ways that are pretty much publication-ready. First, let’s load all the packages we need:

library("brms")
library("tidyverse")

First we need to simulate the experiment data. The simulation below allows for individual variability in the control condition (varying intercept) and also individual variability in sensitivity to my condition effect (varying slope).It then does a quick and dirty plot of the distribution of two conditions.

The simulation code is based on chapter 13 of Solomon Kurtz / McElreath: https://bookdown.org/connect/#/apps/1850/access

pp_n <- 20  #number of participants
trials_n <- 10  #number of trials
sigma_y <- 1 #within-participant variability / measurement error

a       <-  0  # average control condition mean
b       <- .5    # average effect size
sigma_a <- .5    # std dev in intercepts (control means, between participant variability)
sigma_b <-  1  # std dev in slopes (effect size)
rho     <- -.5   # correlation between intercepts and slopes

#combine means and standard deviations for multivariate gaussian sampling
mu <- c(a, b)
sigmas <- c(sigma_a, sigma_b)          # standard deviations
rho    <- matrix(c(1, rho,             # correlation matrix
                   rho, 1), nrow = 2)

# now matrix multiply to get covariance matrix
sigma <- diag(sigmas) %*% rho %*% diag(sigmas)

#let's create our participant level parameters
set.seed(1)  # used to replicate example
vary_effects <- MASS::mvrnorm(pp_n, mu, sigma) #use multivariate gaussian

#rename columns for varying intercepts and varying slopes
vary_effects <- vary_effects %>% as_tibble() %>% rename(pp_a = V1,
                                                        pp_b = V2)

  
#now simulate observations.
#this piping operation is a little involved, but essential it makes a dataframe of with trial_n entries of all  combinations of participant intercepts and slopes. Then adds a binary column to signify whether the condition is the experimental (1) or control (0). It then uses this column to switch 'on' or 'off' the experimental effect to generate a column of participant condition means. Then it uses the specified (and  constant) within-individual variability (or measurement error) to generate observations
set.seed(1)  # used to replicate example
data <- vary_effects %>% 
  mutate(pp = 1:pp_n) %>% 
  expand(nesting(pp, pp_a, pp_b), trial = 1:trials_n) %>% 
  mutate(condition = rep(0:1, times = n() / 2)) %>% 
  mutate(mu = pp_a + pp_b * condition) %>% 
  mutate(measurement = rnorm(n = n(), mean = mu, sd = sigma_y))

head(data)
#quick plot.
ggplot(data, aes(x=measurement, group=condition)) + geom_density()

Fit a brms model

Now we have our data, we can fit a multi-level model. The model below predicts measurement, and allows for varying intercepts (the 1 | pp bit) and varying slopes (the condition | pp). For our simple model the default priors for brms will be sufficient. We then see the summary output of our model. We can have a quick look here to see how similar the fitted values are to our created values in the code above. Compare Intercept with a, condition with b, sd(Intercept) with sigma_a, sd(condition) with sigma_b, and cor(Intercept,condition) with rho.

if (!file.exists("m1.rda")) {
  m1 = brm(measurement ~ 1+ condition + (1 + condition | pp), data = data,
            family = gaussian(),
            iter = 3000, cores = 2, refresh = 500, chains = 1, warmup = 1000)
  save(m1, file ="m1.rda")
} else {load("m1.rda")}

summary(m1)            
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: measurement ~ 1 + condition + (1 + condition | pp) 
##    Data: data (Number of observations: 200) 
## Samples: 1 chains, each with iter = 3000; warmup = 1000; thin = 1;
##          total post-warmup samples = 2000
## 
## Group-Level Effects: 
## ~pp (Number of levels: 20) 
##                          Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## sd(Intercept)                0.25      0.16     0.01     0.60        459
## sd(condition)                1.08      0.24     0.68     1.62        530
## cor(Intercept,condition)    -0.29      0.45    -0.94     0.80        135
##                          Rhat
## sd(Intercept)            1.00
## sd(condition)            1.00
## cor(Intercept,condition) 1.02
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept     0.02      0.12    -0.20     0.25       2245 1.00
## condition     0.62      0.28     0.08     1.17        913 1.00
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma     0.95      0.05     0.85     1.06       1820 1.00
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

Model Diagnostics

You could get lost in model diagnostics, but it is important to check that your model fitted well. Some key things to look for are whether the posterior distribution matches the raw data, whether the Rhat values in the summary are close to 1 (if they are not, the model probably didn’t converge well), and whether the trace plots look like they are exploring the parameters space well (i.e. they look like a noisy signal and do not get ‘stuck’ anywhere). Brms has some really useful functions of checking your model: plot and pp_check. You can also launch shinystan, which is a fantastic interactive tool for investigating your model.

pp_check(m1) #quick posterior predictive check

#plot(m1)
#launch_shinystan(m1)

Dealing with posterior samples

Now we have our data we can plot the differences. Dealing with the posterior samples to get them into the desired format ready for plotting is not straight forward (at least, it was a learning curve for me). It is initially confusing as to whether I wanted to include participant-level parameters in my graphing, or include the measurement error (which would be predicting new unobserved participants). It is worth keeping in mind that you are trying to estimate the population-level parameters (i.e. the difference between population means). The model used in this example is extremely simple so we could calculate the posterior estimate of condition means directly from the model. The intercept is the control condition and the experimental condition is intercept + condition. However, calculating the means by selecting the combination of coefficients you are interested in quickly gets complicated when you have multiple factors with multiple levels, as will be the case in most factorial experiments. My general approach below gets brms to do the hard work by feeding the fitted model an empty data frame of condition names. There is then some slightly convoluted data wrangling to get the samples into a familiar long-format data frame that we can play around with. There is probably a neater way of doing this…

condition_list <- data.frame(condition=c(0,1))  #empty dataframe with a codnition list of all factors and levels. One entry per combination

#put allow_new_levels = TRUE because you have not specified a participant number. So the model calculates the group level parameters.
df_fit <- fitted(m1, newdata=condition_list, re_formula = NA, allow_new_levels = TRUE, summary=FALSE) #estimate populations mean
tdf_fit <- t(df_fit) #transpose
selection_cbind <- cbind(condition_list,tdf_fit) #add the condition values
estimates <- selection_cbind %>% gather(key = "Rep", value = "measurement", -condition) #you now have the posterior estimates in a familiar long format for plotting with ggplot, which two condition estimates per 'replication'

#Resulting dataframe has nrows = Npp * Ncndt * Nsamples.

head(estimates)

Plotting

For inference we are interested in the contrast between the control and the experimental condition. I’ve created two plots that I quite like to display the contrast (I’m happy to hear more!). The first plot graphs the distribution and adds a bar for the 95% hdi interval. These can look great, but sometimes you risk overplotting (too much unnecessary information can make a plot hard to parse) if you have many conditions. The second plot uses to bars and an annotation to convey useful information. Although we only have one condition in this example I’ve graphed the contrast as one level of a factor so that the code will apply to factorial designs. We will use the packages ggridges and HDInterval. The hdi function will allow us to calculate useful metrics for inference, such as the proportion of mass above zero.

#calculate contrast for each replication..
contrast_df <- estimates %>% 
    group_by(Rep) %>% 
  summarise(contrast = measurement[2] - measurement[1])

head(contrast_df)
library(HDInterval)

#create hdis, use the HDIinterval hdi function to calculate values for point ranges for drawing on segments.
plot_hdis <- contrast_df %>% 
  summarise(HDI2.5 = hdi(contrast, credMass = .95)[1],
            HDI97.5 = hdi(contrast, credMass = .95)[2],
            HDI.5 = hdi(contrast, credMass = .99)[1],
            HDI99.5 = hdi(contrast, credMass =  .99)[2],
            AboveZero = mean(contrast > 0),
            BelowZero = mean(contrast < 0),
            Mn = mean(contrast),
            Med = median(contrast))

print(plot_hdis)
## # A tibble: 1 x 8
##   HDI2.5 HDI97.5  HDI.5 HDI99.5 AboveZero BelowZero    Mn   Med
##    <dbl>   <dbl>  <dbl>   <dbl>     <dbl>     <dbl> <dbl> <dbl>
## 1 0.0570    1.14 -0.104    1.30     0.988     0.012 0.623 0.622
#now you have your hdis. We can plot.
#first, let's add some theme adjustments that make the plot look nice.
theme_mole <- theme_classic() +
  theme(strip.background = element_rect(fill=NA,color=NA), 
        strip.text = element_text(face="bold",colour="black",size="12"), 
        axis.title = element_text(face="bold",colour="black",size="12"),
        axis.text.x = element_text(vjust=-.5),
        axis.text.y = element_text(vjust=.5),
        axis.text = element_text(face="plain",colour="black",size="10"),
        legend.text = element_text(face="plain",colour="black",size="10"),
        legend.title = element_text(face="bold",colour="black",size="12"),
        legend.key = element_blank(),
        panel.grid.major.y = element_line(color="grey85",size=.2, linetype = 2))

library(ggridges)

#this is an awkward fudge so that you could use the code below with multiple levels of factors, and you can see how the 'as.numeric
contrast_df$myfactor <- factor('level1') #add a fake factor for plotting
plot_hdis$myfactor <- factor('level1') #add a fake factor for plotting

#These are so much better if they are factors.
ggplot(contrast_df, aes(x=contrast, y=myfactor)) + geom_density_ridges(fill="gainsboro", scale=.96, rel_min_height=.005)  + 
  #theme_transition +
  geom_vline(xintercept=0, linetype="dashed") + theme(panel.grid.major.x = element_line(color="grey40",size=.2, linetype = 2)) +
  geom_segment(data=plot_hdis, aes(x=HDI2.5, y=as.numeric(myfactor), xend=HDI97.5, yend=as.numeric(myfactor)), size=1.5) + #if you have multiple levels, put the y and yend as as.numeric(myfactor)
  geom_point(data=plot_hdis, aes(x=Mn, y=as.numeric(myfactor)), size=3) + 
  #xlim(c(-.4, .6)) +
  ylab("Factor levels") +
  xlab("Contrast") +
  theme_mole

#If you have many levels the can look crowded using ggridges. Simply graphing the hdi segments is a clearer option.

# this plotting code probably looks confusing and backward. Thgis is because it uses coord_flip since I like displaying the distributions horizontally.

ggplot(data=plot_hdis, aes(x=myfactor, y=Mn, ymin=HDI2.5, ymax=HDI97.5)) +
  geom_hline(yintercept=0, lty=2, col="Grey") +  # add a dotted line at x=0 after flip
  geom_pointrange(aes(ymin=HDI.5, ymax=HDI99.5), fatten=.5, col="Grey", size=1.5) + #thinner grey 99% HDIs
  geom_pointrange(fatten = 2, size=1.5) + #black 95% HDI
  xlab("Factor Levels") + ylab("Contrast (95% HDIs)") +
  theme_mole + geom_text(aes(label = ifelse(AboveZero>BelowZero, paste(AboveZero*100, "%>0", sep=""),paste(BelowZero*100, "%<0", sep="") )), vjust=-2, hjust=.4, size=4) +  
  coord_flip() #+# flip coordinates (puts labels on y axis)