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)
```