Hi! I’ve spent over a year on this problem and need help!

I am trying to model RMSE values produced by external algorithms.

I have 4 factors (for simplicity, being modeled as unordered):

size (3 levels)

noise (3 levels)

shape (2 levels)

interps (9 levels; the repeated measure)

Also, I have:

g_subj (unique id)

g_rep (20 levels; unique id; subjects are nested within replicates; each rep is a replication of the entire experiment)

see Data:

_RMSE_data_unord.csv (379.2 KB)

The RMSE values are y and I am using the custom gamma2() distribution from @paul.buerkner notebook to model the heteroscedasticity.

- in fitting the variance for the Gamma, are we estimating different shape parameters (essentially, different gamma distributions)? Does this completely eliminate the variance-mean relation?

The model I would like is:

```
y ~ g_size * g_noise * g_shape * g_interps + (g_size * g_noise * g_shape * g_interps | g_rep) + (g_interps | g_subj)
v ~ g_size * g_noise * g_interps
```

I am starting simpler with:

```
y ~ g_size * g_noise * g_shape * g_interps + (1| g_rep) + (1| g_subj)
v ~ g_size * g_noise * g_interps
```

but keep getting different combinations of max tree depth and divergents for every iter after warmup.

I was able to successfully implement:

```
y ~ g_size * g_noise * g_shape * g_interps
v ~ g_size * g_noise * g_interps
```

I originally started with nlme but it had trouble modeling 0 variances and 0 random effects. I then moved to lme4 but it could not model the heteroscedasticity, so now I’m using brms.

Run time for the models range from 4 hours to 3 days depending on the complexity.

This is the current model that works:

```
fit0g16_4a_wPriors20sub3_95a <- brm(bf(y ~ g_size * g_noise * g_shape * g_interps,
v ~ g_size * g_noise * g_interps),
prior =c(
# y ~ Intercept + b
prior(normal(4, 3),class="Intercept"),
prior(normal(2, 3),class='b'),
prior(normal(0, 3),class="Intercept", dpar = 'v'),
prior(normal(0, 5),class='b', dpar = 'v')),
family = gamma2, stanvars = stanvars_gamma2,
sample_prior = "yes",
seed = 2011631714,
data = t_long_CNSubset20red1_sub3_unord, warmup = 1000, iter = 2000, chains = 4, core = 4,
control = list(adapt_delta = 0.90, max_treedepth = 10))
```

I’ve tried different priors for SD but get different combinations of max tree depth and divergents after warmup for all iters.

This model does not converge and has 4000 max tree depths:

```
fit0g16_4a_wPriors20sub3_27x10b <- brm(bf(y ~ g_size*g_noise*g_shape*g_interps + (1 | g_rep) + (1 | g_subj),
v ~ g_size*g_noise*g_interps),
prior =c(
# y ~ Intercept + b
prior(normal(4, 3),class="Intercept"),
prior(normal(2, 3),class="b"),
# v ~ Intercept + b
prior(normal(0, 4),class="Intercept",dpar="v"),
prior(normal(0, 5),class="b",dpar="v"),
# y ~ sd
prior(normal(0, 0.2),class='sd')),
family = gamma2, stanvars = stanvars_gamma2,
sample_prior = "yes",
seed = 2011631714,
data = t_long_CNSubset20red1_sub3_unord, warmup = 1000, iter = 2000, chains = 4, core = 4,
control = list(adapt_delta = 0.90, max_treedepth = 10))
```

- Would modeling correlations improve convergence?
- Help!