Contrasts in hierarchical models

[Guido Biele]

Hi,

I was wondering about the correct way to calculate posterior distributions of contrasts in a hierarchical model. Here are the two ways I can think of:

Lets assume I have 2 conditions, A1 and AB and want to get the posterior distribution for A1-A2. Then I could

a)    Take the nodes for the group means of A1 and A2 and subtract them. The posterior distribution is then the distribution of this differences calculated for all mcmc samples.

b)   Subtract A1-A2 for each individual on the level of participant-nodes, and calculate the average of this division. The posterior distribution is then the distribution of this average of differences calculated for al mcmc samples.

I tried both, and b) typically shows better results in that the variance of the posterior distribution is smaller. However, I am unsure if this actually is a permissible way.

The reason why I think that b) is actually the correct way is that I have a design where condition is manipulated within participants, and I think that if I just take the group nodes to estimate contrasts, the variance I see comes from between and within participant differences. Put differently, I assume that the variation in the group node for a parameter depends also on between-participant variation in that parameter. However, if I have a contrast for a within-participants design, I am not interested in the between participant variance of a particular parameter, but only in the between participant variance for the difference between 2 parameters.

What do you think about this issue?

Cheers – Guido

[Michael J Frank]

Hi Guido, 

We were just talking about this issue - there is not currently great support for within-subject effects in hddm. I.e one can estimate different drift rates per individual, and obtain group distributions for each of these on their own. If you subtract the group distributions you can obtain a posterior that allows you to make inferences about whether these drifts differ at the population level, but this is not as powerful as a within-subject contrast as you imply.

One way to handle it currently is to use a regression model rather than estimating separate drifts (or other params): 

Set v = beta1 +beta2 * [condition]

if condition is categorical and you just have two conditions, then the group distribution of beta2 regression parameter is equivalent to the group level inference on within subject effects. You could probably also set up different contrast codings for multiple conditions and different betas, as you would in a glm in frequentist stats.

Thomas and Imri can chime in on this and other implementational issues

Michael

Michael J Frank, PhD, Associate 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+...@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.
 
 

[guido biele]

thanks for the quick reply!
I had thought of hddm-regression, but we are already using the stim-coding version of hddm, in which one cannot define a regression model (as far as I can tell).

After a bit more thinking, taking the samples from the subject-nodes seems generally sound to me. maybe taking the averages for each sample is not such a good idea, but one could simply generate the posterior from all samples (I.e. when my chain has 10000 samples and I have 10 participants, I would have a posterior distribution generated from 100000 difference values calculated on the subject level).
I think so because if I set up two models:
v = beta1 * [condition1 | condition2] + beta2 * [condition2] or
v = beta1b * [condition] + beta2b * [condition2]

then the distribution of beta2 should be identical to the distribution of beta2b-beta1b.

actually, this is something I could test with hddm regression for a data set where stim-coding is not needed.
I'll do this as soon as I can and report back.

cheers-guido

--
Sent from my Android phone with K-9 Mail. Please excuse my brevity.

[Michael J Frank]

Hi Guido,  If I understand your idea properly, it is to use a model that estimates separate drift rates each coming from a group level distribution, but then take each sample of v1 and v2 for each participant, compute v1-v2 at each sample, and then plot the histogram of these for all samples and participants to construct a group distribution of v1-v2.

I don't think that this solution will produce a true group-level posterior on the within subject effects. If you or others on this list disagree with the below or have other suggestions, please reply!

First, the group level posterior has a variance that conveys the uncertainty about the mean value of the within-subject effects across all subjects. But the variance of the histogram you would create with the method above would not reflect this uncertainty in the same way, but rather the degree to which there is variance across samples,  mixing different subjects in the analysis rather than having the values of subjects sampled from a group prior. One way to make this more clear is to think about what would happen if you had different trial counts per subject.  Subjects with more trials should have more peaked posteriors whereas those with less trials would have more spread out posteriors, and would thus inflate the variance of the constructed group posterior more than they should. (We are interested in the variance of the mean values governing the within subject effects, not the variance across samples of each subject). 

Second, the same method would imply that you could just do the same thing for individual v's in a standard hddm analysis (e.g. even with a single drift rate). If you  take all the samples of individual v's and make one big histogram of those across all subjects, it will not be the same histogram as that of the group level. 

The within-subject structure should really be included in the model itself. One of the celebrated effects of hierarchical models is that they provide statistical strength so that if subjects are relatively similar to each other on some parameter, we can have a relatively strong prior to make inferences about the parameter value of another subject with less data. So if all subjects have a within-subject difference in v1 vs v2, then a subject with less data probably also does. If we have a group distribution that specifies the effect of condition on drift rate, as in the regression model, then it will have this property. But if we only have group distributions separately for v1 and v2, it will not. To see this more clearly, consider a case where overall drift rate is highly variable across subjects,  but where every single subject has v1 - v2 = 0.1.  

 

Hopefully the regression model should work though, and can be recoded to apply for stimulus coding.. 

Michael

[Guido Biele]

Hi,

I agree with my Michael that my conflation of sampling variance and between subject variance was mistaken. Thanks for the clarification.

Fortunately, I think there is an easy way incorporate stimulus-coding into the regression model. The idea is simply to use sign of a stimulus regressor to code the stimulus type, and to adjust the priors and the definition of z accordingly. More specifically:

for v:
Design matrix: the stimulus column should have a 1 for one stimulus and a -1 one for the other stimulus. (and there should be no intercept)
Prior: The density of the prior should be non-zero only for values larger zero (e.g. create_family_trunc_normal('v', lower=0, upper=1e3, value=0) ?)

