<jupyter_start><jupyter_text>Hierarchical Bayesian Linear Regression using PyStan**<NAME>, 2021** Here we fit response times against conflict, also checking for a context interaction. Further description of the analysis and visualization of the results can be found in the main manuscript<jupyter_code>import numpy as np
import pandas as pd 
import matplotlib.pyplot as plt
import pystan
import glob as glob 
import time as time
import arviz as az
plt.style.use('ara')

# Load participant data
filename = glob.glob('../data/behaviour/data_all_partipants_20210623102746.csv') 
dat = pd.read_csv(filename[0],index_col = 0)<jupyter_output><empty_output><jupyter_text>Specifying the model<jupyter_code>model = '''
data {
  int<lower=1> N;
  vector[N] ishard;
  vector[N] logrt;
  vector[N] conflict;
  int<lower=1> N_subjects;
  int<lower = 1, upper = N> vpn[N];
  int<lower=1> N_rep;
  vector[N_rep] conflict_rep;
}

parameters {
  // hyperparameters
  real mu_intercept;
  real mu_beta_conflict;
  real mu_beta_ishard;
  real mu_beta_interaction;
  real<lower=0> sigma_intercept;
  real<lower=0> sigma_beta_conflict;
  real<lower=0> sigma_ishard;
  real<lower=0> sigma_interaction;
  
  // paramters 
  vector[N_subjects] intercept;
  vector[N_subjects] beta_conflict;
  vector[N_subjects] beta_ishard;
  vector[N_subjects] beta_interaction; 
  real<lower=0>  sigma; //same data level noise for all subjects 
  
}

model {

  // hyper priors 
  mu_intercept ~ normal(0,10); 
  mu_beta_conflict ~ normal(0,10); 
  mu_beta_ishard ~ normal(0,10); 
  mu_beta_interaction ~ normal(0,10); 
  sigma_intercept ~ normal(0,10); 
  sigma_beta_conflict ~ normal(0,10); 
  sigma_ishard ~ normal(0,10); 
  sigma_interaction ~ normal(0,10); 

  
  //priors 
  intercept ~ normal(mu_intercept,sigma_intercept); 
  beta_conflict ~ normal(mu_beta_conflict,sigma_beta_conflict); 
  beta_ishard ~ normal(mu_beta_ishard,sigma_ishard); 
  beta_interaction ~ normal(mu_beta_interaction,sigma_interaction);
  sigma ~ normal(0, 10);

  logrt ~ normal(intercept[vpn] + beta_conflict[vpn].*conflict + beta_ishard[vpn].*ishard + beta_interaction[vpn] .* conflict .* ishard,sigma);
}

generated quantities {

  vector[N_subjects] intercept_rep;
  vector[N_subjects] beta_conflict_rep;
  vector[N_subjects] beta_ishard_rep;
  vector[N_subjects] beta_interaction_rep;
  vector[N_rep] rt_new_easy;
  vector[N_rep] rt_new_hard;
  
  for (n in 1:N_subjects){
    intercept_rep[n] = normal_rng(intercept[n], sigma_intercept);
    beta_conflict_rep[n] = normal_rng(beta_conflict[n], sigma_beta_conflict);
    beta_ishard_rep[n] = normal_rng(beta_ishard[n], sigma_ishard);
    beta_interaction_rep[n] = normal_rng(beta_interaction[n], sigma_interaction);
    }

  for (n in 1:N_rep){
    rt_new_easy[n] = lognormal_rng(intercept_rep[vpn[n]] + conflict_rep[n] * beta_conflict_rep[vpn[n]] + beta_ishard_rep[vpn[n]]*0 + beta_interaction_rep[vpn[n]] * conflict_rep[n] * 0 ,sigma);
    rt_new_hard[n] = lognormal_rng(intercept_rep[vpn[n]] + conflict_rep[n] * beta_conflict_rep[vpn[n]] + beta_ishard_rep[vpn[n]]*1 + beta_interaction_rep[vpn[n]] * conflict_rep[n] * 1 ,sigma);
    } 
}
'''<jupyter_output><empty_output><jupyter_text>Compiling<jupyter_code>sm = pystan.StanModel(model_code=model,verbose = False)<jupyter_output>INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_1989a5f201a87f58e7e2e78ff7aa273a NOW.<jupyter_text>Specifying the data<jupyter_code>n_iter = 1000
n_warmup = 500
n_chains = 4
control_dict = dict(adapt_delta=0.95)

dat['ishard'] = ((dat['reward'] == 2) | (dat['reward'] == 3)) 
idx = (dat['timeout'] == 0)
logrt =  np.log(dat.loc[idx,['reaction_time']].to_numpy().squeeze())
conflict = dat.loc[idx,['conflict']].to_numpy().squeeze()
ishard = dat.loc[idx,['ishard']].to_numpy(dtype='int').squeeze()
N = len(logrt)
vpn = dat.loc[idx,['vpn']].to_numpy().squeeze() - 100
N_subjects = len(np.unique(vpn))
conflict_rep = conflict
N_rep = len(conflict_rep)

dat_dict = {'N':N,         
            'logrt':logrt,
            'conflict':conflict, 
            'ishard':ishard ,
            'N_subjects':N_subjects,
            'vpn':vpn,
            'N_rep':N_rep,
            'conflict_rep':conflict_rep           
            }<jupyter_output><empty_output><jupyter_text>Sampling posterior<jupyter_code>res = sm.sampling(data=dat_dict, iter=n_iter,  warmup=n_warmup, thin=1, chains=n_chains,control=control_dict,seed=101, verbose = False);<jupyter_output>WARNING:pystan:Maximum (flat) parameter count (1000) exceeded: skipping diagnostic tests for n_eff and Rhat.
To run all diagnostics call pystan.check_hmc_diagnostics(fit)
WARNING:pystan:4 of 2000 iterations ended with a divergence (0.2 %).
WARNING:pystan:Try running with adapt_delta larger than 0.95 to remove the divergences.
{'fraction_non_alphanumeric': 0.12555112324165443, 'fraction_numerical': 0.024564350199454126, 'mean_word_length': 5.403225806451613, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 29, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '367102', 'license': 'MIT', 'n_tokens_mistral': 1933, 'n_tokens_neox': 1797, 'n_words': 319, 'path': 'code/fit_RT_model.ipynb', 'repo_name': 'fmott/context_dependent_planning', 'size': '4768'}