Contenido

Introducción a la Regresión Lineal

import holoviews as hv
hv.extension('bokeh')
hv.opts.defaults(hv.opts.Curve(width=500), 
                 hv.opts.Histogram(width=500))

import numpy as np
import scipy.stats
print(f"NumPy versión: {np.__version__}")
print(f"Scipy versión: {scipy.__version__}")

NumPy versión: 1.21.2
Scipy versión: 1.7.3

14. Introducción a la Regresión Lineal

Muchos problemas pueden formularse como “encontrar una relación” entre dos o más variables. Lo anterior puede interpretarse como predecir o explicar una variable dadas las demás.

Algunos ejemplos:

  • Predecir ventas dada el dinero invertido en anuncios (ads)

  • Predecir el tiempo de mañana para Valdivia dadas las condiciones atmosféricas actuales e históricas

  • Predecir la probabilidad de padecer cancer de pulmón dado el número de cigarrillados fumados diariamente

De lo anterior surgen las siguientes preguntas relevantes:

  • ¿Están estas variables relacionadas?

  • ¿Qué tan fuerte o significativa es dicha relación?

  • ¿Cuál es la naturaleza de esta relación?

Si respondemos estas preguntas estaremos más cerca a entender los procesos subyacentes tras los datos.

14.1. Definiendo regresión

Regresión se refiere a una familia de modelos estadísticos cuyo fin es encontrar relaciones entre variables. En general esta relación se modela como una función

\[ g: X \to Y, \]

que mapea una variable a un conjunto de variables:

  • La variable de salida \(Y\) se llama variable dependiente, también se suele llamar respuesta u objetivo.

  • La/s variable/s de entrada \(X\) se llama/n variable/s independiente/s. También se suelen llamar características (features).

  • El mapeo o función \(g\) se llama predictor o regresor.

El objetivo es aprender \(g\) tal que podamos predecir \(Y\) dado \(X\), es decir \(\mathbb{E}[Y|X]\)

Nota

La regresión se puede definir desde una perspectiva estadística como un caso particular de ajuste de modelos (estimación de parámetros). También se considera como parte del paradigma de aprendizaje supervisado junto a los modelos de clasificación.

Los métodos de regresión pueden clasificar como parámetricos o no paramétricos. En esta lección nos enfocaremos en los primeros.

Ver también

Los modelos no paramétricos de regresión no tienen estructura predefinida y por ende son más flexibles pero más complejos de ajustar que sus contrapartes paramétricas. Algunos ejemplos son las splines, SVR (support vector regression) y los procesos gaussianos.

Definamos en detalle un modelo paramétrico de regresión. Sea

  • \(X\) una variable continua D-dimensional

  • \(Y\) una variable continua unidimensional

  • \(\{x_i, y_i\}\) con \(i=1,\ldots,N\) un conjunto de \(N\) observaciones de \(X\) e \(Y\)

  • \(g_\theta\) un modelo con vector de parámetros \(\theta\)

El problema de regresión paramétrica consiste en contrar el mejor valor de \(\theta\) tal que

\[ y_i \approx g_\theta(x_i),\quad i=1,\ldots, N \]

El modelo paramétrico más simple es cuando la relación entre variables es lineal (aditivo y homogeneo). En este caso hablamos de regresión lineal.

Nota

La regresión lineal es lineal en \(\theta\) pero no necesariamente en \(X\)

Por ejemplo, el siguiente modelo con entrada unidimensional

\[ g_\theta \left(x_i \right) = \theta_0 + \theta_1 x_i + \theta_2 x_i^2, \]

es un modelo lineal en \(\theta\) y

\[ g_\theta(x_i) = \theta_0 + \theta_1 \log(x_i), \]

es también lineal en \(\theta\), sin embargo

\[ g_\theta(x_i) = \theta_0 + \log(x_i + \theta_1), \]

no es lineal en \(\theta\).

14.2. El modelo lineal más simple

Si consideramos una variable unidimensional \(x_i \in \mathbb{R}, i=1,\ldots,N\), entonces el modelo lineal más simple sería

\[ g_\theta(x_i) = \theta_0 + \theta_1 x_i, \]

