Contenido

Regresión Lineal Multivariada

import holoviews as hv
hv.extension('bokeh')
hv.opts.defaults(hv.opts.Curve(width=500),
                 hv.opts.Scatter(width=500, size=4),
                 hv.opts.Histogram(width=500),
                 hv.opts.Slope(color='k', alpha=0.5, line_dash='dashed'),
                 hv.opts.HLine(color='k', alpha=0.5, line_dash='dashed'))   

import pandas as pd
import numpy as np
import scipy.stats
import statsmodels.api as sm

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

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

15. Regresión Lineal Multivariada

A continuación se generalizará el modelo de la lección anterior al caso multivariado, es decir cuando se quiere predecir una variable unidimensional \(Y\) a partir de una variable multidimensional \(X\). Podemos interpretar \(X\) como una tabla donde cada columna representa un atributo o característica en particular.

Veremos primero el formalismo matemático del método de Mínimos Cuadrados Ordinarios y luego como implementarlo en Python.

15.1. Derivación matemática del método de mínimos cuadrados

Consideremos el dataset \(\{x_i, y_i\}_{i=1,\ldots,N}\) con observaciones i.i.d. donde \(y_i \in \mathbb{R}\) y \(x_i \in \mathbb{R}^D\), con \(D>1\). Queremos encontrar \(\theta\) tal que:

\[ y_i \approx \theta_0 + \sum_{j=1}^D \theta_j x_{ij}, \quad \forall i \]

Tal como en la lección anterior, empezamos escribiendo la suma de los errores cuadráticos:

\[ \min_\theta L = \sum_{i=1}^N (y_i - \theta_0 - \sum_{j=1}^D \theta_j x_{ij})^2, \]

pero en este caso lo expresaremos en forma matricial:

\[ \min_\theta L = \| Y - X \theta \|^2 = (Y - X \theta)^T (Y - X \theta), \]

donde

\[\begin{split} X = \begin{pmatrix} 1 & x_{11} & x_{12} & \ldots & x_{1D} \\ 1 & x_{21} & x_{22} & \ldots & x_{2D} \\ 1 & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{N1} & x_{N2} & \ldots & x_{ND} \end{pmatrix}, Y = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_N \end{pmatrix}, \theta = \begin{pmatrix} \theta_0 \\ \theta_1 \\ \vdots \\ \theta_D \end{pmatrix} \end{split}\]

Luego podemos derivar e igualar a cero

\[ \frac{dL}{d\theta} = -(Y - X \theta)^T X = -X^T (Y - X \theta) = 0, \]

para obtener las “ecuaciones normales”:

\[ X^T X \theta = X^T Y \]

cuya solución es:

\[ \hat \theta = (X^T X)^{-1} X^T Y \]

que se conoce como el estimador de mínimos cuadrados ordinario (ordinary least squares, OLS) de \(\theta\)

Relación con la pseudo-inversa de Moore-Penrose

La matriz \(X^{\dagger} = (X^T X)^{-1} X^T \) se conoce como la pseudo-inversa izquierda de Moore-Penrose. También existe la pseudo-inversa por la derecha \(X^T (X X^T)^{-1}\). En conjunto representan una generalización de la inversa para matrices no-cuadradas. En particular, si \(X\) es cuadrada e invertible entonces \(X^{\dagger} = (X^T X)^{-1} X^T = X^{-1} (X^T)^{-1} X^T = X^{-1}\).

Advertencia

La solución de mínimos cuadrados sólo es válida si \(A=X^T X\) es invertible (no-singular). Por construcción \(A \in \mathbb{R}^{D\times D}\) es una matriz cuadrada y simétrica. Para que \(A\) sea invertible se requiere que su determinante no sea cero o de forma equivalente que:

  • El rango de \(A\), i.e. el número de filas o columnas LI, sea igual a \(D\)

  • Los valores propios de \(A\) sea mayores que cero

Nota

La solución que encontramos para el caso univariado es un caso particular de la solución de mínimos cuadrados.

Demostración

La solución del caso univariado era:

\[\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}\]

que puede re-escribirse como:

