Modelo predictivo para la selección de medicamentos anti hipertensivos

Barcelona 2019

Científico de datos: Dr Juan Ignacio barrios Arce.

En este primer bloque vamos a importar todas las librerías que necesitamos para procesar nuestros datos. Destacan la librería "pandas" para de análisis de datos, sklearn para el modelo de regresión y las librerías de diagramación o ploteo

import pandas as pd
import numpy as np
from sklearn import linear_model
from sklearn import model_selection
from sklearn.metrics import classification_report
from sklearn.metrics import confusion_matrix
from sklearn.metrics import accuracy_score
import matplotlib.pyplot as plt
import seaborn as sb
%matplotlib inline

El paso siguiente es importar el set de datos en este caso un archivo tipo .csv . Este set de datos será nuestro dataframe que luego con el comando dataframe.head() nos presentará un primer bloque con 5 casos.

dataframe = pd.read_csv(r"antihipertensivos.csv")
dataframe.head(15)

	Medicamento	Edad	Sexo	Etnia	IMC	Diabetes	Dislipidemias	Hist_HTA	Tabaquismo	Alcoholismo	Presion_inicial	Presion_actual	Infartos	ECG
0	2	4	1	2	3	2	2	1	3	1	1	2	2	4
1	3	3	2	3	4	1	3	1	1	5	3	1	2	4
2	4	4	2	2	4	1	1	1	5	1	3	1	2	1
3	3	4	2	1	2	3	1	1	2	1	3	3	2	2
4	4	4	2	2	1	1	2	3	1	5	1	3	2	5
5	1	4	2	2	1	3	1	3	4	4	1	1	2	3
6	2	3	2	2	2	3	1	1	4	1	3	1	3	4
7	1	4	2	2	2	3	2	2	2	5	3	1	3	3
8	4	3	2	2	2	1	1	1	5	3	3	4	2	5
9	5	4	2	1	2	1	1	1	3	5	3	3	3	5
10	2	4	2	3	4	1	2	1	5	4	3	3	2	5
11	4	2	2	2	1	1	3	3	4	1	3	2	3	2
12	3	1	2	1	2	1	1	1	5	4	3	2	2	5
13	4	4	2	2	1	1	2	1	5	1	3	1	1	2
14	2	4	2	3	1	2	2	1	5	5	3	3	2	2

dataframe.tail(5)

	Medicamento	Edad	Sexo	Etnia	IMC	Diabetes	Dislipidemias	Hist_HTA	Tabaquismo	Alcoholismo	Presion_inicial	Presion_actual	Infartos	ECG
49995	3	4	2	1	3	2	3	1	4	1	3	1	2	3
49996	1	3	2	2	2	2	2	3	4	3	3	2	2	2
49997	2	4	2	2	2	1	2	1	2	5	2	4	3	5
49998	3	3	1	1	3	3	2	2	4	4	3	1	3	2
49999	4	4	2	2	1	2	2	3	2	1	3	2	1	4

Ahora necesitamos conocer el tipo de atributos que tenemos. En este caso hay tres atributos tipo int (numéricos ) y 4 atributos tipo Objeto. El comando dataframe.dtypes nos muestra todos los atributos con su tipo.

dataframe.dtypes

Medicamento        int64
Edad               int64
Sexo               int64
Etnia              int64
IMC                int64
Diabetes           int64
Dislipidemias      int64
Hist_HTA           int64
Tabaquismo         int64
Alcoholismo        int64
Presion_inicial    int64
Presion_actual     int64
Infartos           int64
ECG                int64
dtype: object

En esta primera vista de los datos, vemos que NO hay campos tipo Object. Esto es porque a la hora de incluir los datos no dejaron espacios. Estos espacios los tenemos que convertir a valores cero y luego convertir el atributo a tipo numerico "o sea tipo int" para poder operar con ellos. Por último el comando dataframe.describe() nos muestra los estadísticos descriptivos de todos los atributos.

No hay que transformar ningun dato , se adjuntan ejemplos en caso necesario de que los hubiesen

# dataframe=dataframe.replace(r'\s+', 0, regex=True)
# dataframe['Edad']=dataframe['Edad'].astype(int)
# dataframe['Sexo']=dataframe['Sexo'].astype(int)
# dataframe['Etnia']=dataframe['Etnia'].astype(int)
# dataframe['IMC']=dataframe['IMC'].astype(int)

dataframe.describe()

	Medicamento	Edad	Sexo	Etnia	IMC	Diabetes	Dislipidemias	Hist_HTA	Tabaquismo	Alcoholismo	Presion_inicial	Presion_actual	Infartos	ECG
