Ajuste e interpretación de coeficientes asociados a variables continuas

Comenzamos por cargar nuestros datos:

Datos <- read.csv("ejemplo2RLM.csv", header=TRUE)
summary(Datos)
##        y              X1             X2      
##  Min.   : 7.0   Min.   :10.4   Min.   :-0.1  
##  1st Qu.: 8.8   1st Qu.:14.0   1st Qu.: 1.2  
##  Median : 9.4   Median :15.1   Median : 1.8  
##  Mean   : 9.4   Mean   :15.2   Mean   : 1.9  
##  3rd Qu.:10.0   3rd Qu.:16.4   3rd Qu.: 2.5  
##  Max.   :11.8   Max.   :19.4   Max.   : 5.2

vemos un gráfico inicial para ver un primer comportamiento de nuestros datos

library(GGally)
ggpairs(Datos)

vemos que la variable \(x_{1}\) parece describir bien la variable de respuesta \(y\), pero en el caso de la variable \(x_{2}\) notamos una variabilidad alta. En principio, nos gustaría quedarnos con la variable \(x_{1}\), no obstante, debemos ver si la variable \(x_{2}\) está aportando información al modelo.

# Consideremos el modelo:
# E(y; X1, X2) = b0 + b1X1 + b2X2

fit <- lm(y~ X1 + X2, data=Datos)
summary(fit)
## 
## Call:
## lm(formula = y ~ X1 + X2, data = Datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.09365 -0.03304 -0.00622  0.03107  0.10399 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.05431    0.04161    25.3   <2e-16 ***
## X1           0.49667    0.00262   189.4   <2e-16 ***
## X2           0.40119    0.00495    81.0   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.048 on 97 degrees of freedom
## Multiple R-squared:  0.998,  Adjusted R-squared:  0.998 
## F-statistic: 2.05e+04 on 2 and 97 DF,  p-value: <2e-16

De la prueba F notamos que el modelo tiene sentido (\(\beta_{1}\neq 0\) o \(\beta_1\neq 0\)). Ahora, en las pruebas individuales estamos rechazando la prueba t, es decir, estamos encontrando evidencia estadística a favor de que las variables explicativas, ambas, aportan información al modelo.

En resumen, hasta ahora, tenemos que

\[ H_0: \beta_1=0 \ \ y \ \ \beta_2=0 \ \ vs\ \ H_a: \beta_1\neq0 \ \ o \ \ \beta_2\neq0 \]

Ahora vamos a interpretar este modelo:

Notemos además que al incluir dos variables en el modelo se puede tener un mejor ajuste con base en el coeficiente de determinación. Por ejemplo, los ajustes con una sola variable son:

fitred1 <- lm(y ~ X2, data=Datos)
summary(fitred1)
## 
## Call:
## lm(formula = y ~ X2, data = Datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3496 -0.5958 -0.0625  0.5278  2.2163 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   8.6821     0.2011   43.17   <2e-16 ***
## X2            0.3547     0.0947    3.74    3e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.91 on 98 degrees of freedom
## Multiple R-squared:  0.125,  Adjusted R-squared:  0.116 
## F-statistic:   14 on 1 and 98 DF,  p-value: 0.000305
fitred2 <- lm(y ~ X1, data=Datos)
summary(fitred2)
## 
## Call:
## lm(formula = y ~ X1, data = Datos)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -0.793 -0.308 -0.028  0.240  1.368 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.9733     0.3302    5.98  3.7e-08 ***
## X1            0.4861     0.0216   22.51  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.39 on 98 degrees of freedom
## Multiple R-squared:  0.838,  Adjusted R-squared:  0.836 
## F-statistic:  507 on 1 and 98 DF,  p-value: <2e-16

vemos para el caso del modelo reducido \(y=\beta_{0}+\beta_{2}x_{2}+e\) que el coeficiente de determinación bajó drásticamente. Por otro lado, para el modelo \(y=\beta_{0}+\beta_{1}x_{1}+e\), de nuevo, el coeficiente de determinación bajó, pero no está tan alejado del valor inicial obtenido, por lo cual éste sería un modelo bueno, no obstante, el modelo que incluye a ambas variables es mejor.

#Alternativamente
fitred12 <- lm(y~.-X1, data = Datos)
summary(fitred12)
## 
## Call:
## lm(formula = y ~ . - X1, data = Datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3496 -0.5958 -0.0625  0.5278  2.2163 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   8.6821     0.2011   43.17   <2e-16 ***
## X2            0.3547     0.0947    3.74    3e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.91 on 98 degrees of freedom
## Multiple R-squared:  0.125,  Adjusted R-squared:  0.116 
## F-statistic:   14 on 1 and 98 DF,  p-value: 0.000305
fitred22 <- lm(y~.-X2, data = Datos)
summary(fitred12)
## 
## Call:
## lm(formula = y ~ . - X1, data = Datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3496 -0.5958 -0.0625  0.5278  2.2163 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   8.6821     0.2011   43.17   <2e-16 ***
## X2            0.3547     0.0947    3.74    3e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.91 on 98 degrees of freedom
## Multiple R-squared:  0.125,  Adjusted R-squared:  0.116 
## F-statistic:   14 on 1 and 98 DF,  p-value: 0.000305

donde el . indica que estamos incluyendo todas las variables de la BD; - indica que estamos quitando alguna variable, por ejemplo -X1 indica que estamos quitando esa variable del modelo.

Otra alternativa cuando ya teníamos un modelo predefinido:

# Del modelo fit, lo actualizamos almacenandolo en
# fitred13 y despues de considerar todas las variables
# del modelo fit, quitaremos la variable x1
fitred13 = update(fit, ~.-X1)
summary(fitred13)
## 
## Call:
## lm(formula = y ~ X2, data = Datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3496 -0.5958 -0.0625  0.5278  2.2163 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   8.6821     0.2011   43.17   <2e-16 ***
## X2            0.3547     0.0947    3.74    3e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.91 on 98 degrees of freedom
## Multiple R-squared:  0.125,  Adjusted R-squared:  0.116 
## F-statistic:   14 on 1 and 98 DF,  p-value: 0.000305
# Del modelo fit, lo actualizamos almacenandolo en
# fitred23 y despues de considerar todas las variables
# del modelo fit, quitaremos la variable x2
fitred23 = update(fit, ~.-X2)
summary(fitred23)
## 
## Call:
## lm(formula = y ~ X1, data = Datos)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -0.793 -0.308 -0.028  0.240  1.368 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.9733     0.3302    5.98  3.7e-08 ***
## X1            0.4861     0.0216   22.51  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.39 on 98 degrees of freedom
## Multiple R-squared:  0.838,  Adjusted R-squared:  0.836 
## F-statistic:  507 on 1 and 98 DF,  p-value: <2e-16