Regression Notes

Vásquez Guerra Carlos Fernando

08/2021

Regresión lineal simple

NOTA: Este trabajo es una continuación directa del trabajo Statistical Notes que se puede encontrar en el siguiente enlace

Nuestro objetivo en un modelo de regresión será ajustar o asociar una variable, a la que llamaremos variable respuesta o dependediente, con una o más variables distintas llamadas variables predictoras o independientes, con el fin de poder predecir valores futuros de la variable respuesta y analizar los efectos que tendría esta variable con alteraciones en el resto de variables del modelo.

En este caso, el de la regresión lineal simple consideramos dos variables (1 predictora (\(x\)) y otra como respuesta (\(y\))) y establecemos este modelo de la siguiente manera:

\[ y_{i}=\beta_{0}+\beta_{1}x_{i}+\epsilon_{i} \]

Por algunas propiedades interesantes, podemos ver que \(\mathbf{E}[y_{i}]=\beta_{0}+\beta_{1}x_{i}\) por lo que este modelo considera que la media de \(y\) cambia a razón constante mediante los cambios de \(x\). Como estamos obteniendo una recta, \(\beta_{1}\) es la pendiente de dicha recta y es lo que nos indica la importancia de \(x\) sobre \(y\) (que tan drásticos son los cambios de una sobre otra).

\(\beta_0\) nos indica el valor de \(y\) cuando la variable predictora esta ausente. \(\epsilon\) es simplemente un error aleatorio que nos dará ciertos problemas para la validación de este modelo mediante pruebas de hipótesis pero nos otorga la flexibilidad probabilística para la predicción.

Es sencillo realizar una regresión lineal en R, simplemente podemos utilizar el comando base::lm() agregando una fórmula (~) como parámetro. Por ejemplo, consideremos la base de datos datasets::cars donde se tiene la velocidad y la distancia que recorrió un coche para detenerse; estos datos fueron recolectados en 1920.

cars %>% ggplot(aes(x = speed, y = dist)) + labs(x = "Velocidad", y = "Distancia") + geom_point(colour = "purple3") + general_theme

Ahora, realizaremos nuestro modelo lineal simple:

car_regression <- lm(dist~speed, data = cars)
summary(car_regression)

Call:
lm(formula = dist ~ speed, data = cars)

Residuals:
    Min      1Q  Median      3Q     Max 
-29.069  -9.525  -2.272   9.215  43.201 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -17.5791     6.7584  -2.601   0.0123 *  
speed         3.9324     0.4155   9.464 1.49e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 15.38 on 48 degrees of freedom
Multiple R-squared:  0.6511,    Adjusted R-squared:  0.6438 
F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12

Véase que se nos da mucha información:

  • Estadísticas sobre los residuales.
  • Información sobre los coeficientes de la regresión mediante mínimos cuadrados.
  • Información sobre el ajuste del modelo con los datos.

Estos puntos se verán después, algunos de ellos se verán hasta regresión lineal múltiple ya que tienen una interpretación similar, pero por el momento basta con saber que, de acuerdo a los datos obtenidos, nuestro modelo queda expresado de la siguiente manera:

\[ y_i = -17.5791 + 3.9324x_i \]

Que es lo mismo al siguiente modelo:

\[ \mbox{Distancia del auto }i = -17.5791 + 3.9324\times\mbox{Velocidad del auto }i \]

Y ahora podemos generar la siguiente gráfica

cars %>% ggplot(aes(x = speed, y = dist)) + 
  labs(x = "Velocidad", y = "Distancia") + general_theme +
  geom_point(colour = "purple3") + 
  geom_abline(intercept = car_regression$coefficients[1], 
              slope = car_regression$coefficients[2],
              colour = "red")

Hay cosas importantes que ver en este tema y aquí se comenzará con un pequeño análisis de correlación.

Correlación y linearidad

La manera más sencilla de proponer un conjunto de variables predictoras es ver la relación lineal que tienen estas con la variable respuesta. Esto lo podemos obtener mediante diferentes estadísticas, por ejemplo:

  • Coeficiente de correlación de pearson: \(\rho_{xy} = r_{xy} = \frac{\sum_{i = 1}^n(x_i-\bar{x})(y_i-\bar{y})}{\sqrt{\sum_{i = 1^n}(x_i-\bar{x})^2}\sqrt{\sum_{i = 1}^n(y_i-\bar{y})^2}}\)
  • Coeficiente de correlación de spearman: \(r_s = 1- \frac{6\sum d_i^2}{n(n^2-1)}\) donde \(d_i = rango(X_i)-rango(Y_i)\)
  • Coeficiente de correlación de Kendall: \(\tau = \frac{(\mbox{número de pares concordantes})-(\mbox{número de pares discordantes})}{n \choose 2}\)

Cada uno de las estadísticas anteriores tiene un propósito diferente, \(\rho\) mide la relación lineal entre las variables que se están comparando, el coeficiente de correlación de spearman mide relaciones monótonas (una relación lineal lo es) y el coeficiente de Kendall mide la semejanza en el ordenamiento de los datos cuando se clasifican en rangos por cada una de las cantidades. Todas estas se pueden obtener mediante la función base::cor(x, y, method = c("pearson", "kendall", "spearman")).

Es importante no sólo considerar el primero de los coeficientes de correlación mencionados, ya que puede existir una clara relación entre las variables (como en un conjunto de puntos que dibujen una parábola) y el coeficiente tener un valor de 0. Vamos a tomar los datos dplyr::starwars donde se tienen datos descriptivos de los personajes de la saga Star Wars.

starwars

Sólo para estos fines, eliminaremos un outlier sobre el peso y consideraremos sólo esta y la altura para ver si existe algún tipo de relación entre estas variables.

new_starwars <- starwars %>% filter(mass < max(starwars$mass, na.rm = T))

pearson_star <- cor(new_starwars$height, new_starwars$mass)
spearman_star <- cor(new_starwars$height, new_starwars$mass, method = "spearman")
kendall_star <- cor(new_starwars$height, new_starwars$mass, method = "kendall")

new_starwars %>% ggplot(aes(x = height, y = mass)) + 
  geom_point(color = "royalblue4") + labs(x = "Altura (cm)", y = "Peso (kg)") +
  annotate("text", x = 100, y = 150, label = paste("Pearson: ", pearson_star, collapse = "")) +
  annotate("text", x = 100, y = 140, label = paste("Spearman: ", spearman_star, collapse = "")) +
  annotate("text", x = 100, y = 130, label = paste("Kendall: ", kendall_star, collapse = "")) +
  general_theme

Como vemos, el coeficiente de correlación de pearson tiene un valor cercano al de Spearman, aunque este último nos inidica que existe una relación monótona, por lo que podríamos aplicar una transformación para otorgar un comportamiento más cercano al lineal con una función monóntona como lo es el logaritmo.

pearson_star <- cor(new_starwars$height, log(new_starwars$mass))
spearman_star <- cor(new_starwars$height, log(new_starwars$mass), method = "spearman")
kendall_star <- cor(new_starwars$height, log(new_starwars$mass), method = "kendall")

new_starwars %>% ggplot(aes(x = height, y = log(mass))) + 
  geom_point(color = "royalblue4") + labs(x = "Altura (cm)", y = "Peso (kg)") +
  annotate("text", x = 100, y = 4.8, label = paste("Pearson: ", pearson_star, collapse = "")) +
  annotate("text", x = 100, y = 4.6, label = paste("Spearman: ", spearman_star, collapse = "")) +
  annotate("text", x = 100, y = 4.4, label = paste("Kendall: ", kendall_star, collapse = "")) +
  general_theme

Nuestro coeficiente de correlación lineal aumente significativamente y surge la siguiente pregunta: ¿Cuándo utilizar alguna transformación y cuál? En general, cualquier función que sólo cambie la escala y este definida para los datos que se estén considerando se puede utilizar, como un logaritmo (que nos ayuda a reducir dispersión y dar valores pequeños), una raíz cuadrada, una función inversa, etc; aunque hay que tener cuidado. Por ejemplo ¿Qué pasaría con si tengo valores iguales a cero y aplico un logaritmo? ¿La interpretación sobre los coeficientes cambiaría o no? Este último punto se verá después.

  • Se recomienda revisar el libro The Statistical Sleuth, tercera edición donde se muestran algunos patrones comunes que podrían impedir realizar alguna relación lineal entre variables.

Ahora, si bien es razonable pensar que la altura y el peso estén relacionados de manera lineal ¿Podemos confiar siempre en el resultado de una correlación? No, ya que una correlación no implica causalidad (una relación verdadera entre dos variables). Veamos los siguientes ejemplos del siguiente enlace:

drawing
drawing

Es difícil pensar que realmente existe una relación entre dichas variables en cada caso, así como considerar que existe una relación entre la cantidad de IPhones comprados y la cantidad de personas que mueren al caer de las escaleras o pensar que la cantidad de cigueñas esta relacionado con la cantidad de nacimientos en Europa aunque se tenga un \(r = 0.62\) y un \(p-value = 0.008\) en la prueba de hipótesis correspondiente a \(\beta_i\).

Otro punto interesante es, ¿Qué tanto podemos confiar sólo en nuestras estadísticas? Consideremos la siguiente información obtenida del paquete datasauRus.

library(datasauRus)
dino_data <- datasauRus::datasaurus_dozen %>% filter(dataset == "dino")
  datasaurus_dozen %>% 
    group_by(dataset) %>% 
    summarize(
      mean_x    = mean(x),
      mean_y    = mean(y),
      std_dev_x = sd(x),
      std_dev_y = sd(y),
      corr_x_y  = cor(x, y)
    )

Como vemos, las estadísticas son bastante similares entre todos los conjuntos de datos, incluso la correlación; por lo que tal vez pensaríamos en tener datos similares; bueno veámos como se ven algunos de estos datos.

ggplot(datasaurus_dozen %>% filter(dataset %in% c("dino", "star", "h_lines")), aes(x=x, y=y, colour=dataset))+
  geom_point()+
  theme_void()+
  theme(legend.position = "none")+
  general_theme+
  facet_wrap(~dataset, ncol=3)

Aquí se dejan algunos enlaces sobre estos ejemplos en particular:

Y algunos con temas relacionados a la correlación:

Variables categóricas

Se puede trabajar con variables categóricas en un modelo de regresión para determinar la importante que tiene una cierta segregación en nuestra variable a modelar. Si bien se podría remplazar todos los valores categóricos por valores numéricos, esto nos ayudaría sólo con valores ordinales. Otra técnica común es la creación de variables dummys a lo cual se le conoce como One-Hot-Encoding.

Para esta sección se utilizarán algunas variables de los tabulados de la encuesta Encuesta Nacional de Ingresos y Gastos de los Hogares (ENIGH). 2018 Nueva serie, específicamente para el estado de Aguascalientes. Nos interesan ciertas variables que se puedan relacionar con el ingreso corriente de los hogares de dicha región y en esta sección nos interesa determinar la importancia del sexo de los jefes del hogar.

Vamos a aplicar una regresión lineal entre el ingreso corriente y el sexo del jefe de familia para ver como se comporta la función lm().

data_income_sexo  %>% lm(ing_cor~sexo_jefe, .) %>% tidy()

Al parecer trabaja bien dicha función con datos categóricos. Lo que sucedio es que tal función realizo una codificación sobre el factor sexo_jefe, de hecho aplico esta configuración

contrasts(data_income_sexo$sexo_jefe)
       Mujer
Hombre     0
Mujer      1

Es decir, la variable dummy que se esta considerando en el modelo tiene la siguiente codificación:

  • 1 si una persona es mujer
  • 0 si una persona es hombre

Lo que otorgaría una ponderación sobre las variables del modelo cuando estos son diferentes de cero:

  • \(\beta_0 + \beta_1\) si una persona es mujer.
  • \(\beta_0\) si una persona es hombre.

Estos pueden ser interpretados de la siguiente manera:

\(\beta_0 = 60071.50\) es el promedio del ingreso entre los hombres. \(\beta_0 + \beta_1 = 60071.50 - 9822.53 = 50248.97\) es el promedio del ingreso entre las mujeres \(beta_1 = -9822.53\) es el la diferencia promedio en el ingreso entre los hombres y mujeres

El coeficiente negativo para las mujeres en la regresión sólo indica que la categoría “mujer” está asociado con una disminución en el salario (en relación con los hombres).

Si cambiamos la configuración de la variable dummy, la interpretación y los valores “se conservan”

di_nc <- data_income_sexo %>%
  mutate(sexo_jefe = relevel(sexo_jefe, ref = "Mujer"))
contrasts(di_nc$sexo_jefe)
       Hombre
Mujer       0
Hombre      1

