Interpreting the posterior probability of the parameter for a 2 by 2 within-subject design

[Chuan-Peng Hu]

Dear experts,

Many thanks in advance. I am new to HDDM and apologize if my questions are too naive. 

I am trying to apply HDDM to one of my studies, and not sure what I've were right or wrong. So, I will first describe what I have done and then list my questions. Maybe I did it complete wrong ( where I wrote in blue words), please correct me where I was wrong.

What I have done: 

My experiment has a 2 * 2 * 2 within-subject design, but we usually analyze data from two levels of the first variable separately because previous studies have shown it they are quick different. So the data I used to fit HDDM is a 2 (varA) * 2 (varB) within-subject design. I used HDDM in two steps: first model fitting and comparison; second, compare the parameter from different conditions from the best-fitting model.

As for the model comparison, I fitted my data 16 models in which different parameter(s) were set free to vary. For example, in the first model, v, a, z, and t were all set free to vary, while the next model I fixed one of the parameters; for the last one model, I fixed all the parameters. The code for fitting first mode is like this:

m1 = hddm.HDDM(datL_match2,depends_on = {'v':['varA','varB'],'z':['varA','varB'],'a':['varA','varB'],'t':['varA','varB']}, include=['v', 'z', 'a','t'],p_outlier=.05)

m1.find_starting_values()

m1.sample(20000,burn = 2000, dbname='traces.db', db='pickle')

The model comparison showed that the model with all four parameters free to vary was the best one (smallest DIC, PPC plot also looks good).

Then, I used two ways to do the hypothesis testing. One way was using the posterior probability as suggested by the manual: http://ski.clps.brown.edu/hddm_docs/howto.html#hypothesis-testing; another way was extracting the parameter statistics. I used the following code to extract the parameters:

stats_m = m1.gen_stats()

After I got the data of parameters, I carried a 2 by 2 Bayesian repeated measure ANOVA to the mean of each parameter, using JASP. For v, t, and a, the posterior probability and the Bayes Factor results from JASP were consistent (usually, if BF>3, the corresponding posterior probability is > 0.95), however, for z, the BF10 > 3 didn't corresponding to posterior probability > 0.95. Actually, varA1 > varB1 resulted in BF10 =  5.2, but posterior probability is 0.86. 

My Questions:

Q1: Can I set all v, a, t, z free to vary in a model? As I've found in the published papers, people rarely set all four parameters free. For example, in Frank, Gagne, Nyhus, Masters, Wiecki, Cavanagh, & Badre. (2015), only v and a were set free during the model comparison stage.

Q2: could I used the depends_on for my experimental design? I found that in many threads it was clear the "depends_on" is for between-subject design instead of within-subject design. So is it completely unacceptable to use "dependes_on" for a within-subject design?

Q3: If using HDDMRegressor is better, then,  How could do it right? because I found that I can't get the statistics for each parameter after using the HDDMRegressor, as I had reported in this thread:https://groups.google.com/forum/#!topic/hddm-users/taD9oirzx3M 

Q4: Can I use Bayesian repeated measure ANOVA to do the hypothesis testing? I knew that we CAN'T using frequentist ANOVA on the parameters. So can we use Bayesian ANOVA instead?

Q5: If I can use "depends_on" for my data and the Bayesian ANOVA can be used, how to interpret the BF10 > 5 while the posterior probability?  Are they conflict with each other? A related question is, is there a criterion for just a posterior probability is big enough (85%, 90% or 95%?).

Thank again.

Frank, Gagne, Nyhus, Masters, Wiecki, Cavanagh, & Badre. (2015). fMRI and EEG Predictors of Dynamic Decision Parameters during Human Reinforcement Learning. The Journal of Neuroscience, 35(2), 485-494. doi:10.1523/jneurosci.2036-14.2015

[Thomas Wiecki]

Hi,

Seems like you did a comprehensive analysis.

On Sat, Nov 26, 2016 at 9:40 AM, Chuan-Peng Hu <hcp...@gmail.com> wrote:

Q1: Can I set all v, a, t, z free to vary in a model? As I've found in the published papers, people rarely set all four parameters free. For example, in Frank, Gagne, Nyhus, Masters, Wiecki, Cavanagh, & Badre. (2015), only v and a were set free during the model comparison stage.

Yes, we've done that recently, no good reason it's not possible. 

Q2: could I used the depends_on for my experimental design? I found that in many threads it was clear the "depends_on" is for between-subject design instead of within-subject design. So is it completely unacceptable to use "dependes_on" for a within-subject design?

Yes, can't use it. 

