<jupyter_start><jupyter_text>Calculate differences- Calculate differences using different models on the effect of Trauma vs. Sad on hippocampus<jupyter_code>import warnings
warnings.filterwarnings('ignore')
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import pymc3 as pm
import arviz as az

# regular GLMs
import pingouin as pg
import statsmodels.api as sm
import statsmodels.formula.api as smf
## Load amygdala database
df = pd.read_csv('amg_TraumavsNeutral.csv')#('TraumavsNeutral_threeSes.csv')
# create a long format for analysis
df_long = pd.melt(df, id_vars=['scr_id','groupIdx'], value_vars=['hippo1','hippo2','hippo3'])
df_long
## Start with simple GLM models (using statsmodel)
glm = smf.mixedlm('value ~ variable * groupIdx', data=df_long, groups = df_long.scr_id, 
                 missing='drop')
model1 = glm.fit()
print(model1.summary())
# mixed anova
aov = pg.mixed_anova(data=df_long, dv='value', between='groupIdx', within='variable',
                     subject='scr_id', correction=False, effsize="np2")
pg.print_table(aov)
pg.ttest(df.hippo2[df.groupIdx==1], df.hippo2[df.groupIdx==0], paired=False,
         correction=False)<jupyter_output><empty_output><jupyter_text>Build a Bayesian model- Here we would model each subject under each group<jupyter_code>sns.set(font_scale = 1.15)
sns.set_style("white")
g = sns.catplot(y='value',x='variable', hue='groupIdx', data=df_long, kind='point', ci=68,
               hue_order=[1,0])
g.set_xlabels("")
g.set_xticklabels(
    labels=["Before Treatment", "End of Treatment", "30-days"], rotation=30)

g.set_ylabels("Hippocampus (Trauma > Neutral)")
g._legend.set_title("Group")
new_labels = ['Ketamine', 'Midazolam']
for t, l in zip(g._legend.texts, new_labels): t.set_text(l)
g.savefig('hippo_pointGraph.png', dpi=300, bbox_inches='tight')
#df[['hippo1','hippo2','hippo3', 'hippo4','groupIdx']].groupby('groupIdx').describe()
df_long.groupby(['groupIdx','variable']).describe()<jupyter_output><empty_output><jupyter_text>Use time 1 as covariate to remove differences in baseline<jupyter_code>## Load amygdala database
df = pd.read_csv('amg_TraumavsNeutral.csv')#('TraumavsNeutral_threeSes.csv')
# create a long format for analysis
df_long2 = pd.melt(df, id_vars=['scr_id','groupIdx', 'hippo1'], value_vars=['hippo2','hippo3'])
df_long2
# take indexes for building the model
df_long2['scr_id'] = df_long2.scr_id.astype('category')
sub_idx = df_long2.scr_id.cat.codes.values # get county category index values (613 x 1) vec with nums [0-7]
print(sub_idx)
n_subs = len(np.unique(df_long2.scr_id))
print(n_subs)
# set time as index also
#df_long.groupIdx = pd.Categorical(df_long.groupIdx)
df_long2.variable = pd.Categorical(df_long2.variable, 
                    categories=['hippo2','hippo3'])
timeIDX = df_long2.variable.cat.codes.values
print(timeIDX)
df_long2['time'] = timeIDX

with pm.Model() as hierarchical:
      
    a = pm.Normal('a', mu=0, sd=0.1, shape=n_subs)
    cov = pm.Normal('cov', mu=0.05, sd=0.1)
    b = pm.Normal('b', mu=0, sd=0.1, shape=2)
    
    
   # bMed = pm.Normal('bMed', mu=0, sd=0.1, shape=2) #
   # sigma = pm.HalfCauchy('sigma', 1)
    interact = pm.Normal('interaction', mu=0, sd=1, shape=(2,2))
  
    # Model error
    eps = pm.HalfCauchy('eps', 1)

    y_hat = a[sub_idx] + b[timeIDX] + interact[timeIDX,df_long2.groupIdx] + cov*df_long2.hippo1

    # Data likelihood
    likelihood = pm.Normal('likelihood', mu=y_hat,
                           sigma=eps, observed=df_long2.value)
    
    
    # calculate differences and effect size
    diff_of_means1 = pm.Deterministic('differenceGroups end of treatment', interact[0,1] - interact[0,0])
    diff_of_means2 = pm.Deterministic('differenceGroups 30 days', interact[1,1] - interact[1,0])
   
    # effect size
    effect_size1 = pm.Deterministic('effect size2_1', 
                        abs(diff_of_means1) / eps)
    
    effect_size2 = pm.Deterministic('effect size3_1', 
                        abs(diff_of_means2) / eps)
    
