```python
%matplotlib notebook
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import scipy.optimize
```

# Diagnóstico de cancer usando un regresor logístico

Considere el dataset de [diagnóstico de cancer de mama de la Universidad de Wisconsin](https://archive.ics.uci.edu/ml/datasets/Breast+Cancer+Wisconsin+(Diagnostic))


```python
df = pd.read_csv('cancer.csv', index_col=0)
df['diagnosis'] = df['diagnosis'].map({'M': 1, 'B': 0})
display(df.head())
```

El objetivo de esta actividad formativa es entrenar un regresor logístico para clasificar las muestras como benignas y malignas (columna *diagnosis*) en base a las demás columnas disponibles. Las columnas corresponden al valor medio (1), desviación estándar (2) y "peor caso" (3) de 10 descriptores de la forma de los nucleos celulares presentes en la biopsia del paciente. Cada fila en la tabla corresponde a un paciente distinto.

Revisando este ejercicio aprenderas a

1. formular una regresión logística 
1. estimar los mejores parámetros de un regresor logístico en base a datos
1. implementar un regresor logístico con `numpy` y `scipy`
1. analizar los resultados del regresor logístico

## Modelo de regresión logística

El regresor logístico es la extensión del regresor lineal para problemas de clasificación binaria, es decir donde los datos se agrupan en dos clases

Sean $N$ tuplas $\{\vec x_i, y_i\}$ donde $\vec x_i \in \mathbb{R}^D$ y $y_i \in \{0, 1\}$. El modelo de **regresión logística** para estos datos busca ajustar

$$
y_i \approx \mathcal{S} \left(\theta_0 + \sum_{j=1}^D \theta_j x_{ij}\right),
$$

donde $\mathcal{S}(z) = \frac{1}{1+\exp(-z)}$ se conoce como función logística o sigmoide



```python
z = np.arange(-10, 10, step=0.1)
sigmoide = lambda z : 1/(1+np.exp(-z))

fig, ax = plt.subplots(figsize=(7, 3), tight_layout=True)
ax.plot(z, sigmoide(z));
```

## Función de costo para el regresor logístico

Definamos primero el regresor lineal como 

$$
f_\theta(x_i) = \theta_0 + \sum_{j=1}^D \theta_j x_{ij}
$$

en base a esta definición el modelo de regresión logística es entonces $\mathcal{S} (f_\theta(x_i))$

Para optimizarlo usaremos como función de costo el **logaritmo de la verosimilitud** (ver lección estadística inferencial para más detalles)

$$
\hat \theta = \text{arg} \max_\theta \log \mathcal{L}(\theta)
$$

Asumiremos que las observaciones $y_i$ son independientes y que siguen una [distribución de Bernoulli](https://en.wikipedia.org/wiki/Bernoulli_distribution) con parámetro $p=\mathcal{S} (f_\theta(x_i))$. Luego podemos escribir el logaritmo de la verosimilitud como

$$
\begin{align}
\log \mathcal{L}(\theta) &= \sum_{i=1}^N \log p(y_i | \mathcal{S} (f_\theta(x_i)) ) \nonumber \\
&= \sum_{i=1}^N \log \mathcal{S} (f_\theta(x_i))^{y_i} (1-\mathcal{S} (f_\theta(x_i)) )^{1-y_i} \nonumber \\
&= \sum_{i=1}^N y_i \log (\mathcal{S} (f_\theta(x_i)) ) + (1-y_i) \log(1-\mathcal{S} (f_\theta(x_i)) ) \nonumber
\end{align}
$$

donde la primera igualdad viene de reemplazar la función de masa de la distribución de Bernoulli y la segundo de aplicar el logaritmo

Este modelo no admite una solución analítica cerrada, por lo tanto usaremos un método de optimización iterativo que tome el logaritmo de la verosimilitud como función de costo. En particular usaremos un algoritmo basado en gradiente descendente por lo que necesitamos calcular la derivada de la función de costo

La derivada de la log verosimilitud con respecto a $\theta_j$ es

$$
\frac{d}{d\theta_j} \log \mathcal{L}(\theta) = \sum_{i=1}^N \left ( \frac{y_i}{\mathcal{S} (f_\theta(x_i))} - \frac{1-y_i}{1-\mathcal{S} (f_\theta(x_i))} \right) \frac{d \mathcal{S} (f(x_i, \theta))}{d f_\theta(x_i)} \frac{d f_\theta(x_i)}{d \theta_j} 
$$

donde

$$
\frac{d \mathcal{S}(z)}{dz} =  \mathcal{S}(z) ( 1 - \mathcal{S}(z))
$$

entonces

$$
\begin{align}
\frac{d}{d\theta_j} \log \mathcal{L}(\theta) &= \sum_{i=1}^N \left[ y_i - y_i\mathcal{S} (f_\theta(x_i)) - \mathcal{S} (f_\theta(x_i))  + y_i \mathcal{S} (f_\theta(x_i)) \right] \frac{d f_\theta(x_i)}{d \theta_j} \nonumber \\
&= \sum_{i=1}^N \left[ y_i  - \mathcal{S} (f_\theta(x_i)) \right] \frac{d f_\theta(x_i)}{d \theta_j} \nonumber 
\end{align}
$$

La derivada con respecto a los parámetros es 

$$
\frac{d f(x_i, \theta)}{d \theta_0} =  1
$$

y para $j>0$
$$
\frac{d f(x_i, \theta)}{d \theta_j} =  x_{ij}
$$

Con esto tenemos todo lo necesario para implementar el modelo

## Implementación y optimización del modelo 

Primero implementamos el logaritmo de la verosimilitud y su gradiente

Lo hacemos de tal forma que sean compatibles con `scipy.optimize.minimize`


```python
# Modelo
def modelo(theta, X):
    f = theta[0] + np.sum(theta[1:]*X, axis=1)     
    return sigmoide(f)

# Función de costo 
def neglogverosimilitud(theta, *args):
    X, Y = args
    S = modelo(theta, X)
    # Le agregamos un signo menos ya que quiero minimizar en lugar de maximizar
    return -np.sum(Y*np.log(S+1e-10) + (1-Y)*np.log(1-S+1e-10), axis=0)
    # La función logaddexp es más estable numericamente:
    #return -np.sum(-np.logaddexp(0, -f) - (1-y)*f, axis=0)
    
# Gradiente 
def grad_neglogverosimilitud(theta, *args):
    X, Y = args
    N = len(Y)
    S = modelo(theta, X)
    X1 = np.concatenate((np.ones(shape=(N, 1)), X), axis=1)
    e = (Y - S)
    return -np.sum(e[:, np.newaxis]*X1, axis=0)
```

El ajuste del modelo se muestra a continuación. 

Para facilitar la convergencia del modelo estándarizamos los atributos restándole su media y dividiendo por su desviación estándar

Aprovechando que tenemos la información del gradiente usaremos BFGS para optimizar. BFGS es un método quasi-Newton que usa información de las primeras derivadas

Para validar que el modelo usaremos un conjunto *hold-out*. Otras técnicas para validar son *k-fold* y *bootstrap*


```python
# Obtener base de datos
Y = df["diagnosis"].values
X = df.drop(columns=["diagnosis"]).values
# Estandarizar
X_std = (X - np.mean(X, axis=0))/np.std(X, axis=0)

N = len(Y)
D = X_std.shape[1]
train_proportion = 0.75
idx = np.random.permutation(len(X))
train_idx, test_idx = idx[:int(N*train_proportion)], idx[int(N*train_proportion):]

# Solución inicial
theta_init = 0.1*np.random.randn(1+D)
# Usaremos un callback para guardar el mejor modelo de validación
best_theta = np.zeros_like(theta_init)
best_logl = np.inf
def eval_model(theta):  
    global best_theta, best_logl
    logltrain = neglogverosimilitud(theta, *(X_std[train_idx, :], Y[train_idx]))/len(train_idx)
    logltest = neglogverosimilitud(theta, *(X_std[test_idx, :], Y[test_idx]))/len(test_idx)
    print("Train: %0.4f, Test: %0.4f" %(logltrain, logltest))   
    if logltest < best_logl: # Guardar el mejor modelo de test
        best_theta = theta
        best_logl = logltest

# Mejor valor de theta
res = scipy.optimize.minimize(fun=neglogverosimilitud, x0=theta_init, 
                              method='BFGS', jac=grad_neglogverosimilitud, 
                              args=(X_std[train_idx, :], Y[train_idx]),
                              callback=eval_model, tol=1e-1)

print(res.message)
```

## Evaluación del modelo de clasificación

La salida de este clasificador es un número en el rango $[0, 1]$


```python
p = modelo(best_theta, X_std)

fig, ax = plt.subplots(figsize=(7, 3), tight_layout=True)
ax.hist(p[Y==1], label='Ejemplos con etiqueta maligno', bins=10, range=(0, 1), alpha=0.5, color='r')
ax.hist(p[Y==0], label='Ejemplos con etiqueta benigno', bins=10, range=(0, 1), alpha=0.5, color='b')
ax.set_xlabel('Predicción del modelo')
#ax.set_yscale('log')
ax.legend();
```

Para tomar un decisión se debe seleccionar un umbral $\mathcal{T} \in [0,1]$ tal que

$$
d_i = 
\begin{cases} 
\text{maligno/1}, & \text{si } p(y_i|\theta, \vec x_i) < \mathcal{T} \\ 
\text{benigno/0}, & \text{si } p(y_i|\theta, \vec x_i) \geq \mathcal{T}
\end{cases}
$$

Una vez seleccionado el umbral se puede contar la cantidad de 
- True positives (TP): Tumores malignos clasificados como malignos
- True negative (TN): Tumores benignos clasificados como benignos
- False positives (FP): Tumores benignos clasificados como malignos
- False negative (FN): Tumores malignos clasificados como benignos 

Estas métricas son la base para construir una "tabla o matriz de confusión" para el clasificador

|Clasificado como/En realidad era|Positivo/1:|Negativo/0:|
|---|---|---|
|Positivo/1:|TP | FP |
|Negativo/0:| FN | TN |

Por ejemplo si usamos $\mathcal{T} = 0.5$ tendríamos


```python
def matriz_confusion(Y_real, Y_pred):
    TP = sum(Y_real & Y_pred)
    FP = sum(~Y_real & Y_pred)
    FN = sum(Y_real & ~Y_pred)
    TN = sum(~Y_real & ~Y_pred)
    return np.array([[TP, FP], [FN, TN]])


C = matriz_confusion(Y[test_idx], modelo(res.x, X_std[test_idx]) > 0.5)
display(C)
```

En la próxima unidad veremos formas más sofisticadas para evaluar modelos de clasifición


```python

```

120552decc0a8b22918fa22efbe4effb082a8307
ipynb
Jupyter Notebook
actividades/unidad2/cancer_diagnosis.ipynb
magister-informatica-uach/INFO147
3898eb6f589a22beefb5972a0c911bb9dd098c6d
Unlicense
2019-04-12T21:10:39.000Z
2021-10-12T14:30:09.000Z
actividades/unidad2/cancer_diagnosis.ipynb
magister-informatica-uach/INFO147
3898eb6f589a22beefb5972a0c911bb9dd098c6d
Unlicense
actividades/unidad2/cancer_diagnosis.ipynb
magister-informatica-uach/INFO147
3898eb6f589a22beefb5972a0c911bb9dd098c6d
Unlicense
2019-04-12T20:00:09.000Z
2021-06-17T21:48:53.000Z
Qwen/Qwen-72B
1. YES
2. YES