que tiene \(M=2\) parámetros. Este modelo corresponde a una linea en \(\mathbb{R}^2\) y podemos reconocer:

  • \(\theta_0\) como la intercepta

  • \(\theta_1\) como la pendiente

Si consideramos una variable de entrada bi-dimensional \(x_i = (x_{i1}, x_{i2}) \in \mathbb{R}^2, i=1,\ldots,N\) entonces tendríamos

\[ g_\theta(x_i) = \theta_0 + \theta_1 x_{i1} + \theta_2 x_{i2}, \]

que tiene \(M=3\) parámetros. Esto corresponde a un plano en \(\mathbb{R}^3\).

La forma más general asume una variable D-dimensional \(x_i = (x_{i1}, x_{i2}, \ldots, x_{iD}), i=1,\ldots,N\) con

\[ g_\theta(x_i) = \theta_0 + \sum_{j=1}^D \theta_j x_{ij}, \]

que tiene \(M=D+1\) parámetros. Esto corresponde a un hiperplano en \(\mathbb{R}^M\).

14.2.1. Ajustando el modelo lineal simple: Matemáticamente

Asumiendo que tenemos \(\{x_i, y_i\}_{i=1,\ldots,N}\) observaciones iid de las variables unidimensionales \(X\) e \(Y\)

¿Cómo encontramos \(\theta_0\) y \(\theta_1\) tal que \(y_i \approx \theta_0 + \theta_1 x_i, \forall i\)?

Empezaremos escribiendo el error o residuo cuadrático como:

\[ e_i^2 = (y_i - \theta_0 - \theta_1 x_i)^2, \]

Luego podemos ajustar (entrenar) un modelo en base al siguiente criterio:

\[ \min_{\theta} L = \sum_{i=1}^N e_i^2 = \sum_{i=1}^N (y_i - \theta_0 - \theta_1 x_i)^2, \]

que se conoce como la suma de los errores cuadrados.

Si derivamos la expresión con respecto a los parámetros obtenemos:

\[ \hat \theta_0 = \bar y - \hat \theta_1 \bar x, \]

donde \(\bar x = \frac{1}{N} \sum_{i=1}^N x_i\), \(\bar y = \frac{1}{N} \sum_{i=1}^N y_i\) y

\[ \hat \theta_1 = \frac{N\sum_i x_i y_i - (\sum_i x_i)(\sum_i y_i)}{N \sum_i x_i^2 - (\sum_i x_i)^2}, \]
Demostración

Con

\[\begin{split} \begin{align} \frac{dL}{d\theta_0} &= -2 \sum_{i=1}^N (y_i - \theta_0 - \theta_1 x_i) \nonumber \\ &= -2 \sum_{i=1}^N y_i + 2 N\theta_0 + 2 \theta_1 \sum_{i=1}^N x_i = 0 \nonumber \end{align} \end{split}\]

y

\[\begin{split} \begin{align} \frac{dL}{d\theta_1} &= -2 \sum_{i=1}^N (y_i - \theta_0 - \theta_1 x_i) x_i \nonumber \\ &= -2 \sum_{i=1}^N y_i x_i + 2 \theta_0 \sum_{i=1}^N x_i + 2 \theta_1 \sum_{i=1}^N x_i^2 = 0 \nonumber \end{align} \end{split}\]

se tiene un sistema de dos ecuaciones y dos incognitas:

\[\begin{split} \begin{pmatrix} N & \sum_i x_i \\ \sum_i x_i & \sum_i x_i^2\\\end{pmatrix} \begin{pmatrix} \theta_0 \\ \theta_1 \end{pmatrix} = \begin{pmatrix} \sum_i y_i \\ \sum_i x_i y_i \end{pmatrix} \end{split}\]

cuya solución es:

\[\begin{split} \begin{pmatrix} \hat \theta_0 \\ \hat \theta_1 \end{pmatrix} = \frac{1}{N\sum_i x_i^2 - \left(\sum_i x_i\right)^2}\begin{pmatrix} \sum_i x_i^2 & -\sum_i x_i \\ -\sum_i x_i & N\\\end{pmatrix} \begin{pmatrix} \sum_i y_i \\ \sum_i x_i y_i \end{pmatrix} \end{split}\]

