Parcial 2 - Modelo Lineal

Librerías

library(tidyverse)
library(magrittr) # %<>%
library(broom) # augment, glance, tidy
library(GGally) # ggpairs
library(gghighlight) # gghighlight
library(ggrepel) # geom_label_repel

Punto 1

datos1 <- read.table("datos1.txt") %>% as_tibble
glimpse(datos1)

## Observations: 100
## Variables: 5
## $ y  <dbl> 412.8432, 129.4738, 368.4621, 190.7698, 161.2475, 209.6363, 1…
## $ x1 <dbl> 19.71750, 10.83758, 18.73870, 13.29231, 12.22276, 14.01648, 1…
## $ x2 <dbl> 19.23959, 11.01319, 19.53059, 12.03363, 13.91585, 14.87797, 1…
## $ x3 <dbl> -1.10617180, 0.14439070, 2.66525949, -0.04207274, 3.39226155,…
## $ x4 <dbl> 1.4862349, 0.1340185, 0.3929311, 2.4853584, 0.2753532, 0.1549…

ggpairs(datos1, progress = F)

X1 parece estar muy correlacionada con Y, y también X2, pero ya se observa un problema: las covariable X1 y X2 estan muy correlacionadas entre sí! Por otro lado, X3 y X4 parecen no decir nada sobre Y.

Punto a

modelo <- lm(y ~ x1 + x2 + x3 + x4, data = datos1)
summary(modelo)

## 
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4, data = datos1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.597  -6.584  -1.817   5.968  19.321 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -210.0257     4.1889 -50.138  < 2e-16 ***
## x1            28.8426     0.9018  31.982  < 2e-16 ***
## x2             2.0059     0.8560   2.343   0.0212 *  
## x3            -2.5914     0.3830  -6.767 1.08e-09 ***
## x4            -0.3958     0.3907  -1.013   0.3136    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.959 on 95 degrees of freedom
## Multiple R-squared:  0.9931, Adjusted R-squared:  0.9928 
## F-statistic:  3403 on 4 and 95 DF,  p-value: < 2.2e-16

La regresión fue significativa!

Punto b

Graficar \(r_{i}\) (residuos estandarizados) vs. valores predichos. Qué observa?

datos1_ajuste <- augment(modelo, data = datos1)

ggplot(data = datos1_ajuste) +
  aes(x = .fitted, y = .std.resid) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dotted") +
  geom_smooth(se = F, color = "Gray") +
  labs(
    title = "Punto a",
    subtitle = "Residuos est. vs. valores predichos",
    y = "Residuos estandarizados",
    x = "Valores predichos (y_hat)"
  )

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

¿Observa alguna estructura o tendencia?

Sí! Los residuos tienen una estructura muy definida, que parece cuadrática con respecto a los valores predichos.

¿Cómo sería el gráfico ideal?

El gráfico ideal sería de residuos repartidos alrededor del cero, que se dispersan en igual modo (tanto por encima como por debajo del cero) a medida que se avanza en el eje x (de valores predichos). No se parece en nada al que obtuvimos!

¿Cuáles serían las razones para un gráfico como el que obtuvo?

Podría haber colinearidad entre las covariables (lo que sospechábamos antes) o haber otros factores que explican las Y y que no están en nuestro conjunto de datos (es decir, no los tenemos como covariables). También podría ocurrir que la esperanza de las Y no sea lineal, o que los errores no se distribuyan normalmente.

Punto c

Realizar un QQ-plot de los residuos estandarizados.

ggplot(datos1_ajuste) +
  aes(sample = .std.resid) +
  geom_qq() +
  geom_qq_line() +
  labs(
    title = "Punto b: testear normalidad",
    subtitle = "QQ-plot de los residuos estandarizados",
    y = "Cuantiles observados de los residuos estandarizados",
    x = "Cuantiles teóricos de dist. normal"
  )