count	50000.000000	50000.000000	50000.000000	50000.000000	50000.000000	50000.00000	50000.000000	50000.000000	50000.000000	50000.000000	50000.000000	50000.000000	50000.000000	50000.000000
mean	3.254280	3.172120	1.842760	1.792100	1.799600	1.92316	1.918420	1.694700	3.712080	3.086960	2.526720	2.031620	2.177780	2.877700
std	1.391151	0.966951	0.364031	0.707565	1.114444	0.81402	0.811208	0.869567	1.191256	1.569741	0.711995	1.126697	0.582427	1.487127
min	1.000000	1.000000	1.000000	1.000000	0.000000	1.00000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000
25%	2.000000	3.000000	2.000000	1.000000	1.000000	1.00000	1.000000	1.000000	3.000000	1.000000	2.000000	1.000000	2.000000	2.000000
50%	3.000000	3.000000	2.000000	2.000000	2.000000	2.00000	2.000000	1.000000	4.000000	3.000000	3.000000	2.000000	2.000000	2.000000
75%	5.000000	4.000000	2.000000	2.000000	3.000000	3.00000	3.000000	3.000000	5.000000	5.000000	3.000000	3.000000	3.000000	5.000000
max	5.000000	4.000000	2.000000	3.000000	4.000000	3.00000	3.000000	3.000000	5.000000	5.000000	3.000000	4.000000	3.000000	5.000000

El paso siguiente es describir la variable que sera nuestra "etiqueta" o variable predictiva.

En este caso definimos que la variable "Medicamento" será la variable predictiva . El comando dataframe.groupby agrupará las diferentes "clases" de esta variable.

Nuestra variable de salida es el tipo de medicamento que la persona esta tomando , hay cinco posibilidades

(1) enalapril 20

(2) enalapril 20 + hidroclorotizida

(3) amlodipidina 5 mg

(4) ibersartan 150

(5) Ibersaratan 150 + hidroclorotiazida

En esta distribución de pacientes es usual que se presente un desbalance entre las llamadas "clases", el cual no es recomendable y debemos utilizar algún método para RE-balancearlas. Si hubiese alguna desproporción de una clase con respecto a las otra, debemos recordarlo para incluirlo en el modelo mas adelante

print(dataframe.groupby('Medicamento').size())

Medicamento
1     7346
2     8341
3    11828
4     9223
5    13262
dtype: int64

Acá sacamos del juego la variable de salida con el comando dataframe.drop

##  Ahora examinamos el resto de las variables en histogramas con el siguiente comando
dataframe.drop(['Medicamento'],1).hist()
plt.show()

Tambien se pueden hacer otro tipos de gráficos que favorezcan el análisis. En este próximo conjunto de gráficos se realiza un cruce de las variables, "todas contra todas". Nótese que cuando se cruza una variable contra si misma, se produce un pico. El comando utilizado es sb.pairplot(dataframe.dropna(),)

# sb.pairplot(dataframe.dropna(), hue='Medicamento',height=4,vars=["Edad", "Sexo","Etnia","IMC","Diabetes","Dislipidemias", "Hist_HTA","Alcoholismo","Tabaquismo","Presion_inicial","Presion_actual","Infartos","ECG"],kind='reg')

En todo análisis donde se utiliza aprendizaje automático es recomendable comparar el modelo utilizando al menos dos algoritmos.

Tratándose de que estamos usando algoritmos de clasificación usaremos primero el algoritmo de Regresión Logística.

Regresión logística:

En primer término vamos a decirle al algoritmo de RL que nuestra variable predictiva es "Salida", o sea el tipo de antihipertensivo que la persona está tomando. El sistema nos muestra el total de la muestra con los atributos o variables (sin considerar la variable predictiva)

# Hacemos Regresión logistica para tener un punto de comparación
X = np.array(dataframe.drop(['Medicamento'],1))
y = np.array(dataframe['Medicamento'])
X.shape

(50000, 13)

Ahora aplicamos el modelo. Si hubiese desbalance de clases, acá debemos considerarlo y decir que el peso entre clases es de: 1:X

model = linear_model.LogisticRegression(class_weight={1:1})
model.fit(X,y)

C:\ProgramData\Anaconda3\lib\site-packages\sklearn\linear_model\logistic.py:432: FutureWarning: Default solver will be changed to 'lbfgs' in 0.22. Specify a solver to silence this warning.
  FutureWarning)