#     effect_size3 = pm.Deterministic('effect size4_1', 
#                         abs(diff_of_means3) / eps)
    
    ppChecks = pm.sample_prior_predictive(samples = 8000)
    trace = pm.sample(draws=2000, tune=1000, target_accept=0.9, return_inferencedata=True, cores=4)
az.summary(trace, hdi_prob=.9)
priorC = ppChecks['interaction'][:,0,0] - ppChecks['interaction'][:,0,1] 
sns.displot(priorC)
tr = trace.posterior.stack(draws=("chain", "draw"))
sns.histplot(tr.interaction.values[0,0,:,], label='midazolam', color='blue')
sns.histplot(tr.interaction.values[0,1,:], label = 'ketamine', color='red')
plt.legend()
#tr.interaction.values[0,0,:,]
after = (tr.interaction.values[0,1,:,] - tr.interaction.values[0,0,:,])# / np.sqrt(
az.plot_posterior(after, kind='hist', ref_val=0, hdi_prob=.9)
print(np.sum(after>=0) / len(after))
calcBF(y2=priorC, y1=after)
after2 = (tr.interaction.values[1,1,:,] - tr.interaction.values[1,0,:,])
calcBF(y1=priorC, y2=after2)<jupyter_output><empty_output><jupyter_text>Next we calculate the difference between the groups in each time- Althgough the interaction wasn't significant, we still look at this one - as there's a group difference<jupyter_code>with pm.Model() as model_glm:
    pm.GLM.from_formula('hippo1 ~ groupIdx', df)
    traceSes1 = pm.sample(return_inferencedata=True)
az.summary(traceSes1, hdi_prob=.9)
with pm.Model() as model_glm:
    pm.GLM.from_formula('hippo2 ~ groupIdx', df)
    traceSes2 = pm.sample(return_inferencedata=True)
az.summary(traceSes2, hdi_prob=.9)
with pm.Model() as model_glm:
    pm.GLM.from_formula('hippo3 ~ groupIdx', df)
    traceSes3 = pm.sample(return_inferencedata=True)
az.summary(traceSes3, hdi_prob=.9)<jupyter_output>Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [sd, groupIdx, Intercept]<jupyter_text>Use a different way<jupyter_code>## create a function that will simply run the model for each session
def runModel(variable):
    with pm.Model() as hierarchical:
        nu = pm.Gamma('nu',2,0.1)
        a = pm.Normal('a', mu=0, sd=1) 
        bMed = pm.StudentT('bMed', mu=0,nu=nu , shape=2) #

        # Model error
        eps = pm.HalfCauchy('eps', 5)

        y_hat = a + bMed[df.groupIdx]

        # Data likelihood
        likelihood = pm.Normal('likelihood', mu=y_hat,
                               sigma=eps, observed=df[str(variable)])

        trace = pm.sample(draws=2000, tune=1000,target_accept=0.9, return_inferencedata=False, cores=4)
        diff = trace['bMed'][:,1] - trace['bMed'][:,0]
    return trace, diff