Claramente se ve que los residuos (y por ende, sospechamos, los errores del modelo) no siguen una distribución normal!

¿Que test conoce para testear la normalidad de los residuos? Aplicarlo.

Podemos aplicar el test de Shapiro-Wilk, que confirmaría con un p valor bajo nuestra sospecha de no normalidad:

shapiro.test(datos1_ajuste$.std.resid)

## 
##  Shapiro-Wilk normality test
## 
## data:  datos1_ajuste$.std.resid
## W = 0.95138, p-value = 0.001024

Efectivamente, el p valor (0.001024) es significativo, nos permite rechazar la hipótesis nula de normalidad de los residuos. Concluimos que los residuos no se distribuyen normalmente.

Punto d

Realizar un analisis de Box-Cox de los datos, explicar el objetivo y cómo se interpreta el gráfico y el intervalo que ne él se halla.

lambdas <- MASS::boxcox(modelo)

En primer lugar, para transformar las variable Y como sugiere el método de Box-Cox, necesitamos encontrar el lambda cuya log-verosimiltud para nuestros datos sea máxima. El gráfico de arriba, realizado por MASS::boxcox, nos permite ubicar un lambda óptimo dentro de ese rango dibujado, que es un intervalo de confianza para el valor de lambda.

¿A qué conclusiones llega? ¿Aplicaría alguna transformación? Aplicarla y repetir a, b, c.

Tenemos que transformar los valores de Y según la fórmula (dado que el lambda encontrado no es cero):

\[ Y = (Y^{\lambda} - 1) / \lambda \quad \lambda \neq 0 \]

# Acá me quedo con los 3 lambdas que maximizan la log-likelihood según MASS:

as_tibble(lambdas) %>%
  rename(lambda = x, log_likelihood = y) %>%
  top_n(3, log_likelihood) %>%
  arrange(-log_likelihood)

## # A tibble: 3 x 2
##   lambda log_likelihood
##    <dbl>          <dbl>
## 1  0.545           272.
## 2  0.505           272.
## 3  0.586           261.

El lambda más verosímil es 0.5454545. Lo uso para transformar las Y:

lambda <- 0.5454545
datos1 %<>% mutate(y_lambda = (y^lambda - 1) / lambda)
glimpse(datos1)

## Observations: 100
## Variables: 6
## $ y        <dbl> 412.8432, 129.4738, 368.4621, 190.7698, 161.2475, 209.6…
## $ x1       <dbl> 19.71750, 10.83758, 18.73870, 13.29231, 12.22276, 14.01…
## $ x2       <dbl> 19.23959, 11.01319, 19.53059, 12.03363, 13.91585, 14.87…
## $ x3       <dbl> -1.10617180, 0.14439070, 2.66525949, -0.04207274, 3.392…
## $ x4       <dbl> 1.4862349, 0.1340185, 0.3929311, 2.4853584, 0.2753532, …
## $ y_lambda <dbl> 47.14830, 24.18867, 44.20207, 30.31480, 27.49778, 32.01…

Ahora uso el Y transformado para una nueva regresión lineal:

modelo_boxcox <- lm(y_lambda ~ x1 + x2 + x3 + x4, data = datos1)
summary(modelo_boxcox)

## 
## Call:
## lm(formula = y_lambda ~ x1 + x2 + x3 + x4, data = datos1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.26245 -0.08935 -0.00525  0.08720  0.33909 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.790393   0.064356 -58.897  < 2e-16 ***
## x1           2.465616   0.013856 177.952  < 2e-16 ***
## x2           0.109822   0.013151   8.351 5.46e-13 ***
## x3          -0.178643   0.005883 -30.363  < 2e-16 ***
## x4          -0.004142   0.006002  -0.690    0.492    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1223 on 95 degrees of freedom
## Multiple R-squared:  0.9998, Adjusted R-squared:  0.9998 
## F-statistic: 1.002e+05 on 4 and 95 DF,  p-value: < 2.2e-16