C:\ProgramData\Anaconda3\lib\site-packages\sklearn\linear_model\logistic.py:469: FutureWarning: Default multi_class will be changed to 'auto' in 0.22. Specify the multi_class option to silence this warning.
  "this warning.", FutureWarning)

LogisticRegression(C=1.0, class_weight={1: 1}, dual=False, fit_intercept=True,
                   intercept_scaling=1, l1_ratio=None, max_iter=100,
                   multi_class='warn', n_jobs=None, penalty='l2',
                   random_state=None, solver='warn', tol=0.0001, verbose=0,
                   warm_start=False)

Con el comando model.score obtenemos el Indice de predicción de nuestro modelo. En este caso el modelo tiene un 0.26426 de capacidad predictiva lo que lo hace ser un MAL MODELO. Un buen modelo estaría por arriba del 0.8

model.score(X,y)

0.26606

La matriz de confusión, es la matriz que identifica los valores predictivos para cada clase y a partir de los valores observados de esta matriz se obtienen los falsos negativos y los falsos positivos.

predictions = model.predict(X)
print(confusion_matrix(y, predictions))

[[    0     0   192     0  7154]
 [    0     0   277     0  8064]
 [    0     0   359     0 11469]
 [    0     0   268     0  8955]
 [    0     0   318     0 12944]]

Los valores de precisión (precision) y exactitud (accuracy) y la sensibilidad ( Recall) se presentan en el siguiente cuadro de clasificación

print(classification_report(y, predictions))

              precision    recall  f1-score   support

           1       0.00      0.00      0.00      7346
           2       0.00      0.00      0.00      8341
           3       0.25      0.03      0.05     11828
           4       0.00      0.00      0.00      9223
           5       0.27      0.98      0.42     13262

    accuracy                           0.27     50000
   macro avg       0.10      0.20      0.09     50000
weighted avg       0.13      0.27      0.12     50000

C:\ProgramData\Anaconda3\lib\site-packages\sklearn\metrics\classification.py:1437: UndefinedMetricWarning: Precision and F-score are ill-defined and being set to 0.0 in labels with no predicted samples.
  'precision', 'predicted', average, warn_for)

Opción 2 Algoritmo de Arboles Decisorios

Para este algoritmo es necesario importar alguna librerías adicionales

#arbol
# Imports needed for the script
import numpy as np
import pandas as pd
import seaborn as sb
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams['figure.figsize'] = (16, 9)
plt.style.use('ggplot')
from sklearn import tree
from sklearn.metrics import accuracy_score
from sklearn.model_selection import KFold
from sklearn.model_selection import cross_val_score
from IPython.display import Image as PImage
from subprocess import check_call
from PIL import Image, ImageDraw, ImageFont

De igual forma que como hicimos en el algoritmo anterior aca excluiremos la variable predictiva con el comando dataframe.drop

entradas = dataframe.drop(['Medicamento'], axis=1)
entradas.columns

Index(['Edad', 'Sexo', 'Etnia', 'IMC', 'Diabetes', 'Dislipidemias', 'Hist_HTA',
       'Tabaquismo', 'Alcoholismo', 'Presion_inicial', 'Presion_actual',
       'Infartos', 'ECG'],
      dtype='object')

En este tipo de algoritmo conviene conocer de previo la profundidad idónea del arbol

accuracies = list()
max_attributes = 60
depth_range = range(25, max_attributes + 1)

y_train = dataframe['Medicamento']
x_train = dataframe.drop(['Medicamento',], axis=1).values 


# Testearemos la profundidad del arbol 
for depth in depth_range:
    fold_accuracy = []
    tree_model = tree.DecisionTreeClassifier(criterion='entropy',
                                             min_samples_split=2,
                                             min_samples_leaf=2,
                                             max_depth = depth,
                                             class_weight={1:1})

    model = tree_model.fit(X = x_train, 
                               y = y_train) 
    valid_acc = model.score(X = x_train, 
                                y = y_train) 
    
# calculamos la precision con el segmento de validacion
    fold_accuracy.append(valid_acc)

    avg = sum(fold_accuracy)/len(fold_accuracy)
    accuracies.append(avg)
    
