Estimating z_Intercept

[John Clithero]

Hi everyone,

This topic has come up before, but I'm not sure it ever was resolved publicly.

My issue is this: I seem to have trouble estimating the z_Intercept parameter in HDDMRegression models. I can simulate data (intended to be within-subject) with varying starting points (example z_A=0.4, z_B=0.5, z_C=0.6). Regular HDDM is able to recover those parameters (as in the sample). However, when I set up the model in HDDMRegression, it fails to correctly identity the intercept term. In other words, it can identify z_B=0.5 and z_C=0.6, but fails to correctly identify z_A. If anything, it seems like the posterior has a bound that prevents the intercept estimate for z from going below 0.5.

Perhaps I have made an error in how I am setting up the HDDMRegression model? The link function?

I have constructed an example in a short notebook (attached) that hopefully illustrates the issue. It simulates data and then fits it using both HDDM and HDDMRegression. I show how it correctly recovers the z parameters in the first case but not the second. I also included the relevant posterior plots.

Any help would be very much appreciated!

Thanks,

John

[John Clithero]

Hi Thomas et al.,

Any thoughts on the issue I sent out a few weeks ago (below)? I'm very curious to see a response if this is something silly in my code or something in HDDM. If it is an issue then it would seem to affect anyone wanting to use HDDM Regression with the starting point parameter!

Thanks!

John

[Michael J Frank]

Hi John, 

It may be that the prior is biasing the estimate of Z at the intercept to be 0.5. The priors for the intercept in hddmRegressor are centered around the standard priors (which are 0.5 for z). The priors for the slope terms are centered at 0.  So it is recommended to mean center the values of the IV's that are used for the regressors so that any effect of the regressors would still allow for the mean value of the param in the intercept to be centered at the prior.  

(in your case you could check this by using cond B as the intercept and defining A and C relative to that). 

 Michael

Michael J Frank, PhD, Professor
Laboratory for Neural Computation and Cognition
Brown University
http://ski.clps.brown.edu
(401)-863-6872

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

[John Clithero]

Hi Michael,

Thanks, I finally got around to trying this and it does seem to have helped. However, the new posterior plot for z_intercept (attached) still gives me the impression that there is a bound at zero (or z=0.5)?

Cheers,

John

[Michael J Frank]

Hi John, there shouldn't be a bound at 0 in the actual trace, it looks like the minimum value is above zero and maybe this is just a small bias in the estimate of the data (since now the intercept should actually be 0.5 and the value close to 0). Can you check that the condition effect can now appropriately recover the case where generated z=0.4? 

  

Michael 

[John Clithero]

Hi Michael,

I am attaching an updated Python notebook. It looks like HDDMRegressor does a fairly good job of recovering when I switch the intercept to z=0.5 and make the other two conditions 0.4 and 0.6. As you noted it does seem like there might be a slight upward bias in the z_Intercept estimate (even with a pretty big synthetic dataset).

Thanks again!

John

[Michael J Frank]

Hi John, looks good. And actually I'm not sure there is a bias after all: 

note that when you generate data with gen_rand_data using multiple subjects, the parameters that you set (e.g. z_A = 0.5) are actually just specifying the means of a normal distribution, where each simulated subject has parameters sampled from this distribution, but with variance. So for smallish number of subjects your actual generated data may have actually had means that differ slightly from the prescribed mean. You can check by just looking at the values of sim_params to see how close they are to the recovered values (at individual or group level).

Michael

[Michael J Frank]

Hi John and all,

an important point has come to light about estimating regression models on the starting point, which I believe led to the issue you identified in this thread. (Ultimately you corrected it by choosing a different condition for your baseline, but it turns out that there was an actual problem in the first place, in contrast to what I thought.)

The issue is that your regression used the inverse logit transformation. This is meant to constrain z between 0 and 1, which is sensible, and is what we had used in our tutorial example for stimulus coded regression. But since then, HDDM now applies that transform to the prior for z (for all models with z, including the intercept of regression models), which constrains it. Applying the transform again in the link function leads to a bias (because the invlogit prior on an invlogit transformed variable then forces the intercept to be >0.5). 

So, with recent versions of HDDM, one should not use an inverse logit transform on z in the link function for regression models (unless they change the prior in the guts of their version of HDDM).  It is sufficient to use the constrained prior on the intercept with the regular linear link function (any extreme values of z that would go out of the [0,1] bounds on the full regression would get rejected by the sampler anyway). This also makes z regression coefficients more comparable to those for other models parameters (which are all usually linear), and  easier to interpret the coefficients.

I confirmed this was the culprit for your original case that you shared in this thread, by simply refitting the model to data in which the baseline condition has z<0.5 and using the identity link function for z (ie lambda x:x), and it works properly (recovers z intercepts below and above .5). I also confirmed that the original tutorial code with model recovery on stim-coded regression would fail if run as originally specified (due to the altered prior), but that it works properly when the link function is fixed (in that case still needing to apply the 1-z swap) The upshot  is that there is no bug in HDDM itself, but just that this particular tutorial code was outdated - just updated it now on github (the old docs are still on ski, will fix that).   

