Modeling reaction time data with Weibull hierarchical model

[Xiao He]

Hi, 

I've been trying to fit some reaction time data with a Weibull hierarchical model. I tried to follow the Kidney example here. Below is my attempt. A complete example, with the relevant R code and dataset is here. I wonder if I am on the right track or if I am completely off. Also, since I am trying to implement a varying-intercept, varying-slope model. I wonder if the shape parameter should also vary based on condition, subject, and item. Thanks!

data{
 int<lower=1> N;
 int<lower=1> NS;
 int<lower=1> NI;
 int<lower=1, upper=NS> subject[N]; //subject ID
 int<lower=1, upper=NI> item[N];    //item ID
 vector[N] cond1;        //dummy coded variable
 vector[N] cond2;        //dummy coded variable
 vector<lower=0>[N] rt;  //reaction time
 vector[3] zero;
}

parameters{
 vector[3] beta;
 vector[3] u[NS];     //subject random effects
 vector[3] w[NI];     //item random effects
 vector<lower=0>[3] sigma_u;
 vector<lower=0>[3] sigma_w;
 corr_matrix[3] omega_u;
 corr_matrix[3] omega_w;
 real<lower=0> shape;
}

transformed parameters{
 cov_matrix[3] Sigma_u;
 cov_matrix[3] Sigma_w;
 for(r in 1:3){
   for(c in 1:3){
     Sigma_u[r, c] <- sigma_u[r] * sigma_u[c] * omega_u[r, c];
     Sigma_w[r, c] <- sigma_w[r] * sigma_w[c] * omega_w[r, c];
   }
 }
}

model{
 vector<lower=0>[N] scale;
 for(j in 1:NS){
   u[j] ~ multi_normal(zero, Sigma_u);
 }
 for(k in 1:NI){
   w[k] ~ multi_normal(zero, Sigma_w);
 }
 sigma_u ~ gamma(1.5, 1.0E-4);
 sigma_w ~ gamma(1.5, 1.0E-4);
 shape ~ gamma(1, 0.001);
 for(i in 1:N){
   scale[i] <- exp(-((beta[1] + u[subject[i], 1] + w[item[i], 1]) +
                     (beta[2] + u[subject[i], 2] + w[item[i], 2])*cond1[i] +
                     (beta[3] + u[subject[i], 3] + w[item[i], 3])*cond2[i])/shape);
 }

 rt ~ weibull(shape, scale);
}

-Xiao

[Bob Carpenter]

I don't know enough about the specifics of these kinds of models
to comment on the data model.

That model was translated from BUGS and we'd actually recommend
different priors --- there's a discussion in the Stan manual
chapter on regression spread over three sections.

The biggest problem you'll have to worry about is identifiability.
For instance, you have an additive non-identifiability in

> for(i in 1:N){
> scale[i] <- exp(-((beta[1] + u[subject[i], 1] + w[item[i], 1]) +
> (beta[2] + u[subject[i], 2] + w[item[i], 2])*cond1[i] +
> (beta[3] + u[subject[i], 3] + w[item[i], 3])*cond2[i])/shape);
> }

with beta[1], u[,1] and w[,1], which is somewhat mitigated by drawing
u and w from priors centered on zero.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Shravan Vasishth]

Hi Xiao He,

In our research, we use the log-normal in a similar setting; Jeff Rouder suggested to me that this might make more sense generally than the Weibull. See our tutorial specifically for these kinds of problems:

http://www.ling.uni-potsdam.de/~vasishth/statistics/BayesLMMs.html

Also, you can fit your models much faster if you use LKJ priors as described in the manual.  A fully working example is here:

https://github.com/vasishth/StanJAGSexamples/blob/master/FrankEtAlCogSci2015/FTVCogSci.Stan