Q3: If using HDDMRegressor is better, then,  How could do it right? because I found that I can't get the statistics for each parameter after using the HDDMRegressor, as I had reported in this thread:https://groups.google.com/forum/#!topic/hddm-users/taD9oirzx3M 

Has that issue been resolved? 

Q4: Can I use Bayesian repeated measure ANOVA to do the hypothesis testing? I knew that we CAN'T using frequentist ANOVA on the parameters. So can we use Bayesian ANOVA instead?

You can just do 1-to-1 comparisons of your conditions (i.e. their posteriors) to get the same thing. That's what Kruschke proposes at least. I think in his book he has a chapter on Bayesian ANOVA.

Q5: If I can use "depends_on" for my data and the Bayesian ANOVA can be used, how to interpret the BF10 > 5 while the posterior probability?  Are they conflict with each other? A related question is, is there a criterion for just a posterior probability is big enough (85%, 90% or 95%?).

Not sure what you mean with BF10 > 5?

Usually the criterion is 97.5% (or 2.5%) for two-sided and 95% for one-sided. Kruschke also argues for a ROPE: http://doingbayesiandataanalysis.blogspot.de/2013/08/how-much-of-bayesian-posterior.html

[Chuan-Peng Hu]

Hi, Thomas, 

Thank you a lot for your reply!! Really appreciate that.

On Thursday, December 1, 2016 at 4:40:02 PM UTC+8, Thomas wrote:

Hi,

Seems like you did a comprehensive analysis.

On Sat, Nov 26, 2016 at 9:40 AM, Chuan-Peng Hu <hcp...@gmail.com> wrote:

Q1: Can I set all v, a, t, z free to vary in a model? As I've found in the published papers, people rarely set all four parameters free. For example, in Frank, Gagne, Nyhus, Masters, Wiecki, Cavanagh, & Badre. (2015), only v and a were set free during the model comparison stage.

Yes, we've done that recently, no good reason it's not possible. 

So glad that ok to do that!

 

Q2: could I used the depends_on for my experimental design? I found that in many threads it was clear the "depends_on" is for between-subject design instead of within-subject design. So is it completely unacceptable to use "dependes_on" for a within-subject design?

Yes, can't use it. 

Thank you a lot! At least it I know that what I've done was wrong ;-).

 

Q3: If using HDDMRegressor is better, then,  How could do it right? because I found that I can't get the statistics for each parameter after using the HDDMRegressor, as I had reported in this thread:https://groups.google.com/forum/#!topic/hddm-users/taD9oirzx3M 

Has that issue been resolved? 

Partially, I think. But I think I need to do more homework before asking you guys. 

 

Q4: Can I use Bayesian repeated measure ANOVA to do the hypothesis testing? I knew that we CAN'T using frequentist ANOVA on the parameters. So can we use Bayesian ANOVA instead?

You can just do 1-to-1 comparisons of your conditions (i.e. their posteriors) to get the same thing. That's what Kruschke proposes at least. I think in his book he has a chapter on Bayesian ANOVA.

Q5: If I can use "depends_on" for my data and the Bayesian ANOVA can be used, how to interpret the BF10 > 5 while the posterior probability?  Are they conflict with each other? A related question is, is there a criterion for just a posterior probability is big enough (85%, 90% or 95%?).

Not sure what you mean with BF10 > 5?

Sorry for the ambiguity of my question. What I meant is that what if the results of Bayesian ANOVA and the posterior probability seems conflict, for example, using Bayesian t-test for the mean of v under condition A and condition B results in a BF10 > 5, while the posterior probability of v under condition A > condition B is only about 86%. But now, I am not sure whether this is a problem because I am not sure whether there is equivalent between Bayesian factor and posterior probability  

 

Usually the criterion is 97.5% (or 2.5%) for two-sided and 95% for one-sided. Kruschke also argues for a ROPE: http://doingbayesiandataanalysis.blogspot.de/2013/08/how-much-of-bayesian-posterior.html

This link is great! 

Thanks again!

[Micah Johnson]

Hi Thomas and team,

Thank you for making and maintaining such a brilliant tool. 

I've been looking through forum posts and published material but I've been unable to answer some questions, would you mind helping me out? I apologize if I have overlooked something simple. 

I have a completely within-subjects design with factor 1 (Test, 2 levels: base, task), factor 2 (Site, 4 levels: site1, site2, site3, site4), and factor 3 (Timepoint, 4 levels: tpt1, tpt2, tpt3, tpt4). Factors 1 and 2 are fully nested within each other (i.e., base and task have all site levels, and each site has both base and task), but factor 3 of Timepoint codes for the order of experimental sessions which is different for each subject (i.e., for subj1, tpt1 is site1, tpt2 is site4, etc; for subj2, tpt1 is site3, tpt2 is site1, etc). So I've set up the following patsy model and HDDMRegressor model...

dmatrix("C(Test,Treatment('base'))*C(Site, Treatment('site1'))+C(Timepoint, Treatment('tpt1'))", data)

model = hddm.HDDMRegressor(data, {"v ~ C(Test, Treatment('base'))*C(Site, Treatment('site'))+C(Timepoint, Treatment('tpt1'))"}, keep_regressor_trace=True, p_outlier=.08, is_group_model=True)

This then gives me the following columns:

''Intercept''

"C(Test, Treatment('base'))[T.t]"

"C(Site, Treatment('site'))[T.site2]"

"C(Site, Treatment('site'))[T.site3]"

"C(Site, Treatment('site'))[T.site4]"

"C(Timepoint, Treatment('tpt1'))[T.tpt2]"

"C(Timepoint, Treatment('tpt1'))[T.tpt3]"

"C(Timepoint, Treatment('tpt1'))[T.tpt4]"

"C(Test, Treatment('base'))[T.t]:C(Site, Treatment('site'))[T.site2]"

"C(Test, Treatment('base'))[T.t]:C(Site, Treatment('site'))[T.site3]"

"C(Test, Treatment('base'))[T.t]:C(Site, Treatment('site'))[T.site4]"

So here my questions:

1). How do I interpret the different contrast vectors? Do I interpret all of them, even the interaction terms, relative to the single "Intercept" term (v_Intercept)? For example with the simple means reported in the print_stats option (although I know I will be doing Bayesian hypothesis testing on the full posteriors, but just for easy examples)...

If my v_Intercept = 1.24,

- when mean of v_C(Test, Treatment('base'))[T.task] = 0.20, then does that mean that the "task" condition = 1.44? If so, then how would I look at the "base" condition, would it be like this...v_C(Test, Treatment('base'))[T.base]?? (I'm suspicious of that, since not reported in stats output). Or is it actually the case that when v_C(Test, Treatment('base'))[T.task] = 0.20, this means that the "task" condition has drift rate that is 0.20 higher than the "base" condition drift rate? If so, how do I then figure out what is the mean (or full posterior distribution) of the "base" condition?

- similarly, for interaction terms, when the mean drift rate of C(Test, Treatment('base'))[T.t]:C(Site, Treatment('site'))[T.site2] = -0.35, is that relative to the v_Intercept term? If so, do I then interpret it as follows...the mean drift rate for "task" condition of "site1" = (1.24 + -0.35 = 0.89)?? But if it's not relative to v_Intercept, then I have no idea how to interpret this complex interaction term, since it would seem that the mean drift rate of the interaction term would be relative to two different intercepts, the intercept for the first contrast vector, C(Test, Treatment('base')) which would be "base" condition, and then the intercept for the second contrast vector, C(Site, Treatment('site')) which would be "site1". 

2). When I add the Timepoint factor with its contrast vectors (e.g. v_C(Timepoint, Treatment('tpt1'))) as a "covariate" to my regression model, is it correct to interpret this inclusion as analogous to the "controlling for" aspect of standard regression? What I mean is, can I then say, the drift rates for the factor1xfactor2 interaction terms (e.g., C(Test, Treatment('base'))[T.t]:C(Site, Treatment('site'))[T.site2]) are meaningfully interpreted as "independent of" or "controlled for" the effects of timepoints? (much like an ANCOVA, where the effect of the Timepoint covariate needs to be removed)

3). Trying to understand that complex model in relation to an alternative simpler model with a single factor, e.g. Condition, with 8 levels (site1base, site1task, site2base...site4base, site4task), coded with patsy and hddm as: 

dmatrix("C(Condition,Treatment('site1base'))", data) 

model = hddm.HDDMRegressor(data, {"v ~ C(Condition, Treatment('site1base'))"}, keep_regressor_trace=True, p_outlier=.08, is_group_model=True)

Technically, for each individual subject, seems to me that the effect of timepoint would be inherent in this Condition coding, because I know the condition order for each subject...however, since the condition order is different across subjects yet my HDDM model analyses each Condition at a group level, I would lose the ability to interpret these Condition drift rates as "controlling for" the effects of timepoints, is that correct? In which case, these different models (model1: Test x Site + Timepoints...model2: Condition) would not produce equivalent drift rate posteriors for each basic condition (e.g., the "task" condition of "site1"), is that correct? So for my situation, do you recommend I stick with model1 since I want to "remove" or "control for" the effect of time?

