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

Modeling
1

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.

https://doi.org/10.3102/1076998614559412
https://doi.org/10.3102/10769986031002181

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[2]              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[1],cov[1,1]);
    lambda  ~ normal(mu[2],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]);
        }  
     } 
}
'
2

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.

3

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

4

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).

5

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[2] mu_ipar;              // mean vector for phi and lambda
  vector<lower=0>[2] sigma_ipar;  // sd vector for phi and lambda
  corr_matrix[2] rho;             // correlation matrix between phi and lambda
	
  vector[2] 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);
   ipar         ~ multi_normal(mu_ipar,quad_form_diag(rho, sigma_ipar)); 

   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