Véase que el ingreso promedio entre las mujeres sigue siendo 50248.97 y la diferencia promedio en el ingreso entre los hombres y mujeres sigue siendo de 9822.53 pesos, aunque al tomar a los hombres como el factor con valor 1, indica un aumento en la variable respuesta.

di_nc %>% lm(ing_cor~sexo_jefe, .) %>% tidy()

Paquetes y funciones útiles

Correlación

Al utilizar la función cor() con un conjunto de datos numéricos obtenemos una matriz de correlaciones entre las distintas variables; esta representación no es la mejor ayuda visual para identificar correlaciones fuertes, por lo que podemos utilizar un correlograma con dicha matriz. Vamos a tomar los datos de gasto-ingreso vistos anteriormente poniendo especial atención sobre el ingreso corriente de los hogares de dicha región.

Aquí se ve una aplicación de la función corrplot::corrplot()

library(corrplot)
corrplot(cor(data_income), type = "lower", diag = F)

Tenemos otras alternativas que nos permiten tener un comportamiento similar al paquete ggplot, como lo es la función ggcorrplot::ggcorrplot().

library(ggcorrplot)
ggcorrplot(cor(data_income), 
           #method = "circle"  #Método de visualización, "square" por default.
           hc.order = TRUE, #Orden jerárquico sobre los valores obtenidos.
           outline.col = "white", #Color del margen de los cuadrádos o círculos.
           type = "lower" #Elementos a desplegar en la gráfica
           #lab = TRUE # Anotaciones de los valores obtenidos en cada celda de la matriz.
           #colors = c("green", "black", "red") #Colores 
           )

Además de un correlograma, podemos obtener las correlaciones más importantes en otro tipo de gráfico, esto es lo que hace el paquete lares con las funciones corr_cross() y corr_var(). Para más información se sugiere revisar el siguiente post

library(lares)

corr_cross(data_income,
  max_pvalue = 0.05,  # Correlaciones significantes a mostrar (a un nivel del 5%)
  top = 10            # Cantidad de variables a mostrar con mayor significancia
) + 
  
corr_var(data_income,         
  ing_cor,                    # Nombre de la variable a comparar
  top = 5                  # Cantidad de variables a mostrar con mayor significancia
)

Igual podríamos elegir una visualización más tradicional pero con un formato más elegante, como el proporcionado con la función modelsummary::datasummary_correlation()

library(modelsummary)
datasummary_correlation(data_income, 
                        method = "pearspear" # "pearson", "kendall", "spearman", or "pearspear" (Pearson correlations above and Spearman correlations below the diagonal)
                        )
edad_jefe tot_integ percep_ing ing_cor alimentos vestido calzado vivienda limpieza salud transporte personales educacion esparcimiento
edad_jefe 1 -.10 .09 .02 -.10 -.10 -.14 -.05 -.03 .08 -.02 -.02 -.09 -.09
tot_integ -.14 1 .67 .12 .24 .10 .22 .00 .03 -.02 .08 .13 .18 -.02
percep_ing .11 .61 1 .18 .23 .13 .18 .02 .06 .00 .10 .13 .15 .05
ing_cor .01 .30 .40 1 .37 .42 .28 .18 .41 .13 .28 .33 .22 .37
alimentos -.14 .32 .27 .52 1 .30 .32 .23 .40 .10 .34 .34 .28 .35
vestido -.18 .18 .21 .38 .37 1 .60 .14 .37 .10 .24 .35 .26 .34
calzado -.22 .31 .23 .30 .34 .54 1 .11 .30 .09 .22 .29 .30 .29
vivienda -.02 .11 .12 .25 .23 .14 .08 1 .23 .11 .15 .17 .15 .21
limpieza -.10 .19 .19 .48 .46 .38 .32 .23 1 .10 .39 .38 .31 .40
salud -.01 .09 .14 .29 .28 .25 .20 .12 .28 1 .07 .12 .07 .05
transporte -.08 .20 .23 .62 .47 .31 .27 .24 .44 .26 1 .31 .17 .28
personales -.16 .34 .26 .52 .50 .45 .40 .26 .58 .30 .47 1 .22 .30
educacion -.23 .45 .28 .30 .31 .25 .36 .11 .25 .15 .29 .32 1 .31
esparcimiento -.13 .04 .11 .41 .36 .32 .24 .14 .34 .20 .37 .33 .21 1

Regresión

  • Con ggplot podemos añadir de manera rápida una recta de regresión a nuestros datos; véase como se hace uso de la base de datos de las olimpiadas y de la función ggplot::geom_smooth() para agregar una capa con la recta de regresión.
olympics <- read_csv("Olympic/athlete_events.csv") %>% distinct()
olympics %>% filter(!is.na(Medal) & !is.na(Age)) %>% 
  filter(!(Year %in% c(1994, 1998, 2002, 2006, 2010, 2014))) %>% #Eliminamos algunos valores inluyentes
  group_by(Year) %>%
  count(Medal) %>% 
  ggplot(aes(x = Year, y = n)) + 
  geom_point() + 
  general_theme + 
  geom_smooth(formula = y~x, method = "lm") +
  labs(x = "Año", y = "Cantidad de medallas")

TidyModels::Broom

library(broom)
glance(car_regression)
tidy(car_regression)

Ejemplo

Consideremos la base de datos de la encuesta ingreso-gasto del 2018 en Aguascalientes.

library(skimr)
fancy_summary <- skim_with(
  numeric = sfl(
    Min = min,
    Max = max,
    Q1 = ~ quantile(., probs = .25),
    Median = ~quantile(., probs = .50),
    Q3 = ~ quantile(., probs = .75),
    Mean = mean,
    Sd = sd,
    hist = ~ inline_hist(., 5)
  ),
  append = FALSE
)
fancy_summary(data_income)
Data summary
Name data_income
Number of rows 2306
Number of columns 14
_______________________
Column type frequency:
numeric 14
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate Min Max Q1 Median Q3 Mean Sd hist
edad_jefe 0 1 16.00 96.00 37.00 47.00 60.00 48.65 15.63 ▃▇▆▃▁
tot_integ 0 1 1.00 16.00 2.00 4.00 5.00 3.91 1.97 ▇▃▁▁▁
percep_ing 0 1 0.00 12.00 1.00 2.00 3.00 2.47 1.38 ▇▅▁▁▁
ing_cor 0 1 1873.77 2310762.00 27415.18 43175.40 67358.96 57622.26 69677.79 ▇▁▁▁▁
alimentos 0 1 0.00 76602.70 6238.78 9681.26 14149.07 11298.88 7755.85 ▇▂▁▁▁
vestido 0 1 0.00 26421.70 0.00 332.60 1261.93 1044.05 2026.59 ▇▁▁▁▁
calzado 0 1 0.00 11739.12 0.00 273.91 893.86 681.48 1115.71 ▇▁▁▁▁
vivienda 0 1 0.00 70664.50 1040.04 1759.72 3060.00 2661.50 3317.12 ▇▁▁▁▁
limpieza 0 1 0.00 40309.56 720.72 1277.64 2492.03 2196.77 2874.19 ▇▁▁▁▁
salud 0 1 0.00 184800.25 0.00 63.58 665.21 1257.52 6207.59 ▇▁▁▁▁
transporte 0 1 0.00 353941.46 1935.70 4447.32 8922.68 7505.52 12313.27 ▇▁▁▁▁
personales 0 1 0.00 113246.10 917.55 1749.66 3275.63 2787.67 4087.78 ▇▁▁▁▁
educacion 0 1 0.00 160801.44 0.00 0.00 3048.38 2788.35 6848.96 ▇▁▁▁▁
esparcimiento 0 1 0.00 40726.22 0.00 170.43 1043.05 1074.29 2473.48 ▇▁▁▁▁

… …

Aquí consideramos un modelo lineal por cada variable y mostramos el resultado de la regresión entre el ingreso y la edad del jefe de familia.

all_linear_models <- data_income %>% 
  dplyr::select(-ing_cor)  %>% map(~summary(lm(data_income$ing_cor~.x)))
all_linear_models$edad_jefe

Call:
lm(formula = data_income$ing_cor ~ .x)

Residuals:
    Min      1Q  Median      3Q     Max 
 -56861  -30361  -14048    9919 2253282 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 53422.63    4745.03   11.26   <2e-16 ***
.x             86.33      92.86    0.93    0.353    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 69680 on 2304 degrees of freedom
Multiple R-squared:  0.0003749, Adjusted R-squared:  -5.895e-05 
F-statistic: 0.8641 on 1 and 2304 DF,  p-value: 0.3527

Para comparar los distintos modelos, podemos hacerlo bajo diferentes criterios; por ejemplo con el \(R^2\)

all_linear_models %>% map_df(~`[[`(.x, "r.squared")) %>% 
  gather("Variable predictora", "r", 1:13) %>% arrange(desc(r))

O por el p-value de la variable independiente

all_linear_models_tidy <- data_income %>% 
  dplyr::select(-ing_cor) %>% map(~tidy(lm(data_income$ing_cor~.x)))
all_linear_models_tidy %>% map_df(~.x[2,"p.value"]) %>% 
  mutate("Variable predictora" = names(all_linear_models_tidy)) %>% 
  dplyr::select("Variable predictora", "p.value") %>% arrange(p.value)

… …

Finalmente, aquí se dejan más enlaces útiles

Otros puntos importantes

Hasta este momento sólo se ha construido un modelo de manera sencilla sin haber considerado factores importantes como la validación de supuestos (lo cual determinará si nuestro modelo tiene veracidad o no), el uso de otro tipo de transformaciones, la comparación de modelos, etc. Todo lo antes mencionado se verá en la regresión múltiple ya que es completamente equivante el enfoque.

# cars %>% ggplot(aes(x = speed, y = dist)) + 
#   labs(x = "Velocidad", y = "Distancia") + general_theme +
#   geom_point(aes(color = "Observaciones")) + 
#   geom_abline(color = "dodgerblue4", 
#               aes(linetype = "best fit"),
#               intercept = car_regression$coefficients[1], 
#               slope = car_regression$coefficients[2]) +
#   scale_color_manual(values = "cornflowerblue", name = "Observaciones") + 
#   scale_linetype_manual(name = "Regresión", values = "solid")

Regresión lineal múltiple

El mundo es complejo y aveces tratar de modelar el comportamiento de un evento sólo por una variable no es lo más conveniente; así que el camino natural para mejorar nuestro modelo de regresión simple es hacer múltiple agregando más variables. Así ya no estaríamos trabajando, geométricamente hablando, en una recta, si no en un hiperplano.Entonces, suponiendo \(p\) predictores, una regresión múltiple queda expresada matematicamente como sigue:

\[ f(X) = Y = \beta_0 +\beta_1X_1+\beta_2X_2+\dots + \beta_pX_p+\epsilon = \beta_0 + \sum_{i = 1}^pX_i\beta_i + \epsilon = X\beta + \epsilon \]

Donde \(\beta_i\) son los parámetros del modelo representando el efecto promedio en \(Y\) de un incremento de una unidad en \(X_i\), manteniendo todos los otros predictores fijos y \(\epsilon\sim N(0,\sigma^2)\).

De igual manera, la forma habitual de estimar los parámetros es mediante el uso de mínimos cuadrados para minimizar la suma de los residuales al cuadrado:

\[ RSS = \sum_{i = 1}^n(y_i-f(x_i))^2 = \sum_{i = 1}^n(y_i-\hat y_i)^2 = \sum_{i = 1}^n(y_i-\hat\beta_0-\hat\beta_1x_{i1}-\hat\beta_2x_{i2}-\cdots-\hat\beta_px_{ip})^2 = \sum_{i = 1}^n\left(y_i-\beta_0-\sum_{j = 1}^px_{ij}\beta_j\right)^2 \]

Vamos a seguir tomando el ejemplo de los ingresos en el estado de Aguascalientes de acuerdo a la Encuesta Nacional de Ingresos y Gastos de los Hogares (ENIGH). 2018 Nueva serie. Vamos a ajustar un modelo con las variables que tengan la mejor correlación con el ingreso eliminando previamente algunos outliers que se consiguen en ciertas variables.

Para facilitarnos esto, veamos el siguiente gráfico de correlación con la correlación de spearman para ver relaciones monotonas y después diagramas de dispersión junto a la correlación de pearson que se consigue por pares.

(data_income %>% filter(ing_cor < 2000000 & transporte < 350000 & limpieza < 40000 & personales < 113246) %>% 
  cor(method = "spearman") %>% 
  ggcorrplot::ggcorrplot(hc.order = TRUE, outline.col = "white", type = "lower", 
                         tl.cex = 8#, lab = T
                         )+
   general_theme + 
  theme(legend.text = element_text(size = 8, face = "plain"), 
        legend.title = element_text(size = 8, face = "plain"), 
        legend.position = c(0.30, 0.85),
        legend.direction = "horizontal",
        legend.box = "horizontal",
        panel.background = element_blank())) + 