This is from a paper that is close to being accepted, pending minor revisions (Cognitive Science):

 @unpublished{FrankEtAl2014,
      author = {Stefan L. Frank and Thijs Trompenaars and Shravan Vasishth},
      title = {Cross-linguistic differences in processing double-embedded relative clauses: {W}orking-memory constraints or language statistics?},
      year = {2014},
      note = {submitted},
      pdf = {http://www.ling.uni-potsdam.de/~vasishth/pdfs/FrankTrompenaarsVasishthSubmitted.pdf}
    }

[Shravan Vasishth]

Of course, you could equally well log-transform the data and fit a normal model of the type we fit in the Frank et al example ;). Some people might feel more excited if they fit the "raw" data instead of "changing" it through a log-transform :). A reasonable alternative is to use the Box-Cox procedure to find the appropriate variance stabilizing transform and stay with the normal distribution. This is what we usually do.

[Andrew Ellis]

hi Shravan,

On Wednesday, December 31, 2014 9:50:16 AM UTC+1, Shravan Vasishth wrote:

Hi Xiao He,

In our research, we use the log-normal in a similar setting; Jeff Rouder suggested to me that this might make more sense generally than the Weibull.

Do you have any references discussing the log-normal vs. the Weibull distribution? We recently had a paper reviewed and one of the reviewers questioned our use of the lognormal and wanted a Weibull fit instead. 

Have you considered using a shifted log-normal to model response times? I guess this would require a custom probability function in Stan?

cheers,

Andrew

[Bob Carpenter]

> On Dec 31, 2014, at 4:35 PM, Andrew Ellis <a.w....@gmail.com> wrote:
>
> Have you considered using a shifted log-normal to model response times? I guess this would require a custom probability function in Stan?

If you have a constant c and write:

(y - c) ~ logormal(log_mu, sigma);

Then y has a shifted lognormal. It's a linear transform,
so there's no Jacobian adjustment necessary (ignore any warning
messages that might pop up).

- Bob

[Andrew Ellis]

wow, that was fast. Thanks very much...

But what if I need to estimate the shift parameter c, i.e.  c ~ uniform(0, min(y))?

[Bob Carpenter]

> On Dec 31, 2014, at 4:45 PM, Andrew Ellis <a.w....@gmail.com> wrote:
>
> wow, that was fast. Thanks very much...
>
> But what if I need to estimate the shift parameter c, i.e. c ~ uniform(0, min(y))?

If y and c are both parameters, the transform is just

(y, c) -> (y-c, c)

which gives you a lower triangular Jacobian with
a determinant of 1,

d(y-c)/dy dc/dy 1 0
=
d(y-c)/dc dc/dc -1 1

Just make sure to declare c with the relevant constraints

parameters {
...
real<lower=0, upper=min(y)> c;


If only c is a parameter, then the map is just

c -> (y-c)

and d(y-c)/dc = -1, so the absolute value's 1 and you're good to go
with no Jacobian adjustment.

- Bob

>
>
>
> On Wednesday, December 31, 2014 10:41:17 PM UTC+1, Bob Carpenter wrote:
>
> > On Dec 31, 2014, at 4:35 PM, Andrew Ellis <a.w....@gmail.com> wrote:
> >
> > Have you considered using a shifted log-normal to model response times? I guess this would require a custom probability function in Stan?
>
> If you have a constant c and write:
>
> (y - c) ~ logormal(log_mu, sigma);
>
> Then y has a shifted lognormal. It's a linear transform,
> so there's no Jacobian adjustment necessary (ignore any warning
> messages that might pop up).
>
> - Bob
>

[Shravan Vasishth]

Jeff Rouder's spent a lot of time thinking about this, it seems. See his publications for more. I found the attached paper (no title on it) by just googling around. 

 

--
Shravan Vasishth
Professor for Psycholinguistics and Neurolinguistics

Department of Linguistics
University of Potsdam, Germany
http://www.ling.uni-potsdam.de/~vasishth

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/zEKVSv6_K70/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Andrew Ellis]

thanks. In my case, y is the observed response time in a psychology experiment. The shift parameter should model the non-decision time, i.e. time to execute motor response, etc.

cheers,

Andrew

[Andrew Ellis]

thanks, I'll read it.

Andrew