12  Regresión lineal múltiple

Esta sección está en fase de preparación

12.1 Modelos

De manera sucinta, entenderemos aquí por modelo a una representación abstracta, formal, de un fenómeno o de un sistema de interés. Esta representación debe ser suficientemente precisa y manejable para poder cumplir algún propósito. Pero es imposible maximizar el realismo del modelo y que este sea tratable desde un punto de vista práctico. En este contexto, es inevitable citar el aforismo de George Box “Todos los modelos son erróneos, pero algunos son útiles”; los modelos son simplificaciones de la realidad y, por tanto, imperfectos, pero aun así valiosos para entender y predecir fenómenos, debiendo enfocarse en su utilidad práctica más que en su exactitud absoluta.

12.1.1 Los elementos de un modelo formal

Variable respuesta (habitualmente designada por \(y\)) es aquella magnitud cuya variación se intenta justificar en términos de una o más variables explicativas (\(x\) cuando sea una sola, \(x_1, x_2, \dots, x_k\), cuando en total sean \(k\) variables). Formalmente, la dependencia de la variable respuesta respecto a las variables explicativas se expresa \[y=f\left( x_1, x_2, \dots ,x_k \right)\]

en donde \(f()\) representa a algún tipo de función matemática (sin detallar por ahora). En esta función participaran, además de las variables explicativas, ciertas constantes cuyo valor resultará crítico y que son los parámetros de la función \(\theta_1, \dots, \theta_p\).

Un primer aspecto de interés sobre \(f()\) es si se trata de una relación lineal o no. Una función \(f()\) es lineal respecto a los argumentos \(\gamma_1\) y \(\gamma_2\) si, para cualquier constante \(\lambda_1\) y \(\lambda_2\) se da la relación \[f(\lambda_1 \gamma_1 + \lambda_2 \gamma_2) = \lambda_1 f(\gamma_1) + \lambda_2 f(\gamma_2)\]

distinguiendo si lo es respecto a los parámetros (todos aparecen de forma aditiva), respecto a las variables (son estas las que aparecen de forma aditiva) o respecto a ambos (tanto parámetros como variables aparecen de forma aditiva). Un modelo no lineal se dice que es intrínsecamente lineal cuando mediante alguna transformación de las variables (logaritmos, inversos, etc) se consigue un modelo lineal. Por ejemplo, el modelo potencial \(y=\alpha x^\beta\) no es lineal, pero sí que es intrinsecamente lineal, ya que tomando logaritmos en ambos miembros, se obtiene \(\log{y} = \log{\alpha} + \beta \log{x}\), de modo que considerando \(Y=\log{y}\), \(X=\log{x}\) y \(A=\log{\alpha}\), el modelo original se puede escribir como \(Y=A+\beta X\), que es un modelo lineal respecto a parámetros y variables.

A menudo, resulta conveniente indicar de forma explícita los parámetros de los que depende la función. Si los designamos de forma genérica como \(\theta_1, \dots,\theta_k\), expresaremos el modelo como

\[\begin{equation}\tag{1}\label{moddet} y=f\left( x_1, x_2,\dots,x_k; \theta_1, \dots,\theta_k \right) \end{equation}\]

12.1.2 Modelos estadísticos

La expresión \(\ref{moddet}\) es de naturaleza puramente determinística, para los mismos valores de las variables predictivas, se obtiene el mismo valor de la respuesta \(y\). Sin embargo, desde un punto de vista empírico, en donde el modelo se construye sobre un conjunto de observaciones, hace falta incluir el componente aleatorio \(\epsilon\) que representa a aquella parte de la variabilidad en la respuesta que no queda explicada por el cambio en las variables explicativas. \[\begin{equation}\tag{2}\label{modest} y=f\left( x_1, x_2, \dots ,x_k; \theta_1, \dots,\theta_k; \epsilon \right) \end{equation}\] Este componente podrá aparecer también de forma aditiva o no.

En términos de las variables explicativas:

• Todas continuas Modelos de regresión (simple, si solo hay una variable explicativa y múltiple si hay más)
• Todas categóricas Modelos de análisis de la varianza (ANOVA). Se habla de ANOVA de un factor si hay solo hay una variable explicativa, y de ANOVA factorial si hay más de una.
• Continuas y categóricas, modelos de análisis de la covarianza (ANCOVA)