\[\begin{split} \begin{align} \begin{pmatrix} 1 & 1 & \ldots & 1 \\ x_1 & x_2 & \ldots & x_N \end{pmatrix} \begin{pmatrix} 1 & x_1 \\ 1 & x_2 \\ \vdots & \vdots \\ 1 & x_N \end{pmatrix} \begin{pmatrix} \theta_0 \\ \theta_1 \end{pmatrix} &= \begin{pmatrix} 1 & 1 & \ldots & 1 \\ x_1 & x_2 & \ldots & x_N \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_N \end{pmatrix} \nonumber \\ X^T X \theta &= X^T Y \nonumber \end{align} \end{split}\]

15.2. Perspectiva estadística de mínimos cuadrados

Digamos que se tienen \(\{x_i, y_i\}_{i=1,\ldots,N}\) observaciones i.i.d. de una variable dependiente unidimensional y una variable independiente D-dimensional. Asumiremos que las mediciones \(Y\) corresponden a una combinación aditiva de un modelo y ruido con distribución Gaussiana centrado en cero, es decir:

\[\begin{split} \begin{split} y_i &= f_\theta(x_i) + \varepsilon_i \\ &= \theta_0 + \sum_{j=1}^D \theta_j x_{ij} + \varepsilon_i \end{split} \end{split}\]

donde \(\varepsilon_i \sim \mathcal{N}(0, \sigma^2)\). Con lo anterior la log-verosimilitud de \(\theta\) es:

\[\begin{split} \begin{split} \log L(\theta) &= \log \prod_{i=1}^N \mathcal{N}(y_i | f_\theta(x_i), \sigma^2) \\ &= \sum_{i=1}^N \log \mathcal{N}(y_i | f_\theta(x_i), \sigma^2) \\ &= -\frac{N}{2} \log(2\pi \sigma^2) - \frac{1}{2\sigma^2} \sum_{i=1}^N (y_i - f_\theta(x_i))^2\\ &= -\frac{N}{2} \log(2\pi \sigma^2) - \frac{1}{2\sigma^2} (Y-X\theta)^T (Y - X\theta), \end{split} \end{split}\]

cuya solución de máxima verosimilitud se obtiene de:

\[ \max_\theta \log L(\theta) = - \frac{1}{2\sigma^2} (Y-X\theta)^T (Y - X\theta). \]

Nota

Lo anterior es equivalente a:

\[ \min_\theta \log L(\theta) = \frac{1}{2\sigma^2} (Y-X\theta)^T (Y - X\theta). \]

El estimador de máxima verosimilitud de \(\theta\) es:

\[ \hat \theta = (X^T X)^{-1} X^T Y \]

Importante

La solución de mínimos cuadrados es equivalente a la solución de máxima verosimilitud bajo el supuesto de muestras iid y ruido gaussiano.

15.2.1. Propiedades estadísticas de mínimos cuadrados

Algunas propiedades importantes del estimador de máxima verosimilitud (mínimos cuadrados) de \(\theta\):

  • El estimador es insesgado.

Demostración

Definamos \(\varepsilon = (\varepsilon_1, \varepsilon_2, \ldots, \varepsilon_N)\), donde \(\varepsilon \sim \mathcal{N}(0, I \sigma^2) \quad \forall i\)

Luego:

\[\begin{split} \begin{split} \mathbb{E}[\hat \theta] &= \mathbb{E}[(X^T X)^{-1} X^T Y] \\ &= \mathbb{E}[(X^T X)^{-1} X^T (X \theta + \varepsilon)] \\ &= \mathbb{E}[\theta] + (X^T X)^{-1} X^T \mathbb{E}[\varepsilon] \\ & = \mathbb{E}[\theta] \end{split} \end{split}\]

  • La varianza del estimador es \(\sigma^2 (X^T X)^{-1}\).

