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 \]
Si analizamos la prueba t para el coeficiente \(\beta_1\), se rechaza \(H_0\). Aquí el contraste es \(H_0: \beta_1=0 \ \ vs \ \ H_a: \beta_1 \neq 0\) Esto nos indica, que aún considerando a la variable \(x_{2}\) en el modelo, la variable \(X_1\) nos está agregando información para modelar \(\mathbb{E}(Y;X_1,X_2)\). Análogamente para la otra prueba t referente a la variable \(x_{2}\).
Con base en lo anterior, parece que no podríamos reducir el modelo, es decir, todos los coeficientes parecen significativos.
Ahora vamos a interpretar este modelo:
\(R^{2}\), el coeficiente de determinación, nos indica que se está explicando el 99.8% de la variabilidad observada en \(Y\) a través del modelo \(y=\beta_0+\beta_1 X_1+\beta_2X_2+e\). Parece un muy buen ajuste.
Además, \(\beta_1\) se puede interpretar como, condicionado en un valor fijo de \(X_2\), el promedio de la variable \(Y\) AUMENTA en .5 unidades al aumentar en una unidad la variable \(X_1\).
Por otro lado, \(\beta_2\) se puede interpretar como, condicionado en un valor fijo de \(X_1\), el promedio de la variable \(Y\) AUMENTA en .4 unidades al aumentar en una unidad la variable \(X_2\).
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