Capítulo 7 Regresion lineal

En machine learning el objetivo principal no es estimar y hacer inferencia como en la econometría clásica sino hacer predicciones/clasificar. El modelo de regresión lineal es una herramienta útil para predecir cuando la variable de respuesta es cuantitativa.

Modelo lineal simple

\[\begin{equation} \tag{7.1} y = \beta_0 + \beta_1 x_1 + u \end{equation}\]

donde \(u\) es un término de error aleatorio que captura todo lo que no puede representarse con este modelo simple (factores no observables, errores de medición, etc.).

Modelo lineal múltiple

En notación matricial:¹⁴

\[\begin{equation} \tag{7.2} Y = X \beta + u \end{equation}\]

donde la primera columna de \(X\) es la constante.

Método de Mínimos Cuadrados Ordinarios \((MCO)\)

\[\begin{equation} \tag{7.3} \text{RSS} = \sum_{i=1}^{n}(y_{i} - \beta_{0} - \sum_{j=1}^{p} \beta_{j} x_{ij} )^2 \end{equation}\]

\[\begin{equation} \tag{7.4} \hat{\beta} = min \sum_{i=1}^{n}e_i^2 \end{equation}\]

\[\begin{equation} \tag{7.5} \hat{\beta} = (X'X)^{-1}X'Y \end{equation}\]

Figura 7.1: Modelo lineal y MCO

Teorema de Gauss/Markov (TGM): Bajo los supuestos clásicos el estimador de \(MCO\) es el mejor dentro de los lineales e insesgados (MELI).

Definición: el \(R^2\) mide la proporción de la variabilidad de \(Y\) explicada por el modelo. \(R^2 \in [0,1]\).

\[\begin{equation} \tag{7.6} R^2 = 1 - \frac{SCR}{SCT} \end{equation}\]

Dado que la suma de cuadrados totales es igual a la suma de cuadrados explicada y la suma de cuadrados residual (\(SCT = SCE + SCR\)) y que \(SCT\) es una magnitud fija, por lo tanto, \(MCO\) maximiza el \(R^2\).

7.1 Relacion entre estimacion optima y prediccion optima

Dado:

\[\begin{equation} \tag{7.7} Y_i = X_i^{'} \beta + u_i \end{equation}\]

con \(i=1,...,n\)

La predicción de \(Y\) se define como:

\[\begin{equation} \tag{7.8} \hat{Y_i} \equiv X_i^{'} \hat{\beta} \end{equation}\]

Donde \(\hat{Y_i}\) es el predictor (variable aleatoria) y \(\hat{\beta}\) es el estimador (parámetro). Por el TGM:

• \(E(\hat{Y_i}) = X_i^{'} \beta\) (predictor insesgado)

• \(V(\hat{Y_i}) = X_i^{'} V(\hat{\beta})X_i = \sigma^2X_i^{'}(X'X)^{-1}X_i\)

donde \(V(\hat{\beta}) = \sigma^2(X'X)^{-1}\)

Entonces, si el estimador \(\hat{\beta}\) es insegado y de varianza mínima, \(\hat{Y_i}\) es un predictor insegado y de varianza mínima, ambos en la clase de estimadores/predictores lineales e insesgados.

El resultado anterior surge del hecho que predecir requiere estimar (en este caso \(\beta\)). Retomamos el \(EMC\) en el caso de los estimadores:

\[\begin{equation} \tag{7.9} EMC(\hat{\beta}) \equiv E(\hat{\beta} - \beta)^2 \end{equation}\]

El \(EMC\) mide en promedio cuan lejos esta \(\hat{\beta}\) (estimador) de \(\beta\), el parámetro que se quiere estimar.

Recordar que por definición:

• \(V(\hat{\beta}) \equiv E(\hat{\beta} - E(\hat{\beta}))^2\) (dispersión)

• \(Sesgo(\hat{\beta}) \equiv E(\hat{\beta}) - \beta\) (centro)

A partir de (7.9) y las definiciones anteriores reescribimos el \(EMC\) en términos de la descomposición sesgo-varianza:

\[\begin{equation} \tag{7.10} EMC(\hat{\beta}) = Sesgo^2(\hat{\beta}) + V(\hat{\beta}) \end{equation}\]

Es decir, cuán mal estima \(\hat{\beta}\) depende de cuán descentrado está en relación a la verdad (sesgo) más cuán disperso es en relación a su propio centro (varianza).

Para ver la relación entre el error de estimación y el error de predicción se define el error de pronóstico como:

\[\begin{equation} \tag{7.11} Err(\hat{Y}) \equiv E(Y - \hat{Y})^2 \end{equation}\]

Como vimos antes, dado el modelo genérico¹⁵ \(Y = f(X) + u\) donde \(E(u) = 0\) y \(V(u) = \sigma^2\):

• \(f(X)\) es la parte sistemática

• \(u\) la parte no sistemática

Resultado importante en teoría de la predicción: si se quiere predecir una variable aleatoria \(Y\) con una constante \(m\) el mejor predictor es su esperanza,¹⁶ es decir, \(m = E(Y)\).¹⁷

Entonces, si \(E(Y) = f(X)\), \(u\) es no observable y \(f(X)\) es conocida, \(f(X)\) es el mejor predictor porque, como se dijo arriba, el mejor predictor de una variable aleatoria es su esperanza.

En la práctica \(f(X)\) es desconocida y, por lo tanto, se debe estimar \(\hat{f}(X)\).

\[\begin{equation} \tag{7.12} Err(\hat{Y}) = E(Y - \hat{f})^2 \end{equation}\]

\[\begin{equation} \tag{7.13} Err(Y - \hat{f}) = EMC(\hat{f}) + \sigma^2 \end{equation}\]

En términos de la ecuación (6.3) el error de pronóstico es la suma de un error reducible \((EMC)\) y otro irreducible \((\sigma^2)\). Además, notar que el error de predicción y el \(EMC\) difieren solo por una constante \(\sigma^2\).

Nuevamente, puede verse la relación entre predicción y estimación. Predecir correctamente \((Err(Y - \hat{f}))\) requiere estimar \((EMC(\hat{f}))\) correctamente (porque \(\sigma^2\) no se puede controlar). Es decir, tener bajo sesgo y baja varianza. Utilizando la descomposición:

\[\begin{equation} \tag{7.14} Err(Y-\hat{f}) = \underbrace{Sesgo^2(\hat{f}) + V(\hat{f})}_{\text{EMC}} + \sigma^2 \end{equation}\]

En la econometría tradicional si \(\hat{f}\) es insesgado minimizar el error de pronóstico es minimizar la varianza. Machine learning hace uso de que estrategias sesgadas pueden implicar una reducción drástica en la varianza, por lo tanto, puede ser que el mínimo \(EMC\) ocurra para predictores sesgados.

7.2 Aplicacion practica

Para hacer tablas con los resultados de las regresiones se puede utilizar el paquete stargazer. Algunas referencias en este documento y en este blog.

La biblioteca ISLR2 contiene la base de datos de Boston, que registra medv (mediana del valor de las casas en miles de USD) para \(506\) distritos censales en Boston. Buscaremos predecir medv usando \(12\) predictores como rm (número promedio de habitaciones por casa), age (edad promedio de las casas) y lstat (porcentaje de hogares con bajo status socioeconómico).

library(MASS)
library(ISLR2)
head(Boston)

##      crim zn indus chas   nox    rm  age    dis rad tax ptratio  black lstat
## 1 0.00632 18  2.31    0 0.538 6.575 65.2 4.0900   1 296    15.3 396.90  4.98
## 2 0.02731  0  7.07    0 0.469 6.421 78.9 4.9671   2 242    17.8 396.90  9.14
## 3 0.02729  0  7.07    0 0.469 7.185 61.1 4.9671   2 242    17.8 392.83  4.03
## 4 0.03237  0  2.18    0 0.458 6.998 45.8 6.0622   3 222    18.7 394.63  2.94
## 5 0.06905  0  2.18    0 0.458 7.147 54.2 6.0622   3 222    18.7 396.90  5.33
## 6 0.02985  0  2.18    0 0.458 6.430 58.7 6.0622   3 222    18.7 394.12  5.21
##   medv
## 1 24.0
## 2 21.6
## 3 34.7
## 4 33.4
## 5 36.2
## 6 28.7

Comenzaremos usando la función lm() para ajustar un modelo de regresión lineal simple, con medv como respuesta y lstat como predictor.¹⁸

lm.fit1 = lm(medv ~ lstat, data = Boston)
summary(lm.fit1)

## 
## Call:
## lm(formula = medv ~ lstat, data = Boston)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.168  -3.990  -1.318   2.034  24.500 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 34.55384    0.56263   61.41   <2e-16 ***
## lstat       -0.95005    0.03873  -24.53   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.216 on 504 degrees of freedom
## Multiple R-squared:  0.5441, Adjusted R-squared:  0.5432 
## F-statistic: 601.6 on 1 and 504 DF,  p-value: < 2.2e-16

names(lm.fit1)

##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"

attach(Boston)
plot(lstat, medv, cex=.5)
abline(lm.fit1, col = 'red')

A continuación se examinan algunos gráficos de diagnóstico con la función plot(). En general, este comando producirá un gráfico a la vez (presionando Enter se generará el siguiente gráfico). Sin embargo, se pueden ver los cuatro gráficos juntos con las funciones par() y mfrow(), que le dicen a R que divida la pantalla de visualización en paneles separados. Por ejemplo, par(mfrow = c(2, 2)) divide la región del gráfico en una cuadrícula de paneles de \(2 \times 2\).

• Residuals vs. fitted values: residuos vs. valores ajustados se utiliza para detectar patrones de variables omitidas, heterocedasticidad, etc.
• Scale Location: residuos estandarizados vs. valores ajustados. Similar al anterior y se utiliza para detectar patrones en los residuos.
• Normal Q-Q: cuantiles teóricos de la distribución normal estándar vs. cuantiles reales de residuos estandarizados. Se utiliza para evaluar la normalidad de los errores.
• Residuals vs. Leverage: el leverage es una medida de cuán influyente es una observación en el valor de los coeficientes. Este gráfico se utiliza para detectar las (posibles) observaciones influyentes y los valores atípicos al mismo tiempo.

par(mfrow = c(2, 2))
plot(lm.fit1)

Para obtener la predicción y el \(EMC\):

medv_hat = predict(lm.fit1, data = Boston)
emc = mean((medv - medv_hat)^2)
emc

## [1] 38.48297

Regresión lineal múltiple:

lm.fit2 = lm(medv ~ lstat + age, data = Boston)
summary(lm.fit2)

## 
## Call:
## lm(formula = medv ~ lstat + age, data = Boston)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.981  -3.978  -1.283   1.968  23.158 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 33.22276    0.73085  45.458  < 2e-16 ***
## lstat       -1.03207    0.04819 -21.416  < 2e-16 ***
## age          0.03454    0.01223   2.826  0.00491 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.173 on 503 degrees of freedom
## Multiple R-squared:  0.5513, Adjusted R-squared:  0.5495 
## F-statistic:   309 on 2 and 503 DF,  p-value: < 2.2e-16

BONUS:

library(stargazer)
stargazer(lm.fit1, lm.fit2, type = 'latex', 
          keep.stat=c('n','adj.rsq'),
          dep.var.labels   = 'Valor promedio de la vivienda',
          covariate.labels = c('Porcentaje socioeconómico bajo',
                       'Antigüedad de las casas'))

Utilizando todas las variables de la base de datos:

lm.fit = lm(medv ~ ., data = Boston)
summary(lm.fit)

## 
## Call:
## lm(formula = medv ~ ., data = Boston)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.595  -2.730  -0.518   1.777  26.199 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.646e+01  5.103e+00   7.144 3.28e-12 ***
## crim        -1.080e-01  3.286e-02  -3.287 0.001087 ** 
## zn           4.642e-02  1.373e-02   3.382 0.000778 ***
## indus        2.056e-02  6.150e-02   0.334 0.738288    
## chas         2.687e+00  8.616e-01   3.118 0.001925 ** 
## nox         -1.777e+01  3.820e+00  -4.651 4.25e-06 ***
## rm           3.810e+00  4.179e-01   9.116  < 2e-16 ***
## age          6.922e-04  1.321e-02   0.052 0.958229    
## dis         -1.476e+00  1.995e-01  -7.398 6.01e-13 ***
## rad          3.060e-01  6.635e-02   4.613 5.07e-06 ***
## tax         -1.233e-02  3.760e-03  -3.280 0.001112 ** 
## ptratio     -9.527e-01  1.308e-01  -7.283 1.31e-12 ***
## black        9.312e-03  2.686e-03   3.467 0.000573 ***
## lstat       -5.248e-01  5.072e-02 -10.347  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.745 on 492 degrees of freedom
## Multiple R-squared:  0.7406, Adjusted R-squared:  0.7338 
## F-statistic: 108.1 on 13 and 492 DF,  p-value: < 2.2e-16

Utilizando todas las variables de la base de datos menos age:

lm.fit1 = lm(medv ~ . - age, data = Boston)
summary(lm.fit1)

## 
## Call:
## lm(formula = medv ~ . - age, data = Boston)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.6054  -2.7313  -0.5188   1.7601  26.2243 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  36.436927   5.080119   7.172 2.72e-12 ***
## crim         -0.108006   0.032832  -3.290 0.001075 ** 
## zn            0.046334   0.013613   3.404 0.000719 ***
## indus         0.020562   0.061433   0.335 0.737989    
## chas          2.689026   0.859598   3.128 0.001863 ** 
## nox         -17.713540   3.679308  -4.814 1.97e-06 ***
## rm            3.814394   0.408480   9.338  < 2e-16 ***
## dis          -1.478612   0.190611  -7.757 5.03e-14 ***
## rad           0.305786   0.066089   4.627 4.75e-06 ***
## tax          -0.012329   0.003755  -3.283 0.001099 ** 
## ptratio      -0.952211   0.130294  -7.308 1.10e-12 ***
## black         0.009321   0.002678   3.481 0.000544 ***
## lstat        -0.523852   0.047625 -10.999  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.74 on 493 degrees of freedom
## Multiple R-squared:  0.7406, Adjusted R-squared:  0.7343 
## F-statistic: 117.3 on 12 and 493 DF,  p-value: < 2.2e-16

Bibliografia

Hastie, Trevor, Robert Tibshirani, y Jerome Friedman. 2008. The Elements of Statistical Learning. Data Mining, Inference, and Prediction. Springer. https://hastie.su.domains/Papers/ESLII.pdf.

Wooldridge, Jeffrey. 2012. Introductory Econometrics: A Modern Approach. Vol. 5th edition. South-Western College Publishing.

---

14. Para más detalles véase (Wooldridge 2012) Apéndice E (página 807).↩︎

15. No necesariamente lineal.↩︎

16. Cuando lo “mejor” se mide con el \(EMC\). Para más detalles véase (Hastie, Tibshirani, y Friedman 2008).↩︎

17. En el Capítulo 10 veremos que este resultado es importante para la metodología de arboles de decisión.↩︎

18. Para más detalles véase ?stats::lm.↩︎