La regresión multivariable

El modelo con dos variables, por su simpleza, suele ser inadecuado en la práctica. La teoría económica contempla muchísimas variables y relaciones. Por ejemplo, el consumo, además de la renta, dependerá de otras muchas cuestiones, como de los productos sustitutivos o complementarios, el momento del año o incluso la condición social del comprador.

Es importante por tanto ampliar los modelos bivariables a modelos múltiples, en el que se incluyan más variables que, por teoría, deberían ser incluidas para explicar la variable dependiente. Los modelos de regresión múltiple cuentan, por tanto, con al menos tres variables: una dependiente, y dos regresoras:

\[Y_i = \beta_0 + \beta_1X_{1i} +\beta_2X_{2i} + e_i\]

siend \(β_0\) el término del intercepto, que mide el promedio de Y sobre todas las variables excluidas en el modelo (es decir, cuando \(X_i = 0\)). Los coefi cientes β_2 y β_3 son coeficientes de rgresión parcial. De esta forma, el valor esperado de Y, como variable dependiente, estará condicionada a las \(X_i\) de nuestro modelo descrito:

\[E(Y_i|X_{1i},X_{2i}) = β_0 + β_1X_{1i} + β_{2}X_{2i}\] de tal forma que el valor esperado de Y está condicionado a los valores fijos de sus independientes, o dicho de otra forma, el valor esperado de Y es la respuesta media a los valores dados de las regresoras X.

La explicación de los valores obtenidos por los coeficientes de regresión sería el siguiente: \(β_i\) mide el cambio en el valor de la media de Y, E(Y), por unidad de cambio en \(X_i\), siendo constantes el resto de coeficientes de regresión β. Es decir, que \(β_1\) mide el valor medio de Y por unidad de cambio de \(X_1\) cuando el valor de \(X_2\) es constante.

Principio de no multicolinealidad

Por multicolinealidad se entiende la relación lineal exacta entre dos variables. En este caso nos interesa la inexistencia de esta relación entre las variables \(X_1\) y \(X_2\). El motivo es sencillo: imaginemos que tenemos una relación exacta tal que \(X_2 = 2X_1\). De esta forma tendríamos que:

\[E(Y_i|X_{1i},X_{2i}) = β_0 + β_1X_{1i} + β_{2}X_{2i}\\ E(Y_i|X_{1i}) = β_0 + β_1X_{1i} + β_{1i}2X_{1i}\\ E(Y_i|X_{1i}) = β_0 + \alpha_1X_{1i},\ siendo\ \alpha_{1i}=β_{1i}+2β_{2i}\]

Por tanto, la no colinealidad significa que ninguna de las regresoras puede escribirse como combinación lineal exacta de las regresoras restantes en el modelo. Sin embargo, en la práctica para eliminar variables del modelo no es necesario que su relación sea exacta (es decir, de un \(R^2\) igual a 1), sino que con un valor por encima del 0.8 podríamos asumir que la relación entre las variables dependientes es tan elevada que resultaría adecuado eliminar alguna de ellas.

Para ello tendremos que utilizar diferentes métodos según la distribución de nuestras variables sean o no sean normales. Cuando la distribución de las variables es normal (paramétrico), usaremos el conocido como test de Pearson; en caso contrario, usaremos el test no paramétrico de Spearman.

Con ambos tests nos preguntamos si existe una relación (lineal) entre las variables, aunque sin asumir que una depende de la otra. El test de Pearson (r) se calcula como:

\[r=\frac{\sum xy}{\sqrt{\sum x^2\sum y^2}}\] tomando r un valor que puede ser positivo o negativo, siendo −1≤r≤1.

En el caso del test de Spearman se calcula como:

\[r_S=\frac{6\sum d{_i^2}}{n(n^2-1)}\] siendo n la cantidad de sujetos en nuestra muestra, y \(d_i=x_i-y_i\), es decir la diferencia entre los rangos X e Y. De nuevo, los valores de \(r_s\) toman valores entre -1 y +1.

Para saber cuál de los dos métodos elegir, tendremos que conocer la distribución de nuestras variables. Para ello podríamos utilizar gráficos como el histograma, pero también pruebas estadísticas, como el test Saphiro Wilk, para el cual:

\[ H_0=distribución\ normal\\ H_1=distribución\ no\ normal \]

siendo el valor p la probabilidad de que se cumpla nuestra hipótesis alternativa, y por tanto con un p valor inferior al 0.05 podremos confirmar la no normalidad de la distribución, mientras que con un p valor superior al 0.05, diremos que no tenemos evidencias suficientes como para rechazar la hipótesis nula.

Si creamos una base de datos de prueba, tendremos que:

