Regresión lineal¶

La regresión lineal es un método de aprendizaje supervisado, donde dado un elemento de entrada $x \in \mathbb{R}^d$ se busca aprender una función $f(x; w) \in \mathbb{R}$ que estime valores reales para cada entrada $x$. Como su nombre lo dice, esta estimación la hace de manera lineal. Es decir, estima una recta, una función lineal. De tal forma que $f: X \to Y$ es lineal. Si los datos tienen una correlación lineal alta, la regresión lineal será un buen aproximador de estos valores. Podemos ver también al valor $f(x;w)$ como el valor esperado para $x$.

A continuación presentamos un ejemplo simple de regresión lineal con datos de la paquetería sklearn que describen casos médicos y su progresión de la enfermedad de diabetes. El objetivo es estimar la función que dada las medidas de un caso regrese la progresión de la enfermedad.

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_diabetes
from sklearn.model_selection import train_test_split
import pandas as pd

Preparación de los datos¶

Realizaremos una regresión lineal sobre datos de ''Diabetes dataset'' que se encuentran disponibles a partir de la paquetería de sklearn. En este caso, el dataset está conformado por vectores en $\mathbb{R}^{10}$ asociados a un valor en $\mathbb{R}$.

Cada dimensión del vector de entrada $x$ representa una variable que describe al caso en cuestión; algunas de las variables son la edad, el índice de masa corporal, la presión arterial, etc. Estas variables pueden revisarse a partir de la misma paquetería.

En primer lugar, cargamos el dataset. Para visualizar mejor el dataset, utilizamos la paquetería pandas, la cual nos permitirá visualizar datos tabulares de manera eficiente, así como realizar cálculos y estadísticas sobre estos datos.

In [2]:
#Carga dataset
data = load_diabetes()
#Visualiza variables
print(data.feature_names)

#Formato de pandas
table_data = pd.DataFrame(data.data, columns=data.feature_names)
#Agregamos la clase a la tabla
table_data['Y'] = data.target

['age', 'sex', 'bmi', 'bp', 's1', 's2', 's3', 's4', 's5', 's6']

La misma paquetería sklearn nos proporciona información sobre el dataset. Con el comando .DESCR podemos visualizar esta información. Nos muestra el número de instancias/datos que contiene el dataset. Asimismo, despliega la descripción de las variables en los vectores. Asimismo, proporciona otra información relevante para el dataset.

In [3]:
#Información del dataset
print(data.DESCR)

.. _diabetes_dataset:

Diabetes dataset
----------------

Ten baseline variables, age, sex, body mass index, average blood
pressure, and six blood serum measurements were obtained for each of n =
442 diabetes patients, as well as the response of interest, a
quantitative measure of disease progression one year after baseline.

**Data Set Characteristics:**

  :Number of Instances: 442

  :Number of Attributes: First 10 columns are numeric predictive values

  :Target: Column 11 is a quantitative measure of disease progression one year after baseline

  :Attribute Information:
      - age     age in years
      - sex
      - bmi     body mass index
      - bp      average blood pressure
      - s1      tc, total serum cholesterol
      - s2      ldl, low-density lipoproteins
      - s3      hdl, high-density lipoproteins
      - s4      tch, total cholesterol / HDL
      - s5      ltg, possibly log of serum triglycerides level
      - s6      glu, blood sugar level

Note: Each of these 10 feature variables have been mean centered and scaled by the standard deviation times the square root of `n_samples` (i.e. the sum of squares of each column totals 1).

Source URL:
https://www4.stat.ncsu.edu/~boos/var.select/diabetes.html

