# Syntax for Lognormal Response Time Model - how to specify priors for parameters with multivariate normal distribution

Hi,

I am trying to fit the following model in Stan. I have some idea, but I think it is incomplete. I will try to define the model below and then write the incomplete syntax I have now. I appreciate any input and suggestion. These models are described in papers before for reference.

The logarithm of the response times, Y_{ij} is the response time for the ith person on the jth item, are assumed to be normally distributed with a mean of \phi_j*(\lambda_j - \zeta_i) and variance of \sigma_j^2.

\phi_j : item specific time discrimination parameter
\lambda_j : item specific time intensity parameter
\zeta_i : individual specific latent working speed parameter
\sigma_j^2 : item specific measurement error variance

For model identification, suppose we fix the mean and variance of working speed parameter to 0 and 1. So, I assume \zeta_i ~ N(0,1).

In these papers, they also assume that the item specific parameters, \phi and \lambda, are correlated, and a bivariate normal distribution is assumed to describe the relationship.

\left ( \begin{matrix} \phi_j\\\lambda_j \end{matrix} \right ) \sim N\left (\left [ \begin{matrix} \mu_\phi\\ \mu_\lambda \end{matrix} \right ],\left [\begin{matrix} \sigma_\phi^2 & \sigma_{\phi,\lambda} \\ \sigma_{\lambda,\phi} & \sigma_\lambda^2 \end{matrix} \right ] \right )

Below is my effort to fit this model, but I know it is incomplete and probably wrong. I am struggling to translate this model to Stan language. Especially, the covariance matrix among parameters and how to specify priors for them (normal inverse wishart as suggested in one of these papers).

rt <- '

data {
int <lower=0> I;                  //  number of individuals
int <lower=0> J;                  //  number of items
int <lower=0,upper=1> RT[I,J];    //  matrix of response times
}

parameters {
vector[J]              lambda;   //item specific time intensity parameter
vector <lower=0> [J]   phi;      //item speficic time discrimination parameter
vector <lower=0> [J]   sigma;    // item specific residual standard deviation,
vector[I]              zeta;     // person specific latent working speed
vector              mu;       // mean vector for lambda and phi
matrix[2,2]            cov;       // covariance matrix between lambda and phi
}

model {
zeta   ~ normal(0,1);

//mu and cov priors??????

phi     ~ normal(mu,cov[1,1]);
lambda  ~ normal(mu,cov[2,2]);

sigma   ~ inv_gamma(5,5);         // would this be appropriate??

for (i in 1:I){
for (j in 1:J){
real p = phi[j]*(zeta[i]-lambda[j]);
RT[i,j] ~ normal(p,sigma[j]);
}
}
}
'

Check out the manual chapter on regression. It explains our preferred approach to multivariate priors, which doesn’t involve Wisharts, but rather independent priors on the correlations and on the scales.

You need constrained types for things like covariance. This is also explained in the manual. But why define a covariance matrix if you’re just going to ignore it and sample phi and lambda from teh diagonal?

We don’t recommend inverse gammas, but (5, 5) is probably OK.

You can vectorize that RT sampling statement. You probably need to log it first.

Hi，have you solved your problem? Would you mind telling me how to specify priors for parameters with multrivarite normal distribution?

You coul start from the Stan user’s guide page on multivariate outcomes where there are some examples that can get you started.

If you have a specific problem ask again, there will be someone able to help (be patient this days as traffic is slower then usual).

This was a while ago but thanks for asking so we can close this.

I figured out the multivariate distribution piece. For that part, this post was very helpful in addition to manual.

http://stla.github.io/stlapblog/posts/StanLKJprior.html

Below is a code for the model described above. Hope this helps. Note that there is still an issue to be fixed in the code below. \phi (the first column of the ipar matrix) must be constrained to be positive. I am not sure how to do it. Once you do that, the code below should work for the model described above.

Also note that since then I am using a simplified version of this model by fixing all \phi parameters to 1 across items. van der Linden (2016) argues that \phi is unnecessary in this model and it is like overparameterization. If I have time, I will try post a fully producible example with simulated data .

Also, as Bob mentioned, you can vectorize the sampling statement for Y below.

Hope this helps a little bit.

code_rt<-'

data{
int  J;                  //  number of items
int  I;                  //  number of individuals in condition1
real Y[I,J];             //  matrix of log of response times
}

parameters {

vector mu_ipar;              // mean vector for phi and lambda
vector<lower=0> sigma_ipar;  // sd vector for phi and lambda
corr_matrix rho;             // correlation matrix between phi and lambda

vector ipar[J];            // time intensity and time discrimination parameter matrix

vector <lower=0> [J]  sigma;   // item specific residual error,
real mu_sigma;                 // average error
real sigma_sigma;              // standard deviation of error

vector[I] zeta;                    // person specific latent working speed
}

model{

zeta         ~ normal(0,1); // mean and variance are fixed to 0 and 1 for model identification

mu_ipar      ~ normal(0,5);
sigma_ipar   ~ cauchy(0,2.5);
rho          ~ lkj_corr(1);

mu_sigma    ~ normal(0,5);
sigma_sigma  ~ cauchy(0,2.5);
sigma       ~ normal(mu_sigma,sigma_sigma);

for (i in 1:I) {
for (j in 1:J) {

real p = ipar[j,1]*(ipar[j,2]-zeta[i]);
Y[i,j] ~ normal(p,sigma[j]);
}}

}
'


1 Like