My apologies for a lengthy post, I'm just very excited to figure this out so I can proceed. If you don't have much time to respond, my 1) and 2) questions are most important to me. 

Thank you very much in advance,

Micah

[Chuan-Peng Hu]

Hi, Micah,

It seems that we stuck in a similar situation.

I guess that the real difficulty we faced is how to use patsy to get the effect we are interested in. in your case, is the 2 by 4 design, with control of the time point, in my case is the 2 by 2 within-design. Personally, I get used to the ANOVA way of thinking, so a little confused about how to used the regression model to get similar results as in ANOVA. I guess the documents about patsy, which I am reading,would be useful: http://patsy.readthedocs.io/en/latest/API-reference.html#handling-categorical-data .

I tried this:

dmatrix("C(id, Treatment('self'))*C(moral,Treatment('moral'))", datL.head(10))  # id is my first independent variable (IV), with two levels: self v. other; moral is my second IV, with two levels: moral v. immoral.

it returns:

DesignMatrix with shape (10L, 4L)

  Columns:

    ['Intercept',

     "C(id, Treatment('self'))[T.other]",

     "C(moral, Treatment('moral'))[T.immoral]",

     "C(id, Treatment('self'))[T.other]:C(moral, Treatment('moral'))[T.immoral]"]

  Terms:

    'Intercept' (column 0)

    "C(id, Treatment('self'))" (column 1)

    "C(moral, Treatment('moral'))" (column 2)

    "C(id, Treatment('self')):C(moral, Treatment('moral'))" (column 3)

  (to view full data, use np.asarray(this_obj))

My interpretation of these three columns is as following:

"C(id, Treatment('self'))[T.other]"  ---   The relative effect of "other" to "self" in the "id" variable, so it could be considered as the main effect of "id";

"C(moral, Treatment('moral'))[T.immoral]"  --- the main effect of "moral"

"C(id, Treatment('self'))[T.other]:C(moral, Treatment('moral'))[T.immoral]"  --- the interaction between these two variables.

I am not sure whether I am correct or not.

As for you data, because the time point is different for each participant, I guess it is a between-subject variable, so maybe a more complex design matrix is needed.

[Chuan-Peng Hu]

Dear Experts, 

After the above discussion and searching in this group ( this post is very informative), I could run the within model without errors, but I still have not sure about the interpretation of the results.

I used this code to run the model (I have two within independent variables, id: self vs. other; moral: moral vs immoral; I treated both "self" and "moral" as intercept):

m_within_subj_L1_1 = hddm.HDDMRegressor(datL, {"v ~ C(id, Treatment('self'))*C(moral,Treatment('moral'))","t ~ C(id, Treatment('self'))*C(moral,Treatment('moral'))","a ~ C(id, Treatment('self'))*C(moral,Treatment('moral'))","z ~ C(id, Treatment('self'))*C(moral,Treatment('moral'))"}, include = ['t','a', 'v','z'],keep_regressor_trace=True, p_outlier = 0.05)

m_within_subj_L1_1.sample(10000,burn = 1000,dbname='traces_within1_1.db', db='pickle')

And this code to get the parameter v:

v_Intercept,v_id,v_moral,v_id_moral = m_within_subj_L1_1.nodes_db.ix[['v_Intercept', 

                                                                 "v_C(id, Treatment('self'))[T.other]", 

                                                                 "v_C(moral, Treatment('moral'))[T.immoral]", 

                                                                 "v_C(id, Treatment('self'))[T.other]:C(moral, Treatment('moral'))[T.immoral]"],'node']

I plot these parameters, the plot is as following:

Below is my understanding of the design matrix and the interpretation of each parameter output, But I am not sure whether this is correct. So please give me some feedback if you know.

I think that the design matrix is a dummy coding way of 2 by 2 within-subject design, so it's was as:

                         moral       id

moral-self            0            0

immoral-self        1            0

moral-other          0           1

immoral-other      1           1

Therefore, 

v_Intercept -- moral self;  

v_C(id, Treatment('self'))[T.other] --- (moral other + immoral other) - (moral self + immoral self)     # main effect of id

v_C(moral, Treatment('moral'))[T.immoral]  -- (immoral self + immoral other) - (moral self + moral other)    # main effect of moral 

v_C(id, Treatment('self'))[T.other]:C(moral, Treatment('moral'))[T.immoral]  -- (moral-other -moral-self) - (immoral other - immoral self)    # interaction between id and moral

