1 Extensiones de la regresión lineal

En la clase introductoria a la regresión vimos cómo se ajusta un modelo de regresión lineal simple:

\[\begin{equation} y(x) = \beta_0 + \beta_1 x \tag{1.1} \end{equation}\]

donde \(y(x)\) es la variable dependiente, \(x\) es la indepentiente y \(\beta_0\) y \(\beta_1\) son los coeficientes de regresión estimados. Es difícil, sin embargo encontrar un modelo del mundo real, que pueda ser descrito por una ecuación tan sencilla, pues la mayoría de los fenómenos son resultado de muchos factores. Es por eso que la regresión lineal puede ser generalizada para incluir más variables independientes, de modo que se ajusta un modelo de Regresión lineal múltiple:

\[\begin{equation} y(x_1, x_2, \dots, x_n) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots + \beta_n x_n \tag{1.2} \end{equation}\]

Además de la existencia de múltiples factores que afectan las variables dependientes, estas pueden tener relaciones no-lineales con las variables independientes. Aunque la regresión lineal sólo sirve para estimar efectos lineales como \(\beta_i x_i\), no existe ningún tipo de limitación sobre la forma de \(x_i\), lo que permite estimar parámetros para variables independientes transformadas como \(x_i' = x_i^2, x_i'' = x_i^3\). Este tipo de trasformaciones permite estimar los efectos \(\beta_i\) para ecuaciones polinomiales de cualquier grado para representar relaciones que difieren mucho de la línea recta:

\[\begin{equation} y(x_1, x_2) = \beta_0 + \beta_1 x_1 + \beta_1' x_1^2 + \beta_2 x_2 + \beta_2' x_2^2 \tag{1.3} \end{equation}\]

2 Representación visual

Como bien sabemos, la regresión lineal simple representa la relación graficada en el plano cartesiano entre dos variables como una recta con inclicación \(\beta_1\) (Figura 2.1).

\[ y(x) = \beta_0 + \beta_1 x \]

Regresión lineal simple.

Figure 2.1: Regresión lineal simple.

En ocasiones, los datos de una variable \(y\), pueden tomar formas que no son lineales, por lo que se pueden representar con una ecuación polinomial. La figura 2.2, por ejemplo fue generada con la siguiente ecuación:

\[ y(x) = \beta_0 + \beta_1 x + \beta_1' x^2 + \beta_1'' x^3 \]
Regresión lineal polinomial

Figure 2.2: Regresión lineal polinomial

Un ejemplo más generalizable, que aquellos representados por las figuras 2.1 y 2.2, es aquel de variables dependientes que dependen de más de una variable independiente, siendo posible incluso las interacciones entre variables (2.3):

\[ y(x_1, x_2) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_1 x_2 \]

Regresión lineal múltiple, el eje vertical representa la variable de respuesta.

Figure 2.3: Regresión lineal múltiple, el eje vertical representa la variable de respuesta.

3 Ajustando regresiones lineales en R

La función por default para hacer regresión lineal en R es lm. Pongamos entonces en práctica esta función que hemos utilizado anteriormente para ANOVA, en regresión lineal. Para ello trataremos de estimar los coeficientes \(\beta\) que generaron las figuras 2.1, 2.2 y 2.3. El código utilizado para generarlas es:

set.seed(12345) # Semilla
x1 <- rnorm(1000) 
x2 <- rnorm(1000, mean = 2, sd = 2)
y.lin <- x1 * runif(1, 2, 4) + rnorm(1000, mean = 0, sd = 1) # Simple
y.mult <- x1 * runif(1, 2, 4) + x2 * runif(1, 2, 4) + x1 * x2 * runif(1, -4, 4) + rnorm(1000, mean = 0, sd = 3) #multiple
y.poli <- x1 * runif(1, -4, 4)  + x1^2 * runif(1, -4, -2) + x1^3 * runif(1, -4, 4) +  rnorm(1000, mean = 10, sd = 10) # polinomial

lin.rels <- data.frame(x1, x2, y.lin, y.mult, y.poli)

El proceso de selección de un modelo es complicado, sobretodo cuando desconocemos si existen o no asociaciones entre variables contínuas independientes y la dependiente. Por esta razón, es mucho más sencillo comenzar con datos cuyo origen conocemos.

Las variables con asociación estrictamente lineal son y.lin y x1, de modo que el modelo lineal lo ajustamos con:

m.simp <- lm(y.lin ~ x1, lin.rels)

Y para tener acceso a los coeficientes \(\beta_0\), \(\beta_1\) y el de determinación \(R^2\):

summary(m.simp)
## 
## Call:
## lm(formula = y.lin ~ x1, data = lin.rels)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0321 -0.6232  0.0003  0.6693  3.1658 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.003676   0.031428  -0.117    0.907    
## x1           2.524028   0.031450  80.255   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9928 on 998 degrees of freedom
## Multiple R-squared:  0.8658, Adjusted R-squared:  0.8657 
## F-statistic:  6441 on 1 and 998 DF,  p-value: < 2.2e-16

Los coeficientes estimados están en la columna Estimate. La información para \(\beta_0\) en la fila (Intercept) y para \(\beta_1\) en la fila x1. Recordemos que la hipótesis nula para \(\beta_1\) dice que \(y\) no cambia con respecto de \(x\). En este caso, podemos rechazarla, por lo que es muy probable que \(\beta_1 \neq 0\). Podemos ver, al final del summary el coeficiente de determinación Adjusted R-squared: 0.8657, que es bastante alto.

La intepretación de las pruebas de hipótesis es igual en regresión múltiple, así que podemos continuar con este modelo. Como de antemano sabemos que hay una interacción la fórmula del modelo será x1 * x2 = x1 + x2 + x1 : x2:

m.mult <- lm(y.mult ~  x1 * x2, lin.rels)
summary(m.mult)
## 
## Call:
## lm(formula = y.mult ~ x1 * x2, data = lin.rels)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.1055  -1.8437  -0.1218   1.9550   7.4719 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.02932    0.12570   0.233    0.816    
## x1           2.09167    0.12211  17.129   <2e-16 ***
## x2           2.48605    0.04494  55.324   <2e-16 ***
## x1:x2       -0.55690    0.04420 -12.598   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.861 on 996 degrees of freedom
## Multiple R-squared:  0.7766, Adjusted R-squared:  0.7759 
## F-statistic:  1154 on 3 and 996 DF,  p-value: < 2.2e-16

Como es evidente, podemos rechazar todas las hipótessis con respecto a los efectos de x1 y x2 y x1:x2. Para estos datos, \(R^2\) fue un poco más bajo.

Veamos por último la regresión polinomial. La especificación aquí, es un poco diferente a lo que hemos visto antes, porque tenemos que decirle a R que cada término \(x\) del modelo es sobre la misma variable \(x\) especificada sin ninguna potencia con la función I() (por identical):

m.poli <- lm(y.poli ~ x1 + I(x1^2) + I(x1^3), lin.rels)
summary(m.poli)
## 
## Call:
## lm(formula = y.poli ~ x1 + I(x1^2) + I(x1^3), data = lin.rels)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -28.960  -7.037  -0.249   6.491  39.748 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  10.3130     0.3932  26.226   <2e-16 ***
## x1            0.3394     0.5459   0.622    0.534    
## I(x1^2)      -4.0730     0.2312 -17.618   <2e-16 ***
## I(x1^3)      -1.8641     0.1501 -12.415   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.07 on 996 degrees of freedom
## Multiple R-squared:  0.4644, Adjusted R-squared:  0.4627 
## F-statistic: 287.8 on 3 and 996 DF,  p-value: < 2.2e-16

I la interpretación de los coeficientes es la misma que anteriormente. Para cerrar, podemos ver que \(R^2\) es bastante más bajo que para los ejemplos anteriores, y esto se debe a rnorm(1000, mean = 10, sd = 10), lo cual también se ve reflejado en el intercepto del modelo, con una media de 10.3130.

Regresar al índice del curso