import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib.ticker import NullFormatter
import sklearn.metrics as metrics
import sklearn.linear_model as linear_model
import sklearn.model_selection as model_selection
from sklearn.svm import SVC
from sklearn.preprocessing import MinMaxScaler
from sklearn.feature_selection import SelectKBest
from sklearn.pipeline import Pipeline
import seaborn as sns
sns.set('talk', 'whitegrid', 'dark', font_scale=1.5, rc={"lines.linewidth": 2, 'grid.linestyle': '--'})df_features = pd.read_pickle("df_features_estandarizado.pickle")
df_features.head()| agrupacion_feature | media | ... | std | |||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| feature | delta | theta | alpha | beta | gamma | delta_norm | theta_norm | alpha_norm | beta_norm | gamma_norm | ... | alpha | beta | gamma | delta_norm | theta_norm | alpha_norm | beta_norm | gamma_norm | intra | inter | |
| tipo | indice_paciente | |||||||||||||||||||||
| P | 0 | 0.500218 | -0.240771 | -0.688548 | -0.459638 | -0.406289 | 0.837666 | -0.071074 | -0.666111 | -0.396677 | -0.322384 | ... | -0.603870 | 0.195745 | -0.339861 | 0.542809 | 0.114840 | -0.453845 | 0.337590 | -0.055859 | 1.574465 | 0.652673 |
| 1 | 0.970334 | 0.601456 | -0.762397 | -0.775901 | -0.369460 | 1.079809 | 0.547332 | -0.920190 | -1.062210 | -0.398412 | ... | -0.684283 | -0.636885 | -0.386052 | 0.706450 | 1.161500 | -1.161972 | -1.153092 | -0.429586 | -0.158501 | 1.019270 | |
| 2 | 2.070471 | 0.795315 | -0.646674 | -0.835162 | -0.482753 | 1.370391 | -0.193502 | -0.852261 | -1.294888 | -0.517556 | ... | -0.661683 | -0.924978 | -0.559998 | -0.959542 | -0.148300 | -0.983480 | -1.932346 | -0.810451 | 1.055501 | 1.554689 | |
| 3 | 0.734391 | 1.359020 | -0.410337 | -0.592436 | -0.501728 | 0.578035 | 1.356012 | -0.573282 | -0.948825 | -0.521571 | ... | -0.343123 | -0.445636 | -0.519770 | 1.467841 | 1.439780 | -0.350221 | -0.981742 | -0.788827 | 0.124952 | 1.024255 | |
| 4 | 1.537761 | 0.742978 | -0.572930 | -0.462701 | -0.280550 | 1.062737 | 0.100820 | -0.785426 | -0.944991 | -0.380050 | ... | -0.621974 | -0.682284 | -0.354744 | 0.278225 | 0.464447 | -0.914285 | -1.192041 | -0.408850 | -0.482241 | 1.149974 | |
5 rows × 24 columns
Se definió como "positivo" a aquellos individuos con capacidad cognitiva reducida. Como se verá adelante, esta decisión influencia sobre la efectividad de algunos de los métodos con los cuales se experimenta.
# Agrego 1 y 0 para indicar si son o no casos de capacidad cognitiva reducida
df_features["valor_real"] = 0
df_features.loc["P","valor_real"] = 1c_bandas = ["alpha", "beta", "delta", "gamma", "theta"]
c_bandas_norm = [b+"_norm" for b in c_bandas]
c_entropia = ["inter", "intra"]
c_bandas = [(a,b) for a in ["media", "std"]for b in c_bandas]
c_bandas_norm = [(a,b) for a in ["media", "std"]for b in c_bandas_norm]
c_entropia = [(a,b) for a in ["media", "std"]for b in c_entropia]En esta primera sección se busca analizar la capacidad como clasificador de cada feature por separado. Para ello, dada una feature particular, se plantean dos métodos:
- Utilizar únicamente un umbral tal que serán clasificados como casos positivos aquellas instancias cuyos valores para dicha feature sean mayores al umbral en cuestión.
- Entrenar un modelo de regresión logística para cada feature, prediciendo luego la probabilidad de que cada instancia sea o no positiva. Hecho esto se utilizará nuevamente un umbral para definir cuales instancias son positivas.
def threshold_ROC(expected, predicted, feature_name, ax=None, plot_ROC=True, **plot_args):
fpr, tpr, _ = metrics.roc_curve(expected, predicted)
roc_auc = metrics.auc(fpr, tpr)
if plot_ROC:
ax.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--')
ax.plot(fpr, tpr, lw=4, **plot_args)
ax.set_xlim([0.0, 1.0])
ax.set_ylim([0.0, 1.05])
feature_name = feature_name.replace("estándar ", "estándar\n")
ax.set_title('{name}'.format(name=feature_name, AUC=roc_auc))
return roc_auc
def name_of_feature(c):
return " ".join([c[0], *(c[1].split("_"))]).capitalize().replace("norm", "normalizada").replace("Std", "Desviación estándar")
# COMO METO EL POLIMORFISMO!!1!UNO
def nombre_bonito(s):
return s.capitalize().replace("norm", "normalizada").replace("Std", "Desviación estándar").replace("_", " ")def analisis_univariado(df_, df_metr, axes):
for c in (c_bandas + c_bandas_norm + c_entropia):
feature_name = name_of_feature(c)
plot_args = {
"color": 'darkorange',
"ls": "-."
}
auc = threshold_ROC(
df_["valor_real"].values,
df_[c],
feature_name=name_of_feature(c),
ax=axes[c],
**plot_args
)
df_metr.loc[c,"umbral_auc"] = auc
return axesdef analisis_regresion_logistica(df_, df_metr, axes):
for c in (c_bandas + c_bandas_norm + c_entropia):
clf = linear_model.LogisticRegression()
X = df_[c].values.reshape(df_[c].size,1)
Y = df_["valor_real"]
scores = model_selection.cross_val_score(clf, X, Y)
df_metr.loc[c, "regLog_score"] = sum(scores) / len(scores)
clf.fit(X, Y)
predicted = clf.predict(X)
plot_args = {
"color": "red",
"ls": ":"
}
auc = threshold_ROC(
df_["valor_real"].values,
predicted,
feature_name=name_of_feature(c),
ax=axes[c],
**plot_args
)
df_metr.loc[c, "regLog_auc"] = auc
Para comparar inicialmente la efectividad de ambos métodos se crean las curvas ROC correspondientes a cada feature.
plt.close()
axes = dict()
fig = plt.figure(figsize=(25,25))
for son_bandas, l in [(True, c_bandas+c_bandas_norm), (False, c_entropia)]:
for j, c in enumerate(l):
if son_bandas:
i = (0 if c[0] == "media" else 1) + (2 if j >=10 else 0) + 1
j = j % 5
else:
i = 0
j = j if j < 2 else (j+1)
ax = plt.subplot2grid((5,5), (i,j))
axes[c] = ax
if i == 4:
ax.set_xlabel('FPR')
else:
ax.xaxis.set_major_formatter(NullFormatter())
if j == 0:
ax.set_ylabel('TPR')
else:
ax.yaxis.set_major_formatter(NullFormatter())
df_metricas = pd.DataFrame(
index=pd.MultiIndex.from_arrays(list(zip(*sorted(c_bandas + c_bandas_norm + c_entropia)))))
analisis_univariado(df_features, df_metricas, axes)
analisis_regresion_logistica(df_features, df_metricas, axes)
p_reg_log = mpatches.Patch(color='red')
p_umbral = mpatches.Patch(color='darkorange')
fig.legend(
labels=['Umbral', 'Regresión logística'],
handles=[p_umbral, p_reg_log],
loc="upper center",
bbox_to_anchor=(0.5, 0.88)
)
fig.tight_layout()
fig.subplots_adjust(hspace=0.5)
plt.subplots_adjust(top=0.9)
plt.suptitle("Curvas ROC para features individuales")
plt.show()Para complementar las curvas ROC, se muestra a continuación el valor de AUC para cada una de ellas.
c_feature = "Feature"
df_auc = df_metricas.reset_index()
df_auc[c_feature] = (df_auc["level_0"] + " " + df_auc["level_1"]).apply(nombre_bonito)
c_umbral = "Umbral"
c_reg_log = "Regresión logística"
c_feature = "Feature"
df_auc.rename(columns={"umbral_auc": c_umbral, "regLog_auc": c_reg_log}, inplace=True)
df_auc = df_auc[[c_feature, c_umbral, c_reg_log]].sort_values(c_umbral, ascending=False)
f, ax = plt.subplots(1, 1)
args = {
"data": df_auc,
"ax": ax,
"x": c_feature,
"join": False
}
plt.xticks(rotation=90)
sns.pointplot(y=c_umbral, color='blue', **args)
sns.pointplot(y=c_reg_log, color='green', **args)
ax.legend(handles=ax.lines[::len(df_auc)+1], labels=["Umbral","Regresión logística"])
plt.ylabel("AUC")
ttl = plt.title("Comparación de AUC para ambos métodos y distintas features")
ttl.set_position([.5, 1.05])
plt.show()
La primer observación a realizar es como el método obtiene resultados extremos: por un lado se obtienen altos valores de AUC para Media Delta Normalizada o Desviación Estandar Inter, mientras que con otros features el resultado es completamente distinto, como en la Media Intra. Sin embargo ésto no se debe a que no sean separables los valores de ambos grupos: de hecho, en los casos donde el AUC es prácticamente cero, el problema es que la predicción que se está realizando es completamente inversa a lo que debería ser.
La causa es que el método simplemente utiliza un umbral particular, clasificando luego como negativos a aquellos que estén por debajo de él y como positivos a aquellos que estén por arriba. Si para alguna feature los valores difieren entre ambas clases, pero con la particularidad de que las instancias negativas tienen mayores valores, entonces el método será inútil. Una solución sencilla sería agregar como features al negativo de cada una ya existente y elegir para el método umbral a aquella de las dos que sea conveniente. Si se realizara esta modificación el método umbral no tendría casos donde el AUC sea menor a 0.5.
Así como vemos nuestros resultados de comparación de AUC como con las curvas ROC, se puede notar que el modelo de regresión logística se ajusta bien a nuestros datos, siendo el candidato a elección entre ambos métodos.
Se muestra a continuación el resultado obtenido para los modelos individuales de regresión logística.
c_feature = "Feature"
df_reg_log = df_metricas.reset_index()
df_reg_log[c_feature] = (df_reg_log["level_0"] + " " + df_reg_log["level_1"]).apply(nombre_bonito)
c_metrica = "Métrica"
df_reg_log.rename(columns={"regLog_score": c_metrica}, inplace=True)
df_reg_log = df_reg_log[[c_feature, c_metrica]].sort_values(c_metrica, ascending=False)
plt.xticks(rotation=90)
sns.pointplot(x=c_feature, y=c_metrica, data=df_reg_log, join=False)
plt.ylabel(c_metrica)
ttl = plt.title("Valor de {} para distintas features".format(c_metrica.lower()))
ttl.set_position([.5, 1.05])
plt.show()Podemos ver como se podía esperar por el gráfico anterior que los mejores resultados se obtienen con los features Media Delta Normalizada, Media Beta Normalizada y Desviación Estandar Inter. En el caso de este primer feature, también nos condujo a buenos resultados utilizando umbrales, lo cual nos dice que aporta mucha información en una clasificación. Sería de interés comprobar estos resultados con más sets de prueba.
En esta sección se analizan los resultados obtenidos con clasificadores basados en el algoritmo de Support Vector Machine según dos variantes:
- Utilizar todos los features para entrenar el clasificador.
- Aplicar un pipeline basado en un scaler, feature selection y finalmente la clasificación utilizando SVM
X = df_features.drop(("valor_real", ""), axis=1)
Y = df_features["valor_real"]def analisis_modelo(fitted_clf, X, Y, ax=ax, plot_ROC=True):
classes_prob = fitted_clf.predict_proba(X)
idx_class = np.where(fitted_clf.classes_ == 1)[0][0]
predicted = classes_prob[:,idx_class]
scores = model_selection.cross_val_score(fitted_clf, X, Y)
final_score = sum(scores) / len(scores)
plot_args = {
"color": "black",
"ls": ":"
}
auc = threshold_ROC(Y.values, predicted, "SVM", ax=ax, plot_ROC=plot_ROC, **plot_args)
return auc, final_scorePara comparar los resultados de ambos clasificadores se gráfican las correspondientes curvas ROC y se calcula el score obtenido.
f, axes = plt.subplots(1,2)
# SVM sin pipeline
clf_sin_pipeline = SVC(probability=True)
clf_sin_pipeline.fit(X,Y)
auc_sp, score_sp = analisis_modelo(clf_sin_pipeline, X, Y, ax=axes[0])
# SVM con pipeline
steps = [
("preprocesamiento", MinMaxScaler()),
("selector_features", SelectKBest(k=2)),
("svm", SVC(probability=True))
]
clf_con_pipeline = Pipeline(steps)
clf_con_pipeline.fit(X,Y)
auc_cp, score_cp = analisis_modelo(clf_con_pipeline, X, Y, ax=axes[1])
axes[0].set_title("Sin pipeline")
axes[1].set_title("Con pipeline")
plt.show()
for m, a, s in [("SVM sin pipeline", auc_sp, score_sp), ("SVM con pipeline", auc_cp, score_cp)]:
print("Resultados obtenidos para", m)
print("\tAUC:", "{:0.2f}".format(a))
print("\tScore:", "{:0.2f}".format(s))Resultados obtenidos para SVM sin pipeline
AUC: 1.00
Score: 0.94
Resultados obtenidos para SVM con pipeline
AUC: 1.00
Score: 1.00
Se puede ver que aunque el AUC haya arrajodo el mismo valor, el score final obtenido es mejor utilizando pipeline.
Para este pipeline se utilizó sólo el 10% de los features (2 de 20 en este caso), de modo que puede ser de interés observar los resultados obtenidos utilizanod más. El siguiente gráfico muestra cómo afecta la cantidad de features utlizadas en el score del modelo.
cantidades = list(range(1,21))
df_scores = pd.DataFrame(data={"Score":np.zeros(20)}, index=cantidades)
for mi_k in cantidades:
# SVM con pipeline
steps = [
("preprocesamiento", MinMaxScaler()),
("selector_features", SelectKBest(k=mi_k)),
("svm", SVC(probability=True))
]
clf = Pipeline(steps)
clf.fit(X,Y)
_, score = analisis_modelo(clf, X, Y, plot_ROC=False)
df_scores.loc[mi_k] = score
df_scores = df_scores.reset_index().rename(columns={"index": "Cantidad de features"})
sns.pointplot(data=df_scores, x="Cantidad de features", y="Score", join=False)
plt.ylabel("Score obtenido")
plt.title("Resultado según distintas cantidades de features")
plt.show()Pueden encontrarse varios modelos que obtengan buenos datos en base a nuestro set de datos, tanto univariado como multivariado. Teniendo en cuenta:
- que tenemos pocas muestras.
- los resultados obtenidos.
- la recomendaciónde la cátedra N > 2 d de la cátedra con d la cantidad de features y Nla cantidad de muestras necesarias para entrenar nuestro modelo.
se escoge el último clasficador sobre los demás (basado en SVM con pipeline utilizando el 10% de los features).




