I am trying to fit a multivariate linear mixed (MLMM) for the analysis of multiple response |\underbrace{\mathbf{y}}_{qn\times 1} = \boldsymbol{X\beta} + \boldsymbol{\underbrace{Z}_{qn\times qm}\underbrace{\gamma}_{qm\times 1}} + \boldsymbol{\varepsilon},~~ ||\boldsymbol{\gamma}\sim \mathcal{N}(\textbf{0},I_m\otimes\mathbf{G}),~~| ||\boldsymbol{\varepsilon}\sim \mathcal{N}(\textbf{0},I_n\otimes\mathbf{R})
My stan program is too slow. The Gradient evaluation took 0.136 seconds or 1000 transitions using 10 leapfrog steps per transition would take 1360 seconds. That is too long. How can I improve. My r code and stan code are below
q<-3
m<-20
cs=12
n=cs*m
### Design matrix for the fixed part for y
### contains constant of 50 and for simplicity a randomly-drawn
### covariate from a standard Normal distribution.
intercept = 50
cont_var<-rnorm(q*n,0,1)
X <- as.matrix(cbind(intercept,cont_var))
beta=c(1,1)
fixed_outcome = X %*% beta
G_j = matrix(c(1,.51,.34,.51,1,.26,.34,.26,1), ncol = 3)
#diag_mat_tau =0.25*diag(3) diag_mat_tau %*% R %*% diag_mat_tau # VCV
G = kronecker(diag(m),G_j)
random_eff = mvrnorm(1, rep(0,q*m), G)
b = as.matrix(random_eff)
Z <- matrix(0,q*n,q*m)
j=1
vect1 <- numeric(n)
vect2 <- numeric(n)
vec1 <- numeric(q*n)
vec <- numeric(q*n)
vec2 <- numeric(q*n)
while(j <= q*n-1){
# for(j in 1:(response*students-1)){
ind <- sample(q*m,1,rep=F)
ind2 <- sample(q*m,1,rep=F)
Z[j,ind] <- 1
Z[j+1,ind2] <- 1
Z[j+1,ind] <- 0
T[j,1] <- ind
T[j+1,1] <- ind2
vec[j]=ind
vec2[j]=ind2
j=j+2
}
random_outcome = Z%*%b
##Error part
R_i=matrix(c(1.0, 0.7, 0.49,0.7, 1.0, 0.7,0.49, 0.7, 1.0), ncol = 3)
R=kronecker(diag(n),R_i)
epsilons=mvrnorm(1,mu = rep(0, q*n), Sigma =R )
y=fixed_outcome+random_outcome +epsilons
functions {
matrix block_diag(matrix qq, int n) {
int r = 1;
int sz = rows(qq);
int Q = sz * n;
matrix[Q, Q] G;
for (i in 1:Q) {
for (j in 1:Q) {
G[i, j] = 0;
}
}
while (r < Q) {
G[r:(r + sz - 1), r:(r + sz - 1)] = qq;
r += sz;
}
return G;
}
}
data {
int alpha_low; // alpha boundary value
int alpha_up; // alpha boundary value
int<lower=1> n; // number of students
int<lower=1> q; //number of responses
int<lower=1> m; //number of random effects
int<lower=1> qm;
int<lower=1> qn;
int<lower=1> P; //number of parameters
matrix[qn, P] X; // Fixed effects design matrix
vector[qn] y; // response variable
int vec[qn] ;
int vec2[qn] ;
}
parameters {
vector[P] beta;
vector[qm] gamma;
corr_matrix[q] Omega_G;
corr_matrix[q] Omega_R;
vector<lower=0>[q] sigma_G;
vector<lower=0>[q] sigma_R;
real<lower=alpha_low-0.01, upper=alpha_up+0.01> alpha;
}
transformed parameters {
matrix[q,q] G_j; // VCV matrix
matrix[q,q] R_i; // VCV matrix
matrix[qm, qm] G;
matrix[qn, qn] R;
matrix[qn,qm] Z=rep_matrix(0,qn,qm); //subj ranef design matrix
G_j = quad_form_diag(Omega_G, sigma_G);
R_i = quad_form_diag(Omega_R, sigma_R);
G=block_diag(G_j,m);
R=block_diag(R_i,n);
for(j in 1:n){
Z[2*j-1,vec[2*j-1]]= 1;
Z[2*j,vec2[2*j-1]] =1;
Z[2*j,vec[2*j-1]] = alpha;}
}
model {
beta ~ normal(0, 10);
sigma_G ~ student_t(4,0,1);
sigma_R ~ student_t(4,0,1);
Omega_G ~ lkj_corr(1.0);
Omega_R ~ lkj_corr(2.0);
gamma ~ multi_normal(rep_vector(0, qm), G);
y ~ multi_normal(X * beta + Z * gamma, R);
}