Volvemos a graficar los residuos estandarizados contra los y predichos, ahora con el ajuste boxcox:

datos1_ajuste_boxcox <- augment(modelo_boxcox, data = datos1)

ggplot(data = datos1_ajuste_boxcox) +
  aes(x = .fitted, y = .std.resid) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dotted") +
  geom_smooth(se = F, color = "Gray") +
  labs(
    title = "Punto a con ajuste Box Cox",
    subtitle = "Residuos est. vs. valores predichos",
    y = "Residuos estandarizados",
    x = "Valores predichos (y_hat)"
  )

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Efectivamente, la transformación de Box Cox hace que los errores pierdan la estructura antes observada!

ggplot(datos1_ajuste_boxcox) +
  aes(sample = .std.resid) +
  geom_qq() +
  geom_qq_line() +
  labs(
    title = "Punto b con ajuste Box Cox",
    subtitle = "QQ-plot de los residuos estandarizados",
    y = "Cuantiles observados de los residuos estandarizados",
    x = "Cuantiles teóricos de dist. normal"
  )

Los cuantiles observados se ven más cercanos a los esperados, pero tampoco es ideal.

shapiro.test(datos1_ajuste_boxcox$.std.resid)

## 
##  Shapiro-Wilk normality test
## 
## data:  datos1_ajuste_boxcox$.std.resid
## W = 0.97757, p-value = 0.08579

El p valor aumentó y ahora, por suerte, ya no es significativo: no podemos rechazar la hipótesis de normalidad. El ajuste de Box Cox sirvió.

Punto e

Definir el Cp de Mallows, explicar cómo se utilizaría para seleccionar variables.

El Cp de Mallows permite saber si la cantidad de covariables elegidas para el modelo (p) es similar a la cantidad de covariables realmente explicativas de los valores observados. Se define como:

\[ C_{p} = \frac{ \text{RSS}_{p} }{ s^2 } + 2p - n \approx p \]

En la ecuación de arriba, el \(s^2\) es calculado con el conjunto total de variables disponibles (\(K\) variables, digamos, con \(K > p\)). Se espera que Cp sea aproximadamente p cuando las p variables elegidas son las que explican a las Y.

¿Qué variables seleccionaría si se basar a solamente en este criterio?

Vamos a calcular el Cp de Mallos para todos los subconjuntos posibles de variables, lo que implica hacer una regresión lineal para todos esos subconjuntos. El paquete leaps lo hará por nosotros. Nótese que como variable a explicar ahora usamos la Y transformada por Box Cox.

covariables <- select(datos1, x1, x2, x3, x4)

# ASR: all subsets regression
ASR_call <- leaps::leaps(x = covariables, y = datos1$y_lambda, method = "Cp")

ASR_resultado <-
  tibble(
    p = ASR_call$size,
    cp = ASR_call$Cp
  ) %>%
  rowid_to_column("i") %>%
  mutate(
    which = map(i, ~ ASR_call$which[.x, ]),
    names = map(which, ~ c("x1", "x2", "x3", "x4")[.x]),
    names_peek = map_chr(names, ~ str_c(.x, collapse = " + "))
  )

ASR_resultado %>%
  filter(cp < 1e2) %>% # Quito valores de Cp gigantes, muy lejos de p
  select(p, cp, names_peek)

## # A tibble: 4 x 3
##       p    cp names_peek       
##   <dbl> <dbl> <chr>            
## 1     3 70.8  x1 + x3          
## 2     4  3.48 x1 + x2 + x3     
## 3     4 72.7  x1 + x3 + x4     
## 4     5  5    x1 + x2 + x3 + x4

Basado sólo en el Cp de Mallows, nos interesa encontrar el modelo de menos cantidad de covariables cuyo Cp esté “cerca” de p. En este caso, sería el caso del conjunto de variables: X1, X2, X3, cuyo Cp (3.476) se parece a 3.

Punto f

