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)