Sesión 12. Laboratorio Métodos Cuantitativos Aplicados.

Sobre los coeficientes de la regresión lineal univariada y el \(R^2\).

En los siguientes ejercicios trabajaras con la base de datos advertising, publicada por los autores de An Introduction to Statistical Learning.

library(pacman)
p_load(broom, janitor, moderndive, Rlab, tidyverse)

# Base de datos de publicidad. 
publicidad <- 
  read_csv("http://segasi.com.mx/clases/cide/datos/Advertising.csv") %>%
# Eliminamos el indice
  select(-X1) %>% 
# Limpiamos nombres
  clean_names()

# Exploramos la Base de datos
publicidad

## # A tibble: 200 x 4
##       tv radio newspaper sales
##    <dbl> <dbl>     <dbl> <dbl>
##  1 230.   37.8      69.2  22.1
##  2  44.5  39.3      45.1  10.4
##  3  17.2  45.9      69.3   9.3
##  4 152.   41.3      58.5  18.5
##  5 181.   10.8      58.4  12.9
##  6   8.7  48.9      75     7.2
##  7  57.5  32.8      23.5  11.8
##  8 120.   19.6      11.6  13.2
##  9   8.6   2.1       1     4.8
## 10 200.    2.6      21.2  10.6
## # … with 190 more rows

La base de datos incluye 200 renglones (uno por mercado) y cuatro variables: tv, radio, newspaper y sales. Las primeras tres variables describen el presupuesto gastado en publicidad (en miles de dólares) por una compañía en el medio respectivo en cada mercado. La columna sales indica el número de unidades vendidas (en miles) por la compañía en el mercado correspondiente.

A partir de estos datos, resuelva los siguientes ejercicios:

Ejercicio 1

Estima manualmente, como hicimos en clase, \(\beta_0\) y \(\beta_1\) para el siguiente modelo:

\[sales = \beta_0 + \beta_1 * radio + u\]

¿Qué estamos estimando?

R: Con el intercepto, las ventas base que tendría una empresa que no se anuncia en ninguno de los medios de comunicación.

Con el \(\beta_1\), el efecto sobre las ventas que tiene el invertir un dólar en publicidad en la radio.

Recordemos, del laboratorio pasado, la fórmula para sacar \(\beta_0\) y \(\beta_1\):

Para \(\beta_1\):

\[\beta_1 = \dfrac{cov(y,x)}{var(x)}\]

Donde:

\[cov(x,y) = \dfrac{\Sigma(x_i-\hat{x})(y_i-\hat{y})}{n}\]

Y el \(\beta_0\):

\[\beta_0 = mean(y) - \beta_1*mean(x)\]

# Procesamos la base de datos de publicidad. 
publicidad %>% 
# Calculamos la desviacion de la var. radio  
  mutate(dif_radio = radio - mean(radio),
# Calculamos la desviacion de la var. sales         
         dif_sales = sales - mean(sales),
# Calculamos el producto de las diferencias
         producto_dif = dif_radio*dif_sales,
# Calculamos el producto de las diferencias al cuadrado
         dif_radio_cuadrado = dif_radio^2) %>% 
# Obtenemos beta_1 y beta_0 con las formulas.   
  summarise(beta_1 = sum(producto_dif)/sum(dif_radio_cuadrado),
            beta_0 = mean(publicidad$sales) - beta_1*mean(publicidad$radio))

## # A tibble: 1 x 2
##   beta_1 beta_0
##    <dbl>  <dbl>
## 1  0.202   9.31

Ejercicio 2

Para verificar tus resultados, calcula un modelo lineal con la función lm() y guarda el resultado en un objeto llamado modelo. Después usa las funciones summary() y/o tidy() para obtener los coeficientes.

# Generamos el modelo
modelo <- lm(formula = sales ~ radio, 
             data = publicidad)

# Obtenemos mas información del modelo
summary(modelo)

## 
## Call:
## lm(formula = sales ~ radio, data = publicidad)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.7305  -2.1324   0.7707   2.7775   8.1810 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  9.31164    0.56290  16.542   <2e-16 ***
## radio        0.20250    0.02041   9.921   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.275 on 198 degrees of freedom
## Multiple R-squared:  0.332,  Adjusted R-squared:  0.3287 
## F-statistic: 98.42 on 1 and 198 DF,  p-value: < 2.2e-16

