Thank you for the suggestions. After reviewing the suggestions for a similar post I did years ago for a reduced version of this model, and your additional tips, this is what I came up with. I will test it and see how much faster it runs and how much difference it makes in memory usage.
// Stan model syntax combining the Modified DG-IRT for response accuracy and
// DG-LNRT for response time
data{
int <lower=1> J; // number of examinees
int <lower=1> I; // number of items
int <lower=1> n_obs; // number of observations (J x I - missing responses)
array[n_obs] int<lower=1> p_loc; // person indicator
array[n_obs] int<lower=1> i_loc; // item indicator
array[n_obs] real RT; // vector of log of responses
array[n_obs] int<lower=0,upper=1> Y; // vector of item responses
}
parameters {
// Parameters common for both response time and response accuracy component
vector<lower=0,upper=1>[I] pC; // vector of length I for the probability of item compromise status
vector<lower=0,upper=1>[J] pH; // vector of length J for the probability of examinee item peknowledge
// Parameters for the item parameters component
real mu_beta; // mean for time intensity parameters
real mu_alpha; // mean for log of time discrimination parameters
real mu_b; // mean for item difficulty parameters
vector<lower=0>[3] sigma_I; // vector of standard deviations for item parameters
cholesky_factor_corr[3] Lomega_I; // Cholesky factors of 3x3 correlation matrix for item parameters
// Parameters for the person parameters component
real mu_thetat; // mean for theta_t
vector<lower=0>[2] sigma_P; // vector of standard deviations for latent trait parameters
cholesky_factor_corr[2] Lomega_P; // Cholesky factors of 2x2 correlation matrix for person parameters
real<lower=0> mu_delta_tau; // mean for change from taut to tau_c
real<lower=0> mu_delta_theta; // mean for change from thetat to thetac
vector<lower=0>[2] sigma_D; // vector of standard deviations for change in person parameters
cholesky_factor_corr[2] Lomega_D; // Cholesky factors of 2x2 correlation matrix for person parameters
array[I] vector[3] item;
array[J] vector[2] person;
array[J] vector<lower=0>[2] delta;
}
transformed parameters{
vector[2] mu_P = [0,mu_thetat]'; // vector for mean of person parameters
// mu_taut is fixed to 0
vector[3] mu_I = [mu_alpha,mu_beta,mu_b]'; // vector for mean of item parameters
vector[2] mu_D = [mu_delta_tau,mu_delta_theta]'; // vector for mean of change parameters
// The sum of b parameters is fixed to 0 for the response accuracy component
// this is for model identification for the response accuracy component
vector[I] b_star;
for(i in 1:I){
b_star[i]= item[i][3];
}
vector[I] b = b_star;
b[1] = -sum(b_star[2:I]);
}
model{
sigma_P ~ exponential(1);
mu_thetat ~ normal(0,1);
Lomega_P ~ lkj_corr_cholesky(1);
sigma_I ~ exponential(1);
mu_beta ~ normal(4,1);
mu_alpha ~ lognormal(0,0.5);
mu_b ~ normal(0,1);
Lomega_I ~ lkj_corr_cholesky(1);
sigma_D ~ exponential(1);
mu_delta_tau ~ normal(0,1);
mu_delta_theta ~ normal(0,1);
Lomega_D ~ lkj_corr_cholesky(1);
person ~ multi_normal_cholesky(mu_P,diag_pre_multiply(sigma_P, Lomega_P));
item ~ multi_normal_cholesky(mu_I,diag_pre_multiply(sigma_I, Lomega_I));
delta ~ multi_normal_cholesky(mu_D,diag_pre_multiply(sigma_D, Lomega_D));
pC ~ beta(1,1);
pH ~ beta(1,1);
// Joint density of response time and response accuracy
vector[I] log1m_pC = log1m(pC) ;
vector[I] log_pC = log(pC) ;
vector[J] log1m_pH = log1m(pH) ;
vector[J] log_pH = log(pH) ;
vector[n_obs] norm_p_t ;
vector[n_obs] norm_p_c;
vector[n_obs] bern_p_t ;
vector[n_obs] bern_p_c;
for (i in 1:n_obs) {
// item[i_loc[i]][2] is the beta parameter for the (i_loc[i])th item
// item[i_loc[i]][1] is the alpha parameter for the (i_loc[i])th item
// person[p_loc[i]][1] is tau_t for the (p_loc[i])th person
// (person[p_loc[i]][1]+delta[p_loc[i]][1]) is tau_c for the (p_loc[i])th person
real p_taut = item[i_loc[i]][2] - person[p_loc[i]][1];
real p_tauc = item[i_loc[i]][2] - (person[p_loc[i]][1]+delta[p_loc[i]][1]);
real alpha = 1/exp(item[i_loc[i]][1]);
norm_p_t[i] = normal_lpdf(RT[i] | p_taut, alpha);
norm_p_c[i] = normal_lpdf(RT[i] | p_tauc, alpha);
// b[i_loc[i]] is the b parameter for the (i_loc[i])th item
// person[p_loc[i]][2] is theta_t for the (p_loc[i])th person
// (person[p_loc[i]][2]+delta[p_loc[i]][2]) is theta_c for the (p_loc[i])th person
real p_thetat = person[p_loc[i]][2] - b[i_loc[i]];
real p_thetac = (person[p_loc[i]][2]+delta[p_loc[i]][2]) - b[i_loc[i]];
bern_p_t[i] = bernoulli_logit_lpmf(Y[i] | p_thetat);
bern_p_c[i] = bernoulli_logit_lpmf(Y[i] | p_thetac);
}
matrix[n_obs,4] lprt;
lprt[,1] = log1m_pC[i_loc] + log1m_pH[p_loc] + norm_p_t; // T = 0, C=0
lprt[,2] = log1m_pC[i_loc] + log_pH[p_loc] + norm_p_t; // T = 1, C=0
lprt[,3] = log_pC[i_loc] + log1m_pH[p_loc] + norm_p_t; // T = 0, C=1
lprt[,4] = log_pC[i_loc] + log_pH[p_loc] + norm_p_c; // T = 1, C=1
matrix[n_obs,4] lpr;
lpr[,1] = log1m_pC[i_loc] + log1m_pH[p_loc] + bern_p_t; // T = 0, C=0
lpr[,2] = log1m_pC[i_loc] + log_pH[p_loc] + bern_p_t; // T = 1, C=0
lpr[,3] = log_pC[i_loc] + log1m_pH[p_loc] + bern_p_t; // T = 0, C=1
lpr[,4] = log_pC[i_loc] + log_pH[p_loc] + bern_p_c; // T = 1, C=1
vector[n_obs] lps;
for (i in 1:n_obs) {
lps[i] = log_sum_exp(lprt[i,]) + log_sum_exp(lpr[i,]);
}
}