Regresión lineal múltiple

1.0 La regresión lineal múltiple: conceptos básicos

La regresión lineal múltiple es un método estadístico que se utiliza para modelar la relación entre dos o más variables explicativas (o independientes) y una variable respuesta (o dependiente). Se trata de una extensión de la regresión lineal simple, que sólo considera una variable explicativa.

En la regresión lineal múltiple, el modelo se puede expresar como:

\[Y = β_0 + β_1x_1 + β_2x_2 + … + β_nx_n + u \]

Donde:

- \(Y\) es la variable dependiente que estamos tratando de predecir o explicar.

- \(x_1, x_2, …, x_n\) son las variables independientes o explicativas que se utilizan para predecir \(Y\).

- \(β_0\) es el término del intercepto, que representa el valor esperado de Y cuando todas las variables explicativas son cero.

- \(β_1, β_2, …, β_n\) son los coeficientes de las variables explicativas, que representan el cambio esperado en \(Y\) para un cambio unitario en una variable explicativa, manteniendo constantes las demás variables.

- \(u\) es el término de error, que es una variable aleatoria que agrega el efecto de todas las otras variables que influyen en \(Y\) pero que no se incluyen en el modelo.

La regresión lineal múltiple puede ser una herramienta útil para entender la relación entre varias variables y puede ser utilizada en una variedad de contextos, incluyendo las ciencias sociales, la economía, la medicina y muchas otras áreas de la investigación.

En el apartado de regresión lineal simple vimos con la base de datos marketing del paquete de datos datarium el siguiente modelo:

\[ Ventas_i= \beta_0+ \beta_1Yt_i+ u_i \]

En donde \(Ventas_i\) es el nivel de ventas en miles de dólares y \(Yt_i\) es la inversión en publicidad en Youtube. En la base de datos existen otras variables que pueden explicar las ventas, como son la inversión en publicidad en Facebook, o en otros medios. Esto significa que la inversión en publicidad en Youtube no es el único determinante o factor explicativo de las ventas de la empresa. Si bien queda claro que los modelos son una simplificación de la teoría, ya que como menciona el filósofo y epistemólogo de la economía Uskali Mäki, la teoría económica sirve de marco de referencia para la identificación de las variables a introducir en un modelo, pero no todos los elementos que son parte del armazón lógico y sistemático de la teoría son necesarios para el modelamiento. Es por eso que el uso de supuestos con fines heurísticos (en jerga menos de la filosofía de la ciencia, supuestos irreales) tienen un rol fundamental de aislar teóricamente lo que se intenta observar mediante el modelo (Mäki, 2018). Es por esto que la ciencia económica (y me voy un poco más lejos, las ciencias sociales) es una ciencia de modelos, como menciona Dani Rodrik. Hay contextos en donde ciertos modelos no aplican debido a que las condiciones iniciales son distintas y por ende los supuestos deben de cambiar para estudiarlo (Rodrik, 2015).

Volviendo de nuevo al tema, el modelo de ventas sería más realista basándonos en esta ecuación:

\[ Ventas_i = \beta_0+\beta_1Yt_i+\beta_2Fb_i+\beta_3Per+u_i \]

En donde \(Fb_i\) y \(Per_i\) son la inversión en publicidad en Facebook y en periódico respectivamente. Una pregunta válida a hacerse sería: ¿Cambia la interpretación de los coeficientes? Ahora \(\beta_1\) es el efecto de la inversión en publicidad en youtube manteniendo la inversión en facebook y en períodico constante o ceteris paribus.

Si comparamos este modelo con la primera ecuación usada en el caso de la regresión lineal simple pasa algo curioso: ya que omitimos \(Fb\) y \(Per\) entonces \(\beta_1\) es la derivada del valor esperado de las ventas condicionada a la inversión en publicidad en Youtube:

\[ \beta_1= \frac{\delta[Ventas|Yt]}{\delta Yt} \]

Entonces,

\[ = \frac{\delta[Ventas|Yt, Fb, Per]}{\delta Yt} \color{red}{+\frac{\delta[Ventas|Yt, Fb, Per]}{\delta Fb}+ \frac{\delta Ventas|Yt, Fb, Per]}{\delta Per}} \]

El problema fundamental es que si no controlamos por por las otras variables relevantes, estaríamos estimando un efecto en \(\beta_1 Yt\) que no es real. Este es un problema que, de ser importante las demás variables en el modelo y en teoría y en base a la literatura empírica entendemos que sí, nos puede causar un sesgo por variables omitidas (el cual hablaremos más adelante). Esto rompe con el supuesto de ceteris paribus ya que al existir un problema de endogeneidad no podemos asegurar causalidad.

1.1 Algebra de Mínimos Cuadrados en el contexto de la RLM

Volviendo a la ecuación de a principios de la sección 1.0 \(Y = β_0 + β_1x_1 + β_2x_2 + … + βn*Xn + u\), entendemos que para cada observación tenemos:

\[ y_i= \beta_0+\beta_1x_{1,i}+\beta_{2}x_{2,i}+...+\beta_nx_{n,i} \]