# Lo hacemos tidy para sacarle información mas fácilmente
tidy(modelo)

## # A tibble: 2 x 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)    9.31     0.563      16.5  3.56e-39
## 2 radio          0.202    0.0204      9.92 4.35e-19

Ejercicio 3

Usa la función augment() para extraer otro chunk de datos del objeto modelo.

A partir del tibble que obtendrás, haz una gráfica con dos capas. La primera es de puntos, para representar la relación entre sales y radio (acuérdate que las variables deben respetar cierto orden en los ejes).

La segunda capa es una gráfica de líneas, para representar los valores predichos por el modelo para sales y radio.

# Obtenemos los datos aumentados con augment: 
augment(modelo)

## # A tibble: 200 x 9
##    sales radio .fitted .se.fit   .resid    .hat .sigma   .cooksd .std.resid
##    <dbl> <dbl>   <dbl>   <dbl>    <dbl>   <dbl>  <dbl>     <dbl>      <dbl>
##  1  22.1  37.8   17.0    0.424   5.13   0.00982   4.27   7.22e-3     1.21  
##  2  10.4  39.3   17.3    0.446  -6.87   0.0109    4.26   1.43e-2    -1.62  
##  3   9.3  45.9   18.6    0.552  -9.31   0.0167    4.23   4.09e-2    -2.20  
##  4  18.5  41.3   17.7    0.476   0.825  0.0124    4.29   2.37e-4     0.194 
##  5  12.9  10.8   11.5    0.395   1.40   0.00854   4.28   4.67e-4     0.329 
##  6   7.2  48.9   19.2    0.604 -12.0    0.0200    4.20   8.22e-2    -2.84  
##  7  11.8  32.8   16.0    0.360  -4.15   0.00707   4.28   3.39e-3    -0.975 
##  8  13.2  19.6   13.3    0.311  -0.0806 0.00531   4.29   9.52e-7    -0.0189
##  9   4.8   2.1    9.74   0.527  -4.94   0.0152    4.27   1.05e-2    -1.16  
## 10  10.6   2.6    9.84   0.519   0.762  0.0147    4.29   2.41e-4     0.180 
## # … with 190 more rows

# .fitted es la columna que almacena los valores predichos
## por el modelo de regresión lineal para el valor 
## de X correspondiente

# Graficamos el modelo aumentado: 
augment(modelo) %>% 
  ggplot() +
# 1a capa: linea de los valores realoes
  geom_point(aes(x = radio, y = sales)) +
# 2a capa: linea de los valores predichos
  geom_line(aes(x = radio, y = .fitted), 
            color = "dodgerblue")

Ejercicio 4

Usa el método de bootstrap para generar la distribución muestral de \(\hat \beta_0\) y \(\hat \beta_1\). Después gráfica cada distribución (ojo: vas a tener que hacerlo en por lo menos dos chunks de código).

Primero para el intercepto:

# (Intercept) - B0
publicidad %>% 
  # Tomar 1k muestras bootstrap
  rep_sample_n(size = nrow(publicidad), 
               replace = T, 
               reps = 1000) %>% 
  # Calcular el modelo de regresión para cada muestra y tidyear
  do(lm(formula = sales ~ radio, 
        data = .) %>% 
       tidy()) %>% 
  # Desagrupar
  ungroup() %>% 
  # Filtrar para quedarnos con estimaciones del coeficiente de altura (B_1)
  filter(term == "(Intercept)") %>%   
  # Calcular el error estándar de la distribución muestral
  ggplot(aes(estimate)) +
  geom_density() + 
  geom_vline(xintercept = mean(modelo$coefficients[1]), 
             color = "red", 
             linetype = 2) + 
  labs(title = "Estimación del B_0 (Intercepto) mediante Bootstrap")

Después, para el \(\beta_1\)

# radio - B1
publicidad %>% 
  # Tomar 1k muestras bootstrap
  rep_sample_n(size = nrow(publicidad), replace = T, reps = 1000) %>% 
  # Calcular el modelo de regresión para cada muestra y tidyear
  do(lm(formula = sales ~ radio, data = .) %>% 
       tidy()) %>% 
  # Desagrupar
  ungroup() %>% 
  # Filtrar para quedarnos con estimaciones del coeficiente de altura (B 1)
  filter(term == "radio") %>%   
  # Calcular el error estándar de la distribución muestral
  ggplot(aes(estimate)) +
  geom_density() + 
  geom_vline(xintercept = modelo$coefficients[2],
             color = "red", 
             linetype = 2) +
  labs(title = "Estimación del B_1 (radio) mediante Bootstrap")