# Mostramos los resultados obtenidos
df = pd.DataFrame({"Profundidad": depth_range, "Precision promedio": accuracies})
df = df[["Profundidad", "Precision promedio"]]
print(df.to_string(index=False))

 Profundidad  Precision promedio
          25             0.70360
          26             0.70372
          27             0.70368
          28             0.70384
          29             0.70426
          30             0.70374
          31             0.70400
          32             0.70396
          33             0.70336
          34             0.70384
          35             0.70338
          36             0.70410
          37             0.70386
          38             0.70398
          39             0.70410
          40             0.70346
          41             0.70322
          42             0.70462
          43             0.70386
          44             0.70376
          45             0.70406
          46             0.70360
          47             0.70360
          48             0.70408
          49             0.70448
          50             0.70380
          51             0.70368
          52             0.70382
          53             0.70416
          54             0.70368
          55             0.70410
          56             0.70392
          57             0.70396
          58             0.70436
          59             0.70342
          60             0.70398

#  OJO hay que ajustar los hiperparámetros. 
y_train = dataframe['Medicamento']
x_train = dataframe.drop(['Medicamento'], axis=1).values 

# Parametros del arbol de decisión
decision_tree = tree.DecisionTreeClassifier(criterion='entropy',
                                            min_samples_split=2,
                                            min_samples_leaf=2,
                                            max_depth = 42,
                                            class_weight={1:1})
decision_tree.fit(x_train, y_train)

# exportar el modelo a archivo .dot
with open(r"tree1.dot", 'w') as f:
     f = tree.export_graphviz(decision_tree,
                              out_file=f,
                              max_depth = 42,
                              impurity = True,
                              feature_names = list(dataframe.drop(['Medicamento'], axis=1)),
                              class_names = ['1: ENLP', '2: ENLP+HCT', '3: AMLD', '4: IBST+HCT', '5: IBST+HCT'],
                              rounded = True,
                              filled= True )
        
# Convertir el archivo .dot a png para poder visualizarlo
check_call(['dot','-Tpng',r'tree1.dot','-o',r'tree1.png'])
PImage("tree1.png")

Ahora calculamos la precisión del modelo. logré este valor haciendo overfitting , coloqué en el árbol los mínimos de división y de hojas por rama en 2 y 2 , con 42 niveles de profundidad.

decision_tree.score(X,y)

0.70378

Acá la matriz de confusión, que indica para cada una de las clases que tan bueno es el algoritmo prediciendo

y_pred = decision_tree.predict(X)
print(confusion_matrix(y, y_pred))

[[6606  149  223  166  202]
 [ 926 6723  225  176  291]
 [1215 1235 8766  223  389]
 [1009 1067 1311 5522  314]
 [1365 1363 1725 1237 7572]]

print(classification_report(y, y_pred))

              precision    recall  f1-score   support

           1       0.59      0.90      0.72      7346
           2       0.64      0.81      0.71      8341
           3       0.72      0.74      0.73     11828
           4       0.75      0.60      0.67      9223
           5       0.86      0.57      0.69     13262

    accuracy                           0.70     50000
   macro avg       0.71      0.72      0.70     50000
weighted avg       0.73      0.70      0.70     50000

Para concluir presentamos el modelo predictivo propiamente

En este caso se pueden colocar directamente las valores : Por ejemplo paciente de 50 años, masculino, blanco, no obeso , no diabético, sin dislipidemias, con historia famiiar de Hipertensión, con tabaquismo, No alcohólico, con presión inicial de 150/100, y presión actual 130/90, sin IM, y con cambios en el ECG

# paciente de 50 años, 3
# masculino, 2
# blanco, 1
# no obeso 1
# no diabético 2
# sin dislipidemias 2
# con historis  familiar de HTA, 2
# No Tabaquismo, 4
# No alcohólismo, 4
# con presión inicial 150/100, 3
# y presión actual 130/90, 1
# sin Isquemias (IM ) 2
# ECG con crec izq.  2

Se colocan los códigos en la matriz de acuerdo al orden del Indice

entradas.columns

Index(['Edad', 'Sexo', 'Etnia', 'IMC', 'Diabetes', 'Dislipidemias', 'Hist_HTA',
       'Tabaquismo', 'Alcoholismo', 'Presion_inicial', 'Presion_actual',
       'Infartos', 'ECG'],
      dtype='object')

X_test = [[3,2,1,1,2,2,2,4,4,3,1,2,2]]
decision_tree.predict(X_test)

array([2], dtype=int64)

#(1) enalapril 20 
#(2) enalapril 20 + hidroclorotizida
#(3) amlodipidina 5 mg
#(4) ibersartan 150 
#(5) Ibersartan 150 + hidroclorotiazida

### El resultado sugiere que al paciente se le debe suministrar la opción 5:  Ibersartan 150 + hidroclorotiazida