La regresión lineal es un modelo matemático que describe la relación entre varias variables. Es una técnica de modelado estadístico que se emplea para describir una variable de respuesta continua como una función de una o varias variables predictoras. Puede ayudar a comprender y predecir el comportamiento de sistemas complejos o a analizar datos experimentales, financieros y biológicos.
En resumen, la regresión lineal simple examina la relación lineal entre dos variables continuas: una respuesta \((Y)\) y un predictor \((X)\). Cuando las dos variables están relacionadas, es posible predecir un valor de respuesta a partir de un valor predictor con una exactitud mayor que la asociada únicamente a las probabilidad. Aunque hay una serie de técnicas en las ciencias sociales y en otras disciplinas, inclusive algunas más complejas como las de Machine Learning, o en el contexto de la econometría y la modelización cuantitativa en ciencias sociales, que es mi área de experticia en donde he usado estas herramientas y varaciones (Mella Gómez, 2022), la regresión lineal es la columna vertebral de todos.
Podemos expresar la regresión con la siguiente ecuación:
\[Y_i = \beta_0+ \beta_1x_i+u_i; i= 1,2,...,N\]
En donde:
Muestra o datos \([x_i, y_i]^N_{i=1}= [(x_1, y_1), (x_2, y_2),...,(x_N,y_N)]\), muestra de tamaño \(N\).
Modelo lineal \(y= \beta_0+\beta_1x\) :
\(y\) es la variable dependiente, o variable a explicar (outcome variable).
\(x\) es la variable independiente/control/explicativa.
\(\beta_0\) es el intercepto, el valor de y cuand x=0.
\(\beta_1\) es \(\frac{\Delta y}{\Delta x}\), la pendiente. Es el cambio de la variable dependiente ante un incremento de la variable independiente en una unidad.
\(u_i\) es el término del error.
Dependiendo de la estructura de los datos, de la precisión y de el tipo de datos, se requiere estimar los parámetros que son propios del modelo. Uno de los más comunes es el método de mínimos cuadrados ordinarios (OLS en inglés o MCO en español).
La estimación por mínimos cuadrados ordinarios busca de alguna manera construir una recta (es decir, estimar el \(\beta_0\) y \(\beta_1\)) alrededor de la nube de puntos entre las variables \(x_i\) y \(y_i\) que mínimiza los errores al cuadrado.
En teoría lo que hacemos es tomar los residuos \(u_i\) (que no son observables pero que podemos estimar por medio de la ecuación) \(u_i \equiv y_i-\beta_0-\beta_1x_i,\) \(i= 1, 2,...,n\). Luego los elevamos al cuadrado:
\[ \sum^n_i u^2_i= \sum^N_i(y_i-\beta_0-\beta_ix_i)^2 \]
Se obtiene el \(\beta_0\) y a \(\beta_1\) que minimiza la expresión \(\sum^N_i(y_i-\beta_0-\beta_1x_i)^2\)
Una forma de ver la estimación de MCO es que es el resultado de los momentos de una población:
| Momentos en la población | Momentos en la muestra |
|---|---|
| \[ E[u]= E[y-\beta_0-\beta_1x]=0 \] | \[ N^{-1} \sum^N_{i=1}(y_i-\hat{\beta_o}-\hat{\beta_1} x_i)=0 \] |
| \[ E[xu]= E[x(y-\beta_0-\beta_1x)=0] \] | \[ N^{-1} \sum^N_{i=1}x_i(y_i-\hat{\beta_o}-\hat{\beta_1} x_i)=0 \] |
La primera ecuación representa que desde un punto de vista poblacional, la ecuación representa que en promedio, el modelo de regresión lineal no comete errores sistemáticos en sus predicciones. Es decir, los errores son aleatorios y su valor medio es cero. Esta es una suposición clave en el modelo de regresión lineal y es necesaria para obtener estimadores de mínimos cuadrados imparciales y consistentes. Esto se obtiene por medio de derivar con respecto a \(\beta_0\) (veremos esto más adelante).
Sin embargo, aunque el valor esperado (representado por la letra E, que representa la esperanza matemática) de los errores sea cero, esto no significa necesariamente que cada error individual sea cero. En su lugar, significa que los errores son igualmente probables de estar por encima o por debajo de cero, lo que sugiere que no hay sesgo sistemático en las predicciones del modelo.
En la segunda ecuación nos indica una cuestión importante de la estadística: en su momento poblacional la Esperanza de \(xu\) es igual a 0. Esto es el equivalente a decir estadísticamente que \(E[x|u]= Cov(x,u)\). Por construcción la covarianza se puede expresar como:
\[Cov[x,u]= E[x|u]- E[x]E[u]\]
La segunda expresión por definición (en base a la primera ecuación) sabemos que es 0, porque \(E[u]=0\) y multiplicar un escalar (la esperanza matemática de \(x\)) da 0. Esto significa que la esperanza de \(x\) y \(u\) es 0 también. Esto significa que ambas variables no poseen relación o son ortogonales. De esto da cuenta Jeffrey Wooldridge:
Si \(u\) y \(x\) no están correlacionadas, entonces, como variables aleatorias, no están relacionadas linealmente. Suponer que \(u\) y \(x\) no están correlacionadas es un avance para definir el sentido en el que \(u\) y \(x\) estarán relacionanadas[…]. Sin embargo, el avance no es suficiente, ya que la correlación sólo mide dependencia lineal entre \(u\) y \(x\). La correlación tiene una propiedad un poco contraintuitiva: es posible que \(u\) no esté correlacionada con \(x\) y que, sin embargo, esté correlacionada con funciones de x como, por ejemplo, \(x^2\).[…] Esta posibilidad no es aceptable para la mayoría de los propósitos de la regresión, ya que causa problemas para interpretar el modelo y obtener propiedades estadísticas. Un supuesto mejor envuelve el valor esperado de \(u\) dado \(x\). Como \(u\) y \(x\) son variables aleatorias, se puede definir la distribución condicional de \(u\) dado cualquier valor de \(x\). En particular, para cada \(x\), se puede obtener el valor esperado (o promedio) de \(u\) en la porción de la población descrita por el valor de \(x\). El supuesto crucial es que el valor promedio de \(u\) no depende del valor de \(x\). (Wooldridge, 2010, p. 25)
Como vimos anteriormente en el punto 2.1 el cuadro nos refiere a dos ecuaciones que tienen dos incognitas, por lo que es un sistema de ecuaciones que tiene una solución identificable. \(\beta\) es un parámetro poblacional y \(\hat{\beta}\) es un estimador. La diferencia entre ambos es que el primero no lo conocemos (ya que desconocemos a priori la población y el segundo es una variable aleatoria (que depende de cada muestra) Considerando las dos condiciones de primer orden , las derivadas de \(\sum^N_{i=1}(y_i-\beta_o-\beta_1 x_i)=0\) con respecto de \(\beta_0\) y \(\beta_1\).
\[\sum^{N-1}_{i=1}(y_i-\hat{\beta_o}-\hat{\beta_1} x_i)=0\]
\[\sum^{N-1}_{i=1}x_i(y_i-\hat{\beta_o}-\hat{\beta_1} x_i)=0 \]
De la primera ecuación podemos estimar:
\[ \bar{y}= \hat{\beta_0}+ \hat{\beta_1}\bar{x} \] La notación sigma como tal es una forma de representar \(\bar{x}\) ya que la media es una sumatoria de un vector con i observaciones entre el número de observaciones i. Por lo que:
\[\bar{x}= N^-1 \sum^n_{i=1} x_i= N^-1(x_1+x_2+x_3+...x_n)\] Por lo tanto, podemos decir que \(\hat{\beta_0}\) es la media condicional de \(\bar{y}\) cuando \(\hat{\beta_1}\bar{x}= 0\). Entonces,
\[\hat{\beta_0}= \bar{y}-\hat{\beta_1}\bar{x}\]
En la segunda ecuación podemos transformarla en base a un paso algebráico1:\[\sum^{N}_{i=1}x_i[y_i-(\bar{y}-\hat{\beta_1}\bar{x})-\hat{\beta_1} x_i]=0 \]
\[\Rightarrow \sum^{N}_{i=1}x_i(y_i-\bar{y})=\beta_{0}\sum^N_{i=1} x_i(x_i-\bar{x})=0 \]
Finalmente podemos estimar a \(\hat{\beta}_1\) con la siguiente ecuación:
\[\hat{\beta_1}= \frac{\sum^N_{i=1}x_i(y_i-\bar{y})}{\sum^N_{i=1} x_i(x_i-\bar{x})}= \frac{\sum^N_{i=1}(x_i-\bar{x})(y_i - \bar{y})}{\sum^N_{i=1}(x_i-\bar{x})^2}\]
Esto significa que ambos parámetros/estimadores son:
\[ \hat{\beta_1}= \frac{ \sum^n_{i=1}(x_i-\bar{x})(y_i-\bar{y})}{ \sum^n_{i=1} (x_i- \bar{x})^2} \]
Y la constante es:
\[ \hat{\beta_0}= \bar{y}-\bar{x} \frac{\sum^n_{i=1}(x_i-\bar{x})(y_i-\bar{y})}{\sum^n_{i=1}(x_i-\bar{x})^2} \]
Vamos a simular dos muestras de datos \({y_i, x_i}^N_{i=1}\) para el cual \(y_i= \beta_0+\beta_1x_i+u_i\) con \(\beta_0= 2\) y \(\beta_1=4\). Una de las muestra de \(N=1000\) y la otra es de \(N=100\). Por último \(u_i\) y \(x_i\) son \(\sim
i.i.d\) \(\mathcal{N}(0,1)\).
Esto permite observar varias cuestiones que hemos mencionado sobre los
momentos en la muestra. \(\hat{\beta_0}\) y \(\hat{\beta_1}\) siguen una distribución.
Por lo que bajo diferentes recogidas de la misma muestra (que nos
aseguramos de la replicabilidad limitando el seed con el comando
set.seed(1234)), si se cumplen los supuestos que hemos
visto hasta ahora (y que luego abundaremos de manera más a fondo en la
sección posterior) los betas estimados se aproximan a los parámetros
poblaconales como muestra la tabla de regresión.
set.seed(1234)
N <- 1000
beta0 <- 2
beta1 <- 4
u <- rnorm(N, 0, 1)
x <- rnorm(N, 0, 1)
y <- beta0 + beta1 * x + u
modelo <- lm(y ~ x)
## Modelo 2:
set.seed(1234)
N2<-100
u2 <-rnorm(N2, 0, 1)
x2 <- rnorm(N2, 0, 1)
y2<- beta0+beta1*x2 +u2
modelo2<- lm(y2 ~x2)
Los resultados muestran que hay diferencia en en ambos modelos. Aun así, tanto el de 1000 observaciones y el de 100 contienen dentro de su intervalo de confianza (los números entre paréntesis en la tabla 2) al verdadero \(\beta\) que hicimos en la simulación de arriba. Sobre el por qué sucede esto lo detallaremos más en la sección 3 y 4.
#Crear tabla
library(stargazer)
library(htmltools)
stargazer(modelo, modelo2, title="Tabla 2. Resultados de la regresión", align=TRUE, type="html", covariate.labels = c("x","x2", "(Intercepto)"), ci=T, out = "Regresión_lineal.html")
| Dependent variable: | ||
| y | y2 | |
| (1) | (2) | |
| x | 4.058*** | |
| (3.995, 4.121) | ||
| x2 | 3.975*** | |
| (3.783, 4.168) | ||
| (Intercepto) | 1.973*** | 1.844*** |
| (1.911, 2.034) | (1.646, 2.042) | |
| Observations | 1,000 | 100 |
| R2 | 0.941 | 0.944 |
| Adjusted R2 | 0.941 | 0.943 |
| Residual Std. Error | 0.996 (df = 998) | 1.009 (df = 98) |
| F Statistic | 15,954.380*** (df = 1; 998) | 1,636.609*** (df = 1; 98) |
| Note: | p<0.1; p<0.05; p<0.01 | |
El teorema de Gauss-Markov establece las condiciones bajo las cuales el estimador de mínimos cuadrados ordinarios (OLS) es el mejor estimador lineal insesgado (BLUE - Best Linear Unbiased Estimator) de los parámetros de un modelo de regresión lineal. Los supuestos del teorema de Gauss-Markov son los siguientes:
Cuando todos estos supuestos se cumplen, el estimador de mínimos cuadrados ordinarios es el estimador lineal insesgado que tiene la mínima varianza posible, y por lo tanto es el BLUE. Reciéntemente formulaciones más modernas del Teorema de Gauss-Markov según Bruce Hansen ilustran que Mínimos Cuadrados, dentro de la familia de modelos lineales, no es solo es el mejor estimador lineal insesgado (MILE) si se cumplen los supuestos, sino que dentro de modelos no lineales este es el Best Unbiased Estimator (BUE) o en español el Mejor Estimador Insesgado (MEI) (Hansen, 2022).
Partiendo de lo dicho anteriormente si los supuestos 1-4, en el contexto del análisis de regresión para corte transversal, entonces MCO es insesgado, o en definitiva: \(E[\hat{\beta_0}|x]= \beta_0\) y \(E[\hat{\beta_1}|x]=\beta_1\).
Para probar su insesgadez3
\[ E[\hat{\beta_1}= E \left[ \frac{\sum^N_{i=1}(x_i-\bar{x})(y_i-\bar{y})}{\sum^N_{i=1}(x_i-\bar{x})^2} \right]= E \left[\frac{\sum^N_{i=1}(x_i-\bar{x})y_i}{\sum^N_{i=1}(x_i-\bar{x})^2} \right] \]
Por la propiedad de que \(\sum^N_{i=1}(x_i-\bar{x})(y_i-\bar{y})= \sum^N_{i=1}(x_i-\bar{x})y_i\) entonces,
\[…= E\left[\frac{\sum^N_{i=1}(x_i-\bar{x})(\color{red}{\beta_0+\beta_1x_i+u_i})}{\sum^N_{i=1}(x_i-\bar{x})^2}\right]\]
Partiendo del supuesto 1 (lineal en los parámetros) la relación entre \(y\) y \(x\) es lineal, o a través de una función lineal.
\[…=\frac{\sum^N_{i=1}(x_i-\bar{x})(\beta_0+\beta_1x_i+E[u_i])}{\sum^N_{i=1}(x_i-\bar{x})^2}\]
Por el supuesto 4 \(E[u|x]=0 = E[u_i]=0\)
\[…=\frac{\sum^N_{i=1}(x_i-\bar{x})(\beta_0+\beta_1x_i+\color{red}{0})}{\sum^N_{i=1}(x_i-\bar{x})^2}\]
Mediante transformar algebraicamente la ecuación anterior (propiedad distributiva):
\[…=\frac{\sum^N_{i=1}(x_i-\bar{x})\beta_0}{\sum^N_{i=1}(x_i-\bar{x})^2}+\frac{\sum^N_{i=1}(x_i-\bar{x})\beta_1}{\sum^N_{i=1}(x_i-\bar{x})^2}= 0+ \beta_1 \frac{\sum^N_{i=1}(x_i-\bar{x})^2}{\sum^N_{i=1}(x_i-\bar{x})^2}\]
Esto significa que \(E[{\hat{\beta_1}}]=\beta_1\) ya que el primer término que contiene a \(\beta_0\) es igual a 0. Ya que \(\beta_0\) es una constante y las sumas de las \(x\) lo son, por lo que la esperanza de una constante es igual 0, porque no hay variabilidad.
Partiendo de que la primera condición de momento de Mínimos cuadrados es \(\hat{\beta_0}= \bar{y}-\hat{\beta_1}\bar{x}\) aplicamos la esperanzas a ambos lados de la ecuación:
\[ E[\beta_0]= E[\bar{y}]-E[\hat{\beta_1}\bar{x}] \]
La esperanza de \(\bar{y}\) es igual a \(E[\beta_0+\beta_1\bar{x}+\bar{u}]\) esto significa que:
\[ \beta_0+\beta_1\bar{x}+E[\bar{u}]= \beta_0+\beta_1\bar{x} \]
Dado que \(E[\hat{\beta_1\bar{x}}]= E[\beta_1\bar{x}]\) significa que:
\[ E[\hat{\beta_0}]=\beta_0 \]
# Definir el tamaño de la muestra y el número de simulaciones
N <- 1000
M <- 5000
# Definir los verdaderos valores de los coeficientes
b0 <- 2
b1 <- 4
# Inicializar vectores para almacenar los coeficientes estimados
beta0 <- rep(NA, M)
beta1 <- rep(NA, M)
# Loop para realizar M simulaciones
for (j in 1:M) {
# Generar las variables
x <- rnorm(N, 0, 1)
u <- rnorm(N, 0, 1)
y <- b0 + b1*x + u
# Estimar la regresión y guardar los coeficientes
modelo <- lm(y ~ x)
beta0[j] <- coef(modelo)["(Intercept)"]
beta1[j] <- coef(modelo)["x"]
}
# Resumir los resultados
summary(beta0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.896 1.979 2.000 2.000 2.022 2.105
summary(beta1)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3.889 3.978 4.000 4.000 4.021 4.113
Viéndolo desde un punto de vista gráfico, por medio de dos histogramas vemos que aunque hay valores que no son exactamente los verdaderos \(\beta_0\) y \(\beta_1\) pero que la distribución de datos de la simulación se acerca a los parámetros poblacionales que son 2 y 4 respectivamente. Esto prueba lo que dijimos anteriormente, que la distribución de los \(\beta\) estimados contienen al parámetro poblacional.
#Las librerías necesarias:
library(ggplot2)
library(gridExtra)
#Convertir en data frames.
b0<-as.data.frame(beta0)
b1<-as.data.frame(beta1)
beta0<-ggplot(b0,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(b1, aes(y=beta1))+ geom_histogram(fill= "navy")+
geom_hline(yintercept = 4, color= "darkred", linetype= "dashed")+
labs(y=expression(beta[1]), x= "Recuento")+theme_light()
grid.arrange(beta0, beta1, nrow=1)
En un modelo de regresión lineal, el estimador de mínimos cuadrados ordinarios (MCO) es una manera común de estimar los parámetros del modelo. La varianza del estimador MCO para el parámetro \(β\) es dada por:
\[Var(\beta) = σ² * (X'X)^-1\]
Donde:
\(σ²\) es la varianza del término de error.
\(X\) es la matriz que contiene las variables independientes de su modelo de regresión.
\((X'X)^-1\) es la inversa de la matriz \(X\) transpuesta multiplicada por \(X\).
En el contexto de un modelo de una sola variable (regresión lineal simple), algebráicamente podemos expresarla de esta manera más simple:
\[Var[\hat{\beta_1}]=Var \left[\frac{\sigma^2}{\sum^N_{i=1}(x_i-\bar{x})^2} \right]\]
Esto asume que los supuestos del teorema de Gauss-Markov se cumplen, lo que significa que los errores tienen varianza constante (homoscedasticidad) y no están correlacionados entre sí \(Corr(e_{t}, e_{j})=0, ~~ (t \neq j)\) (no autocorrelación), entre otros supuestos.
Por último, es importante tener en cuenta que aunque la fórmula dada asume que conocemos la varianza del término de error \(σ²\), en la práctica este parámetro es generalmente desconocido y debe ser estimada a partir de los datos. La estimación más común para \(σ²\) es la suma de los cuadrados de los residuos dividida por los grados de libertad del modelo (\(n - k - 1\), donde \(n\) es el número de observaciones y \(k\) el número de predictores). Esto significa que conforme más observaciones el denominador \(\sum^N_{i=1}(x_i-\bar{x})^2\) se hace más grande, minimizando a \(\sigma^2\). Prueba de ello es ver la tabla 2 y comparar los modelos en donde la estimación de la varianza en el modelo con 1000 observaciones es más precisa (intervalo de confianza menos amplio).
La varianza del estimador de mínimos cuadrados ordinarios (MCO) puede ser derivada en el caso simple de una regresión lineal con un solo predictor. Para una regresión de la forma:
\[Y_i = β_0 + β_1x_i + u_i\]
Los estimadores de MCO para \(β_0\) y \(β_1\) son:
\[\hat{\beta_0} = \bar{y} - \hat{\beta_1}\bar{x}\] y \[\hat{\beta_1} = \frac{\sum^N_{i=1}(x_i-\bar{x})(y_i-\bar{y})}{\sum^N_{i=1}(x_i-\bar{x})^2}\]
La varianza del estimador \(\hat{\beta_1}\) se deriva así:
\[Var[\hat{\beta_1}]=Var \left[\frac{\sum^N_{i=1}(x_i-\bar{x})(y_i-\bar{y})}{\sum^N_{i=1}(x_i-\bar{x})^2} \right]\]
Supongamos que los errores \(u_i\) tienen varianza constante \(σ²\) (homocedasticidad), entonces podemos reemplazar \(y_i\) por su valor esperado y se obtiene:
\[Var[\hat{\beta_1}]=Var \left[\frac{\sum^N_{i=1}(x_i-\bar{x})^2\sigma^2}{\sum^N_{i=1}[(x_i-\bar{x})^2]^2} \right]\]
Luego simplificamos:
\[Var[\hat{\beta_1}]=Var \left[\frac{\sigma^2}{\sum^N_{i=1}(x_i-\bar{x})^2} \right]\]
Este resultado muestra que la varianza del estimador MCO de β1 disminuye a medida que aumenta la variabilidad de la variable predictora X.
Para el caso de múltiples predictores, la derivación de la varianza de los estimadores MCO se vuelve más complicada y generalmente requiere el uso de álgebra matricial. Sin embargo, el resultado final es similar: la varianza de los estimadores MCO depende inversamente de la variabilidad de los predictores y directamente de la varianza de los errores.
Dato que \(\hat{\beta_0}\) y \(\hat{\beta_1}\) son variables aleatorias, se requiere de su distribución para hacer inferencias, las cuales siguen, si partimos del supuesto de que los errores son independientes de \(x\) y que \(u\) \(\sim \mathcal{N}(0, \sigma^2)\), entonces:
\[ \beta_n \sim \mathcal{N}(\beta_n,Var[{\hat\beta_n}]) \]
Por lo que:
\[ \frac{(\hat{\beta_0}-\beta_0)}{se(\hat{\beta_0})} \]
Y para \(\beta_1\)
\[ \frac{(\hat{\beta_1}-\beta_1)}{se(\hat{\beta_1})} \]
Estas distribuciones de los betas se distribuyen normalmente \(\sim \mathcal{N}(0, 1)\), un ejemplo de ello es tomar el vector de las simulaciones de los coeficientes y aplicarles un test normalidad. En caso de que no se cumpla los supuestos sobre la varianza de los coeficientes es necesario corregir la varianza (en caso de que no se pueda mejorar el modelo incluyendo una variable o mejorando la especificación4), como por ejemplo usar errores estándares robustos (White, 1980) sin tampoco dejar de lado que aplicarlos no resuelve toda una serie de problemas inherentes con validar un modelo (King and Roberts, 2015).
Un aspecto importante a destacar es que \(\sigma^2\) es un valor que no conocemos. De forma que para fines de las ineferencias y contrastes de hipótesis nosotros la estimamos con los datos de la muestra. Cuando contrastamos si la pendiente del modelo es un resultado diferente de 0 \(H_a\) o si realmente fallamos de rechazar a \(H_0\) se requiere formalizar el test de la siguiente manera:
\(H_0\): \(\beta_1= \beta_{1,0}\) y que \(u \sim \mathcal{N}(0,1)\) entonces:
\[\frac{\hat{\beta_1}-\beta_{1,0}}{se(\hat{\beta_1})} \sim t_{n-2}\]
\(se(.)\) son los errores estándares y \(t_{n-2}\) es la distribución t de Student con \(n-2\) grados de libertad.
Al final la \(Var[\hat{\beta_1}]\) requiere que estimemos \(\sigma^2\) y se obtiene de esta manera:
\[\hat{\sigma}^2= \frac{\sum^N_{i=1}\hat{u}^2}{n-2}\]
shapiro.test(b0$beta0)
##
## Shapiro-Wilk normality test
##
## data: b0$beta0
## W = 0.99967, p-value = 0.6225
shapiro.test(b1$beta1)
##
## Shapiro-Wilk normality test
##
## data: b1$beta1
## W = 0.99941, p-value = 0.1068
Imaginemos que trabajamos en una empresa de Marketing y queremos ver
el impacto que tiene inaversión en publicidad en Youtube sobre el número
de ventas de un producto para hacer un reporte para indicar hacia dónde
debe de ir la estrategia de promoción de la misma. La base de datos (se
llama “marketing”) está en el paquete datarium el cual
podemos descargar usando el comando
install.packages("datarium").
library(datarium)
data("marketing")
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
Existe una relación, según el gráfico 2 se observa una relación entre ventas y youtube positiva.
library(ggplot2)
ggplot(marketing, aes(y=sales, x=youtube))+ geom_point(color= "red")+
geom_smooth(method= "lm", linetype= "dashed", color="green", fill= "darkgreen")+labs(y= "Ventas (en miles de USD)", x= "Inversión en publicidad en Youtube (miles de USD)")+ theme_light()
## `geom_smooth()` using formula = 'y ~ x'
modelosales<-lm(sales ~ youtube, marketing)
stargazer(modelosales, title="Tabla 3. Resultados de la regresión", align=TRUE, type="html", covariate.labels = c("Inversión Youtube", "(Intercepto)"), dep.var.labels = c("Ventas"), out = "Regresión_lineal2.html")
| Dependent variable: | |
| Ventas | |
| Inversión Youtube | 0.048*** |
| (0.003) | |
| (Intercepto) | 8.439*** |
| (0.549) | |
| Observations | 200 |
| R2 | 0.612 |
| Adjusted R2 | 0.610 |
| Residual Std. Error | 3.910 (df = 198) |
| F Statistic | 312.145*** (df = 1; 198) |
| Note: | p<0.1; p<0.05; p<0.01 |
Como vimos anteriormente en los puntos 3 y 4 se observó los supuestos del modelo. Esto significa se requiere de hacer una evaluación exhaustiva para validarlo o para corregirle en caso de que alguno de los supuestos falle.
Podemos hacer diagnósticos usando la función plot()
Al observar el modelo, existen ciertos datos atípicos al mirar los residuos, se presenta cierto sesgo y heterocedasticidad (varianza no constante). Es quizás necesario ajustar un término cuadrático para mejorar la especificación del modelo o incluir otra variable (sesgo de variables omitidas) o aplicar alguna transformación logarítmica. También el hecho de que los datos son en el tiempo es posible que exista correlación e incluir rezagos o adelantos permitan mejorar el modelo.
plot(modelosales)
Los estadísticos de diagnósticos podemos estimarlos con el comando
summary pero incluyendo la función influence.measures().
Cada medida de influencia tiene puntos de cortes para así estimar qué
sucede si se eliminan ciertos valores influyentes. Su significado es el
siguiente:
dffit: Es el cambio en los valores ajustados. El punto de corte es \(2\sqrt{(p+1)/n}\). Siendo \(p\) el número de predictores.
dfb: Es el cambio en los coeficientes, el punto de corte es \(2\sqrt{n}\)
cov.r: Es el cambio en la matriz de varianzas y covarianzas (\(Cov(\hat{\beta})\)) el cual el punto de corte es \(1 \pm 3(p+1)/n\)
cook.d: Es la distancia de Cook, el cual su punto de corte es \(4/n\)
hat: Es el estadístico de apalancamiento (Leverage), su punto de corte es \(2(p+1)/n\)
Si calculamos el punto de corte para el modelo es de:
2*sqrt((1+1)/200)
## [1] 0.2
Y el cambio en la covarianza es de
1+3*((1+1)/200)
## [1] 1.03
1-3*((1+1)/200)
## [1] 0.97
summary(influence.measures(modelosales))
## Potentially influential observations of
## lm(formula = sales ~ youtube, data = marketing) :
##
## dfb.1_ dfb.yotb dffit cov.r cook.d hat
## 26 0.11 -0.23 -0.28 0.97_* 0.04 0.01
## 36 0.17 -0.30 -0.35_* 0.97_* 0.06 0.02
## 43 0.01 -0.01 -0.01 1.03_* 0.00 0.02
## 56 0.00 0.10 0.19 0.97_* 0.02 0.01
## 129 -0.04 0.14 0.21 0.97_* 0.02 0.01
## 179 0.15 -0.28 -0.34_* 0.96_* 0.06 0.02
Como vimos en la sección para el análisis de correlación, la potencia estadística toma un valor entre 0 y 1. Un valor normalmente usado (aunque dependerá del diseño del experimento y del efecto a buscar en el diseño) es de .80. Esto significa que si repetimos el experimento en varias ocaciones el experimento podrá detectar el efecto y rechazar la hipótesis nula corrrectamente en un 80 % de las veces (Cohen, 2016). Es posible estimar la potencia de un experimento o de un estudio en base a una serie de parámetros tanto teóricos como empíricos si este usa una regresión.
Antes de proceder al análisis vamos a recordar uno de los estadísticos de bondad de ajuste del modelo de regresión, el \(R^2\):
\[ R^2 = 1 - \frac{SS_{res}}{SS_{tot}} \]
El cual es 1 menos el cociente de la suma de cuadrados del error y la suma de cuadrados total. Una forma de representarlo sería:
\[ R^2 = 1 - \frac{\sum_{i}(y_i - \hat{y_i})^2}{\sum_{i}(y_i - \bar{y})^2} \]
Para poder observar si el modelo es válido se requiere de probar la significación conjunta del mismo. Es decir, que todas las variables (y esto aplica más para modelos de regresión múltiple). Una forma de estimar el efecto es por medio de la prueba F, la cual podríamos definir como:
\[ F2= \frac{R^2}{(1-R^2)} \]
En R con el la librería pwr podemos estimar la potencia
de nuestro modelo:
Al ver los resultados indica que para la muestra que tenemos y la \(k\) cantidad de predictores (1) del modelo y el efecto que queremos medir, los datos que poseemos poseen una potencia de 100 %, por lo que la muestra de 200 observaciones es más que suficiente para observar el efecto que detectamos.
#Paso 1 estimemos F2:
f2<- 0.612/(1-0.612) # 0.612 es el R2 del modelo.
#Cargamos la librería:
library(pwr)
pwr.f2.test(u=1, #Número de predictores
v=198, #Grados de libertad del modelo (200-1-1)
f2=f2, #Tamaño del efecto
sig.level = 0.05, #Nivel de significancia
power= NULL)
##
## Multiple regression power calculation
##
## u = 1
## v = 198
## f2 = 1.57732
## sig.level = 0.05
## power = 1
la prediccion del modelo se hace con la siguiente ecuación:
\[ \hat{y_i}= \hat{\beta_0}+\hat{\beta_i}x_i \]
Doon las estimaciones de los parámetros del modelo.
\(\bar{x}\) y \(\bar{y}\) son las medias de X e Y, respectivamente.
\(n\) es el número de observaciones, representado como \(i=1, 2,…, n\)
Una vez que tienes las estimaciones de los parámetros, puedes hacer predicciones utilizando la fórmula del modelo, reemplazando \(β_0\) y \(β_1\).
La función universal predict() del paquete
stats en R facilita la realización de predicciones con
nuestros modelos, utilizando nuevos valores de nuestra variable
predictora. Para señalar a la función los valores de la variable
predictora en los que deseamos estimar la variable de respuesta, debemos
conformar un dataframe que incorpore los nuevos valores que pretendemos
analizar (new).
predict(model, newdata = new, interval = "prediction")
La salida incluye las siguientes columnas:
fit: los valores de venta estimados para los tres nuevos presupuestos publicitarios.
lwr y upr: los límites de confianza inferior y superior para los valores proyectados, respectivamente.
Puedes visualizar los resultados mediante la herramienta
ggplot2.
Además de las predicciones concretas, podríamos estar interesados en calcular sus intervalos. Tenemos dos alternativas:
El intervalo de confianza, que exhibe la incertidumbre asociada a las estimaciones medias.
El intervalo de predicción, que expresa la incertidumbre en torno a un valor individual de predicción.
Un intervalo de predicción manifiesta la incertidumbre en torno a una sola estimación, mientras que un intervalo de confianza muestra la incertidumbre en torno a los valores de predicción media. Por ende, un intervalo de predicción generalmente será considerablemente más amplio que un intervalo de confianza para el mismo valor.
¿Cuál debemos seleccionar? Dependerá del contexto y del objetivo del análisis. Usualmente, estamos interesados en predicciones específicas individuales, por lo que un intervalo de predicción sería más adecuado. Ten presente que utilizar un intervalo de confianza cuando deberías usar un intervalo de predicción minimizará considerablemente la incertidumbre en un valor de predicción específico (Bruce et al., 2020)
new<-data.frame(youtube= c(100, 120, 300, 350, 400))
predict(modelosales, newdata = new,
interval ="confidence")
## fit lwr upr
## 1 13.19278 12.51317 13.87239
## 2 14.14351 13.52138 14.76563
## 3 22.70010 21.84743 23.55278
## 4 25.07694 24.00676 26.14711
## 5 27.45377 26.14830 28.75923
Es posible hacer este pasaje en base a esto:
\[ \sum_{i=1}^{N} a_i(b_i - \bar{b}) = \sum_{i=1}^{N} b_i(a_i - \bar{a})= \sum_{i=1}^{N} (a_i - \bar{a})(b_i - \bar{b})\]
La suma de cada valor de a menos su promedio es igual a 0 \(\sum a_i= n*\bar{a}\), por definición la \(n * \bar{a}= n \frac{\sum a_i} {n}\). Esto significa que los n se cancelan y solo queda que \(\sum a_i= \bar{a}\). Esto significa que si tengo una expresión que contiene \(a_i\) y \(b_i\) y a sus promedios, si aplico la propiedad distributiva implica que la sumatoria de \(a_i\) por la sumatoria de \((b_i-\bar{b})\) es igual a la ecuación de arriba.↩︎
Esto aplica para el contexto del análisis de series temporales, no para cuando el análisis de regresión es con datos de corte transversal (Wooldridge, 2010).↩︎
Para fines de simplificación vamos a escribir \(E[.]\) como sinónimo de una esperanza condicionada \(E[.|x]\)↩︎
Sobre forma funcional y modelos de regresión múltiple lo retomaremos en siguentes lecturas).↩︎