Ahora utilizaremos un poco de notación matricial para así simplificar esta ecuación por fines de espacio y para mejor comprensión. \(x_{j,i}\) es un vector compuesto por los siguientes elementos, \(j\) que es la(s) variable(s) \((j=0, 1, 2…,K)\) e \(i\) que son las observaciones \((i=0, 1, 2,…,N)\). En un modelo de regresión tenemos \(K+1\) variables explicativas más una constante (con valor 1). Entonces \(x=(1,x_1, x_2,…x_K)\) es un vector de \(1 * (K+1)\) y \(\beta=(\beta_0, \beta_1, \beta_2,…,\beta_K)'\) es un vector \((K+1)* 1\), por ende podemos definir el modelo de la sección 1.1 en:

\[ y= \mathbf{X_i}\beta+u_i \]

Este estimador \(\hat{\beta}_{MCO}\) se define como:

\[ \hat{\beta}_{\text{MCO}} = \underset{\mathbf{b}\in \mathbb{R}^{K+1}}{\text{argmin}}\sum_{i=1}^{N}(y_i - b_0 - b_1x_{1,i} - \ldots - b_kx_{k,i})^2 \]

Para minimizar la suma de los errores al cuadrado se requiere realizar derivadas con respecto a \(\mathbf{b}\), el cual sabemos que es un vector \((K+1)*1\), el cual da las siguientes condiciones:

\[\sum^N_{i=1}x_{i,j}(y_i-\hat{\beta}_0-\hat{\beta_1}x_{1,i}-...-\hat{\beta_k}x_{k,i})=0, j= 0, 1, 2,...,K. \]

Similar a la regresión lineal simple, es posible pensar la regresión lineal múltiple como un método de momentos estadísticos:

Momentos en la población
\[ E[x_j,u]= E[x_j(y-\beta_0-\beta_1x_1-...-\beta_kx_k)=0 \] \[j=0,1,2..., K\](En donde \(x_0=1\) es una constante.
Momentos en la muestra
\[ N^{-1}\sum^N_{i=1}x_{j,i}(y_i-\hat{\beta}_0-\hat{\beta}_1x_{1,i}-...-\hat{\beta}_kx_{k,i})=0\] \(j=0,1, 2,...,K\)

Esto da como resultado un sistema de ecuaciones con \(K+1\) ecuaciones y \(K+1\) parámetros a estimar.

1.1.2 Los parámetros en la regresión lineal múltiple:

La forma de representar en términos matriciales la estimación de los \(\beta\) de una regresión lineal múltiple es:

\[ \hat{\beta}=(\mathbf{X'X})^{-1}\mathbf{X'y} \]

En donde:

• \(\mathbf{X}\) es una matríz de \(N *(K+1)\), del cual \(N\) contiene cada una de las \(i\) observaciones y en cada una de las \(K+1\) columnas las incógnitas a estimar. Recueden que \(x_{0,i}=1\), \(\forall\) \(i=1, 2, …, N\).

• \(\mathbf{y}\) es un vector \(N*1\) que contiene todas las observaciones de la variable dependiente.

• \((\mathbf{X'X})\) es una matriz \((K+1)*(N+1)\) y \('\) es la transpuesta de la matriz, que en este caso es \(\mathbf{X}\).

• \((\mathbf{X'X})^{-1}\) es la inversa de \((\mathbf{X'X})\)

• \((\mathbf{X'X})^{-1} \mathbf{X'y}\) es un vector \(K+1)*1\) en donde \(K\) es el número de variables independientes.

1.2 Minimizando errores: algebra de MCO.

Partiendo de nuevo de esta ecuación:

\[ \hat{\beta}=(\mathbf{X'X})^{-1}\mathbf{X'y} \]

La mínimización de los errores podemos plantearla como:

\[ \hat{\beta}= \underset{\mathbf{b}\in \mathbb{R}^{K+1}}{argmin} {\mathbf{b(b)'u(b)}}=\underset{\mathbf{b}\in \mathbb{R}^{K+1}}{argmin} (\mathbf{y-Xb})'(\mathbf{y-Xb)} \]

\[ = \underset{\mathbf{b}\in \mathbb{R}^{K+1}}{argmin}[\mathbf{y'y}+\mathbf{b'X'X b}-2\mathbf{b'X'y}] \]

Esto significa que \(\mathbf{u(b)}\equiv \mathbf{y-Xb}\). Si calculamos las derivadas con respecto a \(\mathbf{b}\in\mathbb{R}^{K+1}\) la solución termina siendo un vector de \((K+1)*1\) condiciones de primer orden. Esto nos lleva a que \(2\mathbf{X'Xb}-2\mathbf{X'y}=\mathbf{0}_{K+1}\) el cual es un vector (de tamaño \((K+1)*1\) ) que contiene ceros. Por último, es posible decir que \(\mathbf{X'y}= \mathbf{X'X b}\), \(\Rightarrow \hat{\beta}= (\mathbf{X'X})^{-1}\mathbf{X'y}\).

Las condiciones de segundo orden implicaría tomar la segunda derivada de la última ecuación, el cual debe de ser, por construcción una matriz positiva semidefinida, esto significa que todos sus autovalores son no negativos. Por lo que la función de mínimos cuadrados resulta ser convexa, por lo que el coeficiente puede estimarse de manera única y consistente.

2.0 Lo que se estima con una regresión lineal múltiple

Un concepto útil de la estadística que precisamente es la base de la regresión es que esta técnica busca estimar una esperanza condicional, es decir la ecuación de regresión lineal se puede expresar como \(E[\mathbf{y}|\mathbf{x}]=\beta_0+\beta_1x_1+…+\beta_kx_k\). Volviendo de nuevo a lo dicho en el apartado uno condicionar por una variable es hacerla fija (recordar el concepto de ceteris paribus y lo que implica como supuesto). La esperanza condicional de \(y\) dado \(x\) es una función estocástica en la cual la aleatoriedad de \(y\) fue técnicamente resuelta al condicionar, valga la redundancia, por las variables de interés. Ahora, eso supone que “sabemos” la forma funcional de dicha esperanza condicional.

Debido a esta propiedad podemos descomponer en dos partes que son independientes entre sí (ortogonales):

\[ y= E[y|\mathbf{x}]+u \]

En donde:

1. \(E[u\mathbf{x}]=0\)

2. \(E[h(\mathbf{x}u)]=0\) para cualquier función \(h(.)\)

3. Dado que \(E[u|\mathbf{x}]=0\), por la ley de esperanzas iteradas que establece que \(E_u(u)= E_x[E_u(u|x)]\), esto significa que por consecuencia \(E[u]=0\)

4. Además, debido a la misma ley de esperanzas iteradas si \(E[u|\mathbf{x}]=0\), entonces \(E[xu]=0\). Si multiplicamos ambos lados de la igualdad de la primera expresión por \(x\) nos da que:

   \[ \mathbf{x}E(u|\mathbf{x})= \mathbf{x}0 \]

Esto es lo mismo que decir que \(E(xu|x)=0\), por ende \(E[E(\mathbf{x}u|\mathbf{x})]= E(xu)=E(0)=0\). Este supuesto significa que los errores \(\mathbf{u}\) de la regresión múltiple, al igual que en la regresión lineal simple, no guardan relación o son ortogonales a \(\mathbf{x}\).

3.0 Supuestos del teorema de Gauss-Markov

Al igual que en el modelo de regresión lineal simple los supuestos son:

1. Linearidad en los parámetros: El modelo es una combinación lineal de los parámetros desconocidos. Por ejemplo, el modelo de regresión simple sería \(Y = β_0 + β_1x_1+\beta_2x_2+...+ u\).
2. Muestra aleatoria: \([(y_i, x_{1,i},...,x_{k,i})]^N_{i=1}\) es una muestra aleatoria del modelo descrito en el primer supuesto.
3. Esperanza condicional cero: La expectativa condicional del término de error dado los regresores es cero. Esto se escribe matemáticamente como \(E(u|x) = 0\). Esto implica que el error es independiente de los regresores.
4. No multicolinealidad perfecta: Los regresores no son perfectamente multicolineales. En otras palabras, no existe una relación lineal perfecta entre las variables independientes. Esto se expresa como que \((X'X)\) sea una matriz no singular (que el rango de la matriz sea \(K+1\)). Que dicha matriz es singular significa que es el equivalente a dividir entre 0 la covarianza entre las \(x\) y \(y\).
5. Homoscedasticidad: La varianza del término de error es constante dado los regresores. Esto se escribe matemáticamente como \(Var(\mathbf{u}|\mathbf{x}) = σ^2\mathbf{I}_N\). Cuando este supuesto no se cumple, tenemos heteroscedasticidad.
6. No autocorrelación serial: Los términos de error no están autocorrelacionados.

Bajo los supuestos 1-5 (en el contexto de análisis de regresión para corte transversal) el teorema de Gauss-Markov ilustra que Mínimos Cuadrados Ordinarios \((\hat{\beta}_0, \hat{\beta}_1,...,\hat{\beta}_k)\) son los mejores estimadores lineales insesgados (MELI, BLUE).

3.1 Insesgadez de MCO

Si la Esperanza condcional de el vector de \(\hat{\beta_j}\) condicional en las \(\mathbf{X}\) es insesgado se debe a que:

\[E[\hat{\beta}|\mathbf{X}]= E[(\mathbf{X'X)^{-1}X'y|X}]\]

Luego incluimos el supuesto de linealidad en los parámetros, sustituyendo al vector \(\mathbf{y}\) por \(\mathbf{X}\beta+\mathbf{u}\):

\[E[\hat{\beta}|\mathbf{X}]= E[(\mathbf{X'X)^{-1}X'\mathbf{X}(\beta+\mathbf{u})|X}]\]

Luego aplicamos la propiedad distributiva e introducimos el supuesto de que \(E[\mathbf{u|X}]=E[\mathbf{u}]=0\):

\[=E[(\mathbf{X'X})^{-1}(\mathbf{X'X})\beta|\mathbf{X}]+ E[(\mathbf{X'X})^{-1}\mathbf{X'u|X}]\]

Dado que dijimos que la esperanza condicional de \(\mathbf{u}\) es 0, así como la expresión \((X'X)^{-1}(X'X)\) se cancelan, po lo que nos queda que todo es igual a:

\[E[\hat{\beta}|\mathbf{X}]= \beta\]

3.2 Simulación para probar la insesgadez de MCO:

# Cargamos las librerías necesarias
set.seed(1234)
library(boot)

# Definir los parámetros de la simulación
N <- 1000 # tamaño de la muestra
M <- 1000 # número de simulaciones
b0 <- 2  # valor verdadero del intercepto
b1 <- 4  # valor verdadero de beta1
b2 <- 8 #Valor del verdadero de beta2

# Definir la función de simulación
regsim <- function(data, indices) {
  # Generamos la muestra
  x <- rnorm(N, 0, 1)
  x2 <- rnorm(N, 0, 1) # x2 independiente de x
  u <- rnorm(N, 0, 1)
  y <- b0 + b1*x +b2*x2+ u
  
  # Ajustamos el modelo
  model <- lm(y ~ x + x2)
  
  # Devolvemos los coeficientes
  return(coef(model))
}

# Generamos los datos originales
data <- data.frame(x = rnorm(N), x2 = rnorm(N), u = rnorm(N))
data$y <- b0 + b1*data$x +b2*data$x2+ data$u

# Realizamos la simulación
results <- boot(data, regsim, R = M)

# Creamos un dataframe con los coeficientes de cada simulación
coef_df <- as.data.frame(t(results$t))
coeft_df <- as.data.frame(t(coef_df)) #transponer el dataframe:
# Renombrar nombre de los vectores:
colnames(coeft_df) <- c("beta0", "beta1", "beta2")


# Gráficar los resultados:
library(ggplot2)
library(gridExtra)



beta0<-ggplot(coeft_df,aes(y=beta0))+geom_histogram(fill= "blue")+
  geom_hline(yintercept = 2, color= "red", linetype= "dashed")+
  labs(y=expression(beta[0]), x="Recuento")+theme_light()

beta1<-ggplot(coeft_df, aes(y=beta1))+ geom_histogram(fill= "navy")+
  geom_hline(yintercept = 4, color= "darkred", linetype= "dashed")+
  labs(y=expression(beta[1]), x= "Recuento")+theme_light()
beta2<-ggplot(coeft_df, aes(y=beta2))+
  geom_histogram(fill= "lightblue")+
  geom_hline(yintercept = 8,color="orange",linetype= "dashed")+
  labs(y=expression(beta[2]), x="Recuento")+theme_light()


grid.arrange(beta0, beta1, beta2, nrow=1)

4.0 Varianza de Mínimos Cuadrados.

En notación matricial podemos decir que la Varianza de los coeficientes es:

\[ Var(\hat{\beta}|\mathbf{x})= \sigma^2(\mathbf{X'X})^{-1} \]

La prueba de ello es que \(Var(\hat{\beta}|\mathbf{x})= Var[\mathbf{(X'X)^{-1} X'y}|\mathbf{X}]\)

Luego saco la expresión \((\mathbf{X'X)^{-1}X'}\) y luego introduzco post-multiplicando la expresión \(\mathbf{X(X'X)^{-1}}\), quedando así:

\[ = (\mathbf{X'X})^{-1}\mathbf{X} 'Var[\mathbf{y|X}]\mathbf{X}(\mathbf{X'X})^{-1} \]

Ya que \(Var[\mathbf{y|X}]\) es \(\sigma^2\mathbf{I_N}\) ya que es lo mismo que decir que \(Var[\mathbf{X}\beta+\mathbf{u'|X}]\) debido al supuesto de Homocetasticidad.

Por esto significa que:

\[ = (\mathbf{X'X})^{-1}\mathbf{X} '\sigma^2\mathbf{X}(\mathbf{X'X})^{-1} \] Lo cual podemos arreglarle de esta manera:

\[ = \sigma^2 (\mathbf{X'X})^{-1}\color{red}{\mathbf{X'X(X'X)^{-1}}} \]

La múltiplicación de lo que esta en rojo da como resultado la matríz de identidad \(\mathbf{I}\) , por lo que¹:

\[= \sigma^2(\mathbf{X'X})^{-1}\]

4.1 Simulación de la varianza de MCO

# Definimos el tamaño de la muestra
N <- 1000  

# Generamos las variables independientes y la variable dependiente
x1 <- rnorm(N, 0, 1)
x2 <- rnorm(N, 0, 1) + x1
u <- rnorm(N, 0, 1)
y <- 2 + 4*x1 + 8*x2 + u

# Creamos un dataframe con nuestros datos
data <- data.frame(y = y, x1 = x1, x2 = x2)

# Realizamos la regresión
model <- lm(y ~ x1 + x2, data = data)

# Imprimimos los coeficientes
print(coef(model))

## (Intercept)          x1          x2 
##    2.031700    3.996606    7.976602

# Imprimimos la matriz de varianzas y covarianzas de los estimadores
print(vcov(model))

##               (Intercept)            x1            x2
## (Intercept)  1.022507e-03 -8.477019e-06  8.272123e-06
## x1          -8.477019e-06  2.101960e-03 -1.021665e-03
## x2           8.272123e-06 -1.021665e-03  1.060326e-03

5.0 Ejercicio práctico

Para este ejercicio práctico le recomiendo descargar estos paquetes si no los tiene instalado:

list.of.packages <- c("plot3D",

"rsq",

"heplots",

"caret",

"MASS",

"leaps",

"car",

"relaimpo",

"hier.part")

new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]

data("marketing", package = "datarium")
head(marketing)

##   youtube facebook newspaper sales
## 1  276.12    45.36     83.04 26.52
## 2   53.40    47.16     54.12 12.48
## 3   20.64    55.08     83.16 11.16
## 4  181.80    49.56     70.20 22.20
## 5  216.96    12.96     70.08 15.48
## 6   10.44    58.68     90.00  8.64

Para representar gráficamente modelos de RLM podemos usar el paquete plot3D

library(plot3D)
scatter3D(x=marketing$youtube,
          y=marketing$facebook,
          z=marketing$sales,
          colvar = marketing$newspaper)

5.1 Ajustar el modelo de regresión

mod1<-lm(sales ~youtube, data=marketing)
mod2<-lm(sales ~youtube+facebook, data=marketing)
mod4<-lm(sales~newspaper, data=marketing)
mod3<-lm(sales~youtube+facebook+newspaper, data=marketing)

#Cargamos librería:
library(stargazer)
tabla<-stargazer(mod1, mod2, mod4, mod3, title="Tabla 3. Resultados de la regresión", align=TRUE, type="html", covariate.labels = c("Inversión Youtube","Inversión en Facebook", "Inversión en periodico", "(Intercepto)"), dep.var.labels = c("Ventas"), out = "Regresión_m.html")

Tabla 3. Resultados de la regresión

	Dependent variable:
	Ventas
	(1)	(2)	(3)	(4)
Inversión Youtube	0.048***	0.046***		0.046***
	(0.003)	(0.001)		(0.001)
Inversión en Facebook		0.188***		0.189***
		(0.008)		(0.009)
Inversión en periodico			0.055***	-0.001
			(0.017)	(0.006)
(Intercepto)	8.439***	3.505***	14.822***	3.527***
	(0.549)	(0.353)	(0.746)	(0.374)
Observations	200	200	200	200
R2	0.612	0.897	0.052	0.897
Adjusted R2	0.610	0.896	0.047	0.896
Residual Std. Error	3.910 (df = 198)	2.018 (df = 197)	6.111 (df = 198)	2.023 (df = 196)
F Statistic	312.145*** (df = 1; 198)	859.618*** (df = 2; 197)	10.887*** (df = 1; 198)	570.271*** (df = 3; 196)
Note:	p<0.1; p<0.05; p<0.01

La tasa del error que es definida como \(\hat{\sigma}^2/\bar{y}\) es una medida de bondad de ajuste, la cual se calcula con el siguiente comando:

sigma(mod3)/mean(marketing$sales)

## [1] 0.1202004

Esto significa que la tasa de error es de un 12.02 %.

Por otro lado, podemos observar el ajuste parcial de cada predictor (que tanto cada uno contribuye a la explicación de la varianza:

Esto indica que \(R^2_{Yt}= 84.6\%\), \(R^2_{Fb}=70.9\%\) y \(R^2_{Per}=0.16\%\)

library(rsq)
rsq.partial(mod3)

## $adjustment
## [1] FALSE
## 
## $variable
## [1] "youtube"   "facebook"  "newspaper"
## 
## $partial.rsq
## [1] 0.8459610964 0.7097694432 0.0001593014

5.2 Comparación de modelos

Existen tres enfoques principales para la comparación de modelos:

1. Métodos basados en criterios de información.

2. Pruebas F parciales.

3. Validación cruzada.

Los métodos basados en criterios de información evalúan si la adición de nuevos predictores o interacciones mejora el ajuste del modelo utilizando alguna métrica de rendimiento como los criterios de información AIC, BIC o Cp. En R, podemos hacer esto con la función AIC(model_all, model_red).

Las pruebas F parciales comparan modelos anidados para evaluar la significancia de cada predictor. Estas pruebas nos proporcionan una tabla de análisis de varianza (ANOVA) y podemos realizarlas en R con la función anova(model_red, model_all). También existe la opción de utilizar la función Anova() del paquete car. Aquí la hipótesis nula y alternativa son la siguientes:

• \(H_0\) : \(\beta_i=0| \beta_js \neq 0\) Modelo reducido.

• \(H_1\) : \(\beta_i \neq 0|\beta_js \neq0\) Modelo completo.

Como vemos, la inversión en publicidad no genera un efecto significativo en la predicción de las ventas:

library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:boot':
## 
##     logit

anova(mod3, mod2)

## Analysis of Variance Table
## 
## Model 1: sales ~ youtube + facebook + newspaper
## Model 2: sales ~ youtube + facebook
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1    196 801.83                           
## 2    197 801.96 -1  -0.12775 0.0312 0.8599

La validación cruzada permite comparar modelos en función de su precisión en las predicciones.

5.3 Selección de predictores en base a criterios de información:

Cuando tenemos muchas variables explicativas, puede ser difícil explorar todas las combinaciones para comparar modelos. En estos casos, es útil utilizar métodos más eficientes para seleccionar los predictores del modelo. Los dos métodos principales de selección de variables son la selección por pasos (stepwise) y la selección por subconjuntos.

Al final la función de log-verosimilitud del modelo \(ventas= \beta_0+\beta_1Yt+\beta_2 Fb\) posee más parsimonia, ya que el valor de AIC es menor, por lo que pierde menor

AIC(mod3, mod2)

##      df      AIC
## mod3  5 855.2909
## mod2  4 853.3227

La selección por pasos añade o elimina automáticamente cada variable explicativa, paso a paso, en base a un criterio seleccionado, dando como resultado un solo modelo de regresión. Esto se puede hacer en R utilizando la función stepAIC() . Como observamos anteriormente, el modelo 2 de la tabla de regresión resulta ser el seleccionado.

library(MASS)
step.model<-stepAIC(mod3,
                    direction = "backward",
                    trace= "F")
summary(step.model)

## 
## Call:
## lm(formula = sales ~ youtube + facebook, data = marketing)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.5572  -1.0502   0.2906   1.4049   3.3994 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.50532    0.35339   9.919   <2e-16 ***
## youtube      0.04575    0.00139  32.909   <2e-16 ***
## facebook     0.18799    0.00804  23.382   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.018 on 197 degrees of freedom
## Multiple R-squared:  0.8972, Adjusted R-squared:  0.8962 
## F-statistic: 859.6 on 2 and 197 DF,  p-value: < 2.2e-16

La selección por subconjuntos compara todos los modelos posibles para un número específico de predictores y muestra los mejores modelos en cada caso. El resultado es una serie de modelos y sus estadísticas de ajuste, y el investigador debe decidir cuál modelo elegir. Esto se puede hacer en R con la función regsubsets().

library(leaps)
subsets<- regsubsets(as.matrix(marketing[-4]), #Ingresamos la matriz de preodictores
marketing[,"sales"], #La variable a explicar
            nbest=1, #un solo mejor modelo para cada subconjunto.
            nvmax= NULL, #número máximo de variables 
            method= "backward") #Método de selección

#Resultados:
summary(subsets)

## Subset selection object
## 3 Variables  (and intercept)
##           Forced in Forced out
## youtube       FALSE      FALSE
## facebook      FALSE      FALSE
## newspaper     FALSE      FALSE
## 1 subsets of each size up to 3
## Selection Algorithm: backward
##          youtube facebook newspaper
## 1  ( 1 ) "*"     " "      " "      
## 2  ( 1 ) "*"     "*"      " "      
## 3  ( 1 ) "*"     "*"      "*"

Ahora debemos de ver la parsimonia de dichos modelos para cada subcojunto usando summary(subsets)$BIC , es decir, del criterio de información Bayesiano (BIC). El mejor modelo del subconjunto 2 (el que incluye el preictor Youtube y Facebook) resulta ser más parsimonioso.

summary(subsets)$bic

## [1] -178.6890 -439.0879 -433.8214

5.4 Una aclaración importante con métodos de selección de variables.

Es mejor hacer uso de ellos con fines exploratorios, esto es porque en la práctica empírica solo podemos guiarnos de la teoría y de la literatura empírica sobre cómo especificar el modelo. El hacer uso de estás técnicas con fines de predicción sin ningún criterio teórico-empírico-metodológico puede llevarte a crear un modelo quizás parsimonioso pero muy alejado del “true model”, es decir, que ocurra sesgo por variables omitidas. Y en econometría, como en métodos cuantitativos es preferible tener un modelo algo saturado pero correctamente especificado que uno en donde rompamos ese supuesto de ceteris paribus que mencionamos en la sección 1.0. Pero también, colocar un montón de variables puede llevarte a sobreajustar el modelo (overfitting en inglés).

6.0 Qué hacer cuando hay multicolinealidad

Las variables explicativas a incluir en el modelo deben cumplir con las siguientes condiciones:

• No deben ser redundantes. La multicolinealidad ocurre cuando algunas de las variables explicativas del modelo de regresión están correlacionadas. Si el grado de correlación entre las variables es muy alto, puede causar problemas al ajustar el modelo e interpretar los resultados. Esto se debe a que será imposible aislar la relación entre cada variable explicativa y la variable respuesta, porque las variables explicativas tenderán a cambiar de manera conjunta. En el siguiente apartado veremos cómo medimos la redundancia entre las variable explicativas.

• Deben tener una cierta justificación teórica. Basa tu elección de las variables explicativas en los modelos propuestos por la literatura de tu área de estudio, así tendrás menos posibilidades de cometer errores.

• El tamaño de la muestra debe ser grande. Otro requisito que se utiliza como “regla tácita”, es que debemos tener, además de los 30 casos en general, un mínimo de 10 casos por cada variable adicional que incluyamos en el RLM.

• Deben tener una relación lineal con la respuesta. La relación de las variables explicativas con la variable respuesta debe de ser lineal, es decir, proporcional.

Existen varias estrategias para tratar la multicolinealidad:

1. Eliminar variables: Si dos variables están altamente correlacionadas, a menudo tiene sentido eliminar una de ellas del modelo.

2. Combinar variables: Si varias variables están altamente correlacionadas, a veces se pueden combinar en una sola variable que represente la información compartida.

3. Regularización: Los métodos de regularización, como la regresión de crestas (Ridge) y la regresión de lazo (Lasso), pueden ayudar a lidiar con la multicolinealidad al reducir la magnitud de los coeficientes.

4. Análisis de componentes principales (PCA): La PCA es una técnica de reducción de la dimensionalidad que puede usarse para transformar un conjunto de variables correlacionadas en un conjunto más pequeño de variables no correlacionadas.

Cada una de estas técnicas tiene sus propias ventajas y desventajas, por lo que la elección de la estrategia más adecuada dependerá del contexto específico.

Para medir la multicolinealidad entre los predictores de un modelo se evalua usando el factor de inflación de varianza (VIF). Esto mide cuánto se infla la varianza de un coeficiente de regresión debido a la multicolinealidad del modelo. Aunque no existen reglas estadísticas robustas y hay ciertas discusiones alrededor, una sugerencia es que un VIF de 10 es preocupante mientras que un valor de 5 es ya no recomendable. La fórmula para estimar el VIF parte de comprender otra fórmula, la tolerancia del modelo y se representa así:

\[ T_i= \frac{1}{VIF}= 1-R^2_i \]

La tolerancia mide la correlación múltiple de un determinado predictor \(i\) con los demás. Si hay multicolinealidad, por lo general la tolerancia es baja (correlación alta). Entonces el Factor de Inflación de Varianza no es más que un ratio la tolerancia del modelo:

\[ VIF_i= \frac{1}{T_1}= \frac{1}{1-R^2_i} \]

Para poder estimar si el modelo tiene algún problema de multicolinealidad se puede estimar usando el fomando vif() y al parecer el modelo con todas las variables no posee problemas de multicolinealidad.

vif(mod3)

##   youtube  facebook newspaper 
##  1.004611  1.144952  1.145187

7.0 Contribución de cada predictor al \(R^2\)

Si definimos la variabilidad explicada como \(1-\frac{SS_{res}}{SS_{total}}\), se asume que la introducción de predictores contribuirá a explicar cada vez más \(SS_{total}\). Ahora bien, el tamaño de los coeficientes individuales no necesariamente implica que un predictor de los \(K\) predictores restantes dentro de un modelo de regresión lineal múltiple es el que más contribuye a la varianza explicada. Es por ello que es necesario encontrar técnicas que descompongan la varianza o ver los cambios en el \(R^2\) al introducir las variables. Ahora bien, hay otras manera de ver el cambio medio en el \(R^2\) en Rstudio como es el comando crlm() del paquete relaimpo. Aquí observamos que de la varianza total explicada por el modelo (89.72 %) un 65 % pertenece a \(\beta_1 Yt\), un 32 % a \(\beta_2 Fb\) y un 3 % a \(\beta_3 Per\). Ahora bien, para pensar el tamaño del efecto en términos prácticos se requiere de conocimiento técnico del área, o más bien, conocimieno teórico para interpretar adecuadamente su significancia práctica. Por eso los detalle acerca del área de tu investigación y los objetivos que te planteaste, así como la calidad de los datos y su colección te permitirán tomar mejores decisiones y recomendaciones.

library(relaimpo)
crlm<-calc.relimp(mod3, 
                  type=c("lmg"), #Tipo de modelo: Modelo Lineal General= lmg
                  rela= T)
crlm

## Response variable: sales 
## Total response variance: 39.19947 
## Analysis based on 200 observations 
## 
## 3 Regressors: 
## youtube facebook newspaper 
## Proportion of variance explained by model: 89.72%
## Metrics are normalized to sum to 100% (rela=TRUE). 
## 
## Relative importance metrics: 
## 
##                  lmg
## youtube   0.65232505
## facebook  0.32187149
## newspaper 0.02580346
## 
## Average coefficients for different model sizes: 
## 
##                   1X        2Xs          3Xs
## youtube   0.04753664 0.04632801  0.045764645
## facebook  0.20249578 0.19351941  0.188530017
## newspaper 0.05469310 0.02543180 -0.001037493

8.0 Predictores categóricos en modelos de regresión

En un modelo de regresión, es posible incluir predictores categóricos además de las variables numéricas. Esto se hace a través de un proceso llamado “codificación de variables” o “creación de variables dummy”.

Las variables categóricas son aquellas que tienen un número limitado y usualmente fijo de posibles valores, también conocidos como niveles. Por ejemplo, el sexo (hombre, mujer), color de ojos (azul, verde, marrón), tipo de vivienda (casa, apartamento, condominio) son todos ejemplos de variables categóricas.

Incluir estas variables en un modelo de regresión puede aportar información valiosa, ya que permiten que el modelo tome en cuenta los efectos de estas categorías en la variable de respuesta.

Sin embargo, estas variables categóricas deben transformarse en un formato que el modelo de regresión pueda entender. Esto se hace a través de un proceso llamado “codificación de variables” o “creación de variables dummy”.

En la codificación de variables, cada nivel de la variable categórica se convierte en una nueva variable binaria (0/1). Por ejemplo, si tenemos una variable categórica “color de ojos” con niveles “azul”, “verde” y “marrón”, creamos tres nuevas variables binarias: “color de ojos es azul”, “color de ojos es verde” y “color de ojos es marrón”. Cada una de estas nuevas variables toma el valor 1 si el color de ojos es el correspondiente, y 0 en caso contrario.

Es importante tener en cuenta la “trampa de la dummy”, que ocurre cuando las variables dummy son multicolineales, es decir, cuando una variable puede predecirse perfectamente a partir de las otras. Para evitar esta trampa, generalmente se excluye una de las variables dummy, que actúa como la categoría de referencia contra la cual se comparan las demás.

Al incluir variables categóricas en el modelo de regresión, podemos capturar los efectos diferenciados de cada categoría en la variable dependiente, aumentando así la precisión y utilidad del modelo.

8.1 Caso práctico: variables cualitativas

En la base de datos Salaries del paquete de datos carData encontramos una serie de variables categóricas de distintos niveles. Dentro de ella tenemos el rank, la disciplina, años desde que logró su Ph.D, años en el servcio, su sexo y el salario

data("Salaries", package= "carData")

head(Salaries)

##        rank discipline yrs.since.phd yrs.service  sex salary
## 1      Prof          B            19          18 Male 139750
## 2      Prof          B            20          16 Male 173200
## 3  AsstProf          B             4           3 Male  79750
## 4      Prof          B            45          39 Male 115000
## 5      Prof          B            40          41 Male 141500
## 6 AssocProf          B             6           6 Male  97000

Ahora un pequeño recuento de las estadísticas descriptivas de la base de datos:

summary(Salaries)

##         rank     discipline yrs.since.phd    yrs.service        sex     
##  AsstProf : 67   A:181      Min.   : 1.00   Min.   : 0.00   Female: 39  
##  AssocProf: 64   B:216      1st Qu.:12.00   1st Qu.: 7.00   Male  :358  
##  Prof     :266              Median :21.00   Median :16.00               
##                             Mean   :22.31   Mean   :17.61               
##                             3rd Qu.:32.00   3rd Qu.:27.00               
##                             Max.   :56.00   Max.   :60.00               
##      salary      
##  Min.   : 57800  
##  1st Qu.: 91000  
##  Median :107300  
##  Mean   :113706  
##  3rd Qu.:134185  
##  Max.   :231545

En este caso hay algunas variables que son de formato texto que es necesario, para construir el modelo, transformarlas. Para ello tomaremos unas reglas sencillas:

• Toma los valores:

  1. 0= Ausencia de atributo.
  2. 1= Presencia del atributo.

• Evitar la trampa de las variables dummies. Siendo \(\delta_i\) un vector que contiene \(k\) categorías, el número de categorías óptimas en una regresión a colocar de una variable que deseamos introducir sera:

\[\delta_i= Cat_i-1\] Esto es porque si la introducimos todas, digamos en un ejemplo basado en la base de datos Salaries

\[Salario_i= \beta_0+\delta_1 Mujer_i+\delta_2 Hombre_i+\beta_1Esc.Phd_i+\beta_2 Exp.Servicio_i+\gamma_1 Discip_i+u_i\]

Este modelo no será posible porque el intercepto y las dos variables dummies de hombre y mujer poseen multicolinealidad perfecta, por lo que hay que eliminar una.

data(Salaries)
contrasts(Salaries$sex) # con esto recodificamos numéricamente

##        Male
## Female    0
## Male      1

#Si quiero cambiar la categoría de referencia:

Salaries$sex2<-relevel(Salaries$sex, ref="Male")
#ahora aplicamos el comando contrast:
contrasts(Salaries$sex2)

##        Female
## Male        0
## Female      1

8.1.1 Interpretando los coeficientes de un modelo de regresión lineal con variables cualitativas:

m1<-lm(salary~ sex, data= Salaries)
m2<-lm(salary~sex+ yrs.since.phd, data=Salaries)
m3<-lm(salary~sex+yrs.since.phd+ yrs.service+ discipline, data=Salaries)

#Calcula términos cuadrático de los año de servicios para modelar no linealidades:
Salaries$yrs.service2<-(Salaries$yrs.service)^2

m4<-lm(salary~sex+yrs.since.phd+yrs.service+ yrs.service2+discipline, data=Salaries)
m5<- lm(log(salary)~sex+yrs.since.phd+yrs.service+yrs.service2+discipline, data=Salaries)

stargazer(m1, m2, m3, m4, m5, title="Tabla 4. Resultados de la regresión", align=TRUE, type="html", covariate.labels = c("Sexo","# Años Ph.D", "# Años Servicio","# años servicio²", "Disciplina",  "(Intercepto)"), dep.var.labels = c("Salario", "Log(Salario)"), out = "Regresión_m2.html")

Tabla 4. Resultados de la regresión

	Dependent variable:
	Salario				Log(Salario)
	(1)	(2)	(3)	(4)	(5)
Sexo	14,088.010***	7,923.624*	7,545.313*	7,963.036*	0.076**
	(5,064.579)	(4,684.081)	(4,462.598)	(4,250.040)	(0.035)
# Años Ph.D		958.079***	1,804.145***	1,526.075***	0.013***
		(108.319)	(248.867)	(240.904)	(0.002)
# Años Servicio			-770.100***	1,356.337***	0.015***
			(244.106)	(404.401)	(0.003)
# años servicio²				-43.120***	-0.0004***
				(6.710)	(0.0001)
Disciplina			16,325.420***	15,558.920***	0.143***
			(2,708.903)	(2,582.330)	(0.021)
(Intercepto)	101,002.400***	85,181.820***	71,325.790***	60,769.750***	11.119***
	(4,809.386)	(4,748.315)	(4,981.841)	(5,020.367)	(0.042)
Observations	397	397	397	397	397
R2	0.019	0.182	0.263	0.334	0.381
Adjusted R2	0.017	0.178	0.256	0.325	0.373
Residual Std. Error	30,034.610 (df = 395)	27,468.940 (df = 394)	26,128.800 (df = 392)	24,881.350 (df = 391)	0.206 (df = 391)
F Statistic	7.738*** (df = 1; 395)	43.742*** (df = 2; 394)	35.035*** (df = 4; 392)	39.168*** (df = 5; 391)	48.083*** (df = 5; 391)
Note:	p<0.1; p<0.05; p<0.01

Recordar hacer los diagnósticos del modelo para observar qué correcciones hay que hacer

plot(m4) # Existen valores influyentes.

vif(m4) # Hay que centrar la variable ya que VIF >5.

##           sex yrs.since.phd   yrs.service  yrs.service2    discipline 
##      1.026113      6.165100     17.695416     10.316148      1.060760

library(lmtest) 

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

bptest(m4) #La prueba de Breuch-Pagan permite testear si existe Heterocedasticidad. H0: Var(u|x)=0 ; H1: Var(u|x) =/=0. En este caso se rechaza H0

## 
##  studentized Breusch-Pagan test
## 
## data:  m4
## BP = 67.505, df = 5, p-value = 3.385e-13

sigma(m4)/mean(Salaries$salary) # La tasa de error de predicción es  22.9 %

## [1] 0.2188209

##############
# Modelo Log #
##############

plot(m5) # Existen valores influyentes y datos atípicos.

vif(m5) # Hay que centrar la variable ya que VIF >5.

##           sex yrs.since.phd   yrs.service  yrs.service2    discipline 
##      1.026113      6.165100     17.695416     10.316148      1.060760

bptest(m5) # Se rechaza H0

## 
##  studentized Breusch-Pagan test
## 
## data:  m5
## BP = 72.674, df = 5, p-value = 2.843e-14

library(rsq)


sigma(m5)/mean(log(Salaries$salary)) # La tasa de error de predicción es  1.8 %

## [1] 0.01775133

A su vez, la mujer gana en promedio en el modelo (4) 7,963 dólares menos en promedio que los doctorantes hombres, aún ajustando por años desde que consiguió su Ph.D, Años de servicio y los rendimientos decrecientes de los años de servicios y la disciplina a la que pertenece. En promedio las mujeres ganan 105,840 dólares al año frente a 113,803 dólares en promedio que los hombres.

Por otro lado, existen diferencias singificativas en la disciplina ya que los de disciplina A ganan 102,042 frente a 117,601 dólares anuales en promedio de la disciplina B.

library(emmeans) #Descargar usando install.packages("emmeans)
emmeans(m4, "sex")

##  sex    emmean   SE  df lower.CL upper.CL
##  Female 105840 4029 391    97919   113760
##  Male   113803 1323 391   111202   116404
## 
## Results are averaged over the levels of: discipline 
## Confidence level used: 0.95

emmeans(m4, "discipline")

##  discipline emmean   SE  df lower.CL upper.CL
##  A          102042 2507 391    97114   106970
##  B          117601 2450 391   112785   122417
## 
## Results are averaged over the levels of: sex 
## Confidence level used: 0.95

Referencias

Mäki, U. (2018). Rights and wrongs of economic modelling: refining Rodrik. Journal of Economic Methodology, 25(3), 218–236. https://doi.org/10.1080/1350178x.2018.1488475

Rodrik, D. (2015). Economics rules: The rights and wrongs of the dismal science. WW Norton & Company.

---

1. Cabe recordar que si definimos a \(Var(\mathbf{u|X})= \Omega\) podemos decir que:

   \[ \Omega= E[\mathbf{uu'|X}]+E[\mathbf{u|X}]E[\mathbf{u'|X}]= E[\mathbf{uu'|X}] \]
   Esta matriz \(\Omega\) nos permite construir la Varianza de los \(\beta\) en forma de sandwich. Esto se debe a que podemos escribir la varianza como la esperanza a cuadrado de \(\mathbf{u}\) menos el cuadrado de las esperanzas. Esto nos permite re-escribir \(Var[\hat{\beta}|\mathbf{X}])\) en forma de “sandwich” o de una matriz de varianzas y covarianzas.

   \[ = (\mathbf{X'X})^{-1}\mathbf{X' \Omega X}(\mathbf{X'X})^{-1} \]

   Esto es debido a que por construcción si la estimación de la varianza es la varianza de los errores al cuadrado \(E[\mathbf{uu'|X}]\), esto es equivalente a \(\sigma^2 \mathbf{I_N}\).↩︎