#!/usr/bin/env python
# coding: utf-8

# In[1]:


get_ipython().run_line_magic('pylab', 'inline')
import os
import os.path as op
import glob
import numpy as np
import pandas as pd
import seaborn as sns
import statsmodels.api as sm
from statsmodels.formula.api import ols

# some, perhaps, relevant links/sources:
# http://stats.stackexchange.com/questions/24314/when-is-a-repeated-measures-anova-preferred-over-a-mixed-effects-model
# http://nbviewer.ipython.org/urls/umich.box.com/shared/static/6tfc1e0q6jincsv5pgfa.ipynb
# http://jpktd.blogspot.jp/2013/03/multiple-comparison-and-tukey-hsd-or_25.html
# http://statsmodels.sourceforge.net/devel/anova.html

resdir = '/Users/cplanting/ABR-vis/'
data_bera = pd.read_csv(op.join(resdir,'amp.csv'), sep=';',header=0)
print data_bera.head()

#bin data source: http://stackoverflow.com/questions/16947336/binning-a-dataframe-in-pandas-in-python
bins = np.linspace(5, 85, 9)
int_binned = data_bera.groupby(pd.cut(data_bera.intensity, bins))
def int_round(x,mult):
    #round to closest multiple of 10; 5.0 -> 0.0; 5.1 -> 10.0
    return np.round(x / mult) * mult

data_bera['intensity (dB SPL)'] = int_round(data_bera['intensity'],10.0)
g0 = sns.factorplot("intensity (dB SPL)","amplitude","stimulus_type",kind='bar',
                    ci=68,n_boot=50000, estimator=np.mean, 
                    data=tmp, size=6, palette="Set1")


# In[5]:


g = sns.lmplot(x="intensity (dB SPL)", y="amplitude", hue="stimulus_type", x_estimator=np.mean, data=data_bera,
               y_jitter=10,palette="Set1")


# In[16]:


g = sns.lmplot(x="intensity (dB SPL)", y="latency", hue="stimulus_type", x_estimator=np.mean, data=data_bera,
               x_jitter=4,palette="Set1")
level = [10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90]
lat = [8.73,8.42,8.12,7.82,7.52,7.24,6.98,6.74,6.52,6.33,6.16,6.01,5.9,5.8,5.71,5.64,5.58]
level2 = [20,30,40,50,60,70,80,90]
lat2 = [7.51, 7.13, 6.59, 6.17, 6.01, 5.8, 5.62, 5.55]
lat2_sd = [0.42, 0.37, 0.43, 0.27, 0.27, 0.2, 0.12,0.12]
plot(level,lat,'k-',lw=3,label='Norm Click (Prosser, 1988)')
plot(level2,lat2,'g-',lw=3,label='UMCG norm Click')

lev_chirp = [10,20,30,40,60,80]
lat_chirp = [7.01, 6.42, 5.79, 5.17, 3.93, 2.85]
sd_chirp = [0.49, 0.36, 0.35, 0.40, 0.62, 0.84]
errorbar(lev_chirp,lat_chirp,yerr=sd_chirp,color='k',label='Norm Chirp (Elberling, 2010)')
legend()
g.ax.set_ylim([0,12])


# In[11]:


model = sm.MixedLM.from_formula("amplitude ~ intensity * stimulus_type", tmp, groups=tmp["sub_id"])
result = model.fit()
print result.summary()


# In[28]:


model = sm.MixedLM.from_formula("latency ~ intensity * stimulus_type", tmp, groups=tmp["sub_id"])
result = model.fit()
print result.summary()


# In[26]:


model_anova = ols('latency ~ stimulus_type*intensity',data=tmp).fit()
table = sm.stats.anova_lm(model_anova, typ=2) # Type 2 ANOVA DataFrame
print model_anova.summary()
print table
#interestingly, by using plain ANOVA (without repeated measures), the effect dissolve in the between-subject variance! 
#therefor use MixedLM


# In[22]:


spect_file = '/Users/cplanting/ProjectsScience/ABR/raw_data/spectrum_click_chirp_bone.csv'
spectrum_df = pd.read_csv(spect_file,sep=';',header=0,decimal=',')
#spectrum_df.head


# In[23]:


fig, ax = plt.subplots(figsize=(8, 6))
semilogx(spectrum_df.f,spectrum_df.WN,'k--',lw=2,label='White Noise')
semilogx(spectrum_df.f,spectrum_df.click,'b-',lw=3,label='click')
semilogx(spectrum_df.f,spectrum_df.chirp,'r-',lw=3,label='chirp')
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('dB re 1 muN')
plt.legend()


# In[4]:


# now same but with individual point
g = sns.lmplot(x="intensity (dB SPL)", y="amplitude", hue="stimulus_type", order=1,data=data_bera,
               y_jitter=10,palette="Set1")


# In[3]:


g = sns.lmplot(x="intensity (dB SPL)", y="latency", hue="stimulus_type", order=2,data=data_bera,
               y_jitter=1,palette="Set1")

#Prosser Arslan 1998
level = [10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90]
lat = [8.73,8.42,8.12,7.82,7.52,7.24,6.98,6.74,6.52,6.33,6.16,6.01,5.9,5.8,5.71,5.64,5.58]
level2 = [20,30,40,50,60,70,80,90]

#http://www.baaudiology.org/files/9313/8597/9894/R_04_Abeir_Osman_-_Characteristics_of_Auditory_Brainstem_Response.pdf
lat2 = [7.51, 7.13, 6.59, 6.17, 6.01, 5.8, 5.62, 5.55]
lat2_sd = [0.42, 0.37, 0.43, 0.27, 0.27, 0.2, 0.12,0.12]
plot(level,lat,'k-',lw=3,label='Prosser, 1988')