En términos de la variable respuesta

• Respuesta continua, se trata de los modelos de regresión, ANOVA y ANCOVA
• Respuesta categórica, son los modelos de regresión logística binaria y multinomial, en función de que la respuesta sea binaria o no.
• Recuento: Modelos log-lineales, regresión de Poisson
• Tiempo hasta que ocurre un evento, modelos de análisis de supervivencia

Desde un punto de vista metodológico, es un error decir que los datos se ajustan al modelo. Es el modelo el que se debe de ajustar a los datos. En la mayoría de las circunstancias, para un conjunto dado de datos, va a haber un gran número de modelos posibles. Algunos de ellos serán más o menos plausibles. Un primer criterio de elección entre modelos alternativos es que el modelo debe ser minimal en el sentido del principio de parsimonia, según el cual explicaciones simples son preferibles a explicaciones complejas (es la navaja de Ockham). Esto se traduce en que los modelos deben tener el menor número de parámetros posible, resultando los modelos lineales preferibles a los no lineales. Los experimentos que permiten obtener los datos, deben de estar basados en pocos supuestos - Los modelos deben reducirse hasta un mínimo adecuado -

12.2 Análisis preliminar

• Verificar la presencia y codificación de datos faltantes (identificar si los valores cero representan datos faltantes, lo mismo con codificaciones del tipo 999y similares)
• Examen visual de cada variable (plot directamente con cuantitativas representa el valor de la variable en función de su orden en el data.frame)
• Explorar la relación entre variables (tapply, plot, tree, gam)
• Consideraciones metodológicas: - variables que tiene sentido incluir en el modelo - Necesidad de incorporar alguna transformación de la respuesta o de las variables explicativas - Interacciones que se deben incluir - Términos no lineales - Presencia de pseudoreplicación

Una variable debe ser mantenida en un modelo solo si su exclusión supone un aumento significativo en la devianza.

El principio de parsimonia se traduce en que es preferible un modelo con \(n-1\) parámetros que uno con \(n\) parámetros. Uno con \(K-1\) variables explicativas que uno con \(K\) variables. Un modelo lineal a uno no lineal. Un modelo sin interacciones que uno con interacciones.

12.3 Tipos de modelos el número de observaciones y la cantidad de parámetros

• Modelo saturado: presenta un parámetro por cada observación. El ajuste es perfecto. Carece de grados de libertad. Su poder explicativo es nulo
• Modelo maximal: Contiene todos los factores, interacciones y covariables que pueden resultar de interés. Sus grados de libertad son \(g.l.=n-K-1\), siendo \(n\) el tamaño muestral y \(k\) el número de variables explicativas.
• Modelo minimal (adecuado): Es un modelo simplificado con \(1\le K^\prime \le K\) parámetros. Su ajuste es menor que el del modelo maximal, pero no de forma significativa. Sus grados de libertad son \(g.l.=n-K^\prime-1\). Su poder explicativo es \(R^2=SSR/SSY\)
• Modelo nulo: Presenta un único parámetro, la media de la respuesta \(\bar{y}\)}. Su ajuste es nulo \(SSE=SSY\). \(g.l.=n-1\). Carece de poder explicativo.

---

12.4 Comparación entre modelos alternativos

12.4.1 Contraste del cambio en la varianza explicada

El incremento de varianza explicada al añadir más regresores (modelo \(m_2\)) a un modelo base (\(m_1\)) de referencia, viene dado por \[R^{2}_{2}-R^{2}_{1}\]

Para determinar si la ampliación del modelo supone un aumento significativo en la proporción de varianza explicada, se puede usar el estadístico

\[ F_{exp}= \frac{\left( R^{2}_{2}-R^{2}_{1} \right) / \left(k_2 - k_1 \right) } {\left( 1-R^2_2 \right)/ \left(n-k_2 - 1 \right)} \]

en donde \(R^{2}_{1}\) y \(R^{2}_{2}\) son los coeficientes de determinación de los modelos base, con \(k_1\) regresores, y ampliado, con \(k_2\) regresores, respectivamente y \(n\) el tamaño muestral.

