Specifying a two-process LATER model

Modeling
,
1

Hello Stan community!

I am trying to fit a “two-horse” linear approach to threshold with ergodic rate (LATER) model to reaction time data (similar to this paper: Adam, Bays & Husain, Cognitive Neuroscience 2012).

The model posits two recinormal processes racing to generate a response. The anticipatory process ( A ) starts at t=0 and the reactive process ( R ) starts at t=t_0. Whichever process reaches threshold first triggers the response (“OR” rule). The resulting cumulative probability distribution is given by:

Pr(T \leq t) = \Psi(t | \mu_A, \sigma_A) + \Psi(t-t_0 | \mu_R, \sigma_R)-\Psi(t | \mu_A, \sigma_A)\Psi(t-t_0 | \mu_R, \sigma_R)

where \Psi are the recinormal CDF, defined as \Psi(t | \mu_i, \sigma_i) = 1 - \Phi(\frac{1}{t} | \mu, \sigma) based on the normal CDF \Phi.

The likelihood function is obtained by differentiating the CDF (note: it is not a mixture of PDFs):

Pr(t) = \psi(t | \mu_A, \sigma_A) + \psi(t - t_0 | \mu_R, \sigma_R) - \Psi(t | \mu_A, \sigma_A)\psi(t-t_0 | \mu_R, \sigma_R) - \Psi(t-t_0 | \mu_R, \sigma_R)\psi(t | \mu_A, \sigma_A)

where \psi is the recinormal PDF, defined as \psi(t | \mu_i, \sigma_i) = \phi(\frac{1}{t} | \mu, \sigma)/t^2 based on the normal PDF \phi.

I have written the model as follows:

data {  
int<lower=1> P; // number of participants 
int<lower=1> N; // total number of datapoints 
int<lower=1,upper=P> subjID[N]; // person index for all the datapoints
vector[N] G; // t0 (in ms)
vector[N] Y; // all RT data (in ms)
} 

parameters {  
vector<lower=0>[P] mu_anticip;
vector<lower=0>[P] sigma_anticip;
vector<lower=0>[P] mu_react;
vector<lower=0>[P] sigma_react;
real<lower=0> MU_mu_anticip;
real<lower=0> MU_sigma_anticip;
real<lower=0> MU_mu_react;
real<lower=0> MU_sigma_react;
real<lower=0> SIGMA_mu_anticip;
real<lower=0> SIGMA_sigma_anticip;
real<lower=0> SIGMA_mu_react;
real<lower=0> SIGMA_sigma_react;
}

model {
// Priors
mu_anticip ~ normal(MU_mu_anticip,SIGMA_mu_anticip);
sigma_anticip ~ normal(MU_sigma_anticip,SIGMA_sigma_anticip);
mu_react ~ normal(MU_mu_anticip,SIGMA_mu_anticip);
sigma_react ~ normal(MU_sigma_anticip,SIGMA_sigma_anticip);

// Likelihood
for (i in 1:N) {
Y[i] ~ twoHorses_lpdf(G[i],mu_anticip[subjID[i]],sigma_anticip[subjID[i]],mu_react[subjID[i]],sigma_react[subjID[i]]);
}
}

However, I am struggling to define the log probability function twoHorses_lpdf. I think I have managed to write log(\Psi(t | \mu_i, \sigma_i)) and log(\psi(t | \mu_i, \sigma_i)) properly as:

functions {
  real recinormal_lpdf(real y, real mu, real sigma) {
    if (y>0)
      return normal_lpdf(1/y, mu, sigma)-2*log(y);
    else
      return negative_infinity();
  }
  real recinormal_lcdf(real y, real mu, real sigma) {
    if (y>0)
      return  log(1-normal_cdf(1/y, mu, sigma));
    else
      return negative_infinity();
  }
}

But I don’t know how to define the log probability function based on the equations above. Could you provide some guidance on how to do so?

Thanks a lot!

Pierre

2

This is my attempt based on the formula.

for (i in 1:N){
  real loglik [4];
  real diff_t = Y[i] - G[i];
  real loglik_ant = recinormal_lpdf(Y[i] | mu_anticip[subjID[i]], sigma_anticip[subjID[i]]);
  real loglik_react = recinormal_lpdf(diff_t | mu_react[subjID[i]], sigma_react[subjID[i]]);
  real cumul_ant = recinormal_lcdf(Y[i] | mu_anticip[subjID[i]], sigma_anticip[subjID[i]]);
  real cumul_react = recinormal_lcdf(diff_t | mu_react[subjID[i]], sigma_react[subjID[i]]);
  loglik[1] = loglik_ant;
  loglik[2] = loglik_react;
  // product -> sum on the log scale. You can drop the -1 because it is a constant.
  loglik[3] = - 1 + cumul_ant + loglik_react;
  loglik[4] =  -1 + cumul_react + loglik_ant;
  // for numerical stability
  target += log_sum_exp(loglike);
}

This is quite fiddly and I probably have typos in there and might have misunderstood the formula but I hope it helps you in the right direction. You can also probably speed this up with some clever vectorization but I think you want to make sure (with some simulations) that you get the correct answer first.