Note, anyone who has performed a regression on starting points before and used invlogit in the link function, if you are concerned you can check your results and if got even a single z<0.5 (after transformed back), for any subject/condition, your results should be ok (ie, it would mean you did it in an earlier version of HDDM without that prior, so that the invlogit was applied just once). 

h/t: Ian Krajbich and his lab for first noticing the problem with biased estimates- thanks!

Michael

Michael J Frank, PhD | Edgar L. Marston Professor

Director, Carney Center for Computational Brain Science 

Laboratory of Neural Computation and Cognition  

Brown University
website 

Michael 

To unsubscribe from this group and stop receiving emails from it, send an email to hddm-users+...@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to the Google Groups "hddm-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to hddm-users+...@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to the Google Groups "hddm-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to hddm-users+...@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to the Google Groups "hddm-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to hddm-users+...@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to the Google Groups "hddm-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to hddm-users+...@googlegroups.com.

[valer...@gmail.com]

Hi Michael,

I am a bit confused and I just want to make sure I get this bit right.

Note, anyone who has performed a regression on starting points before and used invlogit in the link function, if you are concerned you can check your results and if got even a single z<0.5 (after transformed back), for any subject/condition, your results should be ok (ie, it would mean you did it in an earlier version of HDDM without that prior, so that the invlogit was applied just once).

I used the HDDM (version 0.8.0) and estimated z in a regression formula with the following link function:

def z_link_func( x, data = data ):

        stim = ( dmatrix( "0 + C(s, [[1], [-1]])",

                                     { 's': data.stim.loc[x.index] },

                                     return_type = 'dataframe' )

                     )

        

        return 1 / ( 1 + np.exp( -np.multiply(x.to_frame(), stim) ) )

In the output data frame, one of the variables is z_Intercept_trans, which should contain the estimates of the decision bias after the inverse logit function has been applied.

In more than 95% of the cases, z_Intercept_trans < 0.5. However, I am not sure I performed the sanity check you suggested correctly, as the variable name (i.e., z_Intercept_trans) suggests that the estimates are already transformed.

• Are these data OK?
• If not,  shall I transform the data with  1 / (  1 + exp(-(z_Intercept_trans) )  )  and then do  z_Intercept_trans < 0.5?

If you could elaborate a bit more on this it would be great.

Thanks a lot for your help.

Best wishes,

Valerio

[Michael J Frank]

 Hi Valerio,

z_intercept_trans is indeed the transformed value so you would have to transform it back (in transformed space, you would need to have a value less than 0). But if you are using HDDM 0.8.0 you should not see any value less than 0, and so the result is biased if you used the inverse logit link function as below -  instead  I would rerun it with the identity link function as in the updated tutorial. 

Michael

To view this discussion on the web visit https://groups.google.com/d/msgid/hddm-users/262edf46-575c-4aba-9f5a-04132611f6ean%40googlegroups.com.

[valer...@gmail.com]

Hi Michael,

thank you so much for the clarification.

I was just wondering whether re-running the model is the only way to fix things or if I can manipulate those data somehow.

I'm asking because this model took almost three weeks to fit.

Thanks,

Valerio

[Michael J Frank]

That is a bummer. 3 weeks is a really long time for a single model fit, was this a huge dataset? 

There is no way to manipulate the data after the fact because it used the wrong prior for your link function and that would bias all the estimated intercepts on z to be 0.5 or higher (after transforming back to the 0-1 scale).  That said, depending on the model and the data you might not get results that are all that different,  if (i) the true intercept is relatively unbiased anyway (0.5), (ii)  most of the variance on z is determined by your regressors, or (iii) the fitted posteriors on z do not seem to be hitting the lower bound (0 in your transformed space)  

To view this discussion on the web visit https://groups.google.com/d/msgid/hddm-users/dcc6d18d-f571-4c57-9854-188c823c5f0en%40googlegroups.com.

[valer...@gmail.com]

Hi Michael,

I believe the best thing to do would be to re-run the model, then.

The dataset is not huge (30 participants who did 264 trials each), but the fitting time has always been in the range of weeks for me: I have to run long MCMCs (100K) to reach a decent convergence.

If you have any tricks up the sleeve to reduce the fitting time, they're more than welcome.

I have another question if you don't mind and I think it'd be useful to post it in this thread.

In my regression model, I also want to incorporate stimulus coding for z, thus I'm using the following link function:

def v_link_func( x, data = data ):

        stim = ( dmatrix( "0 + C(s, [[1], [-1]])",

                                     { 's': data.stim.loc[x.index] },

                                     return_type = 'dataframe' )

                     )

        

        return np.multiply(x.to_frame(), stim)

Is there anything I should amend here to match the new HDDM release/feature?

I am assuming this won't be impacted, but I just thought it'd be worth checking before running my model.

Thanks a lot, I look forward to hearing from you.

Best,

V

[Michael J Frank]

Hi Valerio, 