Bajo la hipótesis nula de que los dos modelos explican la misma proporción de la varianza de la respuesta, la distribución de este estadístico es \[F_{exp} \sim F_{k_2 - k_1,\space n-k_2 -1 } \]

12.4.2 Tamaño del efecto \(f^{2}\) de Cohen

El efecto \(f^{2}\) de Cohen en regresión lineal es una medida del tamaño del efecto que cuantifica la importancia de las variables predictoras, calculando la proporción de varianza explicada \((R^{2})\) respecto a la no explicada \((1-R^{2})\)

\[f^{2}=\frac{R^{2}}{1-R^{2}}\]

\(f^{2}\) se interpreta con umbrales para efectos pequeños (\(f^{2}\le0.02\)), medianos (\(f^{2}\le0.15\)) y grandes (\(f^{2} > 0.35\)).

Es una métrica útil en regresión múltiple para evaluar la magnitud del impacto de un conjunto de predictores, y se puede usar como alternativa al coeficiente de determinación \(R^{2}\) para comparar modelos o evaluar el cambio en la varianza explicada. Mientras que \(R^{2}\) indica la proporción total de varianza explicada, \(f^{2}\) se puede hacer relativo a la contribución específica de un factor o un conjunto de predictores. 

\[f^{2}=\frac{R_{2}^{2}-R_{1}^{2}}{1-R_{2}^{2}}\]

12.5 Técnicas de regularización y selección del mejor modelo

cf. Gil (2018) cf. Gil (2018)

test de estabilidad estructural de Chow cf. Rico cf. script R de Rico

12.5.1 Comparación de la variabilidad residual \((SCCR)\)

Consideramos dos modelos, un modelo más general, M₁ con \(p_1\) variables explicativas y un modelo reducido M₂ constituido por un subconjunto \(p_2 (<p_1)\) de los predictores considerados en M₁ (hablamos entonces de modelos anidados). El estadístico

\[\frac{SCCR_1-SCCR_2}{SCCR_1}\]

Si se cumplen las condiciones de regularidad del modelo (normalidad residual, homocedasticidad e independencia) este criterio es equivalente a la razón de verosimilitudes \[\frac{\max_{\boldsymbol{\beta},\sigma \in M_1}L(\boldsymbol{\beta}, \sigma,y)}{\max_{\boldsymbol{\beta},\sigma \in M_2}L(\boldsymbol{\beta}, \sigma,y)}\]

siendo \(L\) la función de verosimilitud. Consecuentemente, se puede construir el test de homogeneidad basado en

\[ F= \frac{\left( SCCR_2-SCCR_1 \right)/(p_1-p_2)}{SCCR_1/(n-p_1)} = \frac{\left( SCCR_2-SCCR_1 \right)/(gl_1-gl_2)}{SCCR_1/(gl_1)}\sim F_{df_1-df_2, df_1}.\]

Si \(F>F^{\alpha}_{df_1-df_2, df_1}\) supone que la diferencia entre \(SCCR_1\) y \(SCCR_2\) es grande y el modelo con más parámetros es superior al modelo reducido. En caso contrario, una diferencia \(SCCR_1-SCCR_2\) pequeña supone que el ajuste de ambos modelos es similar y, en ese caso, es preferible aquel que tiene menos parámetros.

La verosimilitud depende directamente de los datos, de manera que si cambia el conjunto de datos, las verosimilitudes no son comparables. Tampoco es procedente el uso de este método si los modelos no son anidados, es decir, el modelo reducido no está constituido por un subconjunto de regresores del modelo general.

Cálculo práctico en R: anova(modelo1, modelo2)

12.5.2 Comparación a través del coeficiente de determinación \(R^{2}\)

El contraste presentado anteriormente se puede hacer a través de los coeficientes de determinación. El incremento de variabilidad explicada al añadir más regresores (modelo \(m_2\)) a un modelo base (\(m_1\)) de referencia viene dado por \[R^{2}_{2}-R^{2}_{1}\]

De manera que la expresión anterior se transforma en el estadístico

\[ F_{exp}= \frac{\left( R^{2}_{2}-R^{2}_{1} \right) / \left(k_2 - k_1 \right) } {\left( 1-R^2_2 \right)/ \left(n-k_2 - 1 \right)} \]

