Extracting posterior samples of smooth effects in brms

Interfaces brms
1

I have a simulated data with different smooth group effects and I want to examine the posterior distribution of the difference between the smooth effects of group 1 and group 2. Is there a simple way to do this in brms? If not, is there a clean way of manually computing the smooth effects based on the posterior_samples() output?

EDIT: Adding to the second question, is there something like mgcv’s predict(…, type=“lpmatrix”) to obtain a design matrix for new data?

Here’s my code:

## Simulated data
set.seed(11)

N <- 20    # Number of observed sequences
B <- 20    # Length of each sequence

id <- rep(1:N, each=B)           # ID of each sequence
group <- rep(c(1,2),each=N*B/2)  # group membership 

x <- rep(seq(1,10,length=B),N)   # covariate domain

# Create dataframe for brms input
df <-data.frame(x)
df$id <- as.factor(id)
df$group <- as.factor(group)
df$t <- rep(1:B,N)

# Define AR1 correlation matrix Sigma
Sigma <- corAR1(value=0.8,form = ~ 1|id)
Sigma <- Initialize(Sigma, data=df)
Sigma <- corMatrix(Sigma)
sd <- 2    # Noise sd
tau <- 3   # Random intercept sd
theta <- rnorm(N, mean=0, sd=tau)   # Random intercept

# Simulate response (group 2 is negative of group 1 response)
y <- 10*sin(x)/x + sd*c(t(mvrnorm(N,rep(0,B),Sigma$`1`))) + theta[id]
y[group==2] <- -y[group==2]
df$y <- y

plot(x,y, main="Non-linear response with AR1 residuals", col=group)
lines(x[1:B], y[1:B], type='o', col=1)

Then I fit the following model with brms:

m0 <- brm(y ~ s(x, by=group, bs='cr', k=10) + (1|id),
    autocor = cor_ar(formula = ~t|id),
    data= df)
2

I would recommend using the tidybayes package in combination with brms to obtain this most easily.

Alternatively, more manually, specify a new data frame of x and group values yourself, for which you want to make predictions. Pass this to fitted(m0, newdata = <new data frame, summary = FALSE).
Then, on the level of the posterior samples, compute the differences of the two groups for each x. Afterwards, you can summarise, plot etc. this posterior distribution of the differences.

Response to your EDIT: It is not possible to directly extract the design matrix of new data, but you can obtain it indirectly by calling extract_draws with the approporiate arguments and then looking through its structure for elements called sm which contain all spline related stuff.

2 Likes