I need to put a period to this thread. Below I am posting a Stan code chunk for running a multilevel logit model with multivariate priors, hopefully without any issue. Note that I silence the lines containing array, since I haven’t figured it out yet and will be working on it. Will update once I have solid answers.
stanDiscoursePost20220402.R (3.4 KB)
Attached please also find an updated R script file. Thank you all for help!
data {
int<lower=0> N; // num individuals
int<lower=1> K; // num ind predictors
int<lower=1> J; // num groups
int<lower=1> L; // num group predictors
// array[N] int<lower=1, upper=J> jj; // group for individual
int<lower=1, upper=J> jj[N]; // group for individual
matrix[N, K] x; // individual predictors
row_vector[L] u[J]; // group predictors
// vector[N] y; // outcomes
int<lower=0, upper=1> y[N];
// vector<lower=0, upper=1>[N] y;
}
parameters {
corr_matrix[K] Omega; // prior correlation
vector<lower=0>[K] tau; // prior scale
matrix[L, K] gamma; // group coeffs
// array[J] vector[K] beta; // indiv coeffs by group
vector[K] beta[J];
real<lower=0> sigma; // prediction error scale
}
model {
tau ~ cauchy(0, 2.5);
Omega ~ lkj_corr(2);
to_vector(gamma) ~ normal(0, 5);
{
// array[J] row_vector[K] u_gamma;
row_vector[K] u_gamma[J];
for (j in 1:J) {
u_gamma[j] = u[j] * gamma;
}
beta ~ multi_normal(u_gamma, quad_form_diag(Omega, tau));
}
for (n in 1:N) {
real eta;
real prob;
eta <- x[n] * beta[jj[n]];
prob <- 1/(1+exp(-1*eta));
y[n] ~ bernoulli(prob);
// y[n] ~ inv_logit(x[n] * beta[jj[n]]);
// y[n] ~ bernoulli_logit(x[n] * beta[jj[n]]);
}
}