I have a model estimating rates of change in amphetamine use during the first year of treatment among people receiving opioid agonist treatment (e.g. methadone buprenorphine). Been working on it a while now, following the workflow laid out by Aki Vehtari in this case study and learning an enormous amount as I (very slowly) go.

After much looing and k-folding I arrived at a zero-inflated beta-binomial regression model with a quadratic trend term. It is specified in brms like this

```
fit_zibetabinomial_quad4 <- brm(formula = bf(outcome | trials(set) ~ ats_baseline*yearsFromStart + ats_baseline*I(yearsFromStart^2) + (1 | encID)),
data = workingDF,
family = zero_inflated_beta_binomial(),
prior = c(prior(normal(0, 3),
class = Intercept),
prior(normal(0, 4),
class = b),
prior(cauchy(0, 3),
class = sd)),
save_pars = save_pars(all=TRUE),
seed = 1234,
refresh = 0,
cores = 4)
```

`outcome`

is number of days amphetamine use, `set`

is the number of possible days of amphetamine use in each observation period (this is always 28), `ats_baseline`

is the primary predictor, number of days of amphetamine used in the 28 days prior to starting treatment. `yearsFromStart`

is how many years from start of treatment each measurement was made. `I(yearsFromStart^2)`

is the quadratic term.

This is the output from the model

```
Multilevel Hyperparameters:
~encID (Number of levels: 695)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 1.77 0.11 1.55 1.99 1.00 1054 1608
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -4.84 0.17 -5.19 -4.50 1.00 1566 2189
ats_baseline 0.36 0.02 0.32 0.40 1.00 1264 1385
yearsFromStart 2.30 0.62 1.09 3.49 1.00 2462 2817
IyearsFromStartE2 -2.18 0.73 -3.61 -0.75 1.00 2576 2882
ats_baseline:yearsFromStart -0.81 0.09 -0.98 -0.64 1.00 1898 2440
ats_baseline:IyearsFromStartE2 0.66 0.10 0.47 0.86 1.00 2349 2951
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
phi 3.92 0.42 3.16 4.79 1.00 2107 2725
zi 0.10 0.03 0.06 0.16 1.00 2753 2899
```

All ESS and Rhat good. Good diagnostic plots etc.

I wondered if a random-slopes version of the model would work, with the only change being that `(1 | encID)`

in the formula becomes `(yearsFromStart| encID)`

. I ran that model and got the following output

```
~encID (Number of levels: 695)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.69 0.09 0.51 0.88 1.00 2319 2949
sd(yearsFromStart) 6.15 0.73 4.85 7.63 1.00 477 1260
cor(Intercept,yearsFromStart) 0.92 0.07 0.76 1.00 1.06 61 349
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -4.17 0.12 -4.41 -3.94 1.00 2566 3181
ats_baseline 0.31 0.01 0.28 0.34 1.00 2887 3136
yearsFromStart -1.80 0.90 -3.62 -0.02 1.00 1179 2205
IyearsFromStartE2 -4.51 0.85 -6.16 -2.93 1.00 2685 3065
ats_baseline:yearsFromStart -0.71 0.10 -0.91 -0.51 1.00 1608 2293
ats_baseline:IyearsFromStartE2 0.86 0.11 0.64 1.08 1.00 2820 2696
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
phi 5.91 0.68 4.70 7.30 1.00 1312 2542
zi 0.06 0.02 0.02 0.10 1.00 2523 2117
```

Along with the warning `Warning message: Parts of the model have not converged (some Rhats are > 1.05). Be careful when analysing the results! We recommend running more iterations and/or setting stronger priors.`

with the obvious culprit being the Multilevel hyperparameters. I am wondering if anyone has advice on how to specify these hyperparameters more intelligently than I have.