en donde \(R^{2}_{1}\) y \(R^{2}_{2}\) son los coeficientes de determinación de los modelos base, con \(k_1\) regresores, y ampliado, con \(k_2\) regresores, respectivamente y \(n\) el tamaño muestral.

Bajo la hipótesis nula de que los dos modelos explican la misma proporción de la varianza de la respuesta, la distribución de este estadístico es \[F_{exp} \sim F_{k_2 - k_1,\space n-k_2 -1 } \]

Cuestiones a destacar sobre \(R^2\):

• El coeficiente de determinación no es estrictamente sensible al tamaño muestral, pero en la práctica, más datos suelen generar un mejor ajuste.
• El coeficiente de determinación no debe usarse para comparar el ajuste de modelos construidos sobre diferentes transformaciones de los datos, ya que cambia la escala.
• Al aumentar el número de regresores, \(R^2\) aumenta de forma natural, es por esto que se define el coeficiente de determinación ajustado \(\bar{R^2}\), que considera el cociente entre cuadrados medios en lugar de las sumas de cuadrados consideradas por \(R^2\) (Peña, 2002; pp:388): \[\bar{R^2}=1-\frac{\hat\sigma^2_\epsilon}{s^2_y}\]
• No es recomendable comparar directamente \(R^2\) o \(\bar{R^2}\) entre modelos ajustados sobre conjuntos de datos distintos, ya que ambos dependen de la variabilidad total de la variable respuesta, si esta cambia, el denominador de \(R^2\) cambia, de manera que como porcentaje de variabilidad explicada por el modelo no aluden a una misma variabilidad total.

Para comparar modelos ajustados sobre conjuntos distintos de datos se utilizan otros indicadores, como RMSE, MAE y MAPE. También técnicas de validación cruzada, que evaluan el rendimiento predictivo, técnicas bootstrap (evalúa la estabilidad de las métricas) y diagramas de residuos.

12.5.3 Comparación a través del tamaño del efecto \(f^{2}\) de Cohen

El efecto \(f^{2}\) de Cohen en regresión lineal es una medida que cuantifica la importancia de las variables predictoras calculando la proporción de varianza explicada \((R^{2})\) respecto a la no explicada \((1-R^{2})\)

\[f^{2}=\frac{R^{2}}{1-R^{2}}\]

\(f^{2}\) se interpreta con umbrales para efectos pequeños (\(f^{2}\le0.02\)), medianos (\(f^{2}\le0.15\)) y grandes (\(f^{2} > 0.35\)).

Es una métrica útil en regresión múltiple para evaluar la magnitud del impacto de un conjunto de predictores, y se puede usar como alternativa al coeficiente de determinación \(R^{2}\) para comparar modelos o evaluar el cambio en la varianza explicada. Mientras que \(R^{2}\) indica la proporción total de varianza explicada, \(f^{2}\) se puede hacer relativo a la contribución específica de un conjunto de predictores (ajustando los dos modelos anidados correspondientes y determinando sus respectivos coeficientes de determinación). 

\[f^{2}=\frac{R_{2}^{2}-R_{1}^{2}}{1-R_{2}^{2}}\]

12.5.4 Raíz del error cuadrático medio (\(RMSE\))

Es la diferencia promedio entre los valores pronosticados por un modelo y los valores observados de la respuesta, es decir el promedio de la suma de los residuos crudos al cuadrado \[RMSE=\sqrt{ \frac{\sum{\left( y_i-\hat{y}_i \right)^2}}{n}} =\sqrt{\frac{\sum{e_i^2}}{n}} \]

Cálculo en R. Se puede elaborar la función

rmse<-function(modelo) {sqrt(mean(modelo$residuals^2))}

12.5.5 Error absoluto medio (\(MAE\))

Es el promedio de las diferencias absolutas entre los valores pronosticados y los valores observados, es decir, el promedio de los valores absolutos de los residuos crudos \[MAE=\frac{\sum{| y_i-\hat{y_i}|}}{n} = \frac{\sum{|e_i|}}{n} \]

Cálculo en R. Se puede elaborar la función

mae<-function(modelo) {mean(abs(modelo$residuals))}

12.5.6 Porcentaje medio de error absoluto (\(MAPE\))