(data_income %>% filter(ing_cor < 2000000 & transporte < 350000 & limpieza < 40000 & personales < 113246) %>%
   dplyr::select(ing_cor, transporte, alimentos, limpieza, personales) %>%
   GGally::ggscatmat() + 
   general_theme + 
   theme(axis.text = element_text(size = 8), 
         axis.text.x = element_text(angle = 45)))

Por lo que proponemos distintas transformaciones; el resultado de esto lo podemos ver en la siguiente matriz de dispersión

data_income %>% filter(ing_cor < 2000000 & transporte < 350000 & limpieza < 40000 & personales < 113246) %>%
  dplyr::select(ing_cor, transporte, alimentos, limpieza, personales) %>% 
  mutate(ing_cor = log(ing_cor),
         transporte = sqrt(transporte),
         alimentos = sqrt(alimentos),
         limpieza = sqrt(limpieza),
         personales = sqrt(personales)) %>% 
  GGally::ggscatmat() + general_theme

Entonces, ajustando un modelo de regresión con los siguientes datos obtenemos el siguiente modelo junto con la implementación de dicho modelo en R:

mutated_data_income <- data_income %>% filter(ing_cor < 2000000 & transporte < 350000 & limpieza < 40000 & personales < 113246) 
  #dplyr::select(ing_cor, transporte, alimentos, limpieza, personales) 
#%>% dplyr::select(ing_cor_log, transporte_sqrt, alimentos_sqrt, limpieza_sqrt, personales_sqrt)

first_m_model <- mutated_data_income %>% lm(log(ing_cor)~sqrt(transporte) + sqrt(alimentos) + sqrt(limpieza) + sqrt(personales), data = .)
first_m_model %>% summary()

Call:
lm(formula = log(ing_cor) ~ sqrt(transporte) + sqrt(alimentos) + 
    sqrt(limpieza) + sqrt(personales), data = .)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.5269 -0.3109 -0.0043  0.3218  2.6055 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)      9.3424687  0.0351277 265.957  < 2e-16 ***
sqrt(transporte) 0.0053498  0.0002970  18.010  < 2e-16 ***
sqrt(alimentos)  0.0049731  0.0003956  12.571  < 2e-16 ***
sqrt(limpieza)   0.0045630  0.0006180   7.383 2.15e-13 ***
sqrt(personales) 0.0054216  0.0005851   9.266  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5266 on 2297 degrees of freedom
Multiple R-squared:  0.479, Adjusted R-squared:  0.4781 
F-statistic: 528.1 on 4 and 2297 DF,  p-value: < 2.2e-16
library(equatiomatic)
extract_eq(first_m_model)

\[ \operatorname{log(ing\_cor)} = \alpha + \beta_{1}(\operatorname{sqrt(transporte)}) + \beta_{2}(\operatorname{sqrt(alimentos)}) + \beta_{3}(\operatorname{sqrt(limpieza)}) + \beta_{4}(\operatorname{sqrt(personales)}) + \epsilon \]

O mejor dicho

\[ \log(ing\_cor) = \alpha + \beta_{1}\sqrt{transporte} + \beta_{2}\sqrt{alimentos} + \beta_{3}\sqrt{limpieza} + \beta_{4}\sqrt{personales} + \epsilon \]

Como vemos, con todas las variables se rechaza la hipótesis \(H_0: \beta_i = 0\), además de que también se rechaza la hipótesis \(H_0: \beta_1 = \beta_2 = \beta_3 = \beta_4 = 0\) con la prueba ANOVA y tenemos un \(R^2\) ajustado del 0.4781, por lo que el modelo recupera un 47% de la variabilidad de los datos, lo cual no significa que sea un mal modelo. Veamos un poco más a detalle que significa todo esto, además de hacer ver posibles problemas que hagan que nuestro modelo no sea adecuado.

Pruebas de hipótesis y ANOVA e intervalos

Lo primero que nos otorga el summary del modelo anterior es lo siguiente:

broom::tidy(first_m_model)

Con esto obtenemos:

  1. Los coeficientes de nuestro modelo (columna estimate)
  2. El error estandar sobre para parámetro (columna std.error)
  3. El estadístico de la prueba sobre los coeficientes (columna statistic)
  4. El \(p-value\) sobre la prueba mencionada en el anterior punto (columna p.value)

Hay que recordar que el error estandar nos ayuda a crear los intervalos de confianza sobre cada parámetro del modelo. Tomemos como ejemplo el termino \(\beta_1\) asociado a la raíz cuadrada del gasto en transporte. El intervalo de confianza en este caso sería, aproximadamente, el siguiente:

\[ \hat{\beta_i}\pm 2\cdot SE(\hat{\beta_i}) = \left[\hat{\beta_1}-2\cdot SE(\hat{\beta_1}), \hat{\beta_1}+2\cdot SE(\hat{\beta_1})\right] = [0.004755661, 0.005943853] \]