for z:
to make this work I think one needs to slightly reformulate z under the hood, such that the z used in the likelihood function is calculated as z = .5+z*. The regression model would then fit z*.
Design matrix: same as for v
Prior: The density of the prior should be non-zero only for 0 <= z* <.5 ( (e.g. create_family_trunc_normal('v', lower=0, upper=.5, value=0) ?)


What do you think?

Cheers - Guido

(401)-863-6872 <tel:%28401%29-863-6872>

On Wed, Mar 6, 2013 at 8:40 AM, Guido Biele
<g.p....@psykologi.uio.no

<mailto:hddm-users%2Bunsu...@googlegroups.com>.

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




--
Sent from my Android phone with K-9 Mail. Please excuse my brevity.

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

<mailto:hddm-users%2Bunsu...@googlegroups.com>.

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



--
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/groups/opt_out.

--

---

 Guido Biele 
 Email: g.p....@psykologi.uio.no 
 Phone: +47 228 45172 
 Website 


 Visiting Address 
 Psykologisk Institutt  
 Forskningsveien 3 A 
 0373 OSLO 



  
  
  
 

  Mailing Address 
  Psykologisk Institutt 
  Postboks 1094 
  Blindern 0317 OSLO

[Thomas Wiecki]

I think this is a very interesting direction to take. If this allows for within-subject effects and other tricks like stimulus coding many problems that would otherwise require a hand-tuned version would be solved.

I had an idea that we could include patsy to specify the design matrix http://patsy.readthedocs.org/en/latest/quickstart.html. This would then allow for better syntax along the lines of R in specifying a GLM. E.g.:

"y ~ x1 + x2"

So please share your experiences if you explore this further.

[Thomas Wiecki]

OK, I took a stab at this. In our github repo there is a new branch called 'patsy' which requires the kabuki develop branch. This greatly simplifies regression model specification while making it pretty easy to specify very complex models.

You can now specify a model with, for example, a trial-by-trial regression on theta, within-subject categorical effects for stimulus (WW, LL, WL) and interactions effects as follows:

# Define the regression descriptor where we specify
reg_descr = {'func': "theta*C(stim)",
             'outcome': 'a'}
m_reg = hddm.HDDMRegressor(data, 
                           reg_descr, 
                           depends_on={'v': 'stim'}) 
# if you do not want subj dists for the regressors, you can set group_only_regressors=True

Which tells which covariates are being created:

Adding these covariates:
['Intercept', 'C(stim)[T.WL]', 'C(stim)[T.WW]', 'theta', 'theta:C(stim)[T.WL]', 'theta:C(stim)[T.WW]']

Intercept will be threshold during the 'LL' condition. C(stim)[T.WL] will be the offset during the 'WL' condition. Theta is a trial-by-trial measure for which the regression coefficient will be estimated. The other are interaction terms.

For more information on how the syntax work I refer to the patsy docs:

http://patsy.readthedocs.org/en/latest/quickstart.html

And these slides:

http://jseabold.net/presentations/seabold_pydata2012.html#slide11

Let me know whether you get any experiences with this.

Thomas

[Guido Biele]

Hi Thomas,

I tried to install hddm patsy, but the cdfdif_wrapper is making a problem, because the c source-code seems to be missing.

here is the relevant section of log:

"

building 'cdfdif_wrapper' extension

gcc -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -mavx -fPIC -I/cluster/software/VERSIONS/python2-packages-2.7_2/lib/python2.7/site-packages/numpy/core/include -I/cluster/software/VERSIONS/python2-2.7.3/include/python2.7 -c src/cdfdif_wrapper.c -o build/temp.linux-x86_64-2.7/src/cdfdif_wrapper.o

gcc: src/cdfdif_wrapper.c: No such file or directory

gcc: no input files

error: command 'gcc' failed with exit status 1

"

I'm not sure what the problem is. the cdfdif_wrapper is also in the master branch, but it seems that it is not compiled when installing hddm, or I am missing something ...

cheers - Guido

 

<mailto:hddm-users%2Bunsubscribe@googlegroups.com>.

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




--
Sent from my Android phone with K-9 Mail. Please excuse my brevity.

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

<mailto:hddm-users%2Bunsubscribe@googlegroups.com>.

[Thomas Wiecki]

Hi Guido,

Indeed. I added the files. Let me know if it still fails.

Thomas

[Guido Biele]

Hi Thomas,

thanks, it worked now.

cheers - guido

[Thomas Wiecki]

I did some minor tweaks like adding support for multiple regressors and updated the documentation. I think it's pretty close to final so you might want to update and let me know if you agree.

[guido biele]

great!
I'll look at it tonight. one thing I stumbled over earlier today is your sanity check for the depends on. should it be that a parameter or a regressor-variable should not be used in the regression and the depends_on?
relatedly, it seems that the check unnecessarily rejects a model when only a regression but no depends_on is specified.
cheers-guido

[Thomas Wiecki]

Right, that was a bug I also fixed.

But it is true that you can not put any regressors into depends_on. However, you do not have to as you can do everything with constrasts within the regression. Let me know if you have more questions on this.

[Guido Biele]

Hi Thomas,

This looks quite complete, and I got some basic models to work, I hope also with 2 regressors, e.g.:

reg_descr = [{'func': '1 + drug + block', 'outcome': 'a'},{'func': '1 + drug + block', 'outcome': 'v'}]. I am not sure though that this worked as desired, as I get very weird parameters (e.g. an intercept for v of -15, whereas I get a much more reasonable estimates of v if I use only one regression function for v) I’ll try again with a longer chain tomorrow.

Here are 2 things I couldn’t figure out:

(1) include=’z’

when I tried to run it on the parameter z (after using “include=”z”” in the model setup), this did not work. Instead, I get following error message:

(2) stim-coding

When I want to run the model with “stim_coding” over 2 stimulus-types for the parameter v, I would like to multiply the result of func with a 1 for one stimulus type and with -1 for the other. Were the 1 and -1 supposed to be coded in the theta you were mentioning in an earlier email?
( reg_descr = {'func': "theta*C(stim)",'outcome': 'a'} )
or would I need to approach this differently, for instance using the I(…) syntax of patsy to get arithmetic transformations?

For me the following seemed to work, where drug and block are categorical variables and stim codes one stimulus as 1 and the other as -1:

reg_descr = [{'func': '1 + I(drug*stim) + I(block*stim)', 'outcome': 'v'}]

Cheers - Guido

PS: i would seem to me that it wold be useful to use a bit different priors for stim-coding of v (and z, if this would work). however, I also see that the code can get quite complex when this is also implemented...

[Guido Biele]

Hi,

I made a few more checks, which lead me to revise my earlier email.

First, what remains is that HDDMregressor does not seem to work for z. I tried it with and without the “include = (‘z’)” specification, but neither worked.

What I wrote about stim-coding was, however, not so helpful I think.
Here is what I am thinking now about the stim-coding issue:
(As a reminder: I suggest to implement stim-coding for e.g. v by (a) multiplying all rows of the design matrix for one stimulus with a 1 and for the other stimulus with a -1, and (b) specifying a prior for the intercept to be larger than 0)

I think this cannot easily be solved within the patsy-formulas. Rather,  one could patsy let make a design matrix without stim-coding info, and then add this to the patsy-generated  design matrix.
As far as I can tell, this could be implemented in lines 72+ of hddm_regression.py.

Priors for regression-parameters for stim-coding should probably be adjusted in line 196+ of hddm_regression.py.
I would suggest a gamma distribution as the prior for v the intercept of v and a normal distribution for the other v-regressors.

Unfortunately the case is more complicated for z, because z is constrained to [0 1], and it is a product of multiple regressors. To give an example, if the intercept for z has the value .6, then the parameter for the effect of one level of a  condition should be constrained to [-.6 .4].  I’ll try to find out if pymc has a function to implement such a linear constraint.

Cheers - Guido

[Thomas Wiecki]

Hi Guido,

Thanks for the follow-up.

Stim-Coding:

I took a closer look at patsy and I agree that this special coding is not supported out of the box. However, one can create new contrast codings. This should do the trick:

import numpy as np
from patsy import ContrastMatrix
class StimCoding(object):
    def __init__(self, reference=0):
        self.reference = reference

    def code_with_intercept(self, levels):
        assert len(levels) == 2, "Only 2 levels are supported for stimulus coding."
        return ContrastMatrix(np.array([[1, 0], [0, -1]]),
                              ["[My.%s]" % (level,) for level in levels])

    def code_without_intercept(self, levels):
        raise NotImplementedError, "StimCoding can only be used without intercept. Add 0 + C()"

from patsy import dmatrix
data = ['stim1', 'stim1', 'stim2', 'stim2']
dmatrix("0 + C(stims, StimCoding)", {'stims': data})

I think this implements your idea for the stimulus coding for v. Did you have a design matrix idea for z?

However, I can not reproduce your missing include='z' bug. I added a unittest to the current patsy branch -- see if that fails for you.

Thomas

[Thomas Wiecki]

Also, for constraining z. We before used an invlogit transform. I think this can also be done. Patsy also seems to support this. You can for example use 'np.log(theta)' to log transform theta. Not sure yet if this actually helps us though.

[Guido Biele]

Hi Thomas,

I have to thank you!
I'll look at this tonight rather than writing something now that i feel stupid having written in 3 hours...

about stim-coding for z:
If we have to stimulus types stim1 and stim2, then i would use the regression model to estimate a parameter z*, which is then used to calculate
z(stim1) = z* and
z(stim2) = 1-z*

operationally one would again need a design matrix in which one stimulus has a positive sign and the other has a negative sign, and in which z* ( i.e. intercept*beta0 + theta1*theta2 ... ) is limited to [.5 1]. I am not sure how the logit transform helps to enforce this constraint. I'll look into this tonight too.

cheers - guido

Thomas

First, what remains is that HDDMregressor does not seem to work for z. I tried it with and without the �include = (�z�)� specification, but neither worked.

What I wrote about stim-coding was, however, not so helpful I think.
Here is what I am thinking now about the stim-coding issue:
(As a reminder: I suggest to implement stim-coding for e.g. v by (a) multiplying all rows of the design matrix for one stimulus with a 1 and for the other stimulus with a -1, and (b) specifying a prior for the intercept to be larger than 0)

I think this cannot easily be solved within the patsy-formulas. Rather, �one could patsy let make a design matrix without stim-coding info, and then add this to the patsy-generated �design matrix.

As far as I can tell, this could be implemented in lines 72+ of hddm_regression.py.

Priors for regression-parameters for stim-coding should probably be adjusted in line 196+ of hddm_regression.py.
I would suggest a gamma distribution as the prior for v the intercept of v and a normal distribution for the other v-regressors.

Unfortunately the case is more complicated for z, because z is constrained to [0 1], and it is a product of multiple regressors. To give an example, if the intercept for z has the value .6, then the parameter for the effect of one level of a� condition should be constrained to [-.6 .4].� I�ll try to find out if pymc has a function to implement such a linear constraint.

Cheers - Guido

On Monday, March 11, 2013 11:15:53 PM UTC+1, Guido Biele wrote:

Hi Thomas,

This looks quite complete, and I got some basic models to work, I hope also with 2 regressors, e.g.:

reg_descr = [{'func': '1 + drug + block', 'outcome': 'a'},{'func': '1 + drug + block', 'outcome': 'v'}]. I am not sure though that this worked as desired, as I get very weird parameters (e.g. an intercept for v of -15, whereas I get a much more reasonable estimates of v if I use only one regression function for v) I�ll try again with a longer chain tomorrow.

Here are 2 things I couldn�t figure out:

(1) include=�z�

when I tried to run it on the parameter z (after using �include=�z�� in the model setup), this did not work. Instead, I get following error message:

(2) stim-coding

When I want to run the model with �stim_coding� over 2 stimulus-types for the parameter v, I would like to multiply the result of func with a 1 for one stimulus type and with -1 for the other. Were the 1 and -1 supposed to be coded in the theta you were mentioning in an earlier email?

( reg_descr = {'func': "theta*C(stim)",'outcome': 'a'} )

or would I need to approach this differently, for instance using the I(�) syntax of patsy to get arithmetic transformations?

For me the following seemed to work, where drug and block are categorical variables and stim codes one stimulus as 1 and the other as -1:

reg_descr = [{'func': '1 + I(drug*stim) + I(block*stim)', 'outcome': 'v'}]

Cheers - Guido

PS: i would seem to me that it wold be useful to use a bit different priors for stim-coding of v (and z, if this would work). however, I also see that the code can get quite complex when this is also implemented...

On Monday, March 11, 2013 6:49:54 PM UTC+1, Thomas wrote:

I did some minor tweaks like adding support for multiple regressors and updated the documentation. I think it's pretty close to final so you might want to update and let me know if you agree.

On Mon, Mar 11, 2013 at 8:10 AM, Guido Biele <g.p....@psykologi.uio.no> wrote:

Hi Thomas,

thanks, it worked now.

cheers - guido

On Monday, March 11, 2013 1:03:03 PM UTC+1, Thomas wrote:

Hi Guido,

Indeed. I added the files. Let me know if it still fails.

Thomas

On Mon, Mar 11, 2013 at 5:47 AM, Guido Biele <g.p....@psykologi.uio.no> wrote:

Hi Thomas,

I tried to install hddm patsy, but the�cdfdif_wrapper is making a problem, because the c source-code seems to be missing.

here is the relevant section of log:

"

building 'cdfdif_wrapper' extension

gcc -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -mavx -fPIC -I/cluster/software/VERSIONS/python2-packages-2.7_2/lib/python2.7/site-packages/numpy/core/include -I/cluster/software/VERSIONS/python2-2.7.3/include/python2.7 -c src/cdfdif_wrapper.c -o build/temp.linux-x86_64-2.7/src/cdfdif_wrapper.o

gcc: src/cdfdif_wrapper.c: No such file or directory

gcc: no input files

error: command 'gcc' failed with exit status 1

"

I'm not sure what the problem is. the�cdfdif_wrapper is also in the master branch, but it seems that it is not compiled when installing hddm, or I am missing something ...

cheers - Guido

�

On Saturday, March 9, 2013 11:15:09 PM UTC+1, Thomas wrote:

OK, I took a stab at this. In our github repo there is a new branch called 'patsy' which requires the kabuki develop branch. This greatly simplifies regression model specification while making it pretty easy to specify very complex models.

You can now specify a model with, for example, a trial-by-trial regression on theta, within-subject categorical effects for stimulus (WW, LL, WL) and interactions effects as follows:

# Define the regression descriptor where we specify
reg_descr = {'func': "theta*C(stim)",
             'outcome': 'a'}
m_reg = hddm.HDDMRegressor(data, 
                           reg_descr, 
                           depends_on={'v': 'stim'}) 
# if you do not want subj dists for the regressors, you can set group_only_regressors=True

Which tells which covariates are being created:

Adding these covariates:
['Intercept', 'C(stim)[T.WL]', 'C(stim)[T.WW]', 'theta', 'theta:C(stim)[T.WL]', 'theta:C(stim)[T.WW]']

Intercept will be threshold during the 'LL' condition. C(stim)[T.WL] will be the offset during the 'WL' condition. Theta is a trial-by-trial measure for which the regression coefficient will be estimated. The other are interaction terms.

For more information on how the syntax work I refer to the patsy docs:

http://patsy.readthedocs.org/en/latest/quickstart.html

And these slides:

http://jseabold.net/presentations/seabold_pydata2012.html#slide11

Let me know whether you get any experiences with this.

Thomas

On Thu, Mar 7, 2013 at 6:57 PM, Thomas Wiecki <thomas...@gmail.com> wrote:

I think this is a very interesting direction to take. If this allows for within-subject effects and other tricks like stimulus coding many problems that would otherwise require a hand-tuned version would be solved.

I had an idea that we could include patsy to specify the design matrix�http://patsy.readthedocs.org/en/latest/quickstart.html. This would then allow for better syntax along the lines of R in specifying a GLM. E.g.:

Cheers � Guido

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

<mailto:hddm-users%2Bunsu...@googlegroups.com>.

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




--
Sent from my Android phone with K-9 Mail. Please excuse my brevity.

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

<mailto:hddm-users%2Bunsu...@googlegroups.com>.

�
�

--
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/groups/opt_out.

�
�

--
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/groups/opt_out.

�
�

--
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/groups/opt_out.

�
�

--
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/groups/opt_out.

�
�

[Thomas Wiecki]

Hi Guido,

Unfortunately I realized just now that what I proposed is not what we need.

The design matrix should look something like [1, 1, -1, -1] rather than [[1, 1, 0, 0], [0, 0, -1, -1]] which is what the code does I proposed. Not sure yet how this can be done with patsy but maybe I'll ask their mailing list.

Thomas

On Tue, Mar 12, 2013 at 12:51 PM, Guido Biele <g.p....@psykologi.uio.no> wrote:

Hi Thomas,

Thomas

First, what remains is that HDDMregressor does not seem to work for z. I tried it with and without the “include = (‘z’)” specification, but neither worked.

What I wrote about stim-coding was, however, not so helpful I think.
Here is what I am thinking now about the stim-coding issue:
(As a reminder: I suggest to implement stim-coding for e.g. v by (a) multiplying all rows of the design matrix for one stimulus with a 1 and for the other stimulus with a -1, and (b) specifying a prior for the intercept to be larger than 0)

I think this cannot easily be solved within the patsy-formulas. Rather,  one could patsy let make a design matrix without stim-coding info, and then add this to the patsy-generated  design matrix.

As far as I can tell, this could be implemented in lines 72+ of hddm_regression.py.

Priors for regression-parameters for stim-coding should probably be adjusted in line 196+ of hddm_regression.py.
I would suggest a gamma distribution as the prior for v the intercept of v and a normal distribution for the other v-regressors.

Unfortunately the case is more complicated for z, because z is constrained to [0 1], and it is a product of multiple regressors. To give an example, if the intercept for z has the value .6, then the parameter for the effect of one level of a  condition should be constrained to [-.6 .4].  I’ll try to find out if pymc has a function to implement such a linear constraint.

Cheers - Guido

On Monday, March 11, 2013 11:15:53 PM UTC+1, Guido Biele wrote:

Hi Thomas,

This looks quite complete, and I got some basic models to work, I hope also with 2 regressors, e.g.:

reg_descr = [{'func': '1 + drug + block', 'outcome': 'a'},{'func': '1 + drug + block', 'outcome': 'v'}]. I am not sure though that this worked as desired, as I get very weird parameters (e.g. an intercept for v of -15, whereas I get a much more reasonable estimates of v if I use only one regression function for v) I’ll try again with a longer chain tomorrow.

Here are 2 things I couldn’t figure out:

(1) include=’z’

when I tried to run it on the parameter z (after using “include=”z”” in the model setup), this did not work. Instead, I get following error message:

(2) stim-coding

When I want to run the model with “stim_coding” over 2 stimulus-types for the parameter v, I would like to multiply the result of func with a 1 for one stimulus type and with -1 for the other. Were the 1 and -1 supposed to be coded in the theta you were mentioning in an earlier email?

( reg_descr = {'func': "theta*C(stim)",'outcome': 'a'} )

or would I need to approach this differently, for instance using the I(…) syntax of patsy to get arithmetic transformations?

For me the following seemed to work, where drug and block are categorical variables and stim codes one stimulus as 1 and the other as -1:

reg_descr = [{'func': '1 + I(drug*stim) + I(block*stim)', 'outcome': 'v'}]

Cheers - Guido

PS: i would seem to me that it wold be useful to use a bit different priors for stim-coding of v (and z, if this would work). however, I also see that the code can get quite complex when this is also implemented...

On Monday, March 11, 2013 6:49:54 PM UTC+1, Thomas wrote:

I did some minor tweaks like adding support for multiple regressors and updated the documentation. I think it's pretty close to final so you might want to update and let me know if you agree.

On Mon, Mar 11, 2013 at 8:10 AM, Guido Biele <g.p....@psykologi.uio.no> wrote:

Hi Thomas,

thanks, it worked now.

cheers - guido

On Monday, March 11, 2013 1:03:03 PM UTC+1, Thomas wrote:

Hi Guido,

Indeed. I added the files. Let me know if it still fails.

Thomas

On Mon, Mar 11, 2013 at 5:47 AM, Guido Biele <g.p....@psykologi.uio.no> wrote:

Hi Thomas,

I tried to install hddm patsy, but the cdfdif_wrapper is making a problem, because the c source-code seems to be missing.

here is the relevant section of log:

"

building 'cdfdif_wrapper' extension

gcc -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -mavx -fPIC -I/cluster/software/VERSIONS/python2-packages-2.7_2/lib/python2.7/site-packages/numpy/core/include -I/cluster/software/VERSIONS/python2-2.7.3/include/python2.7 -c src/cdfdif_wrapper.c -o build/temp.linux-x86_64-2.7/src/cdfdif_wrapper.o

gcc: src/cdfdif_wrapper.c: No such file or directory

gcc: no input files

error: command 'gcc' failed with exit status 1

"

I'm not sure what the problem is. the cdfdif_wrapper is also in the master branch, but it seems that it is not compiled when installing hddm, or I am missing something ...

cheers - Guido

 

On Saturday, March 9, 2013 11:15:09 PM UTC+1, Thomas wrote:

OK, I took a stab at this. In our github repo there is a new branch called 'patsy' which requires the kabuki develop branch. This greatly simplifies regression model specification while making it pretty easy to specify very complex models.

You can now specify a model with, for example, a trial-by-trial regression on theta, within-subject categorical effects for stimulus (WW, LL, WL) and interactions effects as follows:

# Define the regression descriptor where we specify
reg_descr = {'func': "theta*C(stim)",
             'outcome': 'a'}
m_reg = hddm.HDDMRegressor(data, 
                           reg_descr, 
                           depends_on={'v': 'stim'}) 
# if you do not want subj dists for the regressors, you can set group_only_regressors=True

Which tells which covariates are being created:

Adding these covariates:
['Intercept', 'C(stim)[T.WL]', 'C(stim)[T.WW]', 'theta', 'theta:C(stim)[T.WL]', 'theta:C(stim)[T.WW]']

Intercept will be threshold during the 'LL' condition. C(stim)[T.WL] will be the offset during the 'WL' condition. Theta is a trial-by-trial measure for which the regression coefficient will be estimated. The other are interaction terms.

For more information on how the syntax work I refer to the patsy docs:

http://patsy.readthedocs.org/en/latest/quickstart.html

And these slides:

http://jseabold.net/presentations/seabold_pydata2012.html#slide11

Let me know whether you get any experiences with this.

Thomas

On Thu, Mar 7, 2013 at 6:57 PM, Thomas Wiecki <thomas...@gmail.com> wrote:

I think this is a very interesting direction to take. If this allows for within-subject effects and other tricks like stimulus coding many problems that would otherwise require a hand-tuned version would be solved.

I had an idea that we could include patsy to specify the design matrix http://patsy.readthedocs.org/en/latest/quickstart.html. This would then allow for better syntax along the lines of R in specifying a GLM. E.g.:

Cheers – Guido

hddm-users+...@googlegroups.com

[Guido Biele]

Hi Thomas,

nice that the patsy method works!!

just to make sure that we are on the same page, here an example for a design matrix with stimulus coding for a design with 2 stimulus types and an orthogonal condition with 3 levels:

data = {"s": ["s1", "s1", "s1", "s2", "s2", "s2"],"c":["c1", "c2", "c3", "c1", "c2", "c3"]}

patsy.dmatrix("0 + C(s, [[1], [-1]]) + C(s, [[1], [-1]]):C(c)", data)

which makes the design matrix

array([[ 1.,  0.,  0.],

           [ 1.,  1.,  0.],

           [ 1.,  0.,  1.],

           [-1., -0., -0.],

           [-1., -1., -0.],

           [-1., -0., -1.]])

where column one gives us the effect for condition c1, column two gives us the difference between c1 and c2, and column 3 the difference between c1 and c3.

The nice thing is that this works with current version of HDDMregressor (i.e. the modification you described in your first email today is not required).

I think I still don't understand how the inverse logit as  a prior helps to implement a linear contraint. if we take the design matrix above to model z, then one could do this by defining z as

z = .5 + X*beta (in matrix notation, where beta is a column-vector with parameters and X is the design matrix.)

the contraint would then be -.5<=X*beta<=.5. What I specifically do not understand if how using inverse logit priors for the three paramerters in beta here helps to implement the constraint. but it might well be that I am just not aware of the relevant literature.

cheers - guido

Thomas

<mailto:hddm-users%2Bunsubscribe@googlegroups.com>.

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




--
Sent from my Android phone with K-9 Mail. Please excuse my brevity.

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

<mailto:hddm-users%2Bunsubscribe@googlegroups.com>.

...

[Thomas Wiecki]

On Tue, Mar 12, 2013 at 4:36 PM, Guido Biele <g.p....@psykologi.uio.no> wrote:

Hi Thomas,

nice that the patsy method works!!

just to make sure that we are on the same page, here an example for a design matrix with stimulus coding for a design with 2 stimulus types and an orthogonal condition with 3 levels:

data = {"s": ["s1", "s1", "s1", "s2", "s2", "s2"],"c":["c1", "c2", "c3", "c1", "c2", "c3"]}

patsy.dmatrix("0 + C(s, [[1], [-1]]) + C(s, [[1], [-1]]):C(c)", data)

which makes the design matrix

array([[ 1.,  0.,  0.],

           [ 1.,  1.,  0.],

           [ 1.,  0.,  1.],

           [-1., -0., -0.],

           [-1., -1., -0.],

           [-1., -0., -1.]])

where column one gives us the effect for condition c1, column two gives us the difference between c1 and c2, and column 3 the difference between c1 and c3.

The nice thing is that this works with current version of HDDMregressor (i.e. the modification you described in your first email today is not required).

Great! Seems like patsy is the way to go then. I realized that this can actually also replace the depends_on functionality since you can just do: '0 + C(depends_column)' and get the same thing. Just that now you can just remove the 0 and get within subject estimation.

 

I think I still don't understand how the inverse logit as  a prior helps to implement a linear contraint. if we take the design matrix above to model z, then one could do this by defining z as

z = .5 + X*beta (in matrix notation, where beta is a column-vector with parameters and X is the design matrix.)

the contraint would then be -.5<=X*beta<=.5. What I specifically do not understand if how using inverse logit priors for the three paramerters in beta here helps to implement the constraint. but it might well be that I am just not aware of the relevant literature.

You're right, transforming the regressors will not help here. We would have to logit transform whatever the output of the linear model specified by patsy (X*beta, not just beta). We have done this in HDDMTransformed so it's certainly possible. I think something similar could be done for HDDMRegressor as well. These transformations are a little weird but I don't really see another way.

Thomas

[Guido Biele]

Hi Thomas,

I agree that the patsy approach is very versatile. I would still not dismiss the depends_on version, because is seems to be much faster when I implement the same design (the depends_on was 4.2 times faster in my quick test).

About inclusion of z and bias: I looked at your unittest (I actually didn't know if/how i can run it with a simple command...). The difference between your test and what i was trying to do is that I tried to have the regression to model z, and you included z while using the regression to model v. I can also do that, but as far as i can tell the regression can at the moment not be used to model z.

I now see how you use the inverse logit to transform the regression-output into a-parameter. Sorry for not picking this up earlier, I think this makes sense, even if this makes the interpretation of regression parameters for z a bit more difficult.

cheers - guido

...

[Guido Biele]

a quick addition on how to model z with stim-coding:

here is how I think this could be done:

(a) a variable "pre_z" is modeled, analog to how v is modeled. The prior for regression parameters are uninformative normals centered on 0
(b) z is then defined as z = .5 + .5*inv_logit(pre_z)

I hope this makes sense,

cheers - guido

PS: I can write a short documentation for the within-subjects stim-coding issue for v and z if you like.

...

[Thomas Wiecki]

On Wednesday, March 13, 2013 10:31:24 AM UTC+1, Guido Biele wrote:

Hi Thomas,

I agree that the patsy approach is very versatile. I would still not dismiss the depends_on version, because is seems to be much faster when I implement the same design (the depends_on was 4.2 times faster in my quick test).

That is interesting. It might be a result of the likelihood function for the regression model that does some additional checks.

But yeah, since we already have the tried-and-tested depends_on (and it's faster) there'd have to be some really good reasons to throw it out the window.

About inclusion of z and bias: I looked at your unittest (I actually didn't know if/how i can run it with a simple command...). The difference between your test and what i was trying to do is that I tried to have the regression to model z, and you included z while using the regression to model v. I can also do that, but as far as i can tell the regression can at the moment not be used to model z.

In case you haven't figured it out: nosetest tests/test_models.py:TestHDDMRegression

I'll change the test to do a regression on z. 

I now see how you use the inverse logit to transform the regression-output into a-parameter. Sorry for not picking this up earlier, I think this makes sense, even if this makes the interpretation of regression parameters for z a bit more difficult.

Right. Maybe a prior_descr could be added along the lines of:

prior_descr = {'z_bottom': {'transform': lambda x: invlogit(x)}, 
               'z_Intercept': {'prior_func': pm.Normal, 'prior_params': {'mu': .5, 'tau':1} }

Not sure yet if this can be implemented but it might actually allow to do what Roger Ratcliff has been requesting -- a way to set certain parameters to a fixed value (e.g. st = .1).

PS: I can write a short documentation for the within-subjects stim-coding issue for v and z if you like.

That'd be well appreciated. Maybe first some notes on the designmatrix approach and then some common examples like within-subject effects and trial-by-trial regressions (which are in the current docs but are outdated once we change to patsy). Ideally you could just edit the howto.rst directly and then do a pull request on github. If you just want to write the text that's fine too of course.

Thomas 

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

[Thomas Wiecki]

FYI: I updated the model specification to make it even simpler. Instead of

HDDMRegressor(data, {'func': 'C(cond)', 'outcome': 'v'}})

you now do:

HDDMRegressor(data, 'v ~ C(cond)')

[Guido Biele]

I think the new way to specify a model is nice and quicker to parse for a the human reader ;-)!

When reading through the patsy documentation to prepare a how-to, i found that they a number of nice "built in" function to prepare design matrices for particular contrasts of interest.

I think this is very attractive, as it provides a fast way to set up design matrices in which the regressor coefficients can directly be read as contrats.

For example, lets assume you have a design with condition drug with 3 ordered levels (e.g. antagonist, placebo, agonist). Then

hddm.HDDMregressor(data,'v ~ C(drug,Diff)',data)

would set up a design matrix where the regressors model 

(1) the reference-level (e.g. drug1_antagonist) as the intercept, 

(2) a regressor capturing the contrast drug2_placebo > drug1_antagonist and 

(3)  a regressor capturing the contrast drug3_agonist > drug3_placebo.*

The problem is, that the stimulus-coding scheme we discussed earlier is not easily combined with using these built in designs (i.e. one can do either the one or the other).

However, one could have both if one does not let the user  implement the stim-coding through patsy, but rather adds an additional variable to the input (e..g "stim_coding = False" as default, which the the user can change to "stim_coding = True".

When stim_coding is set to True, these few lines would add  stim-coding to any existing design matrix X that was previously defined with patsy (assuming the stimulus-column is named stim)

def stimcode(X,data):

    ncols = X.shape[1]

    stimdata = {'s':data['stim']}

    stim = (np.asarray(dmatrix('0 + C(s,[[1],[-1]])',stimdata)))

    mstim = np.kron(stim,np.ones((ncols,1))).reshape(-1,ncols)

    return X*mstim

About the regression for z: The new tests for include = 'z' together with the regression for z does not work for me.

Cheers - Guido

*Similarily hddm.HDDMregressor(data,'v ~ C(drug,Diff)',data) implements dummy coding (interecept for the reference level and 2 regressors for drug2 > reference level and drug3 > reference level, respectively) and hddm.HDDMregressor(data,'v ~ C(drug,Sum)') has an intercept for the grand mean and 2 regressors that compare 2 levels against the grant mean.

...

[Thomas Wiecki]

On Thu, Mar 14, 2013 at 4:16 PM, Guido Biele <g.p....@psykologi.uio.no> wrote:

I think the new way to specify a model is nice and quicker to parse for a the human reader ;-)!

Great, I agree! 

When reading through the patsy documentation to prepare a how-to, i found that they a number of nice "built in" function to prepare design matrices for particular contrasts of interest.

I think this is very attractive, as it provides a fast way to set up design matrices in which the regressor coefficients can directly be read as contrats.

For example, lets assume you have a design with condition drug with 3 ordered levels (e.g. antagonist, placebo, agonist). Then

hddm.HDDMregressor(data,'v ~ C(drug,Diff)',data)

would set up a design matrix where the regressors model 

(1) the reference-level (e.g. drug1_antagonist) as the intercept, 

(2) a regressor capturing the contrast drug2_placebo > drug1_antagonist and 

(3)  a regressor capturing the contrast drug3_agonist > drug3_placebo.*

Very cool. 

The problem is, that the stimulus-coding scheme we discussed earlier is not easily combined with using these built in designs (i.e. one can do either the one or the other).

However, one could have both if one does not let the user  implement the stim-coding through patsy, but rather adds an additional variable to the input (e..g "stim_coding = False" as default, which the the user can change to "stim_coding = True".

When stim_coding is set to True, these few lines would add  stim-coding to any existing design matrix X that was previously defined with patsy (assuming the stimulus-column is named stim)

def stimcode(X,data):

    ncols = X.shape[1]

    stimdata = {'s':data['stim']}

    stim = (np.asarray(dmatrix('0 + C(s,[[1],[-1]])',stimdata)))

    mstim = np.kron(stim,np.ones((ncols,1))).reshape(-1,ncols)

    return X*mstim

Interesting, thanks for providing that. Would be quite easy to add it seems. Changing to patsy actually simplified the code a bit so adding a little bit back in might not be too bad.

 

About the regression for z: The new tests for include = 'z' together with the regression for z does not work for me.

Yeah, they are failing for me as well now which allows me to debug the problem which I'll still have to do. Probably won't get to it next week but the week after.

Thomas 

--

[Thomas Wiecki]

Guido:

I think I fixed the z include bug.

I also added the option to supply a link-function in case the parameter is to be transformed, you can now do:

hddm.HDDMRegressor(data, [('z ~ cov', lambda x: 1 / (1 + np.exp(-x)))], include='z')

In words, if you supply a tuple for the model, the second item can be any function that will be used to transform the output.

Let me know of any issues with this or if you can think of a better interface.

Thomas

[guido biele]

thanks Thomas,
I'll try this beginning next week!
guido

[Thomas Wiecki]

OK, in terms of interface, maybe this is more explicit:

{'model': 'z ~ cov', 'link_func': lambda x: ...}

If there is no link func you can just pass the model string without the dict.

[Guido Biele]

Hi Thomas,

I tried to let the lates version of HDDMregressor run and found a small bug, which I could correct myself. the code in lines 151+ of hddm_regression.py does not get the model and the link function from the list/tuple.

I think the following works:

for m in models:

assert m.has_key('model'), 'HDDMRegressor requires a model specification like ''model'': ''v ~ 1 + C(your_variable)''' 

model = m['model']

if m.has_key('link_func'):

link_func = m['link_func']

else:

link_func = lambda x: x

Still, after I have corrected teh code as above, I get the same 

"After 7.000000 retries, still no good fit found."

error message Dunovan describe a few days ago. I tried to remebdy this by installing the patsy version (instead of the development version which has the most recent hddm_regression.py file) and replacing the hddm_regression.py of the patsy version with the one form the development version.

However, that didn't solve the problem as this didn't work with hierachial.py from kabuki.

I guess you can fix that much faster than I could ...

Cheers - Guido

...

[Thomas Wiecki]

Not sure why it didn't work. In any case, I changed the interface so that you can either pass in a string, or a dict that looks like this:

link_func = lambda x: 1 / (1 + np.exp(-x))
m = hddm.HDDMRegressor(data, {'model': 'z ~ cov', 'link_func': link_func}, include='z')

Can you see if that still doesn't work? If so, please send a reproducing example so that I can debug.

Thomas

--

[Guido Biele]

thanks!

that worked. I think the problem was that i used a "stupid" test lambda function (1/x) which generated inadmissible z values.

my hopefully last question for a while ;-)!

when I want to use a lambda function to do the stim-coding*, it would need to access a column of the data file. is there a way to do this?

of that is not possible, I have an idea of how to do this at another place in the code. but I'd prefer the lambda function to keep everything in one place

Cheers - Guido

* e.g., the lambda function for stimulus coding of z would look like this:

link_func = lambda x: .5 + 2 / (1 + np.exp(-(x*stim)))

where I assume that x and stim are arrays.

...

[Thomas Wiecki]

That's tricky, but it might work.

The input to the lambda function (x) should be a dataframe with an index. This index allows you to locate the appropriate data. E.g.

def link_func(x, data=my_data_frame):
  stims = data.stim.ix[x.index]
  return .5 + 2 / (1 + np.exp(-(x*stim)))

One trick here is that data will be bound to the function when you instantiate it so you don't have to pass it but can access it later, similar to a class. You then select the appropriate rows once the function gets called.

Good luck ;),

Thomas

[Guido Biele]

thanks,
I'll try that, and then I'll also use the right function for
stim-coding of z, which should be more like

(1 / (1 + np.exp(-x))) * stim1 + (1 - 1/(1 + np.exp(-x)))*stim2

where stim1, stim2 are indicator variable that are 1 for one variable
and zero for the other respectively.
anyhow, I'll implement this as a lambda or otherwise.

cheers - guido

On Fri Mar 22 14:50:51 2013, Thomas Wiecki wrote:
> That's tricky, but it might work.
>
> The input to the lambda function (x) should be a dataframe with an
> index. This index allows you to locate the appropriate data. E.g.
>
> |def link_func(x, data=my_data_frame):
> stims = data.stim.ix[x.index]

> return .5 +2 / (1 + np.exp(-(x*stim)))|

> One trick here is that data will be bound to the function when you
> instantiate it so you don't have to pass it but can access it later,
> similar to a class. You then select the appropriate rows once the
> function gets called.
>
> Good luck ;),
> Thomas
>
> On Fri, Mar 22, 2013 at 6:30 AM, Guido Biele
> <g.p....@psykologi.uio.no <mailto:g.p....@psykologi.uio.no>> wrote:
>
> thanks!
> that worked. I think the problem was that i used a "stupid" test
> lambda function (1/x) which generated inadmissible z values.
>
> my hopefully last question for a while ;-)!
>
> when I want to use a lambda function to do the stim-coding*, it
> would need to access a column of the data file. is there a way to
> do this?
> of that is not possible, I have an idea of how to do this at
> another place in the code. but I'd prefer the lambda function to
> keep everything in one place
>
> Cheers - Guido
>
>
> * e.g., the lambda function for stimulus coding of z would look
> like this:

