Regresión lineal múltiple
Luis Fernando Apáez Álvarez
30 de Agosto de 2022
Ajuste y pruebas de hipótesis
Comenzamos por cargar los datos y ver algunas estadísticas de nuestros datos:
# Cargamos los datos
library(GGally)
Datos=read.csv("ejemplo1RLM.csv", header=TRUE )
summary(Datos)## y X1 X2
## Min. : 5.8 Min. :10.4 Min. :15.4
## 1st Qu.: 8.0 1st Qu.:14.0 1st Qu.:19.0
## Median : 8.7 Median :15.1 Median :20.1
## Mean : 8.7 Mean :15.2 Mean :20.1
## 3rd Qu.: 9.6 3rd Qu.:16.4 3rd Qu.:21.2
## Max. :12.6 Max. :19.4 Max. :24.6Con base en lo anterior podemos ver que tenemos dos variables explicativas \(x_{1}\) y \(x_{2}\). Nos apoyamos mediante un gráfico para ver el comportamiento de nuestros datos
par(mfrow=c(1,2), mar=c(4.5, 4.5, 2, 2))
pairs(Datos)
notamos que nuestras variables explicativas están altamente correlacionadas. De manera intuitiva nos dice que, al incluir ya una variable al modelo, la otro ya no está aportando información adicional para describir nuestro modelo.
# Alternativamente con el paquete GGally:
ggpairs(Datos)
Un posible modelo es: \(\mathbb{E}(Y; X_{1}, X_{2})=\beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}\). Para ello, realizamos el ajuste del modelo
# Ajuste del modelo
fit <- lm(y ~ X1 + X2, data=Datos)
# Informacion del ajuste
summary(fit)##
## Call:
## lm(formula = y ~ X1 + X2, data = Datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.873 -0.661 -0.124 0.621 2.080
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.7370 1.8599 0.93 0.35
## X1 0.3540 0.3315 1.07 0.29
## X2 0.0794 0.3300 0.24 0.81
##
## Residual standard error: 0.95 on 97 degrees of freedom
## Multiple R-squared: 0.413, Adjusted R-squared: 0.401
## F-statistic: 34.2 on 2 and 97 DF, p-value: 5.83e-12donde en la prueba F se nos marca un \(p-value\) muy pequeño, es decir, encontramos evidencia estadística para asumir que nuestro modelo tiene sentido, esto es, hay evidencia para rechazar la hipótesis nula:
\[ H_0: \beta_1=0\ \ y\ \ \beta_2=0 \ \ \ vs \ \ \ H_a: \beta_1\neq 0\ \ o \ \ \beta_2\neq 0 \]
Considerando también la parte de Coefficients, en las
pruebas t individuales de cada variable, notamos que los \(p-values\) son grandes y tendremos que
realizar alguna consideración sobre nuestras variables explicativas, lo
cual veremos más adelante.
Podemos obtener los valores de las principales estimaciones puntuales:
# \hat{\beta}
coef(fit)## (Intercept) X1 X2
## 1.737 0.354 0.079# \hat{V(\hat{\beta})}
vcov(fit)## (Intercept) X1 X2
## (Intercept) 3.46 0.51 -0.55
## X1 0.51 0.11 -0.11
## X2 -0.55 -0.11 0.11# \hat{\sigma}=\sqrt{SCE/(n-p-1)}
sigma(fit)## [1] 0.95Luego, procederemos a obtener la prueba de la tabla Anova mediante la prueba lineal general:
library(multcomp)
K <- matrix(c(0,1,0,
0,0,1), ncol=3, nrow=2, byrow=TRUE)
m <- c(0,0)
# Utilizamos la prueba lineal general
summary(glht(fit, linfct=K, rhs=m), test=Ftest())##
## General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate
## 1 == 0 0.3540
## 2 == 0 0.0794
##
## Global Test:
## F DF1 DF2 Pr(>F)
## 1 34.2 2 97 5.83e-12donde
\[ H_{0}: \left(\begin{array}{cc}\beta_{1}\\\beta_2\\\end{array}\right)=\left(\begin{array}{cc}0 \\ 0 \end{array}\right) \]
Una vez que se realiza la prueba asociada a la tabla Anova, la
pregunta es si ambas variables aportan información al modelado. Para
esto podemos hacer pruebas individuales sobre cada coeficiente usando
las pruebas t como vimos antes, o bien utilizando la prueba lineal
general. Notemos que las pruebas realizadas en el
summary(fit) no están hechas de manera simultánea, y en él,
para la primera línea, se considera \(x_{1}=0\) por lo que estamos contrastando
el modelo reducido \(\mathbb{E}(Y; X_{1},
X_{2})=\beta_{0}+\beta_{2}X_{2}\) versus el modelo completo \(\mathbb{E}(Y; X_{1},
X_{2})=\beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}\), en otras
palabras estamos constrastando \(H_0:
\beta_1=0 \ \ \ vs\ \ \ H_a: \beta_1\neq 0\), y dado que
obtuvimos un \(p-value\) grande,
entonces no encontramos evidencia suficiente para poder rechazar la
hipótesis nula. De manera análoga con la prueba t para la variable \(x_{2}\).
Luego, podemos realizar una prueba análoga pero ahora utilizando la
paquetería multcom
# H0: b1 = 0 vs Ha: b1 != 0
K <- matrix(c(0,1,0), ncol=3, nrow=1, byrow=TRUE)
m <- c(0)
# Usando prueba F equivalente (no permite pruebas de una cola con alternativa > o <)
summary(glht(fit, linfct=K, rhs=m), test=Ftest())##
## General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate
## 1 == 0 0.354
##
## Global Test:
## F DF1 DF2 Pr(>F)
## 1 1.14 1 97 0.288Ahora bien, la prueba anterior es sólo un caso particular en donde se
comparan modelos anidados, es decir, donde se comparan el modelo
reducido versus el completo. R tiene otras formas de
obtener las pruebas usando el lenguaje de modelos anidados, por ejemplo,
la función anova nos sirve para comparar dos modelos anidados sin
necesidad de usar la definición matricial (K y m).
# Por ejemplo si b1=0, entonces el modelo reducido es
# y = b0 + b2 X2 + e
# Lo ajustamos en R
fitred1 <- lm(y ~ X2, data=Datos)
summary(fitred1)##
## Call:
## lm(formula = y ~ X2, data = Datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.9186 -0.7020 -0.0832 0.5960 2.2209
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.1010 1.0553 0.10 0.92
## X2 0.4273 0.0522 8.19 9.8e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.95 on 98 degrees of freedom
## Multiple R-squared: 0.406, Adjusted R-squared: 0.4
## F-statistic: 67.1 on 1 and 98 DF, p-value: 9.84e-13Ahora comparamos el modelo completo y el reducido:
# H0: El modelo reducido es plausible
# Ha: El modelo completo es plausible
anova(fit, fitred1)## Analysis of Variance Table
##
## Model 1: y ~ X1 + X2
## Model 2: y ~ X2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 97 87.8
## 2 98 88.8 -1 -1.03 1.14 0.29donde no se rechaza la hipótesis nula por lo cual hay evidencia estadística para considerar plausible el modelo reducido.
Algo similar para obtener la prueba de la tabla anova, donde contrastabamos si las betas eran en simultáneo cero o no, consideraremos la comparación entre el modelo reducido \(y=\beta_0+e\) cuando \(H_0: \beta_1=0 \ \ y \ \ \beta_2=0\) versus el modelo completo
# Modelo reducido simple, b0=0=b1
fitred2 <- lm(y ~ 1, data=Datos)
# Realizamos la prueba de hipotesis pertienente
anova(fitred2, fit)## Analysis of Variance Table
##
## Model 1: y ~ 1
## Model 2: y ~ X1 + X2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 99 149.6
## 2 97 87.8 2 61.9 34.2 5.8e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1se rechaza \(H_{0}\), lo que nos dice que es plausible asumir el modelo completo versus el modelo reducido.
Selección entre los modelos.
Aquí casi los tres modelos indican un mismo coeficiente de
determinación, se seleccionaría en general el de menor número de
parámetros y con el mayor coeficiente de determinación, aunque más
adelante se considerarían otros criterios (BIC y
AIC) o bien se realiza un análisis del poder predictivo del
modelo.
summary(fit)##
## Call:
## lm(formula = y ~ X1 + X2, data = Datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.873 -0.661 -0.124 0.621 2.080
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.7370 1.8599 0.93 0.35
## X1 0.3540 0.3315 1.07 0.29
## X2 0.0794 0.3300 0.24 0.81
##
## Residual standard error: 0.95 on 97 degrees of freedom
## Multiple R-squared: 0.413, Adjusted R-squared: 0.401
## F-statistic: 34.2 on 2 and 97 DF, p-value: 5.83e-12fitred1 <- lm(y ~ X2, data=Datos)
summary(fitred1)##
## Call:
## lm(formula = y ~ X2, data = Datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.9186 -0.7020 -0.0832 0.5960 2.2209
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.1010 1.0553 0.10 0.92
## X2 0.4273 0.0522 8.19 9.8e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.95 on 98 degrees of freedom
## Multiple R-squared: 0.406, Adjusted R-squared: 0.4
## F-statistic: 67.1 on 1 and 98 DF, p-value: 9.84e-13fitred2 <- lm(y ~ X1, data=Datos)
summary(fitred2)##
## Call:
## lm(formula = y ~ X1, data = Datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.882 -0.670 -0.145 0.605 2.053
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.1408 0.7968 2.69 0.0085 **
## X1 0.4328 0.0521 8.30 5.7e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.95 on 98 degrees of freedom
## Multiple R-squared: 0.413, Adjusted R-squared: 0.407
## F-statistic: 69 on 1 and 98 DF, p-value: 5.67e-13con base en lo dicho anteriormente elegiremos el segundo modelo: \(y = \beta_{0}+\beta_{1}x_{1}+e\).