14-03-2021 Guido Cattani
from pathlib import Path
import pandas as pd
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.lines as mlines
import matplotlib.patches as mpatches
from scipy.stats import linregress as linregress
# read difference sensation levels BC - AC simulated with a 65 dB ISDS input signal at BCD side (90 degree angle)
def read_diff_sl():
f_in = '/media/guido/LACIE/Cingle_Guido/Master/Diff_SL.xlsx'
p_in = Path(f_in)
df = pd.read_excel(p_in, engine='openpyxl', sheet_name='Diff_SL_65dB_90deg', header=0, nrows=85)
df = df.drop(['Unnamed: 0'], axis=1)
df = df.drop(columns='Study_ID')
df = df.fillna(pd.NA)
return df
# read measured BC direct thresholds
def read_BC_direct_thresholds():
f_in = '/media/guido/LACIE/Cingle_Guido/Master/BC_direct_thresholds.xlsx'
p_in = Path(f_in)
df = pd.read_excel(p_in, engine='openpyxl', sheet_name=0, header=0, nrows=85)
df = df.drop(['Unnamed: 0'], axis=1)
df = df.drop(columns='Study_ID')
df = df.fillna(pd.NA)
return df
def select_bp110(df):
# select BP110 data
is_bp110 = df['Device']=='BP110'
df_bp110 = df[is_bp110]
df_bp110.pop('Device')
return(df_bp110)
def select_bh5(df):
# select BAHA5P data
is_baha5p = df['Device']=='BAHA5P'
df_baha5p = df[is_baha5p]
df_baha5p.pop('Device')
return(df_baha5p)
read_BC_direct_thresholds()
Device | BCd_250 | BCd_500 | BCd_1000 | BCd_1500 | BCd_2000 | BCd_3000 | BCd_4000 | |
---|---|---|---|---|---|---|---|---|
0 | BP110 | 25.0 | 35.0 | 35 | 40.0 | 45 | 40.0 | 35.0 |
1 | BP110 | 25.0 | 35.0 | 30 | 10.0 | 30 | 50.0 | 70.0 |
2 | BP110 | 30.0 | 40.0 | 20 | <NA> | 25 | 45.0 | 70.0 |
3 | BP110 | 20.0 | 20.0 | 20 | <NA> | 30 | <NA> | 40.0 |
4 | BP110 | 15.0 | 10.0 | 10 | <NA> | 20 | <NA> | 25.0 |
... | ... | ... | ... | ... | ... | ... | ... | ... |
80 | BAHA5P | 30.0 | 45.0 | 20 | 40.0 | 35 | 25.0 | 45.0 |
81 | BAHA5P | 20.0 | 40.0 | 30 | 30.0 | 25 | 35.0 | 45.0 |
82 | BAHA5P | 25.0 | 25.0 | 25 | 5.0 | 15 | 35.0 | 50.0 |
83 | BAHA5P | 30.0 | 45.0 | 20 | 35.0 | 50 | 65.0 | 60.0 |
84 | BAHA5P | 30.0 | 35.0 | 30 | 35.0 | 30 | 15.0 | 35.0 |
85 rows × 8 columns
bct = read_BC_direct_thresholds()
dsl = read_diff_sl()
bct1 = bct[['Device', 'BCd_1000', 'BCd_2000', 'BCd_4000']].copy()
dsl1 = dsl[['1000_Hz', '2000_Hz', '4000_Hz']].copy()
combi = pd.concat([bct1, dsl1], axis=1)
combi = combi.dropna()
combi.head()
Device | BCd_1000 | BCd_2000 | BCd_4000 | 1000_Hz | 2000_Hz | 4000_Hz | |
---|---|---|---|---|---|---|---|
0 | BP110 | 35 | 45 | 35.0 | -15 | -15 | -12.0 |
1 | BP110 | 30 | 30 | 70.0 | -9 | 10 | 14.0 |
2 | BP110 | 20 | 25 | 70.0 | 7 | 6 | -16.0 |
3 | BP110 | 20 | 30 | 40.0 | -6 | 12 | 21.0 |
4 | BP110 | 10 | 20 | 25.0 | 6 | 6 | 13.0 |
combi_bh5 = select_bh5(combi)
combi_bp110 = select_bp110(combi)
combi_bh5.head()
BCd_1000 | BCd_2000 | BCd_4000 | 1000_Hz | 2000_Hz | 4000_Hz | |
---|---|---|---|---|---|---|
20 | 20 | 30 | 35.0 | -10 | 10 | 19.0 |
37 | 45 | 15 | 35.0 | -14 | 17 | -1.0 |
47 | 20 | 30 | 25.0 | -10 | 4 | 3.0 |
51 | 15 | 30 | 30.0 | 3 | 4 | 13.0 |
53 | 20 | 35 | 50.0 | -6 | 5 | 15.0 |
for i in range(1, 4):
bcdt = (combi.iloc[:,i]).to_numpy(dtype='float')
dsl = (combi.iloc[:,i+3]).to_numpy(dtype='float')
lrgrs = linregress(bcdt, dsl)
print(2**(i-1),'kHz', lrgrs)
1 kHz LinregressResult(slope=-0.5672327129466105, intercept=7.81782366526287, rvalue=-0.5172168592408936, pvalue=4.7168204555163794e-07, stderr=0.10365292165925495, intercept_stderr=2.7661979723777543) 2 kHz LinregressResult(slope=-0.597643330218195, intercept=23.806184226565353, rvalue=-0.6016700572944104, pvalue=1.422720814193018e-09, stderr=0.08761626934280549, intercept_stderr=2.8976357435539493) 4 kHz LinregressResult(slope=-0.10450529164386031, intercept=14.198209714322626, rvalue=-0.11238133822398975, pvalue=0.3087790689054309, stderr=0.10204158961898098, intercept_stderr=4.47212549752971)
for i in range(3):
bcdt_bp110 = (combi_bp110.iloc[:,i]).to_numpy(dtype='float')
dsl_bp110 = (combi_bp110.iloc[:,i+3]).to_numpy(dtype='float')
lrgrs = linregress(bcdt_bp110, dsl_bp110)
print(2**i,'kHz', lrgrs)
1 kHz LinregressResult(slope=-0.6188358842587421, intercept=12.793207462853314, rvalue=-0.6485940580525047, pvalue=6.190847586083582e-07, stderr=0.10707424555275978, intercept_stderr=2.929292858105807) 2 kHz LinregressResult(slope=-0.6428948360346777, intercept=23.428320140721194, rvalue=-0.6291316057857251, pvalue=1.6734608791547085e-06, stderr=0.11711364948461658, intercept_stderr=3.993937475900152) 4 kHz LinregressResult(slope=-0.07406505479213776, intercept=10.70099147677857, rvalue=-0.08195876480755489, pvalue=0.5797209972779328, stderr=0.13279309010038703, intercept_stderr=6.126707990147849)
for i in range(3):
bcdt_bh5 = (combi_bh5.iloc[:,i]).to_numpy(dtype='float')
dsl_bh5 = (combi_bh5.iloc[:,i+3]).to_numpy(dtype='float')
lrgrs = linregress(bcdt_bh5, dsl_bh5)
print(2**i,'kHz', lrgrs)
1 kHz LinregressResult(slope=-0.5143484626647145, intercept=1.6325036603221061, rvalue=-0.4654542088220822, pvalue=0.004228068769034862, stderr=0.16773342279481795, intercept_stderr=4.321826314145049) 2 kHz LinregressResult(slope=-0.4617137648131267, intercept=22.17411121239745, rvalue=-0.5444279167359414, pvalue=0.0005976444074497801, stderr=0.12199874142408537, intercept_stderr=3.8606171404022795) 4 kHz LinregressResult(slope=-0.0480964898716561, intercept=14.874427091387041, rvalue=-0.05009373294395535, pvalue=0.7717112912092209, stderr=0.16445418786976793, intercept_stderr=6.667158835938594)
bcdt_1k_bh5 = combi_bh5.iloc[:,0].to_numpy(dtype='float')
dsl_1k_bh5 = (combi_bh5.iloc[:,3]).to_numpy(dtype='float')
lrgrs_1k_bh5 = linregress(bcdt_1k_bh5, dsl_1k_bh5)
r2_1k_bh5 = (lrgrs_1k_bh5[2])**2
r2_1k_bh5 = round(r2_1k_bh5, 2)
s_1k = 'r square 1 kHz =' + str(r2_1k_bh5)
print(s_1k)
r square 1 kHz =0.22
plt.scatter(bcdt_1k_bh5, dsl_1k_bh5, c='black')
plt.title("Relation BC direct thresholds and diff SL at 1 kHz for BAHA5P")
plt.xlabel("BC direct thresholds at 1 kHz [dB HL]")
plt.ylabel("Diff. SL BC-AC at 1 kHz [dB]")
slope = lrgrs_1k_bh5[0]
intercept = lrgrs_1k_bh5[1]
plt.plot(bcdt_1k_bh5, slope*bcdt_1k_bh5 + intercept, color='grey')
plt.text(32, -23, s_1k)
plt.show()
bcdt_2k_bh5 = combi_bh5.iloc[:,1].to_numpy(dtype='float')
dsl_2k_bh5 = (combi_bh5.iloc[:,4]).to_numpy(dtype='float')
lrgrs_2k_bh5 = linregress(bcdt_2k_bh5, dsl_2k_bh5)
r2_2k_bh5= (lrgrs_2k_bh5[2])**2
r2_2k_bh5 = round(r2_2k_bh5, 2)
s_2k = 'r square 2 kHz =' + str(r2_2k_bh5)
print(s_2k)
r square 2 kHz =0.3
plt.scatter(bcdt_2k_bh5, dsl_2k_bh5, c='black')
plt.title("Relation BC direct thresholds and diff SL at 2 kHz for BAHA5P")
plt.xlabel("BC direct thresholds at 2 kHz [dB HL]")
plt.ylabel("Diff. SL BC-AC at 2 kHz [dB]")
slope = lrgrs_2k_bh5[0]
intercept = lrgrs_2k_bh5[1]
plt.plot(bcdt_2k_bh5, slope*bcdt_2k_bh5 + intercept, color='grey')
plt.text(32, -1, s_2k)
plt.show()
bcdt_4k_bh5 = combi_bh5.iloc[:,2].to_numpy(dtype='float')
dsl_4k_bh5 = (combi_bh5.iloc[:,5]).to_numpy(dtype='float')
lrgrs_4k_bh5 = linregress(bcdt_4k_bh5, dsl_4k_bh5)
r2_4k_bh5= (lrgrs_4k_bh5[2])**2
r2_4k_bh5 = round(r2_4k_bh5, 2)
s_4k = 'r square 4 kHz =' + str(r2_4k_bh5)
print(s_4k)
r square 4 kHz =0.0
plt.scatter(bcdt_4k_bh5, dsl_4k_bh5, c='black')
plt.title("Relation BC direct thresholds and diff SL at 4 kHz for BAHA5P")
plt.xlabel("BC direct thresholds at 4 kHz [dB HL]")
plt.ylabel("Diff. SL BC-AC at 4 kHz [dB]")
slope = lrgrs_4k_bh5[0]
intercept = lrgrs_4k_bh5[1]
plt.plot(bcdt_4k_bh5, slope*bcdt_4k_bh5 + intercept, color='grey')
plt.text(32, 15, s_4k)
plt.show()