Factor Coefs for Varying Group and Population Errors

I’m trying to fit a model with the group and population errors varying as a function of two factors. Something of the form:

y = a_j + \varepsilon
\varepsilon \sim N(0, \sigma = Xb)
a_j \sim N(0, \sigma_a = Gg)

My problem is that I’m not able to get the factor coefficients for predicting the errors (my bs and gs above).

I’ve simulated an example. An RMarkdown file is given below showing the full code, and two Word Docs showing the output for one example that works and one that does not.

The example works or does not depending on the form of the model matrix I pass. With variables Var1 and Var2 both as two-level factors (levels a & b) the example works if my matrix is of the form (but I don’t get the factor coefficients bs and gs that I want, just the group errors):

# created using interaction(Var1, Var2)
# rows omitted for brevity
aa  ba  ab  bb
 1   0   0   0
 0   1   0   0
 0   0   1   0
 0   0   0   1

But does not work if the matrix is of the form:

# created using model.matrix(lm(Y ~ Var1 + Var2, data = df))
# rows omitted for brevity...
Var1a  Var1b  Var2b
    1      0      0
    0      1      0
    0      0      1

My stan code is:

  int<lower=0> N;                // total number of obs across entire data set
  int<lower=0> J;                // number of groups
  int<lower=0> K;                // number of fixed effects
  vector[N] y;                   // response vector
  matrix[N,J] Z;                 // indicator matrix to identify group
  matrix[N,K] X;                 // indicator matrix to identify Var1 * Var2
  matrix[J,K] G;                 // indicator matrix to identify Var1 * Var2 at group level
  vector[J] a;                   // group center
  vector<lower=0>[K] b;          // coefficients for modeling variation at pop level
  vector<lower=0>[K] g;          // coefficients for modeling variation at group level
transformed parameters{
  vector<lower=0>[N] sig;
  vector[N] aj;
  vector<lower=0>[J] sig_a;

  sig = X*b;
  aj = Z*a;
  sig_a = G*g;
  real mu_a = 0.0;
  a ~ normal(mu_a, sig_a);     // group locations
  b ~ normal(10,10);           // priors for population coefficients
  g ~ normal(25,25);           // priors for group coefficients
  y ~ normal(aj, sig);
generated quantities{
  real y_rep[N];                 // for posterior predictive checks
                                 // simulated data based on model
  y_rep = normal_rng(aj, sig);

Is there something I’m doing wrong that’s not allowing me to get the bs and gs correctly when my X ang G matrices are of the second form shown above?

Equivalently, in brms I can do:

    y ~ 0 + (1 | gr(J, by = interaction(Var1, Var2) )),
    sigma ~ 0 + interaction(Var1, Var2)This text will be hidden


    y ~ 0 + (1 | gr(J, by = interaction(Var1, Var2) )),
    sigma ~ 0 + Var1 + Var2   # I want something like this at both the group and population level

But in having to use interaction I am not able to get the coefficients for the factors.

simulated_stan_forums_model.Rmd (11.6 KB)
correct but not factor coefs.pdf (495.6 KB)
not correct.pdf (496.3 KB)

I found my errors… Forcing b and g to be positive.