ses1trace, ses1Diff = runModel('hippo1')
pm.summary(ses1trace)
sns.distplot(ses1Diff)
print(f'Group difference is {sum(ses1Diff>0)/len(ses1Diff)}')
ses2trace, ses2Diff = runModel('hippo2')
pm.summary(ses2trace)
sns.distplot(ses2Diff)
print(f'Group difference is {sum(ses2Diff>0)/len(ses2Diff)}')
ses3trace, ses3Diff = runModel('hippo3')
pm.summary(ses3trace)
sns.distplot(ses3Diff)
print(f'Group difference is {sum(ses3Diff>0)/len(ses3Diff)}')<jupyter_output>Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [likelihood_missing, eps, bMed, a]<jupyter_text>Plotting<jupyter_code># set variables
def plotDiff(depVar,indepVar, trace):
    sns.set_style("ticks") # set style
    y = str(depVar)
    dfPlot = df
    ci = np.quantile(trace.posterior[str(indepVar)], [.05,.95])
    fig, (ax1, ax2) = plt.subplots(1,2, figsize=(3, 5),gridspec_kw={'width_ratios': [1, .2],
                                                            'wspace':.1})
    g1 = sns.stripplot(y= y, x='group', data=dfPlot, size = 8, ax=ax1)
    sns.boxplot(y= y, x='group', data=dfPlot,  ax=ax1,
                boxprops=dict(alpha=.3))
    g2 = sns.distplot(trace.posterior[str(indepVar)], ax = ax2, vertical=True)
    ax2.vlines(x=0.001,ymin=ci[0], ymax=ci[1], color='black', 
               linewidth = 3, linestyle = "-")

    #g3.set_ylim(-.7, .7)
    #ax1.set_ylim(-.7,.7)
    ax2.set_ylim(g1.get_ylim()) # use first graph's limits to get the relevant for this one
    ax2.yaxis.tick_right()
    ax2.set_xticks([])
    ax2.set_ylabel("Difference between groups", fontsize=14) 
    ax2.yaxis.set_label_position("right")
    ax1.set_ylabel("Hippocampus reactivity to traumatic script", fontsize=12)
    ax1.set_xlabel("Group", fontsize=14)
    fig.savefig( str(depVar) + 'Reactivity.png', dpi=300, bbox_inches='tight')
    return fig
def plotSubfigs(subfigs,subNum, depVar,indepVar, trace, title):
    # gets subfig and the index number of it to plot the same as above
    sns.set_style("ticks") # set style
    y = str(depVar)
    dfPlot = df
    ci = np.quantile(trace.posterior[str(indepVar)], [.05,.95])
    axsLeft = subfigs[subNum].subplots(1,2, gridspec_kw={'width_ratios': [1, .2],'wspace':.01})
    ax1 = axsLeft[0]
    ax2 = axsLeft[1]
    g1 = sns.stripplot(y= y, x='group', data=dfPlot, size = 8, ax=ax1)
    sns.boxplot(y= y, x='group', data=dfPlot,  ax=ax1,
                 boxprops=dict(alpha=.3))
    g2 = sns.distplot(trace.posterior[str(indepVar)], ax = ax2, vertical=True)
    ax2.vlines(x=0.001,ymin=ci[0], ymax=ci[1], color='black', 
                linewidth = 3, linestyle = "-")

    ax2.set_ylim(-1,1)#(g1.get_ylim()) # use first graph's limits to get the relevant for this one
    ax1.set_ylim(-1, 1)
    ax2.yaxis.tick_right()
    ax2.set_xticks([])
    ax2.set_ylabel("Difference between groups", fontsize=14) 
    ax2.yaxis.set_label_position("right")
    ax1.set_ylabel("Hippocampus reactivity to traumatic script", fontsize=12)
    ax1.set_xlabel('')#("Group", fontsize=14)
    ax1.set_title(title, fontsize=15)
    #ax1.title.set_style(font="large")
fig = plt.figure(constrained_layout=True, figsize=(15, 6))
subfigs = fig.subfigures(1, 3, wspace=0.1)
plotSubfigs(subfigs, 0, 'hippo1','groupIdx',traceSes1, "Before Treatment")

plotSubfigs(subfigs, 1, 'hippo2','groupIdx',traceSes2, 'After Treatment')

plotSubfigs(subfigs, 2,'hippo3','groupIdx',traceSes3, '30-Days')

fig.savefig("Hipp_reactivity.png",  dpi=300, bbox_inches='tight')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15338508401297044, 'fraction_numerical': 0.03232779797582785, 'mean_word_length': 4.385185185185185, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 28, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '614665', 'license': 'MIT', 'n_tokens_mistral': 4050, 'n_tokens_neox': 3858, 'n_words': 731, 'path': 'func/Stats_TrtVsNeutralHippo.ipynb', 'repo_name': 'orduek/KPE', 'size': '10177'}