For more information see:
Bradley Efron, Trevor Hastie, Iain Johnstone and Robert Tibshirani (2004) "Least Angle Regression," Annals of Statistics (with discussion), 407-499.
(https://web.stanford.edu/~hastie/Papers/LARS/LeastAngle_2002.pdf)

Estadísticas y visualización de los datos¶

Una parte importante del aprendizaje automático es su relación con el análisis estadístico. Por tanto, antes de pasar a aplicar el algoritmo para estimación de la regresión lineal podemos computar algunas estadísticas. Una de la más importantes en este caso, en tanto modelo lineal, es la correlación lineal. Con la tabla de pandas que hemos creado, podemos definir una correlación entre las variables. Ya que buscamos establecer una regresión lineal, los elementos que muestren una mayor correlación con la variable $Y$ (progresión de la enfermedad) serán los que nos servirán mejor.

La correlación de Pearson (lineal) se estima de la siguiente forma:

$$\rho(X,Y) = \frac{Cov(X,Y)}{\sigma_X \sigma_Y}$$

Donde $Cov(X,Y)$ es la covarianza entre las variables $X$ e $Y$, la cuál se calcula como $Cov(X,Y) = \mathbb{E}\big[ (X-\mu_x)(Y-\mu_y) \big]$ donde $\mu_i$ son las medias de las variables correspondientes. Asimismo, $\sigma_i$ es la desviación estándar de las variables.

In [4]:
#Correlaciones lineales
table_data.corr()
Out[4]:
age sex bmi bp s1 s2 s3 s4 s5 s6 Y
age 1.000000 0.173737 0.185085 0.335428 0.260061 0.219243 -0.075181 0.203841 0.270774 0.301731 0.187889
sex 0.173737 1.000000 0.088161 0.241010 0.035277 0.142637 -0.379090 0.332115 0.149916 0.208133 0.043062
bmi 0.185085 0.088161 1.000000 0.395411 0.249777 0.261170 -0.366811 0.413807 0.446157 0.388680 0.586450
bp 0.335428 0.241010 0.395411 1.000000 0.242464 0.185548 -0.178762 0.257650 0.393480 0.390430 0.441482
s1 0.260061 0.035277 0.249777 0.242464 1.000000 0.896663 0.051519 0.542207 0.515503 0.325717 0.212022
s2 0.219243 0.142637 0.261170 0.185548 0.896663 1.000000 -0.196455 0.659817 0.318357 0.290600 0.174054
s3 -0.075181 -0.379090 -0.366811 -0.178762 0.051519 -0.196455 1.000000 -0.738493 -0.398577 -0.273697 -0.394789
s4 0.203841 0.332115 0.413807 0.257650 0.542207 0.659817 -0.738493 1.000000 0.617859 0.417212 0.430453
s5 0.270774 0.149916 0.446157 0.393480 0.515503 0.318357 -0.398577 0.617859 1.000000 0.464669 0.565883
s6 0.301731 0.208133 0.388680 0.390430 0.325717 0.290600 -0.273697 0.417212 0.464669 1.000000 0.382483
Y 0.187889 0.043062 0.586450 0.441482 0.212022 0.174054 -0.394789 0.430453 0.565883 0.382483 1.000000

Como podemos observar, la variable $BMI$ (índice de masa corporal) muestra una correlación positiva alta con la variable de salida $Y$; por tanto, puede servirnos para realizar la correlación lineal. De igual manera pudimos tomar la variable $s_5$ (el nivel de triglicéridos), ya que también muestra una correlación positiva considerable. Tomaremos, entonces, esta variable como nuestra variable de entrada. Podemos visualizar entonces esta variable con respecto a la salida $Y$ para observar la correlación.

In [5]:
#Valores de X para BMI
X_BMI = table_data[['bmi']].to_numpy()
#Valores de salida Y
Y = data.target

#Visualización
plt.scatter(X_BMI, Y)
plt.title('Dataset completo')
plt.ylabel('Y')
plt.xlabel('BMI')
plt.show()
No description has been provided for this image

Separación de los datos¶

Una parte importante en toda metodología de aprendizaje automático es la separación de los datos en entrenamiento y evaluación (otras veces también se usa evaluación, pero no es el caso). El esquema que usamos para separación de datos es 70-30; es decir, tomamos 70% para entrenamiento y 30% para evaluación.

Usamos la función train_test_split de la paquetería sklearn que separa los datos aleatoriamente. Hacer la separación de manera aleatoria es importante, pues si los datos están ordenados bajo ciertos criterios, estos criterios pueden introducir sesgos en nuestro aprendizaje.

In [6]:
#Separación de los datos
X_train, X_test, Y_train, Y_test = train_test_split(
    X_BMI, Y, test_size=0.3, random_state=123
)
#Número de datos
print('Datos de entrenamiento: {}\nDatos de evaluación: {}'.format(len(X_train), len(X_test)))

Datos de entrenamiento: 309
Datos de evaluación: 133

Hasta ahora estamos considerando regresión lineal unidimensional. La entrada es un real, por lo que podemos visualizarlo. Visualizamos tanto los datos de entrenamiento como los de evaluación.

In [7]:
#Visualización datos de entrenamiento
plt.scatter(X_train,Y_train)
plt.title('Train dataset')
plt.show()

#Visualización datos de evaluación
plt.scatter(X_test,Y_test)
plt.title('Test dataset')
plt.show()
No description has been provided for this image
No description has been provided for this image

Regresión lineal¶

La regresión lineal, como ya lo señalamos, es un modelo de aprendizaje supervisado lineal que estima una función lineal $f: \mathbb{R}^d \to \mathbb{R}$ debe ajustar los datos de entrenamiento, de tal forma que pueda predecir datos nuevo. El modelo de regresión lineal se divide en los siguientes módulos:

  1. Predicción: El modelo de regresión lineal estima un valor real a partir de un elemento de entrada $x$ a partir de la función lineal: $$f(x) = \sum_{i=1}^d w_i x_i + b= wx+b$$ De tal forma que podemos ver los parámetros a aprender como un vector $w \in \mathbb{R}^d$. El valor $b$ conocido como bias o sesgo es un valor de traslación para que la recta o hiperplano que ajusta los datos no siempre sea centrada en el 0.
  2. Aprendizaje: El aprendizaje se determina como un problema de optimización. En tanto que lo que buscamos es que $f(x)$ sea lo más cercano a $y$, la función de riesgo/objetivo que definimos es un error cuadrático: $$R(w) = \frac{1}{2} ||f(X) - Y||^2$$ El mínimo de esta función se puede estimar obteniendo la derivada de la función e igualando a cero. La derivada de la función está dada por: \begin{align} \nabla_w R(w) &= (f(X) - Y)X \\ &= (Xw - Y)X \\ &= X^TXw-X^TY \end{align} De aquí que, cuando $\nabla_w R(w) = 0$, entonces: $$w = (X^T X)^{-1}X^TY$$ Nota: Para simplificar el método, agregamos una columna con 1's a los datos de entrada, de tal forma que $x' = [x;1]$. este 1 servirá para representar el bias, de tal forma que la función $f(x) = wx'$; es decir, el bias será una entrada (la última) del vector w.

La clase que definimos incorpora estas funciones: predicción (predict) y aprendizaje (fit).

In [8]:
class LinearRegression(object):
    """Clase del modelo de regresión lineal"""
    def __init__(self):
        #Parámetros
        self.w = None
        self.b = None
        
    def predict(self, x):
        """Predicción del modelo"""
        f = np.dot(x,self.w) + self.b
        
        return f
        
    def fit(self, x,y):
        """Estimación de los parámetros"""
        #Se concatena una columna de 1's a X
        Xb = np.concatenate((x, np.ones((len(x),1))), axis=1)
        #XX^-1 
        term1 = np.linalg.inv(np.dot(Xb.T,Xb))
        #XY
        term2 = np.dot(Xb.T,y)
        #Vector solución
        w = np.dot(term1,term2)
        #Guarda parámetros
        self.w = w[:-1]
        self.b = w[-1]

Ahora bien, aplicamos el modelo a nuestro problema y lo entrenamos con los datos de entrenamiento para obtener los valores de los parámetros.

In [9]:
#Creación del modelo
model = LinearRegression()
#Entrenamiento del modelo
model.fit(X_train,Y_train)

print('Valor de w: {}\nValor de b: {}'.format(model.w, model.b))

Valor de w: [958.96106016]
Valor de b: 153.34795164225764

La recta que se ajusta los datos está definida por la función $f(x)$ que depende de $w$. Por tanto, podemos visualizar qué tanto esta recta se ajusta a los datos. En este caso, visualizamos tanto con los datos de entrenamiento como con los de evaluación, pero debe recordarse que la evaluación del modelo de regresión sólo debe realizarse sobre datos que no se vieron en el entrenamiento.

In [10]:
#Visualización en datos de entrenamiento
plt.scatter(X_train, Y_train)
plt.plot(X_train, model.predict(X_train), color='r')
plt.title('Regression (train data)')
plt.show()
No description has been provided for this image

Evaluación de los datos¶

Finalmente, podemos ver cómo se comporta el error a través del dataset de entrenamiento y de evaluación. El error que nos interesa reportar (para determinar la capacidad de generalización de nuestro modelo) es el error de evaluación. Ya que estamos usando un modelo de regresión, podemos estimar el error cuadrático y el score $R^2$.

El error cuadrático medio se estima como:

$$MSE = \frac{1}{N}\sum_{i=1}^N (y_i - \hat{y}_i)^2$$

Donde $y_i$ son los valores de regresión reales, mientras que $\hat{y}_i$ son los valores predichos por el modelo. Por otra parte, el score $R^2$ se estima com:

$$R^2 = 1 - \frac{\sum_{i=1}^N (y_i - \hat{y}_i)^2}{\sum_{i=1}^N (y_i - \hat{y}_i)}$$

Recordemos que las métricas de evaluación se hacen sobre los datos de evaluación, no los de entrenamiento. Para evaluar, utilizamos las métricas disponibles en sklearn.

In [11]:
from sklearn.metrics import mean_squared_error, r2_score

#Predicción sobre datos de evaluación
y_pred = model.predict(X_test)
#Obtención de métricas
mse = mean_squared_error(Y_test, y_pred)
r2 = r2_score(Y_test, y_pred)

print('Error cuadrático medio: {}\nScore R2: {}'.format(mse,r2))

Error cuadrático medio: 4113.813814261837
Score R2: 0.3082201463223282

Podemos también visualizar qué tanto se ajusta la recta estimada a los datos de evaluación.

In [12]:
#Visualización en datos de evaluación
plt.scatter(X_test, Y_test)
plt.plot(X_test, model.predict(X_test), color='r')
plt.title('Regression (test data)')
plt.show()
No description has been provided for this image

Regresión lineal multi-variable¶

En el caso anterior sólo hemos tomado la variable $BMI$ como variable de entrada en la regresión (esto permite visualizar la regresión en un plano 2-dimensional). Pero la regresión lineal puede estimarse tomando en cuenta $d$ variables de entrada. Así, podemos tomar las 10 variables de los datos y ver cuáles son las que muestran mayor influencia en la decisión del precio. En este caso, nuestra función de regresión es una función $f:\mathbb{R}^d \to \mathbb{R}$.

En realidad, no hay mayores cambios al modelo, pues este puede trabajar con entradas de diferentes dimensiones. Sólo que ahora los parámetros aprendidos serán 10 (o en general $d$) en lugar de sólo 1 (sin considerar el sesgo $b$). El procedimiento es el mismo. En primer lugar, separamos los datos en entrenamiento y evaluación.

In [13]:
# Separa los datos multi-dimensionales
X_train, X_test, Y_train, Y_test = train_test_split(
    data.data, Y, test_size=0.3, random_state=7772
)

print("Ejemplo de vector de entrada:\n{}".format(X_train[0]))

Ejemplo de vector de entrada:
[ 0.04170844  0.05068012  0.06169621 -0.04009893 -0.01395254  0.00620169
 -0.02867429 -0.00259226 -0.01495969  0.01134862]

Podemos entonces pasar a crear el modelo y a entrenarlo con los datos de entrenamiento:

In [14]:
#Crea el modelo
model_multi = LinearRegression()
#Entrena el modelo
model_multi.fit(X_train, Y_train)

En este caso, los parámetros aprendidos están $w\in \mathbb{R}^d$. Los valores aprendidos para cada parámetro $w_i$ representan la relevancia que la variable correspondiente $X_i$ tiene para la predicción del valor. De hecho, se puede ver una relación con la correlación de Pearson, si bien no son iguales, ambas representan una relación lineal entre la variable de entrada y la variable de salida. Nótese que esto también depende de los datos que se hayan tomado para el conjunto de entrenamiento, ya que eso se hizo de manera aleatoria.

In [15]:
print('Sesgo: {}'.format(model_multi.b))
pd.DataFrame(data=model_multi.w, index=data.feature_names, columns=['wi']).T

Sesgo: 156.14571327447828
Out[15]:
age sex bmi bp s1 s2 s3 s4 s5 s6
wi -1.141976 -218.55122 479.625279 252.966599 -107.772936 -66.490764 -331.936948 12.582841 535.41564 102.508069

Como era esperado, los valores positivos más grandes corresponden a las variable $BMI$ y $s_5$, pues estas variables, como señalamos más arriba, son las que muestran una mayor correlación lineal con la variable de salida $Y$. De las variables con mayor valor negativo tenemos $s_3$ (prevalencia de lipoproteinas de densidad alta), lo que indica que esta variable contribuye a disminuir la variable $Y$, es decir, contribuye a disminuir el progreso de la diabetes.

En el caso multi-dimensional, lo que aprende la regresión lineal es un hiperplano (generalización de recta y plano) que pasa por los puntos. En más de dos dimensiones esto es difícil de visualizar, pues se requiere de un espacio $d+1$ dimensional. Para 2 dimensiones podemos visualizar los datos en un espacio 3D.

Por ejemplo, si visualizamos las variables con los parámetros con el valor positivo más alto, podemos ver que los puntos tienden a estar en un hiperplano que se inclina hacia arriba cuando los valores de las dos variables crecen.

In [16]:
from mpl_toolkits.mplot3d import Axes3D

#Visualización de dos variables
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.view_init(azim=240)
ax.scatter(table_data[['bmi']].to_numpy(), table_data[['s5']].to_numpy(), Y)
plt.title('Variables con correlación positiva')
plt.xlabel('bmi')
plt.ylabel('s5')
ax.set_zlabel('Y')
plt.show()
No description has been provided for this image

Por su parte, si visualizamos las variables que muestran un valor negativo alto; entonces los puntos recaen en un hiperplano que se inclina hacia abajo cuando los valores de las variables crecen.

In [17]:
#Visualización de dos variables
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.view_init(azim=60)
ax.scatter(table_data[['s3']].to_numpy(), table_data[['s1']].to_numpy(), Y)
plt.title('Variables con correlación negativa')
plt.xlabel('s3')
plt.ylabel('s1')
ax.set_zlabel('Y')
plt.show()
No description has been provided for this image

Evaluación¶

Finalmente, pasamos a evaluar el modelo de regresión lineal multi-dimensional usando las misma métricas que para el modelo uni-dimensional; esto es, usamos el error cuadrático medio y el score $R^2$.

In [18]:
#Predicción sobre datos de evaluación
y_pred = model_multi.predict(X_test)

#Métricas
MSE = mean_squared_error(Y_test, y_pred)
R2 = r2_score(Y_test, y_pred)

print('Error cuadrático: {}\nScore R2: {}'.format(MSE, R2))

Error cuadrático: 2976.8765392178
Score R2: 0.5189382469731116