# “Unidentified” hierarchical model in brms/stan - where to go from here?

I am evaluating an intervention in which participants are grouped in teams and each participant fills in a survey before and after the intervention. As such, the data presents a classic multilevel problem in which observations are nested in participants and participants are nested in teams.

I want to test two kinds of hypotheses: (H1) that (a) participants differ in terms of the outcome variable when they come into the intervention and (b) these differences are predicted by certain person-level variables; (H2) that (a) participants differ in the extent to which the intervention works or not and (b) these differences depend on certain person-level variables (the same as in H1).

I tested a series of multilevel models (M1-M4), in each of which Time (0 = before, 1 = after) is the only L1 predictor; the intercept thus represents the outcome variable before the intervention while the coefficient of Time quantifies the change across the intervention. M1 tests H1a as a random-intercept model with Time as a fixed effect. M2 tests H1b by adding a L2 predictors as fixed effects. M3 tests H2a by adding a random slope for Time. M4 tests H2b by adding cross-level interaction terms between L2 predictors and Time.

I estimated all models in brms (Bürkner, 2016), an implementation of Bayesian generalized linear mixed models using Stan. Here's the code for models 4 and 5 where pos0 and neg0 are person-level variables and att is the outcome.

AB4 <- brm(att ~ 1 + Time + pos0 + neg0 + pos0:neg0 + (Time|Person) + (Time|Group), data = AB, iter = 4000, cluster = 4)
AB5 <- brm(att ~ 1 + Time + pos0 + neg0 + pos0:neg0 + Time:pos0 + Time:neg0 + (Time|Person) + (Time|Group), data = AB, iter = 4000, cluster = 4)


While the chains converged for $\sigma_{Group:Intercept}$ and $\sigma_{Group:Time}$, models 4 and 5 do not converge on a unique solution for $\sigma_{Person:Intercept}$, $\sigma_{Person:Time}$, and $\sigma_{Residual}$.

Family: gaussian (identity)
Formula: att ~ 1 + Time + pos0 + neg0 + pos0:neg0 + (Time | Person) + (Time | Group)
Data: AB (Number of observations: 4209)
Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
total post-warmup samples = 8000
WAIC: Not computed

Random Effects:
~Group (Number of levels: 1014)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)           0.25      0.08     0.07     0.38        424 1.01
sd(Time)                0.12      0.09     0.00     0.33        892 1.00
cor(Intercept,Time)    -0.05      0.56    -0.95     0.94       2523 1.00

~Person (Number of levels: 3091)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)           0.92      0.11     0.74     1.16         28 1.14
sd(Time)                0.76      0.26     0.38     1.29         29 1.13
cor(Intercept,Time)     0.32      0.37    -0.27     0.96         35 1.12

Fixed Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept     1.34      0.19     0.98     1.70       2216    1
Time          0.06      0.05    -0.04     0.16       8000    1
pos0          0.59      0.05     0.50     0.69       1698    1
neg0         -1.04      0.08    -1.20    -0.88       2798    1
pos0:neg0     0.15      0.02     0.10     0.19       2724    1

Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma(att)     1.19      0.09     0.96     1.31         26 1.15

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).


Here are the trace and density plots for the parameters in question.

I don't know where to go from here. In particular, I have the following questions:

1. In a maximum likelihood framework, this model would be unidentified as there are fewer bits of information than parameters to be estimated. As far as I understand, it is not that simple for MCMC models. Is there an equivalent rule or guideline for MCMC?
2. Where should I go from here? Is there anything I can do? For our research, it is very important to test both hypotheses. I could leave out the random slope and just estimate the cross-level interaction - but then I wouldn't know whether there was any (or how much) variance to be explained to begin with.

I would be very grateful for any advice and/or literature suggestions.

From your model summary I see that you have 4209 obervations in total and 3091 persons. That is, most persons only have 1 corresponding data point and therefore it will be difficult to estimate the random slope of time for them. This is likely what causes the convergence problems. I am not sure whether this random slope is really required or should rather be dropped from your model.

Anyway, you could try the following priors: Use half normal priors on the SD parameters (instead of half student-t priors) and put a bit more informative priors in the correlations (default is flat over [-1, 1]). Also (as Wayne writes) you could set proper priors on the "fixed effects". I would suggest that you try out the following:

prior = c(set_prior("normal(0, 5)", class = "b"),
set_prior("normal(0, 5)", class = "sd"),  # implicit half-normal prior
set_prior("lkj(2)", class = "cor"))

1. Your convergence problem can with near certainty be solved by simply running the MCMC longer - if the posterior is integrable and the MCMC is set up properly, it must eventually converge.

2. In a Bayesian analysis with proper (integrable) priors, there is no real limit to the number of parameters you can estimate - if you have too little data, you simply end up with something close to your priors.

3. I'm not even sure if it's true that your model is unidentifiable with MLE - do you posit that, or have you tried it out? In a multilevel model, parameters are not entirely free, so you cannot count a full degree of freedom for each data point.

Your number of effective samples is so low for Person -- two orders of magnitude -- that I'm not sure you can iterate your way out of this. Though you could try bumping the number of iterations up to, say, 12000 and see if you make progress.

I'm no expert, but I think the way forward will be to choose some priors instead of using the default uniform (improper) priors. I have no experience with assigning priors to sd, but my guess is that's where you might begin. See prior= in the brm help, and also set_prior (with class="sd").

Discussions like this at Andrew Gelman's blog, How to think about “identifiability” in Bayesian inference?, may be helpful.

I'm not sure, but the decov() prior in the rstanarm package might help. I think for AB4 you can do

library("rstanarm")
AB4 <- stan_glmer(att ~ 1 + Time + pos0 + neg0 + pos0:neg0 + (Time|Person) + (Time|Group),
prior_covariance = decov(),
data = AB, cores = 4)