#Creamos una muestra aleatoria de variables x e y
x1 <- c(10, 10, 10, 10, 12, 12, 12, 12, 12, 12, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 16, 16, 16, 16, 18, 18, 18, 20)
x2 <- c(14, 12, 14, 13, 16, 12, 15, 16, 14, 12, 18, 19, 21, 18, 17, 17, 18, 20, 19, 18, 19, 16, 19, 22, 21, 24, 26, 24, 26, 25)
y <- c(10, 14, 12, 10, 12, 10, 14, 14, 10, 12, 12, 10, 16, 16, 14, 18, 12, 16, 16, 14, 18, 14, 14, 18, 16, 18, 16, 18, 20, 20)

df1 <- data.frame(y, x1, x2)
shapiro.test(df1$x1)
## 
##  Shapiro-Wilk normality test
## 
## data:  df1$x1
## W = 0.92239, p-value = 0.03095
shapiro.test(df1$x2)
## 
##  Shapiro-Wilk normality test
## 
## data:  df1$x2
## W = 0.95204, p-value = 0.1917

Dado que en un ejercicio encontramos no normalidad, y en el segundo sí podríamos asumir la normalidad, trabajaremos con el método de testeo de Spearman, aunque vamos a ver los dos casos, para aprender a resolver el ejercicio.

Lo que buscamos con el test de correlación es resolver la siguiente hipótesis:

\[ H_0 = \rho = 0 \\ H_1 = \rho \neq 0 \]

es decir, nuestra hipótesis nula nos indica que no existen evidencias de relación entre las variables, y la hipótesis alternativa pretende demostrar que existe algún tipo de relación entre las variables. Su grado de relación nos dirá si debemos o no eliminar alguna de ellas.

#Método de correlación paramétrico de Pearson
cor.test(df1$x1, df1$x2, method= c("pearson"))
## 
##  Pearson's product-moment correlation
## 
## data:  df1$x1 and df1$x2
## t = 12.426, df = 28, p-value = 6.511e-13
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8373325 0.9615881
## sample estimates:
##       cor 
## 0.9200536

Como vemos, el p-valor es muy inferior a 0.05. Es decir, podemos garantizar que existe una relación superior a 0 entre las variables, pero además esta relación es de 0.92, es decir, muy muy elevada. Todo indica que debemos eliminar alguna de las dos variables.

#Método de correlación no paramétrica de Spearman
cor.test(df1$x1, df1$x2, method= c("spearman"),exact=FALSE)
## 
##  Spearman's rank correlation rho
## 
## data:  df1$x1 and df1$x2
## S = 339.84, p-value = 3.061e-13
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.9243963

Como vemos, el test de Spearman tiene unos resultados muy muy parecidos a los de Pearson, confirmando los resultados que hemos visto anteriormente. Esta relación la podemos ver también gráficamente:

#install.packages("corrplot")
library(corrplot)
## corrplot 0.84 loaded
COR = cor(df1, method = "pearson")
corrplot(COR, method="number", type="upper")

Como vemos, la relación entre \(X_1\) y \(X_2\) es de 0.92, mientras que la relación entre \(X_1\) e \(Y\) es de 0.75, y entre \(X_2\) e \(Y\) es de 0.77. En este caso, convendría eliminar la variable \(X_1\) de nuestro modelo, ya que \(X_2\) tiene una mayor relación con \(X_2\), y es una variable que explica linealmente el comportamiento de \(X_1\).

Para que entendamos mejor por qué es importante eliminar estas variables, realicemos dos regresiones lineales: una incluyendo las dos variables independientes, y una segunda incluyendo únicamente la segunda.

reg1 = lm(df1$y ~ df1$x1 + df1$x2)
summary(reg1)
## 
## Call:
## lm(formula = df1$y ~ df1$x1 + df1$x2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.8128 -1.3864  0.3817  1.1383  3.9672 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   2.9612     2.1319   1.389   0.1762  
## df1$x1        0.3172     0.3690   0.860   0.3975  
## df1$x2        0.3900     0.2275   1.714   0.0979 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.982 on 27 degrees of freedom
## Multiple R-squared:  0.6066, Adjusted R-squared:  0.5774 
## F-statistic: 20.81 on 2 and 27 DF,  p-value: 3.395e-06
reg2 = lm(df1$y ~ df1$x2)
summary(reg2)
## 
## Call:
## lm(formula = df1$y ~ df1$x2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.9417 -1.1617  0.2034  1.0659  4.1983 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.11197    1.65155   2.490    0.019 *  
## df1$x2       0.56998    0.08872   6.424 5.89e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.972 on 28 degrees of freedom
## Multiple R-squared:  0.5958, Adjusted R-squared:  0.5814 
## F-statistic: 41.27 on 1 and 28 DF,  p-value: 5.89e-07

Como vemos, en el segundo modelo encontramos una mayor explicatividad del \(R^2\) ajustado, es decir, el \(R^2\) que, siguiendo el principio de parsimonia, penaliza a nuestro modelo según aumentamos la complejidad del mismo (aumentamos variables). Pero no solo eso: al eliminar la variable \(x_1\) el coeficiente \(\beta_2\) pasa a explicar la relación conjunta de las dos variables independientes, alcanzando un mayor grado de significatividad en el modelo.