Regresión Lineal Múltiple
Addy Bolivar Cime
5/19/2024
Librerías
# Se llama a la librería readxl para importar datos desde excel
library(readxl)
## Warning: package 'readxl' was built under R version 4.1.2
# Se llama a la librería lattice para gráficas de nubes de puntos
library(lattice)
Importación y visualización de datos
# La instrucción col_names=TRUE sirve para tomar la primera fila como los
nombres de las columnas.
Entrega <- read_excel("Datosentrega.xlsx", col_names=TRUE)
# Se visualizan los datos
View(Entrega)
Estimación por mínimos cuadrados
# Se obtienen los estimadores del modelo de regresión lineal múltiple
modelo.ajust<-lm(Tiempo~Cajas+Distancia, data=Entrega, x=TRUE, y=TRUE)
modelo.ajust
##
## Call:
## lm(formula = Tiempo ~ Cajas + Distancia, data = Entrega, x = TRUE,
## y = TRUE)
##
## Coefficients:
## (Intercept) Cajas Distancia
## 2.34123 1.61591 0.01438
Pruebas individuales de los coeficientes
Se realizan las pruebas t de H0: betaj=0 para cada j.
# Se imprimen los coeficientes estimados, sus pruebas t y el valor del es
tadístico F
summary(modelo.ajust)
##
## Call:
## lm(formula = Tiempo ~ Cajas + Distancia, data = Entrega, x = TRUE,
## y = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.7880 -0.6629 0.4364 1.1566 7.4197
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.341231 1.096730 2.135 0.044170 *
## Cajas 1.615907 0.170735 9.464 3.25e-09 ***
## Distancia 0.014385 0.003613 3.981 0.000631 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.259 on 22 degrees of freedom
## Multiple R-squared: 0.9596, Adjusted R-squared: 0.9559
## F-statistic: 261.2 on 2 and 22 DF, p-value: 4.687e-16
Valores ajustados y residuales
# Valores ajustados (medias estimadas)
Val.ajus<-modelo.ajust$fitted.values
Val.ajus
## 1 2 3 4 5 6 7
8
## 21.708084 10.353615 12.079794 9.955646 14.194398 18.399574 7.155376
16.673395
## 9 10 11 12 13 14 15
16
## 71.820294 19.123587 38.092507 21.593041 12.472991 18.682464 23.328798
29.662928
## 17 18 19 20 21 22 23
24
## 14.913640 15.551379 7.706807 40.887970 20.514179 56.006528 23.357568
24.402854
## 25
## 10.962584
# Residuales
Resid<-modelo.ajust$residuals
Resid
## 1 2 3 4 5 6
7
## -5.0280843 1.1463854 -0.0497937 4.9243539 -0.4443983 -0.2895743 0.8
446235
## 8 9 10 11 12 13
14
## 1.1566049 7.4197062 2.3764129 2.2374930 -0.5930409 1.0270093 1.0
675359
## 15 16 17 18 19 20
21
## 0.6712018 -0.6629284 0.4363603 3.4486213 1.7931935 -5.7879699 -2.6
141789
## 22 23 24 25
## -3.6865279 -4.6075679 -4.5728535 -0.2125839
Diagramas de dispersión
# Matriz de diagramas de dispersión
pairs(~Tiempo+Cajas+Distancia, data=Entrega)
# Diagrama 3D de dispersión (nube de puntos) para el caso de 2 variables
explicativas
cloud(Tiempo~Cajas*Distancia, data=Entrega)
ANOVA
Se realizan pruebas F de H0: betaj=0 para cada j.
anova(modelo.ajust, test="F")
## Analysis of Variance Table
##
## Response: Tiempo
## Df Sum Sq Mean Sq F value Pr(>F)
## Cajas 1 5382.4 5382.4 506.619 < 2.2e-16 ***
## Distancia 1 168.4 168.4 15.851 0.0006312 ***
## Residuals 22 233.7 10.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Se realiza la prueba F de H0: beta1=beta2=…=betak=0.
# Ajuste de un modelo sin variables predictoras (y=b0+epsilon):
fit0 = lm(Tiempo ~ 1, data = Entrega)
anova(fit0, modelo.ajust, test="F")
## Analysis of Variance Table
##
## Model 1: Tiempo ~ 1
## Model 2: Tiempo ~ Cajas + Distancia
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 24 5784.5
## 2 22 233.7 2 5550.8 261.24 4.687e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Valor del punto porcentual de la distribución F con el que
# se compara el valor del estadístico F
qf(0.05, 2,22,lower.tail = FALSE)
## [1] 3.443357
Intervalos de confianza para los coeficientes del modelo
# Intervalos de confianza del 95% para los coeficientes del modelo
confint.lm(modelo.ajust,level=0.95)
## 2.5 % 97.5 %
## (Intercept) 0.066751987 4.61571030
## Cajas 1.261824662 1.96998976
## Distancia 0.006891745 0.02187791
Cálculo de variables de interés
# Matriz X
X=modelo.ajust$x
X
## (Intercept) Cajas Distancia
## 1 1 7 560
## 2 1 3 220
## 3 1 3 340
## 4 1 4 80
## 5 1 6 150
## 6 1 7 330
## 7 1 2 110
## 8 1 7 210
## 9 1 30 1460
## 10 1 5 605
## 11 1 16 688
## 12 1 10 215
## 13 1 4 255
## 14 1 6 462
## 15 1 9 448
## 16 1 10 776
## 17 1 6 200
## 18 1 7 132
## 19 1 3 36
## 20 1 17 770
## 21 1 10 140
## 22 1 26 810
## 23 1 9 450
## 24 1 8 635
## 25 1 4 150
## attr(,"assign")
## [1] 0 1 2
# Matriz C=(X'X)^{-1}
C=solve(t(X)%*%X)
# Tamaño de la muestra
n=nrow(X)
n
## [1] 25
# Suma de cuadrados de residuales
SSres=sum(Resid^2)
SSres
## [1] 233.7317
# Número de coeficientes en el modelo (p=k+1)
p=ncol(X)
p
## [1] 3
# Cuadrado medio residual (estimación de sigma^2)
MSres=SSres/(n-p)
MSres
## [1] 10.62417
# Error estándar de regresión
error.est=sqrt(MSres)
error.est
## [1] 3.259473
# Coeficientes estimados
betas=modelo.ajust$coefficients
betas
## (Intercept) Cajas Distancia
## 2.34123115 1.61590721 0.01438483
Intervalo de confianza para la respuesta media en x0
# Vector x0 de interés
x0=c(1,8,275)
# Media estimada
med.est=t(x0)%*%betas
# Varianza estimada
var.est=MSres*t(x0) %*% C %*% x0
# Nivel de confianza
nivel1=0.95
alfa1=1-nivel1
# Punto porcentual de la t(n-p)
tcuantil=qt(alfa1/2,n-p, lower.tail=FALSE)
# Límite inferior del intervalo
LI=med.est-tcuantil*sqrt(var.est)
# Límite superior del intervalo
LU=med.est+tcuantil*sqrt(var.est)
# Intervalo de confianza
c(LI,LU)
## [1] 17.65390 20.79474
Intervalo de predicción de la respuesta y0 en x0
# Vector x0 de interés
x0=c(1,8,275)
# Media estimada
med.est=t(x0)%*%betas
# Varianza estimada
var.est=MSres*t(x0) %*% C %*% x0
# Nivel de confianza
nivel1=0.95
alfa1=1-nivel1
# Punto porcentual de la t(n-p)
tcuantil=qt(alfa1/2,n-p, lower.tail=FALSE)
# Límite inferior del intervalo
LI=med.est-tcuantil*sqrt(1+var.est)
# Límite superior del intervalo
LU=med.est+tcuantil*sqrt(1+var.est)
# Intervalo de predicción
c(LI,LU)
## [1] 16.62294 21.82569
Gráficas de verificación de los supuestos
# Se obtienen diversas gráficas de verificación de supuestos
plot(modelo.ajust)
# Gráfica de los valores ajustados contra los residuales
plot(Val.ajus, Resid, xlab="Valores ajustados", ylab="Residuales")
# Gráfica cuantil cuantil normal de los residuales
qqnorm(Resid)
qqline(Resid)
# Residuales estandarizados
Resid.est=Resid/error.est
qqnorm(Resid.est)
qqline(Resid.est)
# Gráfica de los residuales contra cada variable regresora
cajas=Entrega$Cajas
distancia=Entrega$Distancia
plot(cajas, Resid, xlab="Cajas", ylab="Residuales")
plot(distancia, Resid, xlab="Distancia", ylab="Residuales")