Ejercicio 5

Usa el método de bootstrap para calcular el error estándar de \(\hat \beta_0\) y \(\hat \beta_1\). Después usa tidy(modelo) para comparar tus resultados con los que calculó lm().

# (Intercept) - B0
publicidad %>% 
  # Tomar 1k muestras bootstrap
  rep_sample_n(size = nrow(publicidad), replace = T, reps = 5000) %>% 
  # Calcular el modelo de regresión para cada muestra y tidyear
  do(lm(formula = sales ~ radio, data = .) %>% 
       tidy()) %>% 
  # Desagrupar
  ungroup() %>% 
  # Filtrar para quedarnos con estimaciones del coeficiente de altura (B 1)
  filter(term == "(Intercept)") %>%   
  # Calcular el error estándar de la distribución muestral
  summarise(error_est = sd(estimate))

## # A tibble: 1 x 1
##   error_est
##       <dbl>
## 1     0.391

# radio - B1
publicidad %>% 
  # Tomar 1k muestras bootstrap
  rep_sample_n(size = nrow(publicidad), replace = T, reps = 5000) %>% 
  # Calcular el modelo de regresión para cada muestra y tidyear
  do(lm(formula = sales ~ radio, data = .) %>% 
       tidy()) %>% 
  # Desagrupar
  ungroup() %>% 
  # Filtrar para quedarnos con estimaciones del coeficiente de altura (B 1)
  filter(term == "radio") %>%   
  # Calcular el error estándar de la distribución muestral
  summarise(error_est = sd(estimate))

## # A tibble: 1 x 1
##   error_est
##       <dbl>
## 1    0.0218

tidy(modelo)

## # A tibble: 2 x 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)    9.31     0.563      16.5  3.56e-39
## 2 radio          0.202    0.0204      9.92 4.35e-19

Los errores son algo distintos. Sin embargo, los datos del error estándar son mejores estimaciones del error estándar que los obtenidos por el modelo, que aplica una aproximación teórica, más que una simulación computacional.

Ejercicio 6

Usa summary(modelo) o glance(modelo) para obtener -entre otras cosas- la \(R^2\) de modelo.

# Generamos nuevamente el modelo
modelo <- lm(formula = sales ~ radio, data = publicidad)
# Glanceamos el modelo
glance(modelo)

## # A tibble: 1 x 11
##   r.squared adj.r.squared sigma statistic  p.value    df logLik   AIC   BIC
##       <dbl>         <dbl> <dbl>     <dbl>    <dbl> <int>  <dbl> <dbl> <dbl>
## 1     0.332         0.329  4.27      98.4 4.35e-19     2  -573. 1153. 1163.
## # … with 2 more variables: deviance <dbl>, df.residual <int>

Ahora échale un ojo a esta ecuación:

\[\small r^2 = \frac{\sum (y_i - \bar y)^2 - \sum (y_i - \hat y)^2}{\sum (y_i - \bar y)^2}\]

¿Cómo llamamos en clase a los tres componentes que están del lado derecho del símbolo de igual?
- El primer elemento del numerador es la varianza total
- En el segundo elemento del numerador está varianza no-explicada
¿Cómo llamamos en clase a al numerador y denominador de esta ecuación?
- El denominador es la varianza total
Usa esta fórmula para calcular a patita la \(R^2\).

# Calculamos el modelo
modelo <- lm(formula = sales ~ radio, data = publicidad)

# Aumentamos el modelo
augment(modelo) %>% 
# Datos menos la media (al cuadrado)
  mutate(e_1 = (sales - mean(sales))^2,
# Datos menos la estimación, al cuadrado         
         e_2 = (sales - .fitted)^2) %>% 
# Suma del primer término  
  summarise(suma_e_1 = sum(e_1),
# Suma del segundo término            
            suma_e_2 = sum(e_2)) %>% 
# Aplicamos la fórmula, una vez que tenemos todos los 
# ingredientes. 
  mutate(r_2 = (suma_e_1 - suma_e_2)/suma_e_1)

## # A tibble: 1 x 3
##   suma_e_1 suma_e_2   r_2
##      <dbl>    <dbl> <dbl>
## 1    5417.    3618. 0.332