@paul.buerkner - here is an update, and unfortunately I am a bit stuck.

I have successfully implemented the random walk prior in Stan (not `brms`

). To do this, I have defined this `rw1()`

function, curtesy of @Bob_Carpenter -

```
real rw1_lpdf(vector alpha, real sigma) {
int N = cols(alpha);
return normal_lpdf(alpha[2:N] - alpha[1:N - 1] | 0, sigma);
}
```

When I set the distribution of the random effects, instead of setting them to be the usual `normal`

, I set them to be `rw1`

. I also add a sum-to-zero constraint.

The problem in trying to implement this in `brms`

is that for the formula `Y ~ (1|agecat)`

, `brms`

automatically sets the random effects as normally distributed, and I can’t figure out how to remove that code. I tried doing what you suggested and setting the standard deviation of that normal distribution to be `constant(10000)`

, and then added the `rw1`

distribution on top of it, but it caused all sorts of sampling problems.

Do you have any thoughts on if there is a way to “remove” the existing normal distribution? Do you know why just setting that normal distribution to be super-flat (large sd) and then adding the second `rw1`

distribution on top didn’t work? I saw that you can use `gr()`

to replace the normal with a t-distribution for the random effects, but it doesn’t allow any other distribution such as `rw1`

.

Here’s some simulated data and my best attempt at the `brms`

code, if you want to give it a shot:

```
# Simulate Data
N <- 100
d <- tibble(agecat = sample.int(8, N, replace = TRUE)) %>%
mutate(mu = sin(pi*agecat/7),
Y = mu + rnorm(N, sd = 1))
# Set up stanvars and priors
sv <- stanvar(scode = 'real rw1_lpdf(vector alpha, real sigma) {
int N = num_elements(alpha);
return normal_lpdf(alpha[2:N] - alpha[1:N - 1] | 0, sigma);
}',
block = "functions") +
stanvar(scode = 'r_1_1 ~ rw1(sd_agecat_rw);', block = "model", position = "end") +
stanvar(scode = "real<lower=0> sd_agecat_rw;", block = "parameters", position = "end") +
stanvar(scode = "sd_agecat_rw ~ normal(0,1);", block = "model", position = "end") +
stanvar(scode = "sum(r_1_1) ~ normal(0, N_1*.001);", block = "model", position = "end")
pr <- prior(constant(10000), class = "sd", group = "agecat")
# Fit model
fit <- brm(Y ~ (1|agecat), data = d, prior = pr, stanvars = sv)
```