100K samples is a *lot* - sometimes might be nice if you want really smooth posteriors or to be more confident that you have converged but are you actually getting R-hats or other convergence metrics that tell you haven't converged after, say, 5K samples?

re:the link function for stimulus coding of z, the updated tutorial code has the appropriate syntax: 

def z_link_func(x, data=mydata):
stim = (np.asarray(dmatrix('0 + C(s, [[0], [1]])',
                           {'s': data.stimulus.loc[x.index]}))
)
# Apply z = (1 - x) to flip them along 0.5
z_flip = stim - x
# The above inverts those values we do not want to flip,
# so invert them back
z_flip[stim == 0] *= -1 

return z_flip

 

To view this discussion on the web visit https://groups.google.com/d/msgid/hddm-users/f6cd8bd9-55e3-42aa-a46d-c868afa9c670n%40googlegroups.com.

[Valerio Villani]

Hi Michael,

thanks a lot for sharing your expertise and the thorough explanation.

I tried chains of different lengths (1K, 10K, 25K, 50K) and ended up with 100K, eventually. I also fiddled with the burn-ins and thinning factors.

Now that I think about it, however, I believe that z was the parameter whose diagnostic plots were problematic!

Maybe using the outdated link functions might have determined this issue to some degree...

I will have to amend the code and re-fit the HDDM following the updated tutorial now. Maybe I can try with 5K chains and see how things go.

I will keep you posted and reply here if anything interesting comes up.

Best wishes,

Valerio

___________________________________

               

___________________________________

You received this message because you are subscribed to a topic in the Google Groups "hddm-users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/hddm-users/k8dUBepPyl8/unsubscribe.
To unsubscribe from this group and all its topics, send an email to hddm-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/hddm-users/CAGYaZ%2Bd3oUX%2BEaAmVeK-4pnm4ucrf_rWgN4_8DE1JsM7soTTPg%40mail.gmail.com.

[valer...@gmail.com]

Hi Michael,

I read the updated tutorial you linked me to, but I noticed it concerns an accuracy-coded model.

I'm not sure this would work for me and last time I might have formulated my question unclearly.

I'm using a stimulus-coded model with HDDMRegression, where I also estimate the decision bias z in different conditions.

I also want the drift rate v to vary by condition while taking into account the two stimuli/boundaries.

How should I go about that?

In my original model, I had used both the z_link_func and the v_link_func:

# DEFINE LINK FUNCTION FOR DECISION BIAS

def z_link_func( x, data = data ):

        stim = ( dmatrix( "0 + C(s, [[1], [-1]])",

                                     { 's': data.stim.loc[x.index] },

                                     return_type = 'dataframe' )

                     )

        

        return 1 / ( 1 + np.exp( -np.multiply(x.to_frame(), stim) ) )

# DEFINE LINK FUNCTION FOR DRIFT RATE

def v_link_func( x, data = data ):

        stim = ( dmatrix( "0 + C(s, [[1], [-1]])",

                                     { 's': data.stim.loc[x.index] },

                                     return_type = 'dataframe' )

                     )

        

        return np.multiply(x.to_frame(), stim)

Now I won't need the z_link_func above, but I don't understand whether I have to use:

• (a) the z_link_func in the updated tutorial for z and a lambda function for v;

OR

• (b) a lambda function for z and the v_link_func above;

OR

• (c) the z_link_func in the updated tutorial for z and the v_link_func above.

If you could point me in the right direction and suggest which tutorial I should focus on, it'd be greatly appreciated.

Thank you for your help and apologies for the long post.

Looking forward to your reply,

Valerio

[Michael J Frank]

Hi Valerio,

The stim vs accuracy coding distinction can sometimes be confusing and I think this might be the issue. 

You would only want to flip the “z” (from z to 1-z for stim1 vs stim2) when using accuracy coding.  eg if the subject has a bias to respond stim1, you estimate a bias toward stim1 by flipping z (so that z=.7 for stim1 and z=.3 for stim2 both represent a bias toward stim1 *if* the response column is accuracy coded).   

Since you are using stimcoding and the response columns are the responses themselves rather than correct/incorrect, then you don’t want to flip z at all, and you just estimate one z. The link function you would use is just lambda x:x and with no flipping, and you get a single estimate z of the bias toward stim1.  You still flip the v, as in your example (b) above. (When using HDDMStimCoding you do this with split_param = v but since you are doing the regression you can do it in the link function as you specified). 

 

I suspect that this might be the reason you were having issues with convergence (which is separate from anything to do with the inverse logit). If you flipped before it would mean that when stim1 is present then the bias is to stim1 but when stim2 is present the bias is to stim2, and that doesn’t make sense (ie the subject can’t know the stim before the trial to set the starting point), and instead it should only be the drift that flips. 

To view this discussion on the web visit https://groups.google.com/d/msgid/hddm-users/cfce2ae6-e741-4d60-9ea7-f4a0c382e836n%40googlegroups.com.

[valer...@gmail.com]

Hi Michael,

thanks for the clarification, I'll follow your advice and use the v_link_func.

I will also try to reduce the number of samples and see if the model converges.

Best regards,

Valerio