Demostración
\[\begin{split} \begin{split} \mathbb{E}[(\hat \theta - \mathbb{E}[\hat\theta])(\hat \theta - \mathbb{E}[\hat\theta])^T] &= \mathbb{E}[((X^T X)^{-1} X^T \varepsilon) ((X^T X)^{-1} X^T \varepsilon)^T] \\ &= (X^T X)^{-1} X^T \mathbb{E}[\varepsilon \varepsilon^T] X ((X^T X)^{-1})^T \\ &= (X^T X)^{-1} X^T \mathbb{E}[(\varepsilon-0) (\varepsilon-0)^T] X (X^T X)^{-1} \\ & =\sigma^2 (X^T X)^{-1} \end{split} \end{split}\]

  • (Teorema de Gauss-Markov) El estimador de \(\theta\) tiene la menor varianza entre todos los estimadores insesgados (Ver Hastie 3.2.2 en la bibliografía del curso).

15.2.2. Inferencia y tests de hipótesis para OLS

Dado que conocemos el valor esperado y la varianza de \(\theta\) podemos utilizar las propiedades de la estimación de máxima verosimilitud para escribir:

\[ \hat \theta \sim \mathcal{N}(\theta, \sigma^2 (X^T X)^{-1}). \]

Por otro lado el estimador de lavarianza es proporcional a:

\[ \hat \sigma^2 \sim \frac{1}{(N-M)}\sigma^2 \chi_{N-M}^2. \]

Con esto tenemos lo necesario para formular un test de hipótesis sobre \(\hat \theta\). Típicamente buscamos evaluar la significancia de nuestro modelo bajo las siguientes hipótesis:

  • (t-test) Uno de los parámetros (pendientes) es cero, es decir:

    \(\mathcal{H}_0: \theta_i = 0\)

    \(\mathcal{H}_A: \theta_i \neq 0\)

  • (f-test) Todos los parámetros son iguales a cero, es decir:

    \(\mathcal{H}_0: \theta_1 = \theta_2 = \ldots = \theta_D = 0\)

    \(\mathcal{H}_A:\) Al menos un parámetro es distinto de cero

  • (ANOVA) un subconjunto de los parámetros son cero, es decir:

    \(\mathcal{H}_0: \theta_i = \theta_j =0 \)

    \(\mathcal{H}_A:\) \(\theta_i \neq 0 \) o \(\theta_j \neq 0 \)

15.3. Mínimos cuadrados ordinarios en Python

A continuación veremos como implementar una regresión lineal y los test estadísticos anteriormente mencionados en lenguaje de programación Python.

Primero consideremos el siguiente dataset de ejemplo, el cual importamos con la librería pandas:

df = pd.read_csv('data/ice_cream.csv', header=0, index_col=0)
df.columns = ['Consumption', 'Income', 'Price', 'Temperature']
display(df.head())
Consumption Income Price Temperature
1 0.386 78 0.270 41
2 0.374 79 0.282 56
3 0.393 81 0.277 63
4 0.425 80 0.280 68
5 0.406 76 0.272 69

donde cada fila corresponde a un día, y las columnas corresponden a:

  • consumo de helados promedio (pintas per capita)

  • ingreso familiar promedio (dolares)

  • temperatura promedio (grados Fahrenheit)

  • precio promedio de los helados (dolares)

Intentaremos responder si:

¿Existe una relación lineal entre el consumo de helados e ingreso, temperatura y precio?

El problema a modelar es el siguiente:

\[ \text{consumption} = \theta_0 + \theta_1 \cdot \text{income} + \theta_2 \cdot \text{price} + \theta_3 \text{temperature} \]

Para ajustar este modelo de regresión lineal utilizaremos la clase OLS del la librería statsmodels.

Primero creamos el objeto OLS entregándole las variables dependiente e independientes. Las variables pueden entregars como DataFrames de pandas o arreglos de numpy. La primera opción es más conveniente pues los resultados entregados por OLS seguiran los nombres y etiquetas del DataFrame.

Y = df["Consumption"]
X = df[["Income", "Price", "Temperature"]]
X = sm.add_constant(X)

model = sm.OLS(Y, X)

Nota

Añadimos un valor constante igual a uno con la función add_constant()sobre variable independiente. De esta forma OLS ajustará el parámetro \(\theta_0\) (intercepta).

El objeto creado tiene un método llamado fit que aplica mínimos cuadrados ordinarios (máxima verosimilitud) para obtener los parámetros del regresor lineal:

result = model.fit()

result.params