Then, to examine the each effect, I compare each effect with 0 (either smaller or greater than 0) to infer whether the effect exist:

print "P(v_self > v_other ) = ", (v_id.trace() < 0).mean()

print "P(v_moral > v_immoral) = ", (v_moral.trace() < 0).mean()

print "P(v_interaction) = ", (v_id_moral.trace() > 0).mean()

Many thanks!

[nitzan shahar]

Chung,

I'm no expert in hddm, but I'm also trying to figure out similar things. If you are trying to figure out the two main effects and the two way interaction  - I think you can code your factors with 0 and 1, then use something like {'v ~ id + self + self*id'}, in HDDMregressor for each parameter. If you do that, I would guess that the intercept will reflect changes for the subjects factor, id and self for the main effects, and self*id for the interaction. 

I have to say that this might get you exactly the same result, but at least for me it's easier to describe it this way.

Keep update on your analysis so we could all learn.

BW,

Nitzan

[Thomas Wiecki]

Yes, C() creates a dummy-coding for you, so it should be equivalent to what Nitzan suggests. The effects are then also relative to the intercept. Your coding is quite complex and I don't have the full context of your study but generally this all makes sense.

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

[Chuan-Peng Hu]

Hi, Nitzan, Thomas,

Thank you all for the reply! 

PS: I will update my analysis, Nitzan.

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

[Nitzan Shahar]

OK - so following this topic, I wanted to please add two questions (apologies for any cross-posting) : 

1. Given the above 2x2 within design, if we ran HDDM.regressor = {'v ~ FactorA + FactorB + FactorA*B}, participants's individual scores will give the intercept per participant. If I understand correctly, these per subject intercept scores reflect differences in drift regardless of the manipulations (right?). Is there a way we can calculate participants scores per condition in this design?

2. Previous posts discussed the fact that depends_on is for between. Yet, if I have a 2x2 within design, and I set one factor as regressor (within) and the other as between (e.g., HDDM.regressor = {'v ~ FactorA', depends_on={'v': ['FactorB]}) - what will be the harm in that? I understand that I will be missing the interaction and that statistical power will be probably reduced, as the design is regarding the levels of FactorB as sampled from different groups. Other than that - is there a strict reason why this should be avoided?

Thank you,

Nitzan

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

[Chuan-Peng Hu]

Hi, Nitzan,

Sorry that I can provide an answer but add my own situation, hope that more experienced people to provide some cue.

On Monday, December 12, 2016 at 11:08:13 PM UTC+8, nitzan shahar wrote:

OK - so following this topic, I wanted to please add two questions (apologies for any cross-posting) : 

1. Given the above 2x2 within design, if we ran HDDM.regressor = {'v ~ FactorA + FactorB + FactorA*B}, participants's individual scores will give the intercept per participant. If I understand correctly, these per subject intercept scores reflect differences in drift regardless of the manipulations (right?). Is there a way we can calculate participants scores per condition in this design?

I have the same question ! It seems that HDDM.regressor only provided the effect ( intercept for each participant), but not the values of the parameter (e.g., v) for each participant. I also wondering whether there is a way to get those score just as using "depends_on".

 

There was a PDF file about the equivalence between 2 by 2 AVOVA and regression, which told how to get the value for each condition by manual, but I am not sure it applied to HDDM.regressor.

 

2. Previous posts discussed the fact that depends_on is for between. Yet, if I have a 2x2 within design, and I set one factor as regressor (within) and the other as between (e.g., HDDM.regressor = {'v ~ FactorA', depends_on={'v': ['FactorB]}) - what will be the harm in that? I understand that I will be missing the interaction and that statistical power will be probably reduced, as the design is regarding the levels of FactorB as sampled from different groups. Other than that - is there a strict reason why this should be avoided?

There was a post about mixed effect design (both between and within) https://groups.google.com/forum/#!searchin/hddm-users/Mixed$20design|sort:relevance/hddm-users/E3VgnJCk-vI/meCWssPNqkIJ .

I also hope to know the answer of this interesting question.

[Annik]

Hello all,

I too am stuck on the same problem - this discussion has been extremely useful! I have a 2x2 within subj design (I also have a 2x3 design but thought I'd tackle the simpler one first).

Has there been any resolution to the questions raised above? In particular, I would like to isolate the drift-rates for each of the four conditions (i.e. combinations of factor levels). For a one way anova I see the logic, but  at the moment I'm finding the patsy coded variables difficult to intepret for a 2x2?

Any help or update on your progress much appreciated.