I would appreciate any help to update my brmsfit object with a modified brms-generated stan model because I want to pass various columns of weights to the likelihood in a way that brms does not support yet probably.

My objective is to obtain the posterior distribution of effects after marginalizing over the distribution of weights as discussed here.

I need to do this in brms or stanarm rather than stan directly because I want to use functions of https://github.com/mjskay/tidybayes that are currently not supported by a stanfit object.

The idea is probably similar to the one discussed at Creating a brmsfit object with a modified brms-generated Stan model but I am not certain how to do that in practice, so any help would be much appreciated.

Additionally, I am not understanding how to implement in practice the second and third lines of this explanation:

**#the brms model to be modified and updated**

```
dt = read.table(header = TRUE, text = "
n r r/n group treat c2 c1 weights
62 3 0.048387097 1 0 0.1438 1.941115288 1.941115288
96 1 0.010416667 1 0 0.237 1.186583128 1.186583128
17 0 0 0 0 0.2774 1.159882668 3.159882668
41 2 0.048780488 1 0 0.2774 1.159882668 3.159882668
212 170 0.801886792 0 0 0.2093 1.133397521 1.133397521
143 21 0.146853147 1 1 0.1206 1.128993008 1.128993008
143 0 0 1 1 0.1707 1.128993008 2.128993008
143 33 0.230769231 0 1 0.0699 1.128993008 1.128993008
73 62 1.260273973 0 1 0.1351 1.121927228 1.121927228
73 17 0.232876712 0 1 0.1206 1.121927228 1.121927228")
m <-brm(r | trials(n) + weights(weights) ~ treat*c2+(1|group),
data=dt, family=binomial(link=logit))
```

**#the modified brms-generated stan model modified_brms.stan**

```
functions {
}
data {
int<lower=1> N; // total number of observations
int Y[N]; // response variable
int trials[N]; // number of trials
real<lower=0> weights[N, 10]; // model weights
int<lower=1> K; // number of population-level effects
matrix[N, K] X; // population-level design matrix
// data for group-level effects of ID 1
int<lower=1> J_1[N];
int<lower=1> N_1;
int<lower=1> M_1;
vector[N] Z_1_1;
int prior_only; // should the likelihood be ignored?
}
transformed data {
int Kc = K - 1;
matrix[N, K - 1] Xc; // centered version of X
vector[K - 1] means_X; // column means of X before centering
for (i in 2:K) {
means_X[i - 1] = mean(X[, i]);
Xc[, i - 1] = X[, i] - means_X[i - 1];
}
}
parameters {
vector[Kc] b; // population-level effects
real temp_Intercept; // temporary intercept
vector<lower=0>[M_1] sd_1; // group-level standard deviations
vector[N_1] z_1[M_1]; // unscaled group-level effects
}
transformed parameters {
// group-level effects
vector[N_1] r_1_1 = sd_1[1] * (z_1[1]);
}
model {
vector[N] mu = temp_Intercept + Xc * b;
for (n in 1:N) {
mu[n] += r_1_1[J_1[n]] * Z_1_1[n];
}
// priors including all constants
target += student_t_lpdf(temp_Intercept | 3, 0, 10);
target += student_t_lpdf(sd_1 | 3, 0, 10)
- 1 * student_t_lccdf(0 | 3, 0, 10);
target += normal_lpdf(z_1[1] | 0, 1);
// likelihood including all constants
if (!prior_only) {
for (n in 1:N)
for (w in 1:10) {
target += weights[n, w] * binomial_logit_lpmf(Y[n] | trials[n], mu[n]);
}
}
}
generated quantities {
// actual population-level intercept
real b_Intercept = temp_Intercept - dot_product(means_X, b);
}
```

**#I am able to do this successfully, but not certain how to proceed from here:**

mod_modified_brms_stan <- stan_model(‘modified_brms.stan’)

This question has also been posted here. Thanks in advance for any help.