import pandas as pd
import numpy as np
Vamos a trabajar con el dataset del titanic: https://www.kaggle.com/c/titanic/overview
titanic = pd.read_csv('titanic/train.csv')
titanic.head()
Información básica de las columnas:
titanic.info()
Seleccionamos dos columnas para estudiar las odds ratio y las proporciones
dt = titanic[['Survived','Sex']]
dt.head()
Vamos a ver el porcentaje de mujeres y hombres:
dt['Sex'].value_counts()
100*dt['Sex'].value_counts()/len(dt['Sex'])
Podemos hacer un bar plot:
import seaborn as sns
sns.countplot(dt['Sex'])
Podemos sacar una tabla como la de antes:
pd.crosstab(index=dt['Survived'],columns=dt['Sex'], margins=True) # margins para sacar los totales
Podemos hacer la tabla con las proporciones/probabilidades de cada clase:
pd.crosstab(index=dt['Survived'],columns=dt['Sex'], margins=True).apply(lambda r: r*100/len(dt))
Igual pero relativo por cada fila:
pd.crosstab(index=dt['Survived'],columns=dt['Sex']).apply(lambda r: r*100/r.sum(), axis=1)
# De los que sobreviven cuántos son hombres
109/342
Relativo por columnas:
pd.crosstab(index=dt['Survived'],columns=dt['Sex']).apply(lambda r: r*100/r.sum(), axis=0)
# De los hombres, cuántos sobreviven
109/577
Podemos hacer otro barplot:
pd.crosstab(index=dt['Survived'],columns=dt['Sex']).apply(lambda r: r*100/r.sum(), axis=0).plot(kind='bar')
Vamos a calcular odds y proporciones:
p_mujer_vive=233/314
p_mujer_vive
p_mujer_muere=1-p_mujer_vive
p_hombre_vive=109/577
p_hombre_muere=1-p_hombre_vive
# Odds mujer
odds_mujer=p_mujer_vive/p_mujer_muere
odds_mujer
# Odds hombre
odds_hombre=p_hombre_vive/p_hombre_muere
odds_hombre
# Odds ratio
odds_ratio=odds_hombre/odds_mujer
odds_ratio
odds_mujer/odds_hombre
Hay funciones que calculan los odds ratio (pero no entraremos demasiado)
import scipy.stats as stats
table=pd.crosstab(index=dt['Survived'],columns=dt['Sex'])
oddsratio, pvalue =stats.fisher_exact(table)
oddsratio
El modelo:¶
dt.head()
Para transformar la variable Sex en unos y ceros usamos get_dummies:
dt=pd.get_dummies(dt)
dt.head()
Como es redundante tiramos una:
dt=dt.drop('Sex_female',axis=1)
dt.head()
Hagamos la regresión logística:
from sklearn.linear_model import LogisticRegression
help(LogisticRegression)
# Definimos el modelo
logreg=LogisticRegression(random_state=0, solver='lbfgs') #Se puede fijar un solver
# parámetros del modelo
logreg.get_params()
# Definimos X e y
X = dt.drop('Survived',axis=1) # Para que me de un dataframe y no una serie
y=dt['Survived']
# Ajustamos
logreg = logreg.fit(X,y)
Coeficientes del modelo: $\beta$ = coef_, $\alpha$ = intercept_.
print('alpha: ', logreg.intercept_)
print('beta: ', logreg.coef_)
Odds de la mujer: $e^\alpha$.
np.exp(logreg.intercept_)
odds_mujer
Odds del hombre: $e^{\alpha+\beta}$
np.exp(logreg.intercept_+logreg.coef_)
odds_hombre
Odds ratio: $e^\beta$
np.exp(logreg.coef_)
odds_ratio
# Predecimos
y_pred = logreg.predict(X)
Hagamos la confusion matrix:
from sklearn.metrics import confusion_matrix
confusion_matrix(y_pred,y)
La precisión del modelo la podemos medir con su accuracy. Se puede obtener con .score, que nos evalúa unos datos sobre un modelo y calcula el error o con accuracy_score que calcula la accuracy entre un vector de y reales otro de predichas
logreg.score(X,y)
from sklearn.metrics import accuracy_score
accuracy_score(y,y_pred)
Hagamos otro ejemplo con más variables:
dt=titanic[['Survived','Sex','Age','Pclass']]
dt.head()
dt.info()
# Me voy a transformar Pclass en string porque es una variable categórica
dt.Pclass=dt.Pclass.apply(lambda x: str(x))
dt.info()
# Me saco los dummies de las variables categóricas
dt=pd.get_dummies(dt)
dt.head()
# Filtramos las clases que sobran
dt=dt.drop('Pclass_3',axis=1)
dt=dt.drop('Sex_female',axis=1)
dt.head()
Voy a tratar con los NaN
np.sum(dt.isnull())
# Me quito los NaN
dt=dt.dropna()
Hacemos el modelo:
# Seleccionamos X e y
X=dt.drop('Survived', axis=1)
y=dt['Survived']
# Definir:
logreg= LogisticRegression(random_state=0,solver='lbfgs')
# Ajustar:
logreg=logreg.fit(X,y)
# Predecimos
y_pred=logreg.predict(X)
y_pred
Me ha predicho la variable target como unos y ceros, pero me puede interesar ver cómo lo predice como probabilidades. (Básicamente fija un umbral de 0,5 y devuelve 1 si la probabilidad es mayor que el umbral y 0 en caso contrario).
probs=logreg.predict_proba(X)
probs
Podemos sacar la accuracy:
logreg.score(X,y)
Coeficientes:
logreg.coef_, logreg.intercept_
Hay varios coeficientes ($\beta_1, \beta_2, \beta_3, \beta_4$), uno por cada variable y el otro es el $\alpha$.
Por último, hallemos la ROC curve:
from sklearn.metrics import roc_curve, auc
fpr, tpr, threshold = roc_curve(y,probs[:,1])
roc_auc = auc(fpr,tpr)
import matplotlib.pyplot as plt
t=np.arange(0,5,0.2)
plt.title('Receiver Operating Characteristic')
plt.plot(fpr,tpr, 'b', label = 'AUC = %0.2f' % roc_auc)
plt.legend(loc = 'lower right')
plt.plot([0,1],[0,1], 'r--')
plt.xlim([0,1])
plt.ylim([0,1])
plt.ylabel('True Positive Rate')
plt.xlabel('False Positive Rate')
plt.show()
