30 January 2021 Guido Cattani
from pathlib import Path
import pandas as pd
import numpy as np
from scipy.stats import mannwhitneyu as mannwhitneyu
from scipy.stats import wilcoxon as wilcoxon
from scipy.stats import rankdata as rankdata
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns
pd.set_option("precision", 3)
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)
def read_SII():
f_in = '/media/guido/LACIE/Cingle_Guido/Master/SII_Sfront65_Nbcd65.xlsx'
p_in = Path(f_in)
df = pd.read_excel(p_in, sheet_name='SII_Sfront65_Nbcd65', header=0, nrows=85)
df = df.drop(['Unnamed: 0', 'Study_ID'], axis=1)
return df
SII = read_SII()
SII
Device | AC_path | BC_path | AC&BC_path | diff_combi_AC | |
---|---|---|---|---|---|
0 | BP110 | 0.651 | 0.114 | 0.651 | 0.000 |
1 | BP110 | 0.602 | 0.130 | 0.567 | -0.035 |
2 | BP110 | 0.646 | 0.123 | 0.613 | -0.033 |
3 | BP110 | 0.643 | 0.115 | 0.636 | -0.007 |
4 | BP110 | 0.651 | 0.114 | 0.651 | 0.000 |
... | ... | ... | ... | ... | ... |
65 | BAHA5P | 0.651 | 0.185 | 0.633 | -0.018 |
66 | BAHA5P | 0.651 | 0.290 | 0.651 | 0.000 |
67 | BAHA5P | 0.651 | 0.265 | 0.606 | -0.046 |
68 | BAHA5P | 0.625 | 0.228 | 0.592 | -0.033 |
69 | BAHA5P | 0.651 | 0.203 | 0.645 | -0.006 |
70 rows × 5 columns
# group data by device type and perform calculation of quantiles 10, 50, 90
dvc = SII.groupby('Device')
quantiles = [0.10, 0.50, 0.90]
q = dvc.quantile(q=quantiles)
q= q.round(decimals=4)
q = q.reset_index()
diq = {0:'BAHA5P P10', 1:'BAHA5P P50', 2:'BAHA5P P90',
3:'BP110 P10', 4:'BP110 P50', 5:'BP110 P90'}
q = q.rename(index=diq)
q.drop(['Device', 'level_1'], axis=1, inplace = True)
q
AC_path | BC_path | AC&BC_path | diff_combi_AC | |
---|---|---|---|---|
BAHA5P P10 | 0.627 | 0.205 | 0.587 | -5.660e-02 |
BAHA5P P50 | 0.651 | 0.247 | 0.631 | -1.850e-02 |
BAHA5P P90 | 0.651 | 0.281 | 0.650 | -9.000e-04 |
BP110 P10 | 0.626 | 0.102 | 0.572 | -5.960e-02 |
BP110 P50 | 0.651 | 0.115 | 0.627 | -6.700e-03 |
BP110 P90 | 0.651 | 0.138 | 0.651 | 0.000e+00 |
SII_bp110 = select_bp110(SII)
SII_bh5 = select_bh5(SII)
len(SII_bp110), len(SII_bh5)
(35, 35)
bh5 = SII_bh5.T.to_numpy()
bp110 = SII_bp110.T.to_numpy()
d1 = dict()
d2 = dict()
for i in range(0, 4):
(stat, pvalue) = mannwhitneyu(bp110[i], bh5[i], use_continuity=False, alternative='two-sided')
d1.update({i : stat})
d2.update({i : pvalue})
mwu = pd.DataFrame.from_dict([d1, d2])
rws = {0: 'Mann-Whitney U statistic', 1: 'p-value (two-sided)'}
clmns = {0 : 'AC_path', 1 : 'BC_path', 2 : 'AC&BC_path', 3 : 'diff_combi_AC'}
mwu.rename(index = rws, columns = clmns, inplace = True)
mwu = mwu.round(decimals=4)
mwu
AC_path | BC_path | AC&BC_path | diff_combi_AC | |
---|---|---|---|---|
Mann-Whitney U statistic | 438.500 | 0.0 | 622.500 | 742.500 |
p-value (two-sided) | 0.014 | 0.0 | 0.906 | 0.126 |
# Wilcoxon signed-rank test
w_bp110 = wilcoxon(bp110[3])
w_bp110
WilcoxonResult(statistic=0.0, pvalue=3.789619441580871e-06)
w_bh5 = wilcoxon(bh5[3])
w_bh5
WilcoxonResult(statistic=0.0, pvalue=7.952884142692727e-07)
# make dataframe with results wilcoxon
dwx = {'BAHA5P P50': [w_bh5[0], w_bh5[1]], 'BP110 P50' : [w_bp110[0], w_bp110[1]]}
dfx = pd.DataFrame.from_dict(dwx, orient = 'index', columns = ['Wilcoxon s-r test diff_combi_AC', 'p-value (2s)'])
dfx
Wilcoxon s-r test diff_combi_AC | p-value (2s) | |
---|---|---|
BAHA5P P50 | 0.0 | 7.953e-07 |
BP110 P50 | 0.0 | 3.790e-06 |
# make dataframe with all results
analysis_output = pd.concat([q, mwu])
analysis_output = pd.concat([analysis_output, dfx], axis = 1)
analysis_output
AC_path | BC_path | AC&BC_path | diff_combi_AC | Wilcoxon s-r test diff_combi_AC | p-value (2s) | |
---|---|---|---|---|---|---|
BAHA5P P10 | 0.627 | 0.205 | 0.587 | -5.660e-02 | NaN | NaN |
BAHA5P P50 | 0.651 | 0.247 | 0.631 | -1.850e-02 | 0.0 | 7.953e-07 |
BAHA5P P90 | 0.651 | 0.281 | 0.650 | -9.000e-04 | NaN | NaN |
BP110 P10 | 0.626 | 0.102 | 0.572 | -5.960e-02 | NaN | NaN |
BP110 P50 | 0.651 | 0.115 | 0.627 | -6.700e-03 | 0.0 | 3.790e-06 |
BP110 P90 | 0.651 | 0.138 | 0.651 | 0.000e+00 | NaN | NaN |
Mann-Whitney U statistic | 438.500 | 0.000 | 622.500 | 7.425e+02 | NaN | NaN |
p-value (two-sided) | 0.014 | 0.000 | 0.906 | 1.262e-01 | NaN | NaN |
# write to xlsx file
analysis_output.to_excel("/media/guido/LACIE/Cingle_Guido/Analysis_results/analysis_SII_Sfront65_Nbcd65.xlsx",
sheet_name='SII_Sfront65_Nbcd65')
# make a figure to plot SII for the 3 paths
ttl = 'BP110: SII for the combination path, air conduction path, bone conduction path, S in front 65 dB, N BCD side 65 dB'
SII_combi = SII_bp110['AC&BC_path']
SII_ac = SII_bp110['AC_path']
SII_bc = SII_bp110['BC_path']
tp = pd.melt(pd.DataFrame( {'AC & BC' : SII_combi, 'AC' : SII_ac, 'BC' : SII_bc }),
var_name = 'Transmission path', value_name = 'Speech Intelligibility Index')
fig, ax = plt.subplots(constrained_layout=True)
fig.set_figheight(6)
fig.set_figwidth(11)
ax = sns.swarmplot(data = tp, x = 'Transmission path', y = 'Speech Intelligibility Index',
hue = 'Transmission path', size=5, palette={'silver', 'grey', 'black'})
ax.set_title(ttl)
plt.legend(bbox_to_anchor=(0.15, 0.25), fontsize='large')
# save figure to file
plt.savefig('/media/guido/LACIE/Cingle_Guido/Analysis_results/BP110_SII_Sfront65_Nbcd65.tiff',
transparent=False, dpi=500, bbox_inches="tight")
plt.show()
# make a figure to plot SII for the 3 paths
ttl = 'BAHA5P: SII for the combination path, air conduction path, bone conduction path, S in front 65 dB, N BCD side 65 dB'
SII_combi = SII_bh5['AC&BC_path']
SII_ac = SII_bh5['AC_path']
SII_bc = SII_bh5['BC_path']
tp = pd.melt(pd.DataFrame( {'AC & BC' : SII_combi, 'AC' : SII_ac, 'BC' : SII_bc }),
var_name = 'Transmission path', value_name = 'Speech Intelligibility Index')
fig, ax = plt.subplots(constrained_layout=True)
fig.set_figheight(6)
fig.set_figwidth(12)
ax = sns.swarmplot(data = tp, x = 'Transmission path', y = 'Speech Intelligibility Index',
hue = 'Transmission path', size=6, palette={'silver', 'grey', 'black'})
ax.set_title(ttl)
plt.legend(bbox_to_anchor=(0.15, 0.25), fontsize='large')
# save figure to file
plt.savefig('/media/guido/LACIE/Cingle_Guido/Analysis_results/BAHA5P_SII_Sfront65_Nbcd65.tiff',
transparent=False, dpi=500, bbox_inches="tight")
plt.show()