donde asumimos que el determinante de la matriz no es cero.

14.2.2. Ajustando el modelo lineal simple: Python

Podemos ajustar un modelo de linea en Python usando la librería scipy.stats con

scipy.stats.linregress(x, # Vector con N valores o matriz de Mx2 (si y no se especifica)
                       y=None, # Vector con N valores
                      )

Esta función retorna un objeto cuyos atributos son:

  • slope: Equivalente a \(\hat \theta_1\)

  • intercept: Equivalente a \(\hat \theta_0\)

  • rvalue: El coeficiente de correlación de Pearson

  • pvalue: Un p-value para la hipótesis nula de que \(\theta_1 =0\)

Más adelante ahondaremos en el coeficiente de correlación y los test de hipótesis para la regresión lineal. De momento pongamos en práctica la función anterior utilizando datos sintéticos:

np.random.seed(12345)
theta, sigma = [0.5 , -1], 0.5
x = np.random.rand(25)*5

def model(x, theta):
    return theta[0] + theta[1]*x

y = model(x, theta) + sigma*np.random.randn(len(x))

Ajustamos los parámetros a los datos con:

res = scipy.stats.linregress(x, y)
hat_theta = np.array([res.intercept, res.slope])

print(f"hat theta0: {hat_theta[0]:0.5f}, hat theta1: {hat_theta[1]:0.5f}")

hat theta0: 0.42758, hat theta1: -0.93605

Y luego utilizamos el modelo ajustado para inferir sobre nuevos valores \(\hat x\). La siguiente gráfica muestra los datos utilizados para entrenar y la inferencia con el modelo ajustado en el rango \([-1,6]\).

hat_x = np.linspace(-1, 6, num=50)
hat_y = model(hat_x, hat_theta)

p_fitted = hv.Curve((hat_x, hat_y), label='Fitted model')
p_data = hv.Scatter((x, y), label='Training data').opts(color='k', size=5)

hv.Overlay([p_data, p_fitted])

Se puede apreciar que el modelo sigue bastante de cerca a los datos.

Para evaluar la calidad del ajuste también podemos graficar los residuos, es decir la distancia entre las muestras y la linea ajustada. También podemos inspeccionar el histograma de los residuos.

Para el ejemplo anterior:

residuals = y - model(x, hat_theta)
bins, edges = np.histogram(residuals, density=True)

p_residuals = hv.Scatter((y, residuals), 'Target variable', 'Residuals').opts(color='k', size=5, width=350)
p_zero = hv.HLine(0).opts(color='k', line_dash='dashed', alpha=0.5)
p_hist = hv.Histogram((edges, bins), kdims='Residuals', vdims='Density').opts(width=350)

hv.Layout([p_residuals * p_zero, p_hist]).cols(2)

En general buscamos residuos que:

  • Estén concentrados en torno a cero

  • No estén correlacionados, es decir que se vean como ruido blanco

Advertencia

La existencia de correlación en los residuos es un signo de que el modelo escogido (en este caso la linea) no era adecuada para modelar la relación entre las variables de interés

14.2.3. Coeficiente de determinación

Podemos medir la fuerza de la relación lineal entre las variables utilizando el coeficiente de determinación o \(r^2\) definido como

\[ r^2 = 1 - \frac{\sum_i (y_i - \hat y_i)^2}{\sum_i (y_i - \bar y_i)^2} \in [0, 1], \]

es decir uno menos la suma de los residuos dividido por la varianza de \(y\). Para la regresión lineal simple el estadístico \(r\) es equivalente al coeficiente de correlación de Pearson.

Interpretación de \(r^2\):

  • Si \(r^2 = 1\), los datos son predichos perfectamente por el modelo. El regresor explica toda la variación en \(y\).

  • Si \(r^2 = 0\), el modelo ajustado es una linea horizontal. El regresor no explica la variación en \(y\).

Advertencia

Si las variables tienen una relación fuerte, pero de naturaleza no-lineal, esta podría no ser detectada por \(r^2\).

anterior

13. Análisis de Asociación de Variables Aleatorias

siguiente

15. Regresión Lineal Multivariada