MAPE (Mean Absolute Percentage Error) expresa el error de forma relativa a la magnitud observada.

\[MAPE=\frac{100}{n}\sum{\frac{| y_i-\hat{y_i}|}{y_i}}= \frac{100}{n}\sum{\frac{|e_i|}{y_i}} \]

Por ejemplo, si \(MAPE=10%\), los pronósticos se desvían, en promedio, un 10% respecto al valor observado. Al ser un indicador adimensional, es útil para comparar modelos en diferentes escalas, pero no es un indicador adecuado cuando las observaciones de la respuesta son próximas a cero.

Cálculo en R. Se puede elaborar la función

mape<-function(modelo) {mean(abs(modelo$residuals)/modelo$model[,1])*100}

El componente modelo$model es un data.frame cuya primera columna es siempre la variable respuesta usada en el ajusta. Si la respuesta se ha transformado, aquí se recoge dicha trasnsformación. Con lo que hay que tener precaución es con la posible presencia de valores faltantes (NA).

12.5.7 Criterio de información de Akaike (\(AIC\))

El criterio de información de Akaike (AIC, Akaike Information Criterion) es una medida basada en la Teoría de la Información que permite comparar modelos estadísticos de regresión. Sirve como criterio para determinar qué modelo, de entre varios posibles, equilibra mejor la bondad del ajuste y la complejidad, ya que es una métrica que penaliza por el número de parámetros. Se define como

\[AIC=-2log(\hat{L})+2p\]

en donde \(\hat{L}\) es la verosimilitud máxima del modelo y \(p\) su número de parámetros (incluyendo la constante). El AIC no tiene una escala absoluta, es decir, no indica la calidad absoluta de un modelo. Se debe interpretar de forma relativa, comparando los valores obtenidos para modelos alternativos ajustados a los mismos datos. Un AIC menor indica un modelo que mejora el equilbrio entre la calidad del ajuste y su complejidad.

Puede favorecer modelos más complejos si el tamaño de muestra es reducido. De aquí que se haya definido una versión corregida AICc para muestras pequeñas:

\[AIC_c = AIC + \frac{2p(p+1)}{n - p - 1}\]

el factor de corrección penaliza más el aumento de complejidad cuando el tamaño de muestra \(n\) es reducido. La regla más aceptada para utilizar esta corrección es que sea \(n/p<40\). En realidad, si el tamaño muestral es grande, la corrección se vuelve insignificante

Obtención en R: AIC(modelo1, modelo2,...)

12.5.8 Criterio de información bayesiano (\(BIC\))

El BIC (Bayesian Information Criterion) es un criterio para comparar modelos alternativos sobre los mismos datos similar al AIC, pero este realiza una penalización más fuerte por el aumento de complejidad

\[BIC = -2 \log(\hat{L}) + p \log(n)\]

en donde \(\hat{L}\) es la verosimilitud máxima del modelo y \(p\) su número de parámetros (siempre incluyendo la constante) y \(n\) el tamaño muestral.

Como ocurre con el AIC, carece de escala absoluta y su interés es como criterio comparativo, un BIC menor supone un mejor modelo. Esta métrica favorece más los modelos más simples que el AIC, por lo tanto es un buen criterio cuando se busca mayor parsimonia, resultando especialmente interesante en contextos en los que el tamaño de muestra es grande y se quiere evitar el sobreajuste. Obtención en R: BIC(modelo1, modelo2,...)

12.5.9 \(Cp\) de Mallows

Es un criterio que compara el ajuste de un modelo reducido -que incluye a un subconjunto de regresores- frente a un modelo completo -que incluye todas las variables disponibles. Se obtiene como:

\[C_p = \frac{SCCR_p}{\hat{\sigma}^2_{\epsilon}} - (n - 2p)\]

en donde \(SCCR_p\) es la suma de cuadrados residuales del modelo con \(p\) parámetros. \(\hat\sigma^2_\epsilon\) es la estimación de la varianza residual por parte del modelo completo, y \(n\) y \(p\) son el tamaño de la muestra y el número de parámetros en el modelo, respectivamente.

La construcción de \(Cp\) se basa en la esperanza del error cuadrático medio del modelo reducido \((ECM_p)\).

