In [1]:

    

import pandas as pd
import matplotlib.pyplot as plt
from patsy import dmatrix
import hddm
print hddm.__version__


    

    




0.5.4

In [2]:

    

data = hddm.load_csv('/Users/doby/Dropbox/current_projects/TBS_to_PFC/Analyses/modeling/data_for_HDDM.csv')


    

In [3]:

    

data.head(10)


    

    
Out[3]:






  

    

      

      
subj_idx
      
tms_site
      
Attended
      
SAT
      
rt
      
response
    
  
  

    

      
0
      
 1
      
 S1
      
 Invalid
      
    Speed
      
 1.151184
      
 1
    
    

      
1
      
 1
      
 S1
      
 Invalid
      
    Speed
      
 1.348829
      
 0
    
    

      
2
      
 1
      
 S1
      
 Invalid
      
    Speed
      
 0.834477
      
 1
    
    

      
3
      
 1
      
 S1
      
 Invalid
      
 Accuracy
      
 0.853608
      
 0
    
    

      
4
      
 1
      
 S1
      
   Valid
      
 Accuracy
      
 0.762255
      
 0
    
    

      
5
      
 1
      
 S1
      
   Valid
      
 Accuracy
      
 0.545990
      
 1
    
    

      
6
      
 1
      
 S1
      
   Valid
      
 Accuracy
      
 0.885909
      
 1
    
    

      
7
      
 1
      
 S1
      
 Invalid
      
    Speed
      
 0.659734
      
 0
    
    

      
8
      
 1
      
 S1
      
 Invalid
      
    Speed
      
 0.437045
      
 0
    
    

      
9
      
 1
      
 S1
      
 Invalid
      
 Accuracy
      
 1.191697
      
 0
    
  

In [4]:

    

data = hddm.utils.flip_errors(data)


    

In [5]:

    

# Instantiate model object passing it our data (no need to call flip_errors() before passing it).
# This will tailor an individual hierarchical DDM around your dataset.
m = hddm.HDDM(data, depends_on={'a': ['tms_site', 'Attended', 'SAT'], ...
                                'v': ['tms_site', 'Attended', 'SAT']}, p_outlier=.05)
# find a good starting point which helps with the convergence.
m.find_starting_values()

# start drawing 7000 samples and discarding 5000 as burn-in
m.sample(2000, burn=20)

#save the model
model.save('no within-subject effects model')


    

    




 [-----------------100%-----------------] 2000 of 2000 complete in 819.0 sec





    
Out[5]:





<pymc.MCMC.MCMC at 0x10cc86a50>

In [40]:

    

m.nodes_db.map


    

    
Out[40]:





a(Invalid.Accuracy.DLPFC)    1.815564
a(Invalid.Accuracy.FEF)      1.832361
a(Invalid.Accuracy.S1)       1.824281
a(Invalid.Accuracy.aPFC)     1.798484
a(Invalid.Speed.DLPFC)       1.630211
a(Invalid.Speed.FEF)         1.657380
a(Invalid.Speed.S1)          1.660465
a(Invalid.Speed.aPFC)        1.642362
a(Valid.Accuracy.DLPFC)      1.829784
a(Valid.Accuracy.FEF)        1.914307
a(Valid.Accuracy.S1)         1.850517
a(Valid.Accuracy.aPFC)       1.851659
a(Valid.Speed.DLPFC)         1.650483
a(Valid.Speed.FEF)           1.682566
a(Valid.Speed.S1)            1.687847
...
wfpt(Valid.Speed.FEF).14.0     NaN
wfpt(Valid.Speed.S1).14.0      NaN
wfpt(Valid.Speed.aPFC).14.0    NaN
wfpt(Valid.Speed.DLPFC).15.0   NaN
wfpt(Valid.Speed.FEF).15.0     NaN
wfpt(Valid.Speed.S1).15.0      NaN
wfpt(Valid.Speed.aPFC).15.0    NaN
wfpt(Valid.Speed.DLPFC).16.0   NaN
wfpt(Valid.Speed.FEF).16.0     NaN
wfpt(Valid.Speed.S1).16.0      NaN
wfpt(Valid.Speed.aPFC).16.0    NaN
wfpt(Valid.Speed.DLPFC).17.0   NaN
wfpt(Valid.Speed.FEF).17.0     NaN
wfpt(Valid.Speed.S1).17.0      NaN
wfpt(Valid.Speed.aPFC).17.0    NaN
Name: map, Length: 904, dtype: float64

In [45]:

    

v_accuracy, v_speed = m.nodes_db.node[['v(Valid.Accuracy.S1)', 'v(Valid.Accuracy.FEF)']]
hddm.analyze.plot_posterior_nodes([v_accuracy, v_speed])
plt.xlabel('drift-rate')
plt.ylabel('Posterior probability')
plt.title('Posterior of drift-rate group means')
plt.savefig('Posteriors.pdf')


    

In [42]:

    

print "P(Speed < Accuracy) = ", (v_speed.trace() < v_accuracy.trace()).mean()


    

    




P(Speed < Accuracy) =  0.513636363636

In [8]:

    

dmatrix("C(tms_site , Treatment('S1'))", data)


    

    
Out[8]:





DesignMatrix with shape (32640, 4)
  Columns:
    ['Intercept',
     "C(tms_site, Treatment('S1'))[T.DLPFC]",
     "C(tms_site, Treatment('S1'))[T.FEF]",
     "C(tms_site, Treatment('S1'))[T.aPFC]"]
  Terms:
    'Intercept' (column 0)
    "C(tms_site, Treatment('S1'))" (columns 1:4)
  (to view full data, use np.asarray(this_obj))

In [53]:

    

# Instantiate model object passing it our data (no need to call flip_errors() before passing it).
# This will tailor an individual hierarchical DDM around your dataset.
m = hddm.HDDMRegressor (data, "v ~ C(tms_site , Treatment('S1'))", "a ~ C(tms_site , Treatment('S1'))", depends_on={'a': ['tms_site', 'Attended', 'SAT'], 'v': ['tms_site', 'Attended', 'SAT']}, p_outlier=.05)
# find a good starting point which helps with the convergence.
m.find_starting_values()
# start drawing 7000 samples and discarding 5000 as burn-in
m.sample(5000, burn=200)


    

    




---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-53-f31de08d455c> in <module>()
      1 # Instantiate model object passing it our data (no need to call flip_errors() before passing it).
      2 # This will tailor an individual hierarchical DDM around your dataset.
----> 3 m = hddm.HDDMRegressor (data, depends_on={'a': ['tms_site', 'Attended', 'SAT'], 'v': ['tms_site', 'Attended', 'SAT']}, p_outlier=.05)
      4 # find a good starting point which helps with the convergence.
      5 m.find_starting_values()

TypeError: __init__() takes at least 3 arguments (2 given)

In [ ]: