Please share your Stan program and accompanying data if possible.
When including Stan code in your post it really helps if you make it as readable as possible by using Stan code chunks (```stan) with clear spacing and indentation. For example, use
//
data {
int<lower=0> N; //number of observations
int T; //number of timepoint
int<lower=0> Ngroups; //number of unique subjects
int subjects[N]; //60
vector[N] treatment;
int familyRichness[N];
real IgA[N];
matrix[N, T] X_t; //intercept per time point
matrix[N, T] X_trt; // treatment effect per time point
matrix[N, T] X_fam; //family richness per timepoint
}
// The parameters accepted by the model. Our model
parameters {
vector[T] mux; //intercept per time point for poisson model
vector[T] muy; //intercept per time point for normal model
vector[T] beta; // treatment effect per time point for normal
vector[T] alpha; // treatment effect for poisson model
vector[2] phi;
vector[T] delta;
real<lower=0> p4;
real<lower=0> sigma;
real<lower=0> sigma_iga; // random intercept variance for iga
real<lower=0> sigma_rich; // random intercept variance for richness
cholesky_factor_corr[2] L_iga_rich;
matrix[2, Ngroups] z_iga_rich; //random effect from N(0, 1) to post multiply
}
transformed parameters {
real<lower=0> lambda[N];
vector[T] phis;
real mu[N];
matrix[2, Ngroups] ran_effs;
vector[2] L_sigma;
vector[T] gamma;
L_sigma[1] = sigma_iga;
L_sigma[2] = sigma_rich;
phis[1] = 0;
phis[2] = 0;
phis[3] = phi[1];
phis[4] = phi[2];
for(i in 1:4)
gamma[i] = phi[i]*delta[i];
ran_effs = diag_pre_multiply(L_sigma, L_iga_rich) * z_iga_rich;
for(i in 1:N) {
lambda[i] = exp(X_t[i, ] * mux + X_trt[i, ] * alpha + ran_effs[2, subjects[i]]);//poisson
mu[i] = X_t[i, ] * muy + X_trt[i] * beta + X_fam[i, ]* gamma + ran_effs[1, subjects[i]]; //normal
}
}
model {
int p1;
int p2;
L_iga_rich ~ lkj_corr_cholesky(2.0);
to_vector(z_iga_rich) ~ normal(0, 1);
sigma ~ cauchy(0, 5);
sigma_iga ~ cauchy(0, 5);
sigma_rich ~ cauchy(0, 5);
mux ~ cauchy(0, 10);
muy ~ cauchy(0, 10);
alpha ~ normal(0, 1);
beta ~ normal(0, 1);
delta ~ normal(0, 1);
phi[1] ~ bernoulli(0.5);
phi[2] ~ bernoulli(p4);
p4 ~ uniform(0,1);
for(i in 1:N) {
familyRichness[i] ~ poisson(lambda[i]);
IgA[i] ~ normal(mu[i], sigma);
}
}
generated quantities {
real rho;
matrix[2, 2] Sigma;
//compute sigma
Sigma = diag_pre_multiply(L_sigma, L_iga_rich * L_iga_rich');
//compute rho
rho = Sigma[1, 2] / (pow(Sigma[1, 1], 0.5) * pow(Sigma[2, 2], 0.5));
}
X_{ij} \sim dpois(\lambda_{ij}), log(\lambda_{ij}) = \mu_X[j] + b1_i + \alpha[j] Z_i and
Y_{ij} \sim dnorm(\mu_{ij},\sigma), \mu_{ij} = \mu_Y[j] + b2_i + \beta[j] Z_i + \phi[j]*\delta[j] X_{ij}.
The \phi[j]'s are the challenge here. I fixed its first two components to be 0, while the last two components are to be simulated from the Bernoulli distribution. I understand that stan does not allow for declaration of discrete parameters and one has to marginalize over discrete parameters, but how do I do this in my case? Thanks in anticipation.