I am working on double hierarchical generalized model with data simulated.
The model I set was
y ~ normal(mean, variance)
mean = Xbeta + Zu1
where X and Z are design matrices, beta is coefficients for fixed effects and u1 is coefficients for random effects.
Similarly,
log(variance) = phi + Z*u2
phi is just a intercept term and u2 is coefficients for random effects. Before defining priors for parameters, I used the idea from Gelman (2008) that u1 and u2 can be
u1[i] = dummy1 * Xi1[i]
u2[i] = dummy2 * Xi2[i]
Each element of Xi1 has normal(0, var = tauXi1), each element of X2 has normal(0, var = 1/tauXi2). Indeed, each of u1 and u2 have normal(0, tauMu) and normal(0, tauSigma) where tauMu = absolute(dummy1)/sqrt(tauXi1) and tauSigma = absolute(dummy2)/sqrt(tauSigma).
I am sampling tauXi1 and tauXi2 from gamma(1.5, 37.5) and dummy1 and dummy2 from normal(0,1). For other parameters, I used non-informative priors normal(0, 100^2).
Here is my stan code:
data {
int<lower=0> N;
int<lower=0> Nx1;
int<lower=0> Nx2;
int<lower=0> Nz;
matrix[N,Nx1] x1;
matrix[N,Nx2] x2;
matrix<lower=0,upper=1>[N,Nz] z;
real y[N];
vector[Nz] mvnMu;
matrix<lower=0>[Nz,Nz] mvnDiag;
}
parameters {
vector[Nz] u1;
vector[Nz] u2;
vector[Nx1] beta;
vector[Nx2] phi;
real dummy1;
real dummy2;
real<lower=0> invtauXi1Sq;
real<lower=0> invtauXi2Sq;
}
transformed parameters {
vector[N] mu;
vector[N] logSigma;
vector[N] sigma;
real<lower=0> tauMu;
real<lower=0> tauSigma;
matrix<lower=0>[Nz,Nz] mvnSigma1;
matrix<lower=0>[Nz,Nz] mvnSigma2;
mu = multiply(x1, beta) + multiply(z, u1);
logSigma = multiply(x2, phi) + multiply(z, u2);
sigma = exp(logSigma);
tauMu = abs(dummy1)/sqrt(invtauXi1Sq);
tauSigma = abs(dummy2)/sqrt(invtauXi2Sq);
mvnSigma1 = multiply(pow(tauMu,2), mvnDiag);
mvnSigma2 = multiply(pow(tauSigma,2), mvnDiag);
}
model {
for (i in 1:N) {
y ~ normal(mu[i], sigma[i]);
}
beta ~ normal(0, 100);
phi ~ normal(0, 100);
dummy1 ~ normal(0,1);
dummy2 ~ normal(0,1);
invtauXi1Sq ~ gamma(1.5, 37.5);
invtauXi2Sq ~ gamma(1.5, 37.5);
u1 ~ multi_normal(mvnMu, mvnSigma1);
u2 ~ multi_normal(mvnMu, mvnSigma2);
}
mvnDiag is just an identity matrix and mvnMu is vector of zero.
This stan code takes too long to complete! Can someone help me to improve running time?