> link_func = lambdax: .5+ 2/ (1+ np.exp(-(x*stim)))

>
> where I assume that x and stim are arrays.
>
> On Friday, March 22, 2013 2:10:46 PM UTC+1, Thomas wrote:
>
> Not sure why it didn't work. In any case, I changed the
> interface so that you can either pass in a string, or a dict
> that looks like this:
>

> |link_func =lambda x:1 / (1 + np.exp(-x))

> ncols = X.shape[__1]

> stimdata = {'s':data['stim']}
> stim = (np.asarray(dmatrix('0 +
> C(s,[[1],[-1]])',stimdata)))
> mstim =

> np.kron(stim,np.ones((ncols,1)____)).reshape(-1,ncols)

> tests/test_models.py:__TestHDDMR____egression

> from patsyimport ContrastMatrix

> class StimCoding(object):
> def __init__(self, reference=0):
> self.reference = reference
>
> def code_with_intercept(self, levels):

> assert len(levels) ==2,"Only 2 levels are supported for stimulus coding."
> return ContrastMatrix(np.array([[1,0], [0, -1]]),
> ["[My.%s]" % (level,)for levelin levels])

>
> def code_without_intercept(self, levels):
> raise NotImplementedError,"StimCoding can only be used without intercept. Add 0 + C()"
>

> from patsyimport dmatrix

> -I/cluster/software/VERSIONS/__p________ython2-packages-2.7_2/lib/__pyth________on2.7/site-packages/numpy/__core________/include
> -I/cluster/software/VERSIONS/__p________ython2-2.7.3/include/python2.__7
> -c
> src/cdfdif_wrapper.c
> -o
> build/temp.linux-x86_64-2.7/__sr________c/cdfdif_wrapper.o

> http://patsy.readthedocs.org/__e________n/latest/quickstart.html

> <http://patsy.readthedocs.org/en/latest/quickstart.html>
>
> And
> these
> slides:

> http://jseabold.net/__presentati________ons/seabold___pydata2012.html#__sl______ide11

> http://patsy.__readthedoc________s.org/en/latest/__quickstart.__htm______l
> <http://patsy.readthedocs.org/en/latest/quickstart.html>.

> create_family_trunc_normal('v'__________,lower=0,upper=1e3,value=0)

> create_family_trunc_normal('v'__________,lower=0,upper=.5,value=0)

> <mailto:hddm-users%2Bunsu...@googlegroups.com>.

> For more options, visit https://groups.google.com/groups/opt_out.
>
>
>
> --
> 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/groups/opt_out.
>
>

--
------------------------------------------------------------------------

Guido Biele
Email: g.p....@psykologi.uio.no
Phone: +47 228 45172

Website <https://sites.google.com/a/neuro-cognition.org/guido/home>

[Guido Biele]

Hi Thomas,

I have implemented the method you suggested and it works.

for the benefit of others who want to test similar things, here what is I did:

I created test-data (with hddm.generate) of an experiment with 2 stimulus types and 3 conditions, in which each condition influences drift rate and bias.

Because i want interpretable v and z parameters, I set the hddm analysis up with accuracy coded data.

z and v were modeled with an intercept regressor (modeling condition 1) + 2 dummy regressors for conditions 2 and 3, respectively. In addition I used a link function for the parameter z, which insures that the effective z is "z" for stimulus type A and "1-z" for stimulus type B.

I implemented the model in the development branch. the one thing I had to adjust is line 77 in hddm_regrssion.py.

the original code is

predictor = link_func((design_matrix * params).sum(axis=1))

which I changed to 

predictor = link_func(pd.DataFrame((design_matrix * params).sum(axis=1), index=data.index)

because when implementing the stimulus coding via a link function the predictor has to be a pandas data frame.

For the experiment described above , the link function for stimulus-coding of z is 

def z_link_func(x, data=mydata):

    stim = (np.asarray(dmatrix('0 + C(s,[[1],[-1]])',{'s':data.stim.ix[x.index]})))    

    return 1 / (1 + np.exp(-(x * stim)))

the regression models for z and v are 

z_reg = {'model': 'z ~ 1 + C(condition)', 'link_func': z_link_func}

v_reg = {'model': 'v ~ 1 + C(condition)', 'link_func': lambda x : x}

which one can further use to set up the model

reg_descr = [z_reg, v_reg]

m_reg = hddm.HDDMRegressor( mydata, reg_descr,include='z')

when looking at the results (after only 15000 samples for a data set of 9000 trials) the mean z values were closer to the true values than the mean v values.

[keep link functions and regression model in mind when calculating v or z from the regression results ...]

one thing I was wondering about when setting this up if the way the link function is set up now means that the function dmatrix will be called for each mcmc sample. if so, I think a better way would be to use dmatrix only once outside the link function and to submit the array with the stim coding directly to the link_function (if this is possible, I'll try it later).

thanks again to thomas for making this possible.

I'll write up a short documentation on using HDDMRegressor with patsy (with and without stimcoding ...) in the next days.

cheers-guido

...

[Thomas Wiecki]

On Tue, Mar 26, 2013 at 2:48 PM, Guido Biele <g.p....@psykologi.uio.no> wrote:

Hi Thomas,

I have implemented the method you suggested and it works.

for the benefit of others who want to test similar things, here what is I did:

I created test-data (with hddm.generate) of an experiment with 2 stimulus types and 3 conditions, in which each condition influences drift rate and bias.

Because i want interpretable v and z parameters, I set the hddm analysis up with accuracy coded data.

z and v were modeled with an intercept regressor (modeling condition 1) + 2 dummy regressors for conditions 2 and 3, respectively. In addition I used a link function for the parameter z, which insures that the effective z is "z" for stimulus type A and "1-z" for stimulus type B.

I implemented the model in the development branch. the one thing I had to adjust is line 77 in hddm_regrssion.py.

the original code is

predictor = link_func((design_matrix * params).sum(axis=1))

which I changed to 

predictor = link_func(pd.DataFrame((design_matrix * params).sum(axis=1), index=data.index)

because when implementing the stimulus coding via a link function the predictor has to be a pandas data frame.

Thanks, I made that change and committed. 

For the experiment described above , the link function for stimulus-coding of z is 

def z_link_func(x, data=mydata):

    stim = (np.asarray(dmatrix('0 + C(s,[[1],[-1]])',{'s':data.stim.ix[x.index]})))    

    return 1 / (1 + np.exp(-(x * stim)))

the regression models for z and v are 

z_reg = {'model': 'z ~ 1 + C(condition)', 'link_func': z_link_func}

v_reg = {'model': 'v ~ 1 + C(condition)', 'link_func': lambda x : x}

which one can further use to set up the model

reg_descr = [z_reg, v_reg]

m_reg = hddm.HDDMRegressor( mydata, reg_descr,include='z')

when looking at the results (after only 15000 samples for a data set of 9000 trials) the mean z values were closer to the true values than the mean v values.

[keep link functions and regression model in mind when calculating v or z from the regression results ...]

That's quite an elaborate model so it's fantastic that it actually worked. If you have the code we could include that as an example. Do you want to comment and do a PR on that?

one thing I was wondering about when setting this up if the way the link function is set up now means that the function dmatrix will be called for each mcmc sample. if so, I think a better way would be to use dmatrix only once outside the link function and to submit the array with the stim coding directly to the link_function (if this is possible, I'll try it later).

Fortunately not. The reason is that here:

 def func(args, design_matrix=dmatrix(reg['model'], data=data), link_func=reg['link_func']):

we are computing the dmatrix once when the function is defined, not when it gets called. This is a neat way to make variables persistent and accessible to functions. I always cursed about this as beginners often stumble over this when creating mutable default kwargs but I came to appreciate it ;).

thanks again to thomas for making this possible.

I'll write up a short documentation on using HDDMRegressor with patsy (with and without stimcoding ...) in the next days.

Cool, here is the github issue we created for that:

https://github.com/hddm-devs/hddm/issues/26

Do you have an account?

Thomas 

--