const          0.197315
Income         0.003308
Price         -1.044414
Temperature    0.003458
dtype: float64

Para evaluar la calidad podemos visualizar

  • el consumo predicho contra el valor real

  • el residuo contra el valor real

Para predecir utilizamos el método predict del modelo (especificando los parámetros) o del resultado

hatY = result.predict(X)

p1 = hv.Scatter((Y, hatY), 'Real', 'Predicted').opts(width=330) * hv.Slope(slope=1, y_intercept=0)
p2 = hv.Scatter((Y, Y - hatY), 'Real', 'Residuals').opts(width=330) * hv.HLine(0)
hv.Layout([p1, p2]).cols(2)

También podemos visualizar el residuo en función de las variables independientes:

p = []
for var_name in ["Income", "Price", "Temperature"]:
    p.append(hv.Scatter((df[var_name].values, Y - hatY), var_name, 'Residuals').opts(width=330) * hv.HLine(0))
hv.Layout(p).cols(3).opts(hv.opts.Scatter(width=280, height=250))

Aparentemente el valor predicho sigue el valor real bastante bien y no hay correlación en los residuos. Sin embargo quedan algunas preguntas importantes:

  • ¿Qué tan significativa es la contribución de cada variable independiente a la predicción?

  • ¿Cómo podemos medir la calidad del modelo ajustado de forma cuantitativa?


Para esto podemos utilizar los test de hipótesis antes mencionados.

El atributo summary de la estructura retornada por el método fit de OLS contiene:

  • R-squared: el estadístico \(r^2\) del modelo

  • F-statistic: el estadístico f del modelo y su p-value

  • Una tabla con los errores estándar de los parámetros del modelo, estadístico t, intervalos de confianza y p-value para cada uno

Veamos lo que resulta en este caso:

result.summary(alpha=0.05)
OLS Regression Results
Dep. Variable: Consumption R-squared: 0.719
Model: OLS Adj. R-squared: 0.687
Method: Least Squares F-statistic: 22.17
Date: Tue, 08 Nov 2022 Prob (F-statistic): 2.45e-07
Time: 14:58:52 Log-Likelihood: 58.619
No. Observations: 30 AIC: -109.2
Df Residuals: 26 BIC: -103.6
Df Model: 3
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
const 0.1973 0.270 0.730 0.472 -0.358 0.753
Income 0.0033 0.001 2.824 0.009 0.001 0.006
Price -1.0444 0.834 -1.252 0.222 -2.759 0.671
Temperature 0.0035 0.000 7.762 0.000 0.003 0.004
Omnibus: 0.565 Durbin-Watson: 1.021
Prob(Omnibus): 0.754 Jarque-Bera (JB): 0.047
Skew: 0.038 Prob(JB): 0.977
Kurtosis: 3.179 Cond. No. 1.27e+04


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.27e+04. This might indicate that there are
strong multicollinearity or other numerical problems.

Nota

  • En base al estadístico f, se rechaza la hipótesis de que todos los parámetros son nulos.

  • En base al estadístico t, se rechanza la hipótesis de que el parámetro asociado a temperatura e ingreso son nulos. No hay suficiente evidencia para rechazar en el caso del precio y la constante.

Importante

Sólo podemos confiar en los tests si nuestros supuestos eran correctos. En este caso supusimos:

  • La relación entre \(X\) e \(Y\) es lineal

  • Los errores siguen una distribución normal con covarianza \(I\sigma^2\)

Sugerencias:

  1. Si hay problemas numéricos en el ajuste, normalizar los atributos puede ayudar.

  2. Se puede verificar la normalidad de los residuos con los tests no parámetricos que vimos en lecciones anteriores. Si hay outliers deberían removerse antes de hacer el ajuste con OLS.

  3. Si existe correlación en los residuos entonces la relación no es lineal o faltan variables independientes.

  4. Si la varianza de los datos no es constante entonces se debe utilizar otro modelo: Mínimos Cuadrados Ponderados (Weighted Least Squares, WLS). WLS está implementado en statsmodels.

anterior

14. Introducción a la Regresión Lineal

siguiente

16. Regresión lineal con funciones base y Regularización