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 b
s and g
s 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 b
s and g
s 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:
data{
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
}
parameters{
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;
}
model{
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 b
s and g
s correctly when my X
ang G
matrices are of the second form shown above?
Equivalently, in brms I can do:
bf(
y ~ 0 + (1 | gr(J, by = interaction(Var1, Var2) )),
sigma ~ 0 + interaction(Var1, Var2)This text will be hidden
)
or
bf(
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)