\[E[ECM_p] = (n - p)\,\sigma_\epsilon^2 + \text{sesgo}^2\]

El término \(sesgo\) aparece cuando faltan regresores relevantes en el modelo, de manera que si el modelo reducido incluye a todas las variables realmente explicativas, el sesgo es prácticamente nulo y \(Cp\approx p\). Sin embargo, cuando se omiten regresores importantes, el sesgo aumenta, lo que incrementa el \(ECM\) y, por tanto, el \(Cp\). Así \(Cp>p\) es indicativo de que el modelo no captura toda la estructura de la relación. Por otra parte, un valor \(Cp<p\) es un síntoma de sobreajuste, es decir, de que el modelo tiene demasiadas variables. En resumen:

• Si \(C_p \approx p\), el modelo tiene poco sesgo y es adecuado.
• Si \(C_p \gg p\), el modelo está subajustado (faltan variables).
• Si \(C_p < p\), el modelos puede estar sobreajustado (demasiadas variables).

\(Cp\) puede servir como criterio comparativo a la hora de seleccionar el mejor subconjunto de predictores, se trata de obtener el modelo con \(C_p\) cercano a \(p\) y lo más parsimonioso posible. Por tanto, el propósito de \(Cp\) es similar al de otras métricas como AIC y BIC, pero este índice está basado en la comparación con el modelo completo y la estimación del error cuadrático medio.

El coeficiente \(Cp\) puede resultar negativo cuando ocurra \[\frac{SCCR_p}{\hat{\sigma}^2} < (n - 2p)\]

Típicamente, esto sucede cuando el modelo reducido presenta sobreajuste, es decir, incluye demasiadas variables explicativas para un tamaño de muestra reducido. La estimación de la varianza residual, \(\hat\sigma^2_\epsilon\), del modelo completo es grande en relación a la estimación del modelo reducido y la penalización \((n-2p)\) resulta excesiva. Se debe destacar que en ningún caso se persigue que el modelo de lugar a un valor \(Cp\) negativo, como decimos, esto es una indicación de sobreajuste. Reiterando que el óptimo es \(Cp \approx p\)

Obtención en R. El lenguaje base no presenta implementado el cálculo de esta métrica. Hay varios paquetes que sí permiten obtenerla, por ejemplo olsrr: olsrr::ols_mallows_cp(modelo)

NotaComparativa entre AIC, AICc, BIC y \(Cp\) de Mallows

Criterio	¿Qué mide?	¿Cuándo usarlo?	Ventajas	Limitaciones
AIC	Equilibrio entre la calidad de ajuste y la complejidad	Muestras grandes \((n/p\ge 40)\)	Fácil de calcular, ampliamente usado	Puede favorecer modelos complejos en muestras pequeñas
AICc	Calidad del ajuste con penalización extra por complejidad	Muestras pequeñas \((n/p<40)\)	Corrige sesgo del AIC en muestras pequeñas	como AIC
BIC	Ajuste penalizado por tamaño de muestra	Cuando se busca parsimonia, especialmente con \(n\) grande	Penaliza más la complejidad que AIC, favorece modelos simples	Puede seleccionar modelos demasiado simples si \(n\) es pequeño
\(Cp\) de Mallows	Sesgo y varianza respecto al modelo completo	Selección de variables en regresión lineal	Intuitivo: \(Cp \approx p\) indica buen modelo	Requiere el modelo completo y la estimación \(\hat{\sigma}^2_\epsilon\)

12.6 Regresión robusta de Huber

El estimador de Huber es un método robusto que combina características de la regresión lineal clásica y la regresión robusta, reduciendo la influencia de valores atípicos. Sin embargo, la distribución exacta de los coeficientes bajo este estimador no sigue la misma teoría clásica que el estimador de mínimos cuadrados ordinarios (OLS).

12.6.1 Desventajas

12.6.1.1 Interpretación menos directa

• Los coeficientes estimados en regresión robusta se obtienen mediante ponderaciones adaptativas (M-estimadores), no por mínimos cuadrados. Esto hace que:

  • La interpretación sea menos intuitiva para quienes están acostumbrados a OLS.
  • Las pruebas clásicas (F-test, ANOVA) no sean aplicables directamente.