Regresión lineal múltiple

Especialización en Estadística Aplicada

Roberto Trespalacios

El vector de valores ajustados \( \boldsymbol{\hat{y}} \) en un modelo de regresión lineal se puede expresar como

\[ \begin{align*} \boldsymbol{\hat{y}} &= \boldsymbol{X \hat\beta} \\ &= \boldsymbol{X}(\boldsymbol{X}'\boldsymbol{X})^{-1} \boldsymbol{X}'\boldsymbol{y}\\ &=\boldsymbol{H} \boldsymbol{y} \end{align*} \]

La matriz \( \boldsymbol{H}_{n \times n} = \boldsymbol{X}(\boldsymbol{X}'\boldsymbol{X})^{-1} \boldsymbol{X}' \) es llamada regularmente, la matriz sombrero (matriz “hat” en inglés). Ella, “coloca” el sombrero sobre los valores de la variable \( \boldsymbol{y} \) y los “convierte” en valores predichos.

\[ \begin{align*} \hat{\boldsymbol{\varepsilon}} & =\boldsymbol{y} - \hat{\boldsymbol{y}} \\ & = \boldsymbol{y} -\boldsymbol{X}\hat{\boldsymbol{\beta}} \\ & = \boldsymbol{y} -\boldsymbol{Hy} \\ & =(\boldsymbol{I} - \boldsymbol{H})\boldsymbol{y} \end{align*} \]

El estimador MC (LS “least squares”) para \( \boldsymbol{\beta} \), tiene las siguientes propiedades:

• \( \hat{\boldsymbol{\beta}} \) es insesgado, o sea \( E(\hat{\boldsymbol{\beta}}) = \boldsymbol{\beta} \).

• \( Var(\hat{\boldsymbol{\beta}})=\sigma^2(\boldsymbol{X}'\boldsymbol{X})^{-1} \).

• Si no se asume normalidad de los errores, entonces el estimador de mínimos cuadrados MC, \( \hat{\boldsymbol{\beta}} \) es el mejor estimador dentro de los estimadores lineales insesgados de \( \boldsymbol{\beta} \).

• Si se asume normalidad de los errores, entonces \( \hat{\boldsymbol{\beta}} \) es el mejor estimador entre todos los estimadores insesgados de \( \boldsymbol{\beta} \).

Recordemos que:

\[ s^2 = \displaystyle \frac{SSE}{n-p-1} = \frac{\sum_{i=1}^n e^2_i}{n-p-1}= \frac{\boldsymbol{e}'\boldsymbol{e}}{n-p-1}= \frac{(\boldsymbol{y} - \boldsymbol{X} \boldsymbol{\beta})' (\boldsymbol{y} - \boldsymbol{X}\beta) }{n-p-1} \]

Pero \[ \begin{align*} SSE & = (\boldsymbol{y} - \boldsymbol{X} \boldsymbol{\beta})' (\boldsymbol{y} - \boldsymbol{X}\beta) \\ & = (\boldsymbol{y} - \boldsymbol{X}(\boldsymbol{X}'\boldsymbol{X})^{-1} \boldsymbol{X}' \boldsymbol{y})' (\boldsymbol{y} - \boldsymbol{X}(\boldsymbol{X}'\boldsymbol{X})^{-1} \boldsymbol{X}' \boldsymbol{y})\\ & = \boldsymbol{y}' (\boldsymbol{I} - \boldsymbol{H})' (\boldsymbol{I} -\boldsymbol{H})\boldsymbol{y}\\ & = \boldsymbol{y}' (\boldsymbol{I} -\boldsymbol{H})\boldsymbol{y} \end{align*} \]

donde \( \boldsymbol{H}=\boldsymbol{X}(\boldsymbol{X}'\boldsymbol{X})^{-1} \boldsymbol{X}' \) es la matriz “Hat” (\( \boldsymbol{H}^2=\boldsymbol{H} \) y \( \boldsymbol{H}'=\boldsymbol{H} \)).

Por lo tanto, la varianza estimada de los errores puede ser escrita como:

\[ s^2 = \displaystyle \frac{SSE}{n-p-1} = \frac{\boldsymbol{y}'(\boldsymbol{I} - \boldsymbol{X}(\boldsymbol{X}'\boldsymbol{X})^{-1} \boldsymbol{X}')\boldsymbol{y}}{n-p-1}= \frac{\boldsymbol{y}'(\boldsymbol{I} -\boldsymbol{H}) \boldsymbol{y}}{n-p-1} \]

Regordemos que el modelo matricial de regresión múltiple \( \boldsymbol{y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon} \), podemos escribirlo para cada \( y_i \), así,

\[ y_i=\beta_{0}+\beta_{1}x_{i1}+\cdots +\beta_{ip}x_{ip}+\varepsilon_i \]

Observaciones

• El estimado del coeficiente de regresión lineal \( \beta_j \), para \( j=1,2,\dots, p \), se representará por \( \hat{\beta}_j \).

• \( \hat{\beta}_j \) indica el cambio promedio en la variable de respuesta \( Y \) cuando la variable predictora \( X_j \) cambia en una unidad adicional asumiendo que las otras variables predictoras permanecen constantes.

Ejemplo 3

Sean las variables predictoras altura(cm) y edad(años) vs. la variable respuesta peso (kg). Interprete el siguiente modelo. \[ \widehat{peso}=-0.7 + 1.3 \times altura + 0.4 \times edad \]

Dentro de los métodos de inferencia, se encuentran:

• Pruebas de hipótesis eintervalos de confianza acerca de los coeficientes del modelo de regresión poblacional.

• Intervalos de confianza de las predicciones que se hacen con el modelo.

Observación

Suponemos que \( \varepsilon_i \stackrel{i.i.d}{\sim} N(0,\sigma^2 \boldsymbol{I}_n) \) equivalentemente, que \( \boldsymbol{y} \stackrel{i.i.d}{\sim} N(\boldsymbol{X}\boldsymbol{\beta}, \sigma^2 \boldsymbol{I}_n) \)

La decisión de aceptar o rechazar \( H_0 \) se va a tomar en base al estadístico \( F \), que se obtiene a partir de este análisis de la varianza:

\[ \frac{MSR}{MSE}=\frac{\frac{SSR}{p}}{\frac{SSE}{n-p-1}} \]

Este estadístico tiene una distribución \( F_{(p,n-p-1)} \) (bajo \( H_0 \)) y por lo tanto, la regla de decisión es de la siguiente forma: Rechazamos \( H_0 \), al nivel de significación \( \alpha \), cuando

\[ F=\frac{\frac{SSR}{p}}{\frac{SSE}{n-p-1}}>F_{(\alpha, gl_{num}=p, gl_{den}=n-p-1)} \]

Eliminación backward (hacia atrás)

El procedimiento comienza construyendo el modelo con todas las predictoras y en cada paso se elimina una variable. La secuencia del procedimiento es la siguiente. Se define un nivel de significación fijo \( \alpha \).

• El modelo inicial contiene todos los potenciales predictores (que denominado \( p \)).
• Si todas las variables producen una contribución parcial significativa (es decir, un estadístico \( t \) con \( p-valor < \alpha \)) entonces el modelo completo es el modelo final.
• De otro modo, se elimina la variable que tenga la menor contribución parcial (es decir, mayor \( p-valor \) del estadístico \( t \)) cuando las demás están en el modelo.
• Se ajusta el nuevo modelo con (\( p-1 \)) predictores y se repiten los pasos 2 y hasta que todas las variables en el modelo tengan un coeficiente estimado cuyo \( p-valor \) asociado al estadístico t sea menor a \( \alpha \)

Observación

Si hay una alta multicolinealidad en el conjunto de los \( p \) predictores, este procedimiento no es muy recomendable.

Selección forward (incorporando variables)

• Aquí en el modelo inicial se considera una regresión lineal simple que incluye a la variable predictora que da la correlación más alta con la variable de respuesta.
• Se incluye una segunda variable en el modelo, que es aquella variable dentro de las no incluidas aún, que da el \( p-value \) más bajo para la prueba \( t \) o el valor de la prueba de \( t \) más grande en valor absoluto.
• Se siguen incluyendo variables, notando que una vez que ésta es incluida ya no puede ser sacada del modelo.
  
• El proceso termina cuando los \( p-values \) para la prueba \( t \) de todas las variables que aún no han sido incluidas son mayores que 0.05 ó la prueba de \( t \) es menor que 2 para dichas variables. Si se usa la prueba de \( F \), entonces el proceso termina cuando todas las \( F \) son menores que 4.

Observación

Si se usa un punto de corte pequeño (digamos \( \alpha < 0.01 \)) se pueden perder covariables importantes. Si se usa un punto de corte grande (\( \alpha < 0.20 \)) probablemente, el modelo contendrá más variables. Una vez que el procedimiento finaliza, no todas las variables en el modelo necesariamente tendrán coeficientes parciales significativos.

Selección stepwise (paso a paso)

Es una modificación del procedimiento forward que elimina una variable en el modelo si ésta pierde significación cuando se agregan otras variables.

• La aproximación es la misma que la selección forward excepto que a cada paso, después de incorporar una variable, el procedimiento elimina del modelo las variables que ya no tienen contribución parcial significativa.

• Una variable que entró en el modelo en una etapa, puede eventualmente, ser eliminada en un paso posterior.

• En este caso será necesario definir un punto de corte para que ingrese una variable \( \alpha_I \) y otro para eliminarla del modelo \( \alpha_E \).

• Uno puede desear ser menos exigente (mayor \( p-valor \)) en el punto de corte para que una variable salga del modelo una vez que ingresó, o usar el mismo valor para ambos. Este procedimiento, en general produce modelos con menos variables que la selección forward.

Del análisis preliminar se pueden extraer las siguientes conclusiones:

• Las variables que tienen una mayor relación lineal con la esperanza de vida son: asesinatos (\( r= -0.78 \)), analfabetismo (\( r= -0.59 \)) y universitarios (\( r= 0.58 \)).

• Asesinatos y analfabetismo están medianamente correlacionados (\( r = 0.7 \)) por lo que posiblemente no sea útil introducir ambos predictores en el modelo.

• Las variables habitantes, área y densidad poblacional muestran una distribución exponencial, una transformación logarítmica posiblemente haría más normal su distribución.

En este caso se van a emplear la estrategia de stepwise mixto. El valor matemático empleado para determinar la calidad del modelo va a ser AIC (Akaike information criterion).

Observación 1

Puede usar en criterio BIC (Bayesian information criterion).

Observación 2

La función step se usa para los diferentes métodos de selección del mejor modelo (both, backward, forward). En el caso de “both”, será el metodo stepwise (paso a paso).

step(object = mod1, direction = "both", trace = 1)

Para los datos mtcars, de R,

1. Encuentra las correlaciones entre las variables. Concluya.

2. Replique lo hecho en la selección de variables (usar cualquiera de los métodos).

3. Ajuste el mejor modelo encontrado en la selección de variables. Concluya usando el \( R^2 \) y el \( R^2_a \); así como los \( p-values \) para las diferentes hipótesis sobre los \( \beta \).

4. Construya intervalos de confianza para los coeficientes \( \beta \). Interprete.

5. Interprete los coeficientes del mejor modelo ajustado y escribalo.

Recordemos que definimos el modelo de regresión clásico, como la relación entre la variable dependiente \( y_i \) y el vector \( \boldsymbol{x}_i \) de \( p \) variables predictoras. Tomando la j-esima potencia de cada una (\( j = 1,2,\dots,p \)),

\[ {\displaystyle y_{i}=\beta_{0}1+\beta_{1}x_{i1}+\beta_{1}x^2_{i2}\cdots +\beta_{p}x^p_{ip}+\varepsilon_{i},\qquad i=1,\ldots ,n,} \]

Estas \( n \) ecuaciones, pueden ser escritas de forma matricial, de la siguiente manera:

\[ \boldsymbol{y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}, \]

La ecuación toma la forma expandida: \[ \underbrace{\begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}}_{\boldsymbol{y}} = \underbrace{\begin{bmatrix} 1&x_{11}&x^2_{12} &\cdots &x^p_{1p}\\ 1&x_{21}&x^2_{22} & \cdots &x^p_{2p}\\ \vdots &\vdots &\ddots &\vdots \\ 1&x_{n1}&x^2_{n2}&\cdots &x^p_{np} \end{bmatrix}}_{\boldsymbol{X}} \underbrace{\begin{bmatrix} \beta_{0}\\ \beta_{1}\\ \beta_{2}\\ \vdots \\ \beta_{p} \end{bmatrix}}_{\boldsymbol{\beta}} + \underbrace{\begin{bmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_n \end{bmatrix}}_{\boldsymbol{\varepsilon}} \]

Para el modelo

\[ {\displaystyle y_{i}=\beta_{0}1+\beta_{1}x_{i1}+\beta_{1}x^2_{i1}\cdots +\beta_{p}x^2_{ip}+\varepsilon_{i},\qquad i=1,\ldots ,n,} \]

Se deben satisfacer los siguientes supuestos:

• \( E(y_i|\boldsymbol{x}_i)=\boldsymbol{x}_i'\boldsymbol{\beta} \)
• \( Var(y_i|\boldsymbol{x}_i)=\sigma^2 \)
• \( y_i|\boldsymbol{x}_i \stackrel{ind}{\sim} N(\boldsymbol{x}'_i\boldsymbol{\beta}, \sigma^2) \).

En R se pueden generar modelos de regresión polinómica de diferentes formas:

• Identificando cada elemento del polinomio (usando I()):

mod <- lm(formula = y ~ x + I(x^2) + I(x^3), data = datos)

Observación

El uso de I() es necesario ya que el símbolo \( \hat{} \) tiene otra función dentro de las fórmula de R.

• Con la función poly():

mod <- lm(formula = y ~ poly(x, 3), data = datos)

Cada una de las pendientes de un modelo de regresión lineal múltiple se define del siguiente modo:

• Si el resto de variables se mantienen constantes, por cada unidad que aumenta el predictor en cuestión, la variable \( Y \) varía en promedio tantas unidades como indica la pendiente.
• En el caso del predictor volume, si el resto de variables no varían, por cada unidad de volumen que aumenta el libro el peso se incrementa en promedio 0.71795 unidades.

Observaciones sobre el modelo

• Cuando un predictor es cualitativo, uno de sus niveles se considera de referencia (el que no aparece en la tabla de resultados) y se le asigna el valor de 0.
• El valor de la pendiente de cada nivel de un predictor cualitativo se define como el promedio de unidades que dicho nivel está por encima o debajo del nivel de referencia.
• Para el predictor tipo_tapas, el nivel de referencia es tapas blandas por lo que si el libro tiene este tipo de tapas se le da a la variable el valor 0 y si es de tapas duras el valor 1.
• Acorde al modelo generado, los libros de tapa dura son en promedio 184.04727 unidades de peso superiores a los de tapa blanda.

El modelo final tendrá la forma siguiente:

\[ \widehat{Peso \ libro} = 13.91557 + 0.71795 \ volumen + 184.04727 \ tipo\_tapas \]

• El modelo es capaz de explicar el \( 92.75\% \) de la variabilidad observada en el peso de los libros (\( R^2= 0.9275 \)).

• El valor de \( R^2_a \) es muy alto y cercano al \( R^2 \) lo que indica que el modelo contiene predictores útiles.

• El test F muestra un p-value de \( 1.455\times 10^{-07} \) por lo que el modelo en conjunto es significativo. Esto se corrobora con el p-value de cada predictor, en ambos casos significativo.