Lo que significa que, hay aproximadamente un 95% de probabilidad de que el intervalo \([0.004755661, 0.005943853]\) contiene el verdadero valor de \(\beta_1\). También podríamos haber utilizado el hecho de que suponemos en la construcción de intervalos de confianza que todos los coeficientes se distribuyen de manera normal (\(\hat{\beta_i} = N(\beta_i, \sigma^2C_{(i+1)(i+1)})\) donde \(C_{(i+1)(i+1)})\) es el i-ésimo coeficiente de la diagonal de la matriz \((X'X)^{-1}\)) y haber calculado los intervalos como:

\[ \hat{\beta_i}\pm 1.96\cdot SE(\hat{\beta_i}) = \left[\hat{\beta_1}-2\cdot SE(\hat{\beta_1}), \hat{\beta_1}+1.96\cdot SE(\hat{\beta_1})\right] = [0.004755661, 0.005943853] \]

Aunque para ser más precisos, los intervalos de confianza, de manera general, están determinados por la siguiente ecuación:

\[ \hat{\beta_i}\pm t_{n-k-1}^{\alpha/2} \sqrt{\hat{\sigma}^2C_{(i+1)(i+1)}} \]

Para nuestro caso podemos calcularlos de la siguiente manera:

sfm <- summary(first_m_model)
interval_confidence_firstM <- matrix(
  c(sfm$coefficients[,1] - qt(0.975, df = sfm$df[2]) * sfm$coefficients[, 2],
  sfm$coefficients[,1] + qt(0.975, df = sfm$df[2]) * sfm$coefficients[, 2]),
  ncol = 2
)
row.names(interval_confidence_firstM) <- row.names(sfm$coefficients)
colnames(interval_confidence_firstM) <- c("lower", "upper")
interval_confidence_firstM
                       lower       upper
(Intercept)      9.273583316 9.411354045
sqrt(transporte) 0.004767247 0.005932267
sqrt(alimentos)  0.004197295 0.005748825
sqrt(limpieza)   0.003351097 0.005774906
sqrt(personales) 0.004274251 0.006569003

Los cuales se pudieron haber obtenido con la función stats::confint.lm()

confint(first_m_model)
                       2.5 %      97.5 %
(Intercept)      9.273583316 9.411354045
sqrt(transporte) 0.004767247 0.005932267
sqrt(alimentos)  0.004197295 0.005748825
sqrt(limpieza)   0.003351097 0.005774906
sqrt(personales) 0.004274251 0.006569003

Una manera elegante que podemos utilizar para visualizar la anterior información no las proporciona la función modelsummary::modelplot()

modelplot(first_m_model) +  general_theme + theme(panel.background = element_blank())

Por la escala que tenemos en nuestro intercepto, tenemos poca apreciación del resto de coeficientes

modelplot(first_m_model, coef_omit = "Intercept") + general_theme + theme(panel.background = element_blank())

Retomando el ejemplo, podemos decir que en la ausencia de cualquiera de nuestras variables, el ingreso, en promedio, estará entre $10652.86 (=exp(9.273583316)) y 12226.41 pesos mexicanos. Y, por ejemplo, con un aumento de $1,000 pesos en el gasto de alimentos, el ingreso debe aumentar, en promedio para los ciudadanos de Aguascalientes, entre $1000.018 y $1000.033 pesos mexicanos.

Respecto a los valores obtenidos en la columna statistic, estos corresponden a un estadístico \(t\) para determinar si hay una relación entre la variable asociada a dicho parámetro y la variable a predecir. Para tales fines nuestro estadístico mide el número de desviaciones estandar que nuestro coeficientes lejanas desde el 0, es decir:

\[ t = \frac{\hat{\beta_i}-0}{SE(\hat{\beta_i})} \]

Y la prueba de hipótesis en la que se utiliza dicho estadístico queda determinada de la siguiente manera:

\[ \begin{array}{c} H_0: \mbox{No hay alguna relación entre }X \mbox{ y }Y \equiv \beta_i = 0\\ H_a: \mbox{Hay alguna relación entre }X \mbox{ y }Y \equiv \beta_i \neq 0 \end{array} \]

Entonces lo que buscamos es que \(\hat{\beta_i}\) este lo más alejado del 0, es decir que si \(SE(\hat{\beta_i})\) es pequeño, \(\hat{\beta_i}\) puede ser pequeño y si \(SE(hat{\beta_i})\) es grande, entonces \(\hat{\beta_i}\) debe ser lo suficientemente grande en valor absoluto. En la siguiente gráfica podemos ver que todas nuestras variables son relevantes para modelar el ingreso (gracias a que los \(p-values\) indican de la información no es compatible con la hipótesis nula), así como una comparación entre los valores de cada uno de los coeficientes.

GGally::ggcoef_model(first_m_model) + general_theme

Por lo que, al considerar que la raíz cuadrada es una función monótona y creciente, al igual que el logaritmo en el rango que estamos considerando (\(\mathbb{R^{+}}\)), los cambios en los gastos personales, en promedio, son lo que más puede afectar en el ingreso si los demás gastos permanecen constantes; de hecho por cada $1,000 más en los gastos personales, se necesitará entre $1000.018 y $1000.043 pesos mexicanos que aumente el ingreso.

Lo último que se visualiza en un summary de un modelo lineal son algunas estadísticas sobre el rendimiento del modelo e información sobre la prueba ANOVA, por el momento sólo se analizará esto último.

glance(first_m_model) %>% dplyr::select(statistic, p.value, df)

La prueba ANOVA (Análisis de Varianza) es un prueba que nos ayudará a determinar la significancia del modelo, es decir que se desea probar si la ecuación de regresión no explica una proporción considerable de la variabilidad en la variable respuesta, contra la hipótesis alternativa de que sí la explica, es decir:

\[ \begin{array}{c} H_0 : \beta_1 = \beta_2 = \cdots = \beta_p = 0\\ vs\\ H_1 : \beta_i \neq 0 \mbox{ p.a }i, i\in\{1,2,\dots, p\} \end{array} \]

El estadístico para esta prueba hace entender el nombre de la prueba y para obtenerlo necesitamos diferentes estadísticas. De acuerdo a la bibliografía se pueden encontrar diferentes nomenclaturas, aquí se colocan algunas ejemplos:

  • \(TSS = SC_{total} = \sum(y_i-\bar{y_i})^2\): Total sum of squares o la suma de cuadrados total.
  • \(SC_{reg} = SSM = SSR = \sum(\hat{y}_i-\bar{y})^2\): La suma de cuadrados de la regresión o del modelo.
  • \(RSS = SC_{error} = SSE = \sum(y_i-\hat{y_i})^2\): La suma de cuadrados del error o de los residuales

La segunda estadística la podemos interpretar como la cantidad de varianza explicada por la regresión, por el modelo o por las variables; mientras la tercera se puede interpretar como la cantidad de varianza que no se explica por la regresión.

  • \(CM_{reg} = SC_{reg}/p; MSM = MSR = SSM/p\): Cuadrado medio de la regresión o Mean Square of model/regression
  • \(CM_{error} = SC_{error}/n-p-1; MSE = SSE/n-p-1\): Cuadrado medio del error

Y con esto ya podemos obtener nuestro estadístico \(F\):

  • \(F = \frac{CM_{reg}}{CM_{error}}\)
  • \(F = \frac{MSM}{MSE}\)

\[ F = \frac{MSR}{MSE} = \frac{CM_{reg}}{CM_{error}} = \frac{(TSS-RSS)/p}{RSS/(n-p-1)} \]

Este estadístico tiene una distribución \(F_{p, n-p-1}\) y rechazaremos la hipótesis nula cuando \(F>F^{(\alpha)}_{k, n-k-1}\). Si los supuestos del modelo de regresión lineal son validos, se puede ver que \(\mathbb{E}[RSS/(n-p-1)] = \sigma^2\), es decir la \(\sigma^2\) estimada para todo el modelo, y bajo la hipótesis nula \(\mathbb{E}[(TSS-RSS)/p] = \sigma^2\). Por lo que cuando NO hay relación entre la variable dependiente y los predictores, se esperaría que el estadístico \(F\) sea cercano a 1.

Si \(H_a\) es cierta \(\mathbb{E}[(TSS-RSS)/p] > \sigma^2\), por lo que tendríamos valores más grandes que 1 en el estadístico.

Entonces podemos interpretar este estadístico como la razón entre la variabilidad explicada por los regresores entre la variabilidad no explicada por el modelo, ponderando con los respectivos grados de libertad que contiene cada estadístico. Por lo que buscamos que nuestro modelo contenga la mayor cantidad de información proporcionada por los datos. Sólo para aclarar, \(TSS-RSS = SC_{total}-SC_{error} = SC_{reg}\) y esto es gracias a la siguiente igualdad

\[ SC_{total} = SC_{reg} + SC_{error} \]

La cual proviene de la ley total de la varianza donde los sumandos igual pueden ser interpretados como la cantidad de varianza explicada y no explicada:

\[ Var(Y) = \mathbb{E}[Var(Y|X)] + Var(\mathbb{E}[Y|X]) \]

Esto puede ser entendido de manera sencilla con la siguiente gráfica

drawing


Para resumir la información de dicha prueba tenemos la tabla ANOVA:

\[ \begin{array}{|c| c| c| c| c|} \hline &Grados\ de\ libertad & Suma\ de\ Cuadrados & Cuadrado\ Medio & Prueba\ F \\ \hline \hline Regresión & k & SSR = \hat{\beta}'X'\underline{Y}-n\overline{y}^2 & MSR = \frac{SC_{reg}}{k} & \frac{CM_{reg}}{CM_{error}} \\ \hline Error & n-k-1& SSE = \underline{Y}'(I-H)\underline{Y} & MSE = \frac{SC_{error}}{n-k-1} & -\\ \hline Total & n-1 & TSS = \underline{Y}'\underline{Y}-n\overline{y}^2 & - & - \\ \hline \end{array} \]

Y en R podemos utilizar la siguiente función con nuestro modelo previamente ajustado

anova(first_m_model)
m1 <- mutated_data_income %>% lm(ing_cor_log~transporte_sqrt + alimentos_sqrt +limpieza_sqrt + personales_sqrt, data = .)
m2 <-  mutated_data_income %>% lm(ing_cor_log~ alimentos_sqrt +limpieza_sqrt + personales_sqrt, data = .)
anova(m2, m1)

Aquí se dejan algunos enlaces que pueden ayudar a aclarar todos los conceptos

Interacción y selección de variables

Entre todos los factores que tenemos que considerar en nuestro modelo lineal múltiple, tenemos dos que pueden causar graves problemas: Aditividad y linearidad entre los predictores y la variable respuesta. El concepto (y supuesto) de Linearidad es claro para este momento (1 cambio en X implica un cambio constante en Y), por otro lado hemos utilizado la aditividad en nuestras interpretaciones de los coeficientes.

La aditividad significa que el efecto que tiene una variable predictora sobre el modelo es independiente de otra regresora sobre la variable respuesta lo cual podría no suceder. Véamos un ejemplo con la misma base de los ingresos en Aguascalientes. Para esto vamos a considerar un modelo más grande, tomando en cuenta las variables de calzado y vestido. Primero veamos la siguiente gráfica en la cual se plasman los datos ya transformados, considerando los filtros anteriores, agregando la significancia de la variable regresora sobre modelos lineales simples

m_calzado <- lm(log(ing_cor)~sqrt(calzado) , data = mutated_data_income) %>% summary()
m_vestido <- lm(log(ing_cor)~sqrt(vestido) , data = mutated_data_income) %>% summary()

(mutated_data_income %>%
  ggplot(aes(x = sqrt(calzado), y = log(ing_cor))) +
  geom_point() + geom_smooth(method = "lm") + 
  annotate("text", x = 70, y = 9, family="Dosis",
           label = paste0("p-value: ", m_calzado$coefficients[2,4])) + 
  general_theme) +
(mutated_data_income %>% 
  ggplot(aes(x = sqrt(vestido), y = log(ing_cor))) +
  geom_point() + geom_smooth(method = "lm") + 
  labs(y = NULL) + 
  annotate("text", x = 100, y = 9, family="Dosis",
           label = paste0("p-value: ", m_vestido$coefficients[2,4])) + 
  general_theme)

Como vemos, ambas variables variables tiene un p-value muy pequeño para no considerarlas importantes en nuestro modelo, veamos si mejora nuestro modelo con dichas variables y de paso veamos el uso de la función stats::update()

update(first_m_model, . ~ .+ sqrt(calzado) + sqrt(vestido), data = mutated_data_income) %>% summary()

Call:
lm(formula = log(ing_cor) ~ sqrt(transporte) + sqrt(alimentos) + 
    sqrt(limpieza) + sqrt(personales) + sqrt(calzado) + sqrt(vestido), 
    data = mutated_data_income)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.4248 -0.3146 -0.0112  0.3175  2.6405 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)       9.3719803  0.0352492 265.878  < 2e-16 ***
sqrt(transporte)  0.0052574  0.0002955  17.789  < 2e-16 ***
sqrt(alimentos)   0.0046875  0.0003982  11.773  < 2e-16 ***
sqrt(limpieza)    0.0040612  0.0006197   6.554 6.90e-11 ***
sqrt(personales)  0.0045262  0.0006040   7.493 9.52e-14 ***
sqrt(calzado)    -0.0001498  0.0007305  -0.205    0.838    
sqrt(vestido)     0.0032308  0.0006085   5.309 1.21e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.523 on 2295 degrees of freedom
Multiple R-squared:  0.4867,    Adjusted R-squared:  0.4854 
F-statistic: 362.7 on 6 and 2295 DF,  p-value: < 2.2e-16

Comparado con el anterior modelo, nuestra \(R^2\) ajustada aumento y sigue siendo válida la prueba ANOVA, pero ¿Qué notamos ahora en nuestro modelo?

Al considerar la variable \(\sqrt{calzado}\), esta ya no es importante para el modelo de acuerdo a la prueba \(H_0: \beta_i = 0\). ¿Esto por qué sucede?

Al parecer el efecto que tenía la raíz cuadrado del calzado sobre el logaritmo de los ingresos ya no es relevante y esto es porque de alguna manera esta relación queda explicada por la otra variable que agregamos, es decir: \(\sqrt{vestido}\). Lo cual tiene sentido ya que si se hace un gasto en el vestido, es muy probable que sea haga un gasto en el calzado y si se hace un gasto en calzado, seguramente los encuestados consideraron esto en la variable vestido. De hecho, veamos que tan correlacionadas están estas variables

mutated_data_income %>% dplyr::select(vestido, calzado, ing_cor) %>% cor() %>% ggcorrplot(type = "lower", lab = T) + general_theme + 
  theme(legend.text = element_text(size = 8, face = "plain"), 
        legend.title = element_text(size = 8, face = "plain"), 
        legend.direction = "vertical",legend.position = "right",
        panel.background = element_blank())

Entoces, al tener que los aumentos en el calzado aumentaran los del vestido y veceversa, esto viola el supuesto de aditividad en el modelo lineal. Para solucionar esto agregamos una interacción al modelo con estas variables, es decir que agregaremos el término \(\beta_iX_{calzado}X_{vestido}\) al modelo, entonces ahora nuestro modelo sería el siguiente:

\[ \begin{split} \log(ing\_cor) = & \beta_{0} + \beta_{1}\sqrt{transporte} + \beta_{2}\sqrt{alimentos} + \beta_{3}\sqrt{limpieza} + \beta_{4}\sqrt{personales} + \beta_5\sqrt{calzado} + \beta_6\sqrt{vestido} + \beta_7\sqrt{calzado}\times \sqrt{vestido} + \epsilon\\ & \beta_{0} + \beta_{1}\sqrt{transporte} + \beta_{2}\sqrt{alimentos} + \beta_{3}\sqrt{limpieza} + \beta_{4}\sqrt{personales} + \widetilde{\beta_5}\sqrt{calzado} + \beta_6\sqrt{vestido} + \epsilon \end{split} \]

En la última expresión suponemos que \(\widetilde{\beta_5} = \beta_5+\beta_7\sqrt{X_{vestido}}\) y así podemos seguir viendo nuestro modelo “sin interacciones,” sólo que ahora un cambio en el calzado tendrá un efecto en el vestido y veceversa

Para agregar una interacción entre variables, utilizamos var1:var2 en la formula. Veamos el summary que obtenemos de nuestro nuevo modelo

second_m_model <- update(first_m_model, . ~ .+ sqrt(calzado) + sqrt(vestido) + sqrt(calzado):sqrt(vestido),
       data = mutated_data_income) 
second_m_model %>% summary()

Call:
lm(formula = log(ing_cor) ~ sqrt(transporte) + sqrt(alimentos) + 
    sqrt(limpieza) + sqrt(personales) + sqrt(calzado) + sqrt(vestido) + 
    sqrt(calzado):sqrt(vestido), data = mutated_data_income)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.3987 -0.3134 -0.0097  0.3153  2.6469 

Coefficients:
                              Estimate Std. Error t value Pr(>|t|)    
(Intercept)                  9.351e+00  3.592e-02 260.359  < 2e-16 ***
sqrt(transporte)             5.253e-03  2.951e-04  17.802  < 2e-16 ***
sqrt(alimentos)              4.568e-03  3.997e-04  11.431  < 2e-16 ***
sqrt(limpieza)               4.136e-03  6.192e-04   6.680 3.00e-11 ***
sqrt(personales)             4.578e-03  6.033e-04   7.588 4.68e-14 ***
sqrt(calzado)                1.467e-03  9.206e-04   1.594  0.11108    
sqrt(vestido)                4.718e-03  7.974e-04   5.916 3.79e-09 ***
sqrt(calzado):sqrt(vestido) -5.275e-05  1.832e-05  -2.879  0.00403 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5221 on 2294 degrees of freedom
Multiple R-squared:  0.4886,    Adjusted R-squared:  0.487 
F-statistic: 313.1 on 7 and 2294 DF,  p-value: < 2.2e-16

Es evidente, bajo una significancia del 5% que la ponderación de la interacción es estadísticamente diferente de cero, dándole un peso relevante a nuestro modelo. Con esto podemos decir que por cada $100 pesos de aumento en el calzado, tendremos un aumento de \((\beta_5+\beta_7\sqrt{X_{vestido}})\times 10 = (0.001467 + 0.004718\times X_{vestido})\times 10 = 0.01467 + 0.04718\times X_{vestido}\) unidades sobre el logaritmo de los ingresos.

Véase que seguimos viendo que \(\sqrt{calzado}\) sigue sin ser relevante.entonces ¿Por qué no consideramos un modelo sin esta variable?

Para evitar inconsistencias en los modelos, si se realiza una interacción, en el modelo se debe tener presente las variables de la interacción, a esto se le llama principio jerárquico, independientemente si el coeficiente de alguna de las variables es estadísticamente 0 o no. ¿Qué pasaría por ejemplo si omitimos el calzado en nuestro modelo? en tal caso, si hay un aumento en calzado pero no hay presencial del gasto en vestido \(\beta_7\sqrt{calzado\times vestido} = 0\) cuando debe ser así ya que realmente hubo un aumento en al menos una de las variables. Además de cumplir ahora el supuesto de aditividad, véase que nuestro modelo mejoro en comparación del anterior.

second_m_model %>% 
  glance() %>% gather("Estadística", "Nuevo modelo") %>% 
  inner_join(
    first_m_model %>% glance() %>%  gather("Estadística", "Anterior modelo"),
    by = "Estadística"
  )

Con este nuevo modelo aumento el \(R^2\) ajustada, el \(p-value\) de la prueba ANOVA indica que nuestras variables son significativas y las estadísticas AIC y BIC disminuyeron.

La interacción con datos categóricos cambia un poco en el sentido de que puede ser más fácil identificar dichas interacciones.

Considerando el AIC, por ejemplo, podemos utilizar algunas técnicas para elegir el mejor modelo, de acuerdo a una estadística como esta, omitiendo diferentes variables al modelo o agregándolas. Estas técnicas son llamadas Forward Selection, Backward Selection y Mixed Selection en los cuales simplemente se comienza con un modelo base y se van agregando o eliminando variables de acuerdo a la significancia que se vaya obteniendo con los coeficientes o con otro criterio como el AIC.

Veamos como aplicarlo a nuestro modelo

best_second_model <- second_m_model %>% MASS::stepAIC(direction = "backward")
Start:  AIC=-2983.78
log(ing_cor) ~ sqrt(transporte) + sqrt(alimentos) + sqrt(limpieza) + 
    sqrt(personales) + sqrt(calzado) + sqrt(vestido) + sqrt(calzado):sqrt(vestido)

                              Df Sum of Sq    RSS     AIC
<none>                                     625.42 -2983.8
- sqrt(calzado):sqrt(vestido)  1     2.260 627.68 -2977.5
- sqrt(limpieza)               1    12.164 637.58 -2941.4
- sqrt(personales)             1    15.699 641.11 -2928.7
- sqrt(alimentos)              1    35.621 661.04 -2858.3
- sqrt(transporte)             1    86.397 711.81 -2687.9
summary(best_second_model)

Call:
lm(formula = log(ing_cor) ~ sqrt(transporte) + sqrt(alimentos) + 
    sqrt(limpieza) + sqrt(personales) + sqrt(calzado) + sqrt(vestido) + 
    sqrt(calzado):sqrt(vestido), data = mutated_data_income)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.3987 -0.3134 -0.0097  0.3153  2.6469 

Coefficients:
                              Estimate Std. Error t value Pr(>|t|)    
(Intercept)                  9.351e+00  3.592e-02 260.359  < 2e-16 ***
sqrt(transporte)             5.253e-03  2.951e-04  17.802  < 2e-16 ***
sqrt(alimentos)              4.568e-03  3.997e-04  11.431  < 2e-16 ***
sqrt(limpieza)               4.136e-03  6.192e-04   6.680 3.00e-11 ***
sqrt(personales)             4.578e-03  6.033e-04   7.588 4.68e-14 ***
sqrt(calzado)                1.467e-03  9.206e-04   1.594  0.11108    
sqrt(vestido)                4.718e-03  7.974e-04   5.916 3.79e-09 ***
sqrt(calzado):sqrt(vestido) -5.275e-05  1.832e-05  -2.879  0.00403 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5221 on 2294 degrees of freedom
Multiple R-squared:  0.4886,    Adjusted R-squared:  0.487 
F-statistic: 313.1 on 7 and 2294 DF,  p-value: < 2.2e-16

En este caso tenemos un bueno modelo ya que algorítmicamente si eliminamos alguna de las variables, perderíamos información con el subsequente modelo (esto mediante el AIC). Sólo para ejemplificar, veáse que pasaría con nuestros datos originales

library(MASS)
model_zero <- lm(ing_cor~., data = mutated_data_income) %>% stepAIC(direction = "backward")
Start:  AIC=48539.95
ing_cor ~ edad_jefe + tot_integ + percep_ing + alimentos + vestido + 
    calzado + vivienda + limpieza + salud + transporte + personales + 
    educacion + esparcimiento

                Df  Sum of Sq        RSS   AIC
- tot_integ      1 1.8169e+09 3.2703e+12 48539
<none>                        3.2685e+12 48540
- calzado        1 4.8737e+09 3.2734e+12 48541
- educacion      1 8.1884e+09 3.2767e+12 48544
- vestido        1 8.6011e+09 3.2771e+12 48544
- vivienda       1 1.0226e+10 3.2787e+12 48545
- edad_jefe      1 2.4454e+10 3.2930e+12 48555
- salud          1 3.5137e+10 3.3037e+12 48563
- percep_ing     1 5.0741e+10 3.3193e+12 48573
- personales     1 8.5543e+10 3.3541e+12 48597
- limpieza       1 8.8024e+10 3.3565e+12 48599
- esparcimiento  1 8.9726e+10 3.3582e+12 48600
- transporte     1 1.5210e+11 3.4206e+12 48643
- alimentos      1 1.7074e+11 3.4393e+12 48655

Step:  AIC=48539.23
ing_cor ~ edad_jefe + percep_ing + alimentos + vestido + calzado + 
    vivienda + limpieza + salud + transporte + personales + educacion + 
    esparcimiento

                Df  Sum of Sq        RSS   AIC
<none>                        3.2703e+12 48539
- calzado        1 5.6047e+09 3.2759e+12 48541
- educacion      1 7.5215e+09 3.2779e+12 48543
- vestido        1 9.2882e+09 3.2796e+12 48544
- vivienda       1 1.0559e+10 3.2809e+12 48545
- edad_jefe      1 2.8199e+10 3.2985e+12 48557
- salud          1 3.5594e+10 3.3059e+12 48562
- percep_ing     1 6.8639e+10 3.3390e+12 48585
- personales     1 8.4152e+10 3.3545e+12 48596
- limpieza       1 8.9924e+10 3.3603e+12 48600
- esparcimiento  1 9.5439e+10 3.3658e+12 48603
- transporte     1 1.5196e+11 3.4223e+12 48642
- alimentos      1 1.6893e+11 3.4393e+12 48653

Véase que con todos los datos y sin ninguno tipo de limpieza el AIC determinado por el segundo modelo (se elimino la variable tot_integ) es mejor. Un punto importante que hay que remarcar es que este modelo no significa que sea el mejor que podamos obtener, de hecho nuestro último modelo (second_model) es hasta ahora el mejor.

Por razones que se explorarán posteriormente, vamos a eliminar la interacción y las variables calzada y vestido ya que el modelo es mejor sin estas variables (es mejor en los supuestos, lo cual nos importará más que obtener el mejor AIC o mejor \(R^2\)).

Supuestos y problemas potenciales

Antes de pasar con los supuestos del modelo lineal, véamos un poco más a detalle las transformaciones más comunes para obtener un comportamiento lineal entre nuestros predictores y la variable respuesta.

+Transformaciones

Hasta el momento se han mencionado algunas transformaciones que nos pueden ayudar a conseguir linealidad entre dos variables, por ejemplo cuando deseamos reducir la variabilidad de la variable a predecir o cambiar la escala de la variable predictora:

  • \(log()\): Función monótona creciente definida para los positivos e indefinida para el 0.
  • \(sqrt()\): Función monótona creciente definida para los positivos y el 0.
  • \(1/x\): Función monótona creciente definida para los positivos e indefinida para el 0.

De acuerdo al comportamiento de los datos podemos aplicar cualquiera de estas variables siempre que este bien definida la función. ¿Qué podemos hacer cuando tenemos ceros? ¿Qué otras propuestas comunes existen?

Para el primer punto hay que pensar que los ceros pueden tener diferentes significados en nuestras variables, en una pregunta que se realizo en StackExchange: How should I transform non-negative data including zeros? se porponen diferentes puntos de vista validos que hay que pensar sobre los ceros en nuestros datos:

  • Truncamiento: Puede ser que un modelaje para los datos sea un modelo de mezclas, un modelo de supervivencia, etc.
  • Valores perdidos: Es natural que esto puedan ser valores perdidos y no un valor obtenido de la variable, por lo que se puede imputar información o eliminarla si se considera apropiado.
  • Punto cero natural: Por ejemplo con niveles de ingreso (puede tratarse de una situación de desemepleo) y en tal caso una transformación sería lo indicado.
  • Sensibilidad del instrumento de medición: Una solución podría ser agregar un pequeño valor a la información.

Sea cual sea el caso, podemos eliminarlos si no perdemos una gran cantidad de información o usar dichos datos como un segregador de la información. También podemos utilizar una transformación donde no existan problemas en este valor específico y una de las más populares es una de las parametrizaciones de la familia de transformaciones Box-Cox:

\[ g_{\lambda}(y) = \Bigg\{\begin{array}{rl}\frac{y^{\lambda}-1}{\lambda} & \lambda \neq 0\\ \log(y) & \lambda = 0\end{array} \]

Donde la función inversa es:

\[ y = \Bigg\{\begin{array}{cc}\exp\left(\frac{\log(1+\lambda g_{\lambda}(y))}{\lambda}\right) & \lambda \neq 0 \\ \exp(g_{\lambda}(y) & \lambda = 0\end{array} \]

Y para diferentes valores de \(\lambda\) tenemos diversas funciones asociadas:

  • \(\lambda = 1\): No se requiere alguna transformación ya que se producen resultados idénticos a los originales y no se alteraría la variabilidad
  • \(\lambda = 0.50\): Raíz cuadrada
  • \(\lambda = 0.33\): Raíz cúbica
  • \(\lambda = 0.25\): Raíz cuarta
  • \(\lambda = 0\): Logaritmo natural
  • \(\lambda = -0.50\): raíz cuadrada recíproca
  • \(\lambda = -1\): Transformación inversa

La obtención de \(\lambda\) es mediante el método de máxima verosimilitud para obtener el valor más probable que haya dado lugar a la distribución observada de los residuos del modelo

\[ \log\mathbb{L} = -\frac{n}{2}\log(2\pi) - n\log \sigma - \frac{1}{2\sigma^2}\sum_{i = 1}^n\left[g_{\lambda}(y)-(\beta_0+\beta X)\right]^2 + (\lambda-1)\sum_{i = 1}^n \log Y_i \]

Este es obtenido de manera numérica y en general se elige un valor en el intervalo de confianza para dicho parámetro

\[ \left\{\lambda: \mathbb{L}(\lambda) > \mathbb{L}(\hat{\lambda}-\frac{1}{2}\chi_{1,\alpha}^2)\right\} \]

Tal información la podemos obtener de la función MASS::boxcox() la cual proporciona un gráfica donde podemos apreciar el intervalo de confianza y la verosimilitud para diferentes valores de \(\lambda\). Como argumento le pasamos a la función el modelo de regresión que hemos creado.

Es común que esta transformación es aplicada a la variable respuesta, ya que al hacer esto generalmente se reduce la no normalidad de los residuales, al igual que la no linearidad (relacionado con la independencia de estos) y se mejoran problemas de homocedasticidad.

Con nuestros datos, se decidió eliminar las observaciones donde se tengan ceros en las variables predictoras (ya que estos representan menos del 10% de información y, se puede comprobar, dejando dichas observaciones se obtienen diversos problemas con las pruebas de hipótesis y el análisis de valores influyentes). También se cambiaron algunas transformaciones ya que, al no tener valores que afecten el logaritmo se mejoró la relación lineal de algunos predictores con la variable respuesta

library(MASS)
filter_data <- data_income_sexo %>% 
  filter(ing_cor < 2000000 & transporte < 350000 & limpieza < 40000 & personales < 113246) %>% 
  filter(transporte != 0 & limpieza != 0 & personales != 0 & alimentos != 0) %>% 
  dplyr::select(ing_cor, transporte, alimentos, limpieza, personales)

third_model <- lm(ing_cor~log(transporte) + sqrt(alimentos) + log(limpieza) +
                    log(personales), data = filter_data)

boxCox_third_model <- boxcox(third_model, 
                             plotit = TRUE, #Deseamos ver la gráfica
                             lambda = seq(-0.2, 0.2, by = 0.1) # Con esto modificamos el rango
                             )

Entonces no estábamos mal al proponer al logaritmo en el ingreso como una transformación que ayudaría a nuestro modelo lineal. Podemos obtener el valor exacto en el que se maximiza la verosimilitud

(lambdaBCox_tM <- boxCox_third_model$x[which.max(boxCox_third_model$y)])
[1] 0.02222222

Es decir, esto nos sugiere hacer un modelo con la siguiente variable repuesta

\[ \frac{y^{\lambda}-1}{\lambda} = \frac{y^{0.022}-1}{0.022} \]

Y con esto tendríamos el siguiente summary del correspondiente modelo:

third_model_BX <- lm(((ing_cor^lambdaBCox_tM-1)/lambdaBCox_tM) ~ 
                  log(transporte) + sqrt(alimentos) + log(limpieza) + log(personales), 
                data = filter_data)
third_model_BX %>% summary()

Call:
lm(formula = ((ing_cor^lambdaBCox_tM - 1)/lambdaBCox_tM) ~ log(transporte) + 
    sqrt(alimentos) + log(limpieza) + log(personales), data = filter_data)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.3384 -0.4054 -0.0169  0.3843  3.3192 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)     6.9124443  0.1361592  50.767   <2e-16 ***
log(transporte) 0.2564575  0.0142296  18.023   <2e-16 ***
sqrt(alimentos) 0.0060496  0.0005121  11.813   <2e-16 ***
log(limpieza)   0.1503380  0.0180963   8.308   <2e-16 ***
log(personales) 0.1762014  0.0186708   9.437   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6418 on 2139 degrees of freedom
Multiple R-squared:  0.4886,    Adjusted R-squared:  0.4876 
F-statistic: 510.9 on 4 and 2139 DF,  p-value: < 2.2e-16

¿Cómo se compara con nuestro anterior mejor modelo?

second_m_model %>% 
  glance() %>% gather("Estadística", "Viejo modelo") %>% 
  inner_join(
    third_model_BX %>% glance() %>%  gather("Estadística", "Nuevo modelo") %>%
      dplyr::mutate_if(is.numeric, round, digits = 8),
    by = "Estadística"
  )

Si bien parece ser mejor el nuevo modelo por centésimas en el \(R^2\) ajustado, otras estadísticas como el logLik, IAC y BIC indican que nuestro modelo no mejoro y pero es relevante ver si mejora este modelo omitiendo outliers y valores influyentes.

Otras propuestas son las transformación Box-Cox que permiten un desplazamiento en los datos y la transformación inversa del seno hiperbólico.

\[ g_{\lambda_1, \lambda_2}(y) = \Bigg\{\begin{array}{rl}\frac{(y+\lambda_2)^{\lambda_1}-1}{\lambda_1} & \lambda_1 \neq 0\\ \log(y + \lambda_2) & \lambda_1 = 0\end{array} \]

\[ \begin{array}{lr} f(y, \theta) = \sinh^{-1}(\theta y)/\theta = \log\left(\theta y + (\theta^2y^2+1)^{1/2}\right)/\theta & ;\theta>0 \end{array} \]

Aquí se dejan otros enlaces interesantes sobre la transformación Box-Cox

Interpretaciones

Es importante hacer ver que al modificar alguna de las variables con alguna transformación, ya sea la variable respuesta o alguna de las variables independientes, también se modifica la interpretación de todo. Generalmente podemos aplicar transformaciones inversas para dar una interpretación de nuestras variables originales, hay que tener en cuenta ciertos detalles con algunas funciones.

  • Logaritmo
  • Raíz cuadrada

Validación

Para que nuestros resultados tengan validez, necesitamos cumplir una serie de supuestos y es recomendable verificar algunos problemas potenciales. De lo principal que debemos verificar en cada modelo es lo siguiente:

  1. Linearidad: La respuesta puede ser escrita como una combinación lineal de los predictores: \(\mathbb{E}[Y|X_1 = x_1, \dots, X_k = x_k] = \beta_0+\beta_1x_1+\cdots+\beta_kx_k\).

  2. Homocedasticidad: La varianza de los errores es la misma en cualquier conjunto de los valores predichos: \(\mathbb{V}ar(\epsilon_i) = \sigma^2\), con \(\sigma^2\) constante para \(i = 1, \dots, n\).

  3. Normalidad: La distribución de los errores deben seguir una distribución normal: \(\epsilon_i \sim \mathbb{N}(0,\sigma^2)\) para \(i = 1, \dots, n\).

  4. Independencia: Los errores son independientes: \(\epsilon_1, \dots, \epsilon_n\) son independientes o \(\mathbb{E}[\epsilon_i\epsilon_j] = 0, i\neq j\), es decir, que los errores no estan correlacionados ya que asumimos que son normales.

  5. Aditividad: El efecto que tiene una variable predictora sobre el modelo es independiente de otra regresora sobre la variable respuesta. En caso de tener problemas con esto, podríamos tener problemas de multicolinealidad.

Este último es considerado como supuesto del modelo en los libros [3] y [5] pero es común que sólo los primeros 4 sean considerados los supuestos del modelo. Si bien es cierto que es ideal que todos nuestros supuestos se cumplan, podemos dar cierta tolerancia a ellos por su Robustez. El libro [7] enlista los anteriores supuestos en orden de importancia (de manera descendente) y el libro [4] se dan algunos puntos importantes sobre esto.

  1. Validez: La información debe ser adecuada para responder tu pregunta objetivo

  2. Aditividad y linearidad: El supuesto matemático más importante del modelo de regresión es que su componente determinista es una función lineal de los predictores separados. Esto puede ser corregido con diferentes transformaciones.

Este supuesto puede ser violado cuando la recta de ajuste entre los datos no es recta (una curva por ejemplo) o cuando la cantidad y fuerza de los valores atípicos no permiten que el modelo sea adecuado. Ambas violaciones pueden hacer que las estimaciones de mínimos cuadrados den respuestas engañosas a las preguntas de interés. Las medias y las predicciones estimadas pueden estar sesgadas (subestiman o sobrestiman sistemáticamente la cantidad prevista) y las pruebas y los intervalos de confianza pueden reflejar de forma inexacta la incertidumbre

  1. Independencia de los errores: La falta de independencia no provoca sesgos en las estimaciones de mínimos cuadrados de los coeficientes, pero los errores estándar se ven seriamente afectados. Este problema puede ser dado, por ejemplo, cuando los individuos en un estudio pertenecen a un mismo grupo o comparten características similares que no se consideraron como que tenga la misma dieta, sean de la misma familia o hayan sido expuestos a las mismas condiciones ambientales.

  2. Homocedasticidad: Si la varianza de los errores de regresión es desigual, la estimación se realiza de manera más eficiente utilizando mínimos cuadrados ponderados, donde cada punto se pondera inversamente proporcional a su varianza. En la mayoría de los casos, sin embargo, este problema es menor. La varianza desigual no afecta el aspecto más importante de un modelo de regresión, que es la forma del predictor \(X\beta\).

Las consecuencias de violar este supuesto son las mismas que las del análisis de varianza unidireccional. Aunque las estimaciones de mínimos cuadrados siguen siendo insesgadas incluso si la varianza no es constante, los errores estándar describen de manera inexacta la incertidumbre en las estimaciones. Las pruebas y los intervalos de confianza pueden ser engañosos.

  1. Normalidad: El supuesto de regresión que generalmente es menos importante es que los errores se distribuyen normalmente. De hecho, para estimar la línea de regresión (en comparación con la predicción de puntos de datos individuales), el supuesto de normalidad apenas tiene importancia.

Las estimaciones de los coeficientes y sus errores estándar son robustas a distribuciones no normales. Aunque las pruebas y los intervalos de confianza se originan en distribuciones normales, las consecuencias de violar este supuesto suelen ser menores. La única situación de gran preocupación es cuando las distribuciones tienen colas largas (existen valores atípicos) y los tamaños de muestra son de moderados a pequeños. Si se utilizan intervalos de predicción, por otro lado, las desviaciones de la normalidad se vuelven importantes. Esto se debe a que los intervalos de predicción se basan directamente en la normalidad de las distribuciones de la población, mientras que las pruebas y los intervalos de confianza se basan en la normalidad de las distribuciones muestrales de las estimaciones (que pueden ser aproximadamente normales incluso cuando las distribuciones de la población no lo son).

Por lo que en este punto, si tenemos una buena cantidad de datos, podemos ser flexibles con este último supuesto.

La manera más sencilla de ver la mayoría de estos supuestos es mediante una gráfica de residuales vs los valores ajustados donde, para la linealidad deseamos ver que los residuales tengan un comportamiento aleatorio con una varianza constante y media cero, lo cual implicitamente también ayuda con la homocedasticidad y normalidad. Esta gráfica y otras no las proporciona el comando plot() para un modelo lineal

par(mfrow = c(2, 2))
third_model_BX %>% plot()

O bien con la función ggfortify::autoplot()

library(ggfortify)
third_model_BX %>% autoplot()

Para más detalles de estas gráficas véase la sección de Análisis de valores influyentes. Finalmente aplicamos diversas pruebas para probar cada uno de los supuestos:

  • Normalidad:
    • Todas las pruebas que se vieron anteriormente.
third_model_BX %>% residuals() %>% shapiro.test()

    Shapiro-Wilk normality test

data:  .
W = 0.98185, p-value = 6.856e-16
  • Homocedasticidad:
    • Prueba de Breusch-Pagan: lmtest::bptest()
    • Prueba de White: lmtest::bptest()
library(lmtest)
third_model_BX %>% bptest()

    studentized Breusch-Pagan test

data:  .
BP = 6.2836, df = 4, p-value = 0.1789
  • Autocorrelación (independencia):
    • Durbin Watson: lmtest::dwtest()
    • Breusch-Godfrey Test: lmtest::bgtest()
dwtest(third_model_BX)

    Durbin-Watson test

data:  third_model_BX
DW = 1.8388, p-value = 9.232e-05
alternative hypothesis: true autocorrelation is greater than 0
bgtest(third_model_BX)

    Breusch-Godfrey test for serial correlation of order up to 1

data:  third_model_BX
LM test = 13.716, df = 1, p-value = 0.0002126

Una manera sencilla para detecetar la colinearidad es ver un correlograma y tomar aquellas variables que tengan fuerte relación entre ellas. Así sabiendo que una variable explica a otra, podríamos omitirla del modelo y determinar si esto mejora nuestro rendimiento, o bien crear una nueva variable. Otra solución en caso de colinearidad podría ser una interacción como ya se vio previamente.

El problema con lo anterior es que no todos los problemas de colinearidad pueden ser detectados mediante un correlograma, ya que puede existe una relación entre tres o más variables sin que sea visible tomándolo por pares. A lo anterior es lo que llamaremos multicolinearidad

Lo recomendable, y más sencillo, es visualizar el factor de inflación de varianza: VIF el cual es la razón de la varianza de \(\beta_j\) cuando

\[ VIF\left(\hat{\beta}_i\right) = \frac{1}{1-R^2_{X_j|X_{-j}}} \]

Donde \(R^2_{X_j|X_{-j}}\) es la \(R^2\) de la regresión de \(X_j\) a todos los demás predictores; es decir, la proporción de variabilidad en \(X_j\) que es explicada por su relación lineal a las otras variables del modelo. Cuando \(X_j\) puede ser bien explicada por la adicción de las demas variables; es decir cuando \(R^2_{X_j|X_{-j}}\) es cercano a 1, entonces el valor \(VIF\) será grande e indicaría problas de colinearidad.

Buscamos valores pequeños y lo más pequeño que se puede obtener para cada variable es un \(VIF = 1\) lo cual indica una completa ausencia de colinearidad. Por regla de dedo, un \(VIF\) mayor a 5 o 10 indica problemas de colinearidad. En R podemos obtener el \(VIF\) para cada variable con la función faraway::vif()

library(faraway)
third_model_BX %>% vif()
log(transporte) sqrt(alimentos)   log(limpieza) log(personales) 
       1.409547        1.469732        1.560714        1.649865

No tenemos problemas de colinearidad.

¿Qué sucedió con nuestro segundo modelo?

second_m_model %>% residuals() %>% shapiro.test()

    Shapiro-Wilk normality test

data:  .
W = 0.98339, p-value = 9.574e-16
second_m_model %>% bptest()

    studentized Breusch-Pagan test

data:  .
BP = 51.355, df = 7, p-value = 7.823e-09
second_m_model %>% bgtest()

    Breusch-Godfrey test for serial correlation of order up to 1

data:  .
LM test = 14.285, df = 1, p-value = 0.0001571
second_m_model %>% vif()
           sqrt(transporte)             sqrt(alimentos) 
                   1.475002                    1.546820 
             sqrt(limpieza)            sqrt(personales) 
                   1.631803                    1.816380 
              sqrt(calzado)               sqrt(vestido) 
                   2.500522                    2.971420 
sqrt(calzado):sqrt(vestido) 
                   4.387703

Invalidamos el supuesto de homocedasticidad que ya habíamos conseguido con nuestro tercer modelo, por lo que preferiremos el tercer modelo desde ahora.

Predicción

Uno de los puntos más objetivos principales de crear un modelo de regresión es predecir un valor futuro de acuerdo a diferentes valores de nuestras variables predictoras. Esto es sencillo teniendo ya el modelo. Vamos a suponer que nuestro último modelo es adecuado (ya que aún no hemos tratado el tema de validación de supuestos, los cuales en predicción son vitales para su credibilidad)

En R podemos utilizar el comando predict() para obtener fácilmente las predicciones de nueva información sobre un modelo. En esta caso vamos a realizar un ejemplo con regresión simple ya que es análogo para el modelo de regresión múltiple.

new_data <-  data.frame(vestido = 100)
simple_model <- lm(log(ing_cor)~sqrt(vestido), data = mutated_data_income %>% filter(vestido<7500 & vestido>0))
simple_model %>% predict(new_data)
       1 
10.42712

Véase que R se encarga de las transformaciones en las variables predictoras (en este caso el modelo queda expresado como \(y = 10.26570 + 0.01614\times \sqrt{X_{vestido}}\))

  • \(10.26570 + 0.01614\times100 = 11.8797\)
  • \(10.26570 + 0.01614\times\sqrt{100} = 10.42712\)

Entonces, en este caso podemos decir que con un gasto de $100 mexicanos en vestido (con esta base de datos), un ciudadano de Aguascalientes, en promedio tiene que ganar \(e^{10.26570}\times e^{(0.01614\times\sqrt(100)} = exp(10.42712) = \$33762.97\) pesos mexicanos.

La misma función predict() nos permite obtener los intervalos de confianza para la respuesta media y los intervalos de predicción con el argumento interval = c("confidence", "prediction")

¿Cuál es la diferencia entre estos?

Como sabemos, nuestro modelo sólo es una representación de la realidad ya que mediante algún método obtenemos estimaciones sobre los parámetro (los \(\beta_i\)s); es decir, sólo estamos obteniendo una estimación del verdadero plano poblacional de regresión

\[ f(X) = \beta_0 + \beta_1X_1 + \cdots + \beta_{p}X_p \]

con nuestro plano obtenido por el método de mínimos cuadrados

\[ \hat{Y} = \hat{\beta}_0 + \hat{\beta_1}X_1 + \cdots + \hat{\beta}_{p}X_p \]

en el cual controlamos un error que se puede reducir. Sobre dicha respuesta promedio es donde tenemos nuestros intervalos de confianza. Por ejemplo, si nosotros recolectamos distintas bases de datos con las variables que tenemos en nuestro último modelo del ingreso en Aguascalientes, digamos 100, y obtenemos el intervalo de confianza sobre el promedio del gasto en cada una de esas bases (con las variables explicativas fijas), 95% de dichos intervalos contendrán el verdadero valor promedio del gasto (la respuesta media en este caso). Recordemos que este valor promedio es lo que tenemos con la evaluación de un vector X conocido.

En el caso del intervalo de predicción, estamos haciendo inferencia sobre el error que no controlamos en el modelo (\(\epsilon\)), por lo que es común que el intervalo de predicción sea más amplio que los intervalos de confianza al incorporar esta incertidumbre. Si bien, al igual que en el intervalo de confianza estaríamos dando un intervalo para un valor de la respuesta, en este caso es de manera específica y no sobre una aproximación hacia la población.

Entonces, sí se desea dar un intervalo para un punto en particular, utilizaremos los intervalos de predicción y si deseamos hablar de los valores medios de predicción, usaríamos un intervalo de confianza.

new_data <-  data.frame(vestido = 100, 200, 150)
simple_model %>% predict(new_data, interval = "confidence")
       fit      lwr      upr
1 10.42712 10.37179 10.48246
simple_model %>% predict(new_data, interval = "prediction")
       fit      lwr      upr
1 10.42712 9.206192 11.64806

En este caso de una regresión lineal simple podemos visualizar mejor dichos intervalos

pred_int <- predict(simple_model, interval = "prediction")

mutated_data_income %>% filter(vestido<7500 & vestido>0) %>% bind_cols(pred_int %>% as_tibble()) %>% 
  ggplot(aes(x = sqrt(vestido), y = log(ing_cor))) + 
  geom_point() + 
  geom_smooth(method = "lm") + #Intervalos de confianza
  geom_line(aes(y = lwr), color = "red", linetype = "dashed") + #intervalos de predicción
  geom_line(aes(y = upr), color = "red", linetype = "dashed") + 
  general_theme

Otros puntos importantes de regresión lineal

  • Si se específica el tipo de regresión, ¿Hay otras técnicas de modelaje?

¡Claro! existen muchísimas técnicas de modelaje, en cuanto a la regresión tenemos otros tipos además de la lineal, como la regresión logística, poisson, polinómica, splines, etc.

  • ¿El modelo es lineal en que sentido?

Es lineal en sus coeficientes, es decir que los siguientes modelos también son lineales:

\[ \begin{array}{c} \mu\left\{Y|X_1, X_2\right\} = \beta_0+\beta_1X_1+\beta_2X_2\\ \mu\left\{Y|X_1\right\} = \beta_0+\beta_1X_1+\beta_2X_1^2\\ \mu\left\{Y|X_1, X_2\right\} = \beta_0+\beta_1X_1+\beta_2X_1^2+\beta_3X_1X_2\\ \mu\left\{Y|X_1, X_2\right\} = \beta_0+\beta_1\log(X_1)+\beta_2\log(X_2) \end{array} \]

  • No hay que olvidar toda la teoría que esta por detrás de tal modelo, como la utilización del procedimiento de mínimos cuadrados para estimar los parámetros del modelo.

  • ¿Es suficiente realizar una regresión con todos los datos “limpios?”

En algunos casos no. Piensen en la paradoja de Simpson, la cual menciona que una tendencia que aparece en varios grupos de datos puede desaparece o invertirse cuando se combinan los grupos. Bajo este enfoque, siempre es bueno segregar la información si se considera que dicha segregación mejorará los resultados del modelo.

  • ¿Hay otros métodos para crear un modelo de regresión sin utilizar el método de mínimos cuadrados?

SÍ.

Análisis de valores influyentes

Outliers

Como ya se ha mencionado, ay algunos datos que nos pueden generar problemas, y por lo mismo previamente se han eliminado algunos datos, como ciertos outliers y ceros en nuestras variables.

Un outlier en nuestro modelo de regresión será aquel que sea lejano a su valor predicho. Para ver esto gráficamente vamos considerar nuestro ejemplo de regresión lineal simple visto en la sección de predicción (\(\log(ing\_cor) = 10.265703 + 0.016142\times \sqrt{vestido}\))

library(cowplot)
m1 <- mutated_data_income %>% filter(vestido<7500 & vestido>0) %>% 
  ggplot(aes(x = sqrt(vestido) , y = log(ing_cor))) +
  geom_point() + #Hasta aquí los datos
  geom_abline(intercept = 10.265703, slope = 0.016142, color = "red3") + 
  geom_point(data = tibble(x = 35.440514, y = 7.535708), aes(x = x, y = y), color = "red") + 
  lims(y = c(7, 14)) +
  general_theme 
m2 <- mutated_data_income %>% filter(vestido<7500 & vestido>0) %>% 
  filter(vestido!=1256.03) %>% 
  ggplot(aes(x = sqrt(vestido) , y = log(ing_cor))) +
  labs(y = "") + 
  lims(y = c(7, 14)) + 
  geom_point() + #Hasta aquí los datos
  geom_abline(intercept = 10.267257, slope = 0.016167, color = "red3") + 
  general_theme
plot_grid(m1, m2, ncol = 2, align = "h")

Le punto rojo es nuestro outlier en nuestro modelo inicial y en la misma subgráfica se aprecia la linea de regresión que se ajusta a tales datos. La segunda gráfica nos muestra la linea de regresión que resultad al omitir dicha observación. Como bien se ve, no se afecta realmente el comportamiento de la linea de regresión, pero veamos como cambiaron en otro aspecto los modelos

data_vestido <- mutated_data_income %>% filter(vestido<7500 & vestido>0)
m1 <- data_vestido %>% lm(log(ing_cor)~sqrt(vestido), data = .) 
m2 <- data_vestido %>% 
  filter(vestido!=1256.03) %>% lm(log(ing_cor)~sqrt(vestido), data = .)
l1 <- m1 %>% summary()
l2 <- m2 %>% summary()
cat(paste0("RSE del Modelo con outlier: ",l1$r.squared, "\n"))
RSE del Modelo con outlier: 0.154956818679954
cat(paste0("RSE del Modelo sin outlier: ",l2$r.squared))
RSE del Modelo sin outlier: 0.15805697613396

Si bien el impacto en este caso no se ve tan drástico, este lo puede ser en otros modelos y hay que recordar que el RSE es utilizado para calcular los intervalos de confianza y los \(p-values\), por lo que un cambio tan fuerte para un sólo punto puede tener graves problemas de interpretación.

Para identificar outliers, podemos visualizar dos gráficas: residuales vs valores ajustados y residuales estudentizados vs valores ajustados

res_1 <- (tibble(x = m1$fitted.values, y = residuals(m1)) %>% 
    ggplot(aes(x = x, y = y)) + geom_point(alpha = 0.5) + 
    general_theme + labs(x = "Valores ajustados", y = "Residuales")) 
res_2 <- (tibble(x = m1$fitted.values, y = rstudent(m1)) %>% 
   ggplot(aes(x = x, y = y)) + geom_point(alpha = 0.5) + 
     geom_hline(yintercept = c(-3, 3), alpha = 0.5)+ 
   general_theme + labs(x = "Valores ajustados", y = "Residuales Estudentizados"))
res_1 + res_2

Aveces es difícil determinar si una observación debe ser considerado como outlier con la primera gráfica, pero en la segunda es más sencillo. De acuerdo al libro [3], observaciones con residuales estudentizados que son más grandes que 3, en valor absoluto, pueden ser considerados como posibles outliers. Lo cual tiene sentido, ya que al “normalizar” los residuales así, se esta haciendo que estos sigan una distribución \(t-student\) con \(n-k-1\) grados de libertad, en nuestro caso \(df = 1398-1-1 = 1396\)

library(latex2exp)
tibble(x= c(-4, 4)) %>% 
  ggplot(aes(x = x)) + stat_function(fun = ~dt(.x, 1396)) +
  labs(y = "Densidad") + 
  ggtitle("t-student con 1396 grados de libertad") + 
  general_theme

Entonces, de acuerdo a lo anterior tenemos, en este caso, 9 outliers.

tibble(x = data_vestido$vestido , y = data_vestido$ing_cor, 
       fitted_values = m1$fitted.values, rStudent = rstudent(m1)) %>% 
  mutate(index = row_number()) %>% filter(rStudent< -3 | rStudent>3)

Veamos que sucede con nuestro último modelo

filter_data_ToOutliers <- filter_data %>% 
  mutate(fitted_values = third_model_BX$fitted.values, 
         rStandar = residuals(third_model_BX),
         rStudent = rstudent(third_model_BX),
         fitted_values = third_model_BX$fitted.values)
r_third_model1 <- filter_data_ToOutliers %>% 
  ggplot(aes(y = rStandar, x = fitted_values)) + 
  labs(x = "Valores ajustados", y = "Residuales") + 
  geom_point(alpha = 0.5) + general_theme

r_third_model2 <-filter_data_ToOutliers %>% 
  ggplot(aes(y = rStudent, x = fitted_values)) + 
  geom_hline(yintercept = c(-3, 3), alpha = 0.5) + 
  labs(x = "Valores ajustados", y = "Residuales estudentizados") +
  geom_point(alpha = 0.5) + general_theme

r_third_model1 + r_third_model2

Es decir, tenemos 20 outliers

filter_data_ToOutliers %>% filter(rStudent< -3 | rStudent>3)

Si dicho outlier fue generado por un error de recolección de información, podemos eliminar la observación sin problemas, aunque hay que considerar que si no es el caso, entonces nuestro modelo simplemente no es lo suficientemente bueno para predecir valores que pueden suceder. En este caso vamos a eliminarlos considerando un limite en 3.5.

filter_data_ToOutliers %>% 
  ggplot(aes(y = rStudent, x = fitted_values)) + 
  geom_hline(yintercept = c(-3.5, 3.5), alpha = 0.5) + 
  labs(x = "Valores ajustados", y = "Residuales estudentizados") +
  geom_point(alpha = 0.5, color = "#447e9e") + 
  geom_text( 
    data = filter_data_ToOutliers %>% 
      mutate(index = row_number()) %>% 
      dplyr::filter(rStudent< -3.5 | rStudent>3.5),
    aes(label=index), nudge_x = 0.2
  )+ general_theme

fourth_model_1 <- filter_data_ToOutliers %>% filter(rStudent>=-3.5 & rStudent<=3.5) %>%
  lm(((ing_cor^lambdaBCox_tM-1)/lambdaBCox_tM) ~ 
       log(transporte) + sqrt(alimentos) + log(limpieza) + log(personales), 
     data = .)
fourth_model_1 %>% summary()

Call:
lm(formula = ((ing_cor^lambdaBCox_tM - 1)/lambdaBCox_tM) ~ log(transporte) + 
    sqrt(alimentos) + log(limpieza) + log(personales), data = .)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.1395 -0.4038 -0.0162  0.3866  2.2298 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)     6.8496905  0.1292413  52.999   <2e-16 ***
log(transporte) 0.2514986  0.0134651  18.678   <2e-16 ***
sqrt(alimentos) 0.0057975  0.0004841  11.975   <2e-16 ***
log(limpieza)   0.1593438  0.0171812   9.274   <2e-16 ***
log(personales) 0.1849702  0.0177541  10.418   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6049 on 2127 degrees of freedom
Multiple R-squared:  0.5193,    Adjusted R-squared:  0.5184 
F-statistic: 574.5 on 4 and 2127 DF,  p-value: < 2.2e-16
fourth_model_1 %>% autoplot() + general_theme

Comparison_models <- third_model_BX %>% glance() %>%  
  gather("Estadística", "Tercel Modelo") %>%
  dplyr::mutate_if(is.numeric, round, digits = 8) %>% 
  inner_join(
    fourth_model_1 %>% glance() %>%  gather("Estadística", "Cuarto Modelo_1") %>%
      dplyr::mutate_if(is.numeric, round, digits = 8),
    by = "Estadística"
  )
Comparison_models
tests_linearM <- function(model){
  shapiro <- model %>% residuals() %>% shapiro.test()
  breusch_pagan <- model %>% bptest()
  durbin_watson <- model %>% dwtest()
  breusch_godfrey <- model %>% bgtest()
  p_values <- tibble("Estadística" = c("Shapiro-Wils", "Breusch-Pagan",
                                       "Durbin-Watson", "Breusch_Godfrey"), 
                     "P-value" = c(shapiro$p.value, breusch_pagan$p.value,
                                   durbin_watson$p.value, breusch_godfrey$p.value)) %>% 
    dplyr::mutate_if(is.numeric, round, digits = 8)
  return(list("p-values" = p_values, "vif" = model %>% vif()))
}
fourth_model_1 %>% tests_linearM()
$`p-values`
# A tibble: 4 × 2
  Estadística      `P-value`
  <chr>                <dbl>
1 Shapiro-Wils    0.0152    
2 Breusch-Pagan   0.137     
3 Durbin-Watson   0.00000008
4 Breusch_Godfrey 0.00000019

$vif
log(transporte) sqrt(alimentos)   log(limpieza) log(personales) 
       1.412725        1.476318        1.559839        1.655387

Véase que seguimos “manteniendo” el supuesto de homocedasticidad y obtuvimos un mejor p-value para la prueba de normalidad aunque no mejoramos nuestro modelo para las pruebas de independencia, además de que no tenemos problemas de multicolinearidad.

¿Que hubiera sucedido si hubiéramos aplicado estos filtros a los datos del segundo modelo? ¿Hubiéramos conseguido un mejor modelo?

Valores con alto “apalancamiento”

Mientras que en los outliers nos enfocamos a los valores inusuales de la variable respuesta para un \(x_i\), una observación con alto apalancamiento (leverage) tiene un valor inusual para \(x_i\). El leverage de una observación es una medida de la distancia entre los valores de sus variables explicativas y el promedio de los valores de las variables explicativas en todo el conjunto de datos. Observaciones con alto apalancamiento pueden tener una alta influencia en el modelo.

Cuando sólo tenemos una variable explicativa en nuestro modelo, la estadística de apalancamiento queda definida de la siguiente manera:

\[ h_i = \frac{1}{(n-1)}\left[\frac{x_i-\bar{x}}{s_x}\right]^2 + \frac{1}{n} = \frac{(x_i-\bar{x})^2}{\sum_{i^{'} = 1}^n (x^{'}-\bar{x})^2} + \frac{1}{n} \]

En el caso de regresión lineal simple es sencillo detectar este tipo de observaciones, ya que su valor predicho estará fuera del rango normal de las demás predicciones.

Con la ecuación anterior podemos ver que \(h_i\) incrementa con la distancia entre \(x_i\) y \(\bar{x}\). Esta estadística siempre estará entre \(1/n\) y 1 y el leverage promedio para todas las observaciones será igual a \((p+1)/n\) con \(p\) el número de coeficientes, por lo que si una observación tiene un leverage más grande que \((p+1)/n\) podríamos decir que se tiene un alto apalancamiento aunque algunos estadísticos prefieren utilizar como umbral \(2*p/n\).

Ya habíamos retirado algunos outliers, vamos a ver si estas observaciones también son tienen un alto leverage. En R, podemos utilizar la función stats::hatvalues()

filter_data_ToLeverage <- filter_data_ToOutliers %>% 
  mutate(leverage = third_model_BX %>% hatvalues(),index = row_number())

(Leverage_analysis <- filter_data_ToLeverage %>% 
  ggplot(aes(x = index, y = leverage)) + 
  geom_point(alpha = 0.5, color = "#447e9e") + 
  geom_hline(yintercept = (5+1)/third_model_BX$fitted.values %>% length(), 
             color = "blue3", alpha = 0.5) + 
  geom_hline(yintercept = 2*5/third_model_BX$fitted.values %>% length(), 
             color = "red3", alpha = 0.5) +
  geom_hline(yintercept = 0.016, color = "black", alpha = 0.5)+
  geom_text( 
    data = filter_data_ToLeverage %>% dplyr::filter(leverage >=0.016), 
    aes(label=index), nudge_x = 100
  )+ general_theme)

Las lineas de colores tan sólo son distintos umbrales o cortes para establecer que observación tiene un alto apalancamiento. En este caso no tenemos que algún outlier también tiene un alto apalancamiento, lo cual es lo ideal.

fourth_model_2 <- filter_data_ToLeverage %>% 
  filter(rStudent>=-3.5 & rStudent<=3.5) %>%
  filter(leverage< 0.016) %>% 
  lm(((ing_cor^lambdaBCox_tM-1)/lambdaBCox_tM) ~ 
       log(transporte) + sqrt(alimentos) + log(limpieza) + log(personales), 
     data = .)
fourth_model_2 %>% summary()

Call:
lm(formula = ((ing_cor^lambdaBCox_tM - 1)/lambdaBCox_tM) ~ log(transporte) + 
    sqrt(alimentos) + log(limpieza) + log(personales), data = .)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.14019 -0.40307 -0.01616  0.38674  2.22978 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)     6.8470768  0.1293183  52.947   <2e-16 ***
log(transporte) 0.2510460  0.0134731  18.633   <2e-16 ***
sqrt(alimentos) 0.0058011  0.0004842  11.981   <2e-16 ***
log(limpieza)   0.1602029  0.0174279   9.192   <2e-16 ***
log(personales) 0.1849316  0.0179337  10.312   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6049 on 2125 degrees of freedom
Multiple R-squared:  0.5196,    Adjusted R-squared:  0.5187 
F-statistic: 574.6 on 4 and 2125 DF,  p-value: < 2.2e-16
fourth_model_2 %>% autoplot() + general_theme

Comparison_models <- Comparison_models %>% 
  inner_join(
    fourth_model_2 %>% glance() %>%  gather("Estadística", "Cuarto Modelo_2") %>%
      dplyr::mutate_if(is.numeric, round, digits = 8),
    by = "Estadística"
  )
Comparison_models
fourth_model_2 %>% tests_linearM()
$`p-values`
# A tibble: 4 × 2
  Estadística      `P-value`
  <chr>                <dbl>
1 Shapiro-Wils    0.0148    
2 Breusch-Pagan   0.136     
3 Durbin-Watson   0.00000007
4 Breusch_Godfrey 0.00000017

$vif
log(transporte) sqrt(alimentos)   log(limpieza) log(personales) 
       1.413477        1.476237        1.583429        1.680891

Véase que al retirar estas observaciones el modelo no mejoró aunque no empeoro drásticamente.

Otra estadística para determinar la influencia que tiene una observación en el modelo, es la estadística de Cook, la cual mide la influencia general de la observación, es decir que muestra el efecto que tiene omitir tal observación en el modelo de regresión.

\[ D_i = \sum_{j = 1}^n\frac{(\hat{y}_{j(i)}-\hat{y}_j)^2}{p\hat{\sigma}^2} \]

donde \(\hat{y}_j\) es el \(j-ésimo\) valor en el ajuste usando todos las observaciones y \(\hat{y}_{j(i)}\) es el valor de la regresión omitiendo la observación \(i\); \(p\) es el número de coeficientes (se incluye \(\beta_0\)) y \(\sigma\) la varianza estimada desde el modelo. El calculo de esta estadística es costoso ya que se tienen que hacer diferentes modelos de regresión, por lo que se prefiere utilizar una equivalencia

\[ D_i = \frac{1}{p}(r_i)^2\left(\frac{h_i}{1-h_i}\right) \]

donde \(r_i\) es el residual estandarizado aunque algunos autores prefieren los residuales estudentizados.

El umbral que se establezca también dependerá de la información aunque algunos autores sugieren valores cercanos o más grandes que 1 son valores con gran influencia

filter_data_ToCook <- filter_data_ToLeverage %>% 
  mutate(Cook = third_model_BX %>% cooks.distance())

(Cook_analysis <- filter_data_ToCook %>% 
  ggplot(aes(x = index, y = Cook)) + 
  geom_point(alpha = 0.5, color = "#447e9e") + 
  geom_hline(yintercept = 0.03, color = "black", alpha = 0.5)+
  geom_text( 
    data = filter_data_ToCook %>% dplyr::filter(Cook >= 0.03), 
    aes(label=index), nudge_x = 100
  )+ general_theme)

Es preferible estudiar los valores influyentes y outliers en conjunto, para determinar aquellos valores que tengan más problemas en conjunto

Finalmente, vamos a retirar aquellas observaciones que causan mayores problemas y veamos como se comporta nuestro modelo

fourth_model_3 <- filter_data_ToCook %>% 
  filter(!(index %in% c(881, 190, 1858, 2023, 703, 515))) %>%
  lm(((ing_cor^lambdaBCox_tM-1)/lambdaBCox_tM) ~ 
       log(transporte) + sqrt(alimentos) + log(limpieza) + log(personales), 
     data = .)
fourth_model_3 %>% summary()

Call:
lm(formula = ((ing_cor^lambdaBCox_tM - 1)/lambdaBCox_tM) ~ log(transporte) + 
    sqrt(alimentos) + log(limpieza) + log(personales), data = .)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2355 -0.4063 -0.0165  0.3863  3.3151 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)     6.8390940  0.1334738  51.239   <2e-16 ***
log(transporte) 0.2571418  0.0139177  18.476   <2e-16 ***
sqrt(alimentos) 0.0058523  0.0005006  11.690   <2e-16 ***
log(limpieza)   0.1543125  0.0178434   8.648   <2e-16 ***
log(personales) 0.1843599  0.0184677   9.983   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6263 on 2133 degrees of freedom
Multiple R-squared:  0.5043,    Adjusted R-squared:  0.5034 
F-statistic: 542.5 on 4 and 2133 DF,  p-value: < 2.2e-16
fourth_model_3 %>% autoplot() + general_theme

Comparison_models <- Comparison_models %>% 
  inner_join(
    fourth_model_3 %>% glance() %>%  gather("Estadística", "Cuarto Modelo_3") %>%
      dplyr::mutate_if(is.numeric, round, digits = 8),
    by = "Estadística"
  )
Comparison_models
fourth_model_3 %>% tests_linearM()
$`p-values`
# A tibble: 4 × 2
  Estadística     `P-value`
  <chr>               <dbl>
1 Shapiro-Wils    0        
2 Breusch-Pagan   0.513    
3 Durbin-Watson   0.0000346
4 Breusch_Godfrey 0.0000818

$vif
log(transporte) sqrt(alimentos)   log(limpieza) log(personales) 
       1.412271        1.474882        1.573000        1.670096

Con nuestro último modelo ya omitimos valores influyentes, retiramos outliers, tampoco estamos considerando algunos ceros y aplicamos diversas transformaciones, tanto para la variable dependiente y las variables independientes. Con este modelo logramos un modelo con homocedasticidad, al estar considerando una cantidad importante de datos podemos omitir el supuesto de normalidad aunque ay que tener precaución con las predicciones que se lleguen a dar. Otro punto a considerar es que el supuesto de independencia no se logro mejorar y para solucionar dicho problema tendríamos que considerar otro tipo de modelo.

Aquí se dejan más enlaces:

ANOVA vs Regresión lineal

aov(ing_cor_log~transporte_sqrt + alimentos_sqrt +limpieza_sqrt + personales_sqrt, data = mutated_data_income) %>% summary()
anova(aov(ing_cor_log~transporte_sqrt + alimentos_sqrt +limpieza_sqrt + personales_sqrt, data = mutated_data_income))

La tabla ANOVA no sólo se puede utilizar en el ámbito de la regresión lineal. Si no que se puede utilizar como una prueba paramétrica (por lo que tiene una gran potencia en comparación de las pruebas no paramétricas) para determinar diferencias entre poblaciones.

Aquí se deja un conjunto de enlaces:

Librerías para EDA

[1] Stephanie Glen. “Hypothesis Testing” From StatisticsHowTo.com: Elementary Statistics for the rest of us! https://www.statisticshowto.com/probability-and-statistics/hypothesis-testing/

[2] Sumpter, David. Fútbol y matemáticas. Ariel, 2016.

[3] James, G., Witten, D., Hastie, T., & Tibshirani, R. (2013). An introduction to statistical learning (Vol. 112, p. 18). New York: springer.

[4] Ramsey, F., & Schafer, D. (2012). The statistical sleuth: a course in methods of data analysis. Cengage Learning.

[5] Friedman, J., Hastie, T., & Tibshirani, R. (2001). The elements of statistical learning (Vol. 1, No. 10). New York: Springer series in statistics.

[6] Susan, Milton, Jesse C. Arnold, and Jorge Luis Blanco y Correa Magallanes. Probabilidad y estadística: con aplicaciones para ingeniería y ciencias computacionales. 2004.

[7] Gelman, A., & Hill, J. (2006). Data analysis using regression and multilevel/hierarchical models. Cambridge university press.

[8] Taddy, M. (2019). Business data science: Combining machine learning and economics to optimize, automate, and accelerate business decisions. McGraw Hill Professional.

 

A work by Carlos Vásquez

carlosfvasquez@ciencias.unam.mx