Hide code cell source 
import holoviews as hv
hv.extension('bokeh')
hv.opts.defaults(hv.opts.Curve(width=500), 
                 hv.opts.Histogram(width=500))
from bokeh.plotting import show
import numpy as np
import scipy.stats

import statsmodels.api as sm
from scipy.stats import shapiro

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

NumPy versión: 1.21.5
Scipy versión: 1.9.1

15. 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.

15.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 \]
\[g(x) = \mathbb{E}[Y|X=x] = \int y f(y|x) dy \]

Esta función mapea un conjunto de variables a una variable:

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

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

El objetivo es estimar \(g\) a patir de datos de la forma:

\[(Y_1, X_1), ..., (Y_N, X_N) \sim F_{X, Y}\]

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.

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 (\(X\) puede ser una variable discreta también)

  • \(Y\) una variable continua unidimensional

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

  • \(g_b\) un modelo con vector de parámetros \(b\)

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

\[ y_i \approx g_b(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

Una función es lineal si tiene estas propiedades:

  • Aditividad: \(f(X_1 + X_2) = f(X_1) + f(X_2)\)

  • Homogeneidad: \(f(aX) = af(X)\)

Advertencia

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

Por ejemplo, el siguiente modelo con entrada unidimensional

\[ g_b \left(x_i \right) = b_0 + b_1 x_i + b_2 x_i^2, \]

es un modelo lineal en \(b\) y

\[ g_b(x_i) = b_0 + b_1 \log(x_i), \]

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

\[ g_b(x_i) = b_0 + \log(x_i + b_1), \]

no es lineal en \(b\).

15.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_b(x_i) = b_0 + b_1 x_i, \]

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

  • \(b_0\) como la intercepta (intercept)

  • \(b_1\) como la pendiente (slope)

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_b(x_i) = b_0 + b_1 x_{i1} + b_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_b(x_i) = b_0 + \sum_{j=1}^D b_j x_{ij}, \]

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

15.2.1. Ajustando el modelo lineal simple: Matemáticamente#

Asumiendo que tenemos \(\{x_i, y_i\}_{i=1,\ldots,N}\) observaciones de las variables unidimensionales \(X\) e \(Y\), y \(X\) se considera dado, ¿cómo encontramos \(b_0\) y \(b_1\) tal que

\[y_i \approx b_0 + b_1 x_i, \forall i \qquad?\]

El valor predicho por el modelo se denomina \(\hat y_i\), y:

\[\hat y_i = g_b(x_i)\]

Asumimos un término de error o residuo que cumple:

\[y_i = \hat y_i + e_i = g_b(x_i) + e_i = b_0 + b_1 x_i + e_i, \forall i\]

donde los \(e_i\) son i.i.d. Cada \(e_i \sim \cal N(0, \sigma^2)\), \(e_i\) es independiente de \(x_i\), y los \(e_i\) son independientes entre sí mismos.

Así que el error o residuo para la i-ésima observaicón es:

\[e_i = y_i - g_b(x_i) = y_i - \hat y_i = y_i - b_0 - b_1 x_i\]

El error or residuo cuadrático para la i-ésima observaicón es:

\[ e_i^2 = (y_i - b_0 - b_1 x_i)^2, \]

La suma de los errores cuadrados para todas observaicones es:

\[ L = \sum_{i=1}^N e_i^2 = \sum_{i=1}^N (y_i - b_0 - b_1 x_i)^2, \]

Luego podemos ajustar (entrenar) un modelo en base al siguiente criterio de minimizar la suma de cuadrados residuales:

\[ \min_{b} L = \min_{b} \sum_{i=1}^N (y_i - b_0 - b_1 x_i)^2, \]

Nota

Desde un punto de vista estadístico, este criterio es razonable si las observaciones de entrenamiento (\(x_i\), \(y_i\)) representan extracciones aleatorias independientes de su población. Pero, aunque las \(x_i\) no se hayan extraído aleatoriamente, el criterio sigue siendo válido si las \(y_i\) son condicionalmente independientes dadas las entradas \(x_i\).

Los \(y_i\) tienen que ser i.i.d.? Los \(x_i\) tienen que ser i.i.d.?

Cómo obtiene los valores de \(b_0\) y \(b_1\)?

Demostración

Si derivamos la expresión \(L\) con respecto a \(b_0\) primero, obtenemos:

\[ \hat b_0 = \bar y - \hat b_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 derivamos la expresión \(L\) con respecto a \(b_1\), obtenemos:

\[ \hat b_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} = \displaystyle \frac{\overline{xy}-\bar{x}\bar{y}}{\overline{x^2}-\bar{x}^2} = \frac{\sum_i (x_i-\bar{x})(y_i-\bar{y})}{\sum_i (x_i-\bar{x})^2}= \frac{\frac{1}{n-1}{\sum_i (x_i-\bar{x})(y_i-\bar{y})}}{\frac{1}{n-1}{\sum_i (x_i-\bar{x})^2}} = \frac{s_{xy}}{s_x^2} \]

donde \(s_{xy}\) es la covarianza muestral insesgada (unbiased sample covariance) de \(x\) y \(y\), y \(s_x^2\) es la varianza muestral insesgada de x.

Puedes resolver las derivaciones arriba?

Hecho interestante Cuando tenemos un modelo de sólo intersección \(b_0\) sin pendientes, \(b_0\) es igual a la media \(\bar{y}\).

15.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 b_1\)

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

  • rvalue: El coeficiente de correlación de Pearson

  • pvalue: Un p-value para la hipótesis nula de que \(b_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)
b, sigma = [0.5 , -1], 0.5
x = np.random.rand(25)*5

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

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

Ajustamos los parámetros a los datos con:

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

print(f"hat b0: {hat_b[0]:0.5f}, hat b1: {hat_b[1]:0.5f}")

hat b0: 0.42758, hat b1: -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_b)
Hide code cell source 
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)
plot = hv.Overlay([p_data, p_fitted])
show(hv.render(plot))

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:

hat_y = model(x, hat_b)
residuals = y - hat_y
bins, edges = np.histogram(residuals, density=True)
Hide code cell source 
p_residuals = hv.Scatter((hat_y, residuals), 'Predicted y', '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)

plot = hv.Layout([p_residuals * p_zero, p_hist]).cols(2)
show(hv.render(plot))

print(shapiro(residuals))

ShapiroResult(statistic=0.9754845499992371, pvalue=0.7837038040161133)

En general buscamos residuos que:

  • Estén concentrados en torno a cero (media 0 en general y en cada \(\hat y_i\))

  • No estén correlacionados con \(\hat y_i\), tengan mismos varianzas en cada \(\hat y_i\)

  • Tenga una distribución normal

Advertencia

Si no se cumplen, 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. Veamos cómo abordarlo en el próximo capítulo.

15.2.3. Coeficiente de determinación#

Podemos medir la fuerza de la relación lineal o el “fitness” del modelo por el coeficiente de determinación o \(R^2\) (o \(r^2\)):

\[ R^2 = 1 - \frac{\sum_i (y_i - \hat y_i)^2}{\sum_i (y_i - \bar y_i)^2} = \frac{\sum_i (\hat y_i^2 - \bar 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 en la muestra (sample correlation) entre \(x\) y \(y\). (Si te interesa saber por qué, tratas de hacer la demostración!)

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

  • Decimos \(100R^2\)% de la variación en \(Y\) se explica por la variación en el predictor \(X\).

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

  • Si \(R^2 = 0\), el modelo ajustado es una linea horizontal. El predictor 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\).

El modelo tiene \(R^2\):

res.rvalue**2

0.8900544184564642