Considerar el procedimiento Forward para proponer un modelo con dos variables predictoras. Justificar.

La justificación del procedimiento Forward para elegir variables, en lugar de simplemente calcular la regresión para todos los subconjuntos posibles, es que el espacio de conjuntos de variables a recorrer es mucho menor, en particular cuando la cantidad total de covariables disponibles es grande.

forw <- leaps::regsubsets(y_lambda ~ x1 + x2 + x3 + x4, data = datos1, method = "forward")
forw_result <- summary(forw)

forw_result

## Subset selection object
## Call: regsubsets.formula(y_lambda ~ x1 + x2 + x3 + x4, data = datos1, 
##     method = "forward")
## 4 Variables  (and intercept)
##    Forced in Forced out
## x1     FALSE      FALSE
## x2     FALSE      FALSE
## x3     FALSE      FALSE
## x4     FALSE      FALSE
## 1 subsets of each size up to 4
## Selection Algorithm: forward
##          x1  x2  x3  x4 
## 1  ( 1 ) "*" " " " " " "
## 2  ( 1 ) "*" " " "*" " "
## 3  ( 1 ) "*" "*" "*" " "
## 4  ( 1 ) "*" "*" "*" "*"

El método Forward nos sugiere que las variables X1 y X3 son el mejor subconjunto de dos variables predictoras. Veamos cómo es la regresión al elegir esas dos variables:

modelo_forward_2_variables <- lm(y_lambda ~ x1 + x3, data = datos1) 
summary(modelo_forward_2_variables)

## 
## Call:
## lm(formula = y_lambda ~ x1 + x3, data = datos1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.44449 -0.11791 -0.00658  0.08705  0.45987 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.785618   0.082395  -45.95   <2e-16 ***
## x1           2.576163   0.005320  484.27   <2e-16 ***
## x3          -0.182628   0.007639  -23.91   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1594 on 97 degrees of freedom
## Multiple R-squared:  0.9996, Adjusted R-squared:  0.9996 
## F-statistic: 1.18e+05 on 2 and 97 DF,  p-value: < 2.2e-16

Se ve bastante bien, al menos en significación de las variables, de la regresión en general, y del R cuadrado ajustado.

Punto g

Analizar la matrix de correlación de las covariables.

correlacion <- datos1 %>% select(-y) %>% cor() %>% round(4)
correlacion

##               x1      x2      x3      x4 y_lambda
## x1        1.0000  0.9548 -0.0276 -0.0211   0.9986
## x2        0.9548  1.0000 -0.0495  0.0140   0.9585
## x3       -0.0276 -0.0495  1.0000  0.0337  -0.0768
## x4       -0.0211  0.0140  0.0337  1.0000  -0.0222
## y_lambda  0.9986  0.9585 -0.0768 -0.0222   1.0000

X1 y X2 están muy correlacionadas, como se vio antes. La consecuencia de esto será una alta variabilidad en los betas estimados al cambiar de conjunto de datos, y la solución sería un ajuste robusto para paliar ese problema. Pero si nos mantenemos en el marco de ajuste lineal clásico, las variables X1 y X3 elegidas por Forward parecen ser una buena elección. Perdemos la variable X4, pero su correlación con las Y es muy chica. Me inquieta que el Cp de Mallows me sugeriría, por sí solo, elegir modelos que tienen a ambas X1 y X2, pero la alta correlación entre ellas me hace preferir el modelo sugerido por Forward. El modelo quedaría así:

\[ y_{i}* = \beta_{0} + \beta_{1} x_{1} + \beta_{3} x_{3} + \varepsilon_{i} \\ \frac{ y_{i}^{\lambda} - 1 }{ \lambda } = \beta_{0} + \beta_{1} x_{1} + \beta_{3} x_{3} + \varepsilon_{i} \\ \text{con } \lambda = 0.545 \]

modelo_final <- lm(y_lambda ~ x1 + x3, data = datos1)
summary(modelo_final)

## 
## Call:
## lm(formula = y_lambda ~ x1 + x3, data = datos1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.44449 -0.11791 -0.00658  0.08705  0.45987 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.785618   0.082395  -45.95   <2e-16 ***
## x1           2.576163   0.005320  484.27   <2e-16 ***
## x3          -0.182628   0.007639  -23.91   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1594 on 97 degrees of freedom
## Multiple R-squared:  0.9996, Adjusted R-squared:  0.9996 
## F-statistic: 1.18e+05 on 2 and 97 DF,  p-value: < 2.2e-16

Punto 2

datos2 <- read.table("datos2.txt") %>% as_tibble

# Aplico la transformación Box Cox con el lambda elegido arriba:
datos2 %<>% mutate(y_lambda = (y^lambda - 1) / lambda) %>% select(-y)

# Agrego número de observación
datos2 %<>% rowid_to_column("obs_id")

glimpse(datos2)

## Observations: 200
## Variables: 6
## $ obs_id   <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, …
## $ x1       <dbl> 19.71750, 10.83758, 18.73870, 13.29231, 15.22276, 14.01…
## $ x2       <dbl> 19.23959, 11.01319, 19.53059, 12.03363, 13.91585, 14.87…
## $ x3       <dbl> -1.10617180, 0.14439070, 2.66525949, -0.04207274, 3.392…
## $ x4       <dbl> 1.4862349, 0.1340185, 0.3929311, 2.4853584, 0.2753532, …
## $ y_lambda <dbl> 47.148297, 24.188665, 44.202074, 30.314797, 8.829080, 3…

ggpairs(select(datos2, -obs_id), progress=F)

Estimar por mínimos cuadrados los parámetros del modelo seleccionado en 1.g con estos datos:

modelo_2 <- lm(y_lambda ~ x1 + x3, data = datos2)
summary(modelo_2)

## 
## Call:
## lm(formula = y_lambda ~ x1 + x3, data = datos2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -39.990  -0.283   0.682   1.860   4.044 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.9062     1.9804   0.963   0.3370    
## x1            2.1759     0.1270  17.137   <2e-16 ***
## x3           -0.4365     0.1728  -2.526   0.0123 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.184 on 197 degrees of freedom
## Multiple R-squared:  0.6033, Adjusted R-squared:  0.5993 
## F-statistic: 149.8 on 2 and 197 DF,  p-value: < 2.2e-16

Realizar un gráfico de residuos estandarizados versus valores predichos:

datos2_ajuste <- augment(modelo_2, data = datos2)

residuos_top <- top_n(datos2_ajuste, 6, abs(.std.resid))

ggplot(data = datos2_ajuste) +
  aes(x = .fitted, y = .std.resid) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dotted") +
  geom_smooth(se = F, color = "Gray") +
  labs(
    title = "Punto 2",
    subtitle = "Residuos est. vs. valores predichos",
    y = "Residuos estandarizados",
    x = "Valores predichos (y_hat)"
  ) +
  geom_label_repel(
    data = residuos_top,
    mapping = aes(label = obs_id)
  )

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Parece haber outliers que se llevan mucho residuo y (suponemos) distorsionan los betas estimados.

¿Qué observaciones corresponden a candidatos a outliers según el análisis de estos residuos?

Las observaciones 80, 16, 88, 5, 21 y probablemente la 106 también.

¿Qué medidas propone calcular para decidir si son realmente outliers?

Podemos calcular la distancia de Cook y el leverage de cada observación. El método augment(modelo) del paquete broom lo hace por nosotros! Voy a reproducir el “Plot L-R” de leverages contra residuos que hicimos en algún práctico, y sumarle la distancia de Cook como variable codificada en el tamaño y color de las etiquetas:

Primero tenemos que calcular los \(a_{i}^2\) definidos así:

\[ a_{i} = \frac{ e_{i} }{ \sqrt{e^T e}} \]

Y luego graficar \(a_{i}^2\) vs. \(p_{ii}\). En el mismo gráfico voy a incluir el “umbral” de leverage problemático, \(2p/n\), y con el umbral de Di de Cook problemático: el cuantil 0.5 de una distribución F con (p, n-p) grados de libertad.

norma_del_vector_residuos <- sqrt(sum(modelo_2$residuals^2))

p <- modelo_2$rank
n <- nrow(datos2)
pii_umbral <- 2 * p / n
Di_umbral <- qf(p = 0.5, df1 = p, df2 = n - p)

datos2_ajuste %>%
  mutate(a = .resid / norma_del_vector_residuos) %>%
  ggplot() +
  aes(x = a^2, y = .hat, color = .cooksd > Di_umbral) +
  geom_hline(yintercept = pii_umbral) +
  geom_label(aes(label = obs_id), alpha = .5) +
  labs(title = "Plot L-R",
       subtitle = "Leverage vs. Residuos") +
  scale_color_gradient(low = "Gray", high = "Red") +
  guides(size = F) +
  scale_color_manual(values = c("Black", "Red"))

## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Hay mucha información en el plot. Por un lado, sólo la observación 80 es un outlier basado en su valor de D de Cook. Por otro lado, muchas mas observaciones de las que creíamos tienen un leverage mayor al del umbral. Con todo, esto no implica por sí solo que sean outliers. La combinación de alto residuo transformado (a^2) y alto leverage se aplica en especial a las observaciones 16 y 80. Las observaciones 21, 5 y 88, en cambio, tienen residuo alto pero poco leverage, de modo que pueden ser quitadas pero los betas no cambiarían tanto. En general, los candidatos a outliers que habíamos señalado (106, 80, 16, 88, 21, 5) son extremos en algún sentido en este “plot L-R”.

Comparar con el ajuste robusto:

modelo_robusto <- MASS::rlm(y_lambda ~ x1 + x3, data = datos2)
datos2_ajuste_robusto <- augment(modelo_robusto)

comparacion_OLS_Rob <- bind_rows(
    datos2_ajuste %>% mutate(metodo = "OLS"),
    datos2_ajuste_robusto %>% mutate(metodo = "robusto")
  ) %>%
  mutate(weight = rep(modelo_robusto$w, times = 2))

outliers_candidatos <- c(80, 16, 88, 5, 21, 106)

tibble(obs = datos2$obs_id, peso = modelo_robusto$w) %>%
  filter(obs %in% outliers_candidatos)

## # A tibble: 6 x 2
##     obs    peso
##   <int>   <dbl>
## 1     5 0.00637
## 2    16 0.00425
## 3    21 0.00715
## 4    80 0.00377
## 5    88 0.00535
## 6   106 0.0127

Podemos ver que el ajuste robusto efectivamente les da un peso bastante bajo a los candidatos a outliers. Mi conclusión sería quitarlas a todas y correr de nuevo la regresión, en particular si un ajuste robusto les pone un peso tan bajo.

Comparamos finalmente los parámetros estimados con ajuste de mínimos cuadrados y con ajuste robusto:

message("Ajuste OLS:")

## Ajuste OLS:

modelo_2$coefficients

## (Intercept)          x1          x3 
##   1.9061778   2.1759092  -0.4365215

message("\nAjuste Robusto:")

## 
## Ajuste Robusto:

modelo_robusto$coefficients

## (Intercept)          x1          x3 
##  -3.8125825   2.5756972  -0.1744662

El ajuste robusto le da un valor similar a los betas de ambas variables, X1 y X3, pero cambia el intercept. El cambio podría leerse como que, en el ajuste robusto, la variable X1 adquirió más relevancia que la variable X3, pues los cambios de los parámetros beta_1 y beta_3 ocurrieron en sentido opuesto: beta_1 creció y beta_3 decreció (en valor absoluto).