Examen Final Modelos Lineales 1
Kevin Gargate Osorio
2024-07-13
# Realizar la prueba de Goldfeld-Quandt
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(sandwich)
library(lmtest)Pregunta 1
data = read.csv("procesadores_rendimiento.csv")
data## Tiempo Procesador Carga
## 1 2.05 P1 L1
## 2 1.99 P1 L1
## 3 2.06 P1 L1
## 4 2.15 P1 L1
## 5 1.98 P1 L1
## 6 1.98 P1 L1
## 7 2.16 P1 L2
## 8 2.08 P1 L2
## 9 1.95 P1 L2
## 10 2.05 P1 L2
## 11 1.95 P1 L2
## 12 1.95 P1 L2
## 13 2.02 P1 L3
## 14 1.81 P1 L3
## 15 1.83 P1 L3
## 16 1.94 P1 L3
## 17 1.90 P1 L3
## 18 2.03 P1 L3
## 19 2.41 P2 L1
## 20 2.36 P2 L1
## 21 2.65 P2 L1
## 22 2.48 P2 L1
## 23 2.51 P2 L1
## 24 2.36 P2 L1
## 25 2.45 P2 L2
## 26 2.51 P2 L2
## 27 2.38 P2 L2
## 28 2.54 P2 L2
## 29 2.44 P2 L2
## 30 2.47 P2 L2
## 31 2.44 P2 L3
## 32 2.69 P2 L3
## 33 2.50 P2 L3
## 34 2.39 P2 L3
## 35 2.58 P2 L3
## 36 2.38 P2 L3
## 37 2.02 P3 L1
## 38 1.80 P3 L1
## 39 1.87 P3 L1
## 40 2.02 P3 L1
## 41 2.07 P3 L1
## 42 2.02 P3 L1
## 43 1.99 P3 L2
## 44 1.97 P3 L2
## 45 1.85 P3 L2
## 46 1.93 P3 L2
## 47 1.95 P3 L2
## 48 2.11 P3 L2
## 49 2.03 P3 L3
## 50 1.82 P3 L3
## 51 2.03 P3 L3
## 52 1.96 P3 L3
## 53 1.93 P3 L3
## 54 2.06 P3 L3# Convertir a factores
data$Procesador <- as.factor(data$Procesador)
data$Carga <- as.factor(data$Carga)Procesador: Hipótesis nula (procesador): No hay diferencias significativas en los tiempos de procesamiento entre los diferentes procesadores.
Hipótesis alternativa (procesador): Hay al menos una diferencia significativa en los tiempos de procesamiento entre los diferentes procesadores
Carga:
Hipótesis nula: No hay diferencias significativas en los tiempos de procesamiento entre los diferentes niveles de carga
Hipótesis alternativa: Hay al menos una diferencia significativa en los tiempos de procesamiento entre los diferentes niveles de carga
Interacción Procesador y Carga:
Hipótesis nula: No hay interacción significativa entre el procesador y la carga de trabajo en los tiempos de procesamiento
Hipótesis alternativa: Hay una interacción significativa entre el procesador y la carga de trabajo en los tiempos de procesamiento.
# ANOVA
anova_model <- aov(Tiempo ~ Procesador * Carga, data = data)
summary(anova_model)## Df Sum Sq Mean Sq F value Pr(>F)
## Procesador 2 2.9294 1.4647 170.810 <2e-16 ***
## Carga 2 0.0064 0.0032 0.376 0.689
## Procesador:Carga 4 0.0448 0.0112 1.305 0.283
## Residuals 45 0.3859 0.0086
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1i)Procesador: El efecto del procesador es altamente significativo. Esto significa que hay diferencias significativas entre los tiempos de procesamiento de los diferentes procesadores.
ii)Carga: El efecto de la carga de trabajo no es significativo. Esto significa que no hay diferencias significativas entre los tiempos de procesamiento de los diferentes niveles de carga de trabajo.
iii)Interacción: La interacción entre el procesador y la carga de trabajo no es significativa. Esto significa que el efecto del procesador no depende del nivel de carga de trabajo y viceversa.
Dado que solo el procesador tiene un efecto significativo en el tiempo de procesamiento, debemos centrarnos en seleccionar el procesador con el mejor rendimiento (menor tiempo de procesamiento). No es necesario considerar la carga de trabajo, ya que no tiene un efecto significativo.
tukey <- TukeyHSD(anova_model, "Procesador")
print(tukey)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Tiempo ~ Procesador * Carga, data = data)
##
## $Procesador
## diff lwr upr p adj
## P2-P1 0.4811111 0.40630038 0.55592184 0.0000000
## P3-P1 -0.0250000 -0.09981073 0.04981073 0.6988825
## P3-P2 -0.5061111 -0.58092184 -0.43130038 0.0000000P2 - P1: La diferencia media es de 0.481 segundos, con un intervalo de confianza al 95% entre 0.406 y 0.556 segundos. La p-valor ajustada es < 0.0001, lo que indica que P2 es significativamente más lento que P1.
P3 - P1: La diferencia media es de -0.025 segundos, con un intervalo de confianza al 95% entre -0.100 y 0.050 segundos. La p-valor ajustada es 0.6989, lo que indica que no hay una diferencia significativa entre P3 y P1
P3 - P2: La diferencia media es de -0.506 segundos, con un intervalo de confianza al 95% entre -0.581 y -0.431 segundos. La p-valor ajustada es < 0.0001, lo que indica que P3 es significativamente más rápido que P2
Entonces:
Dado que los resultados muestran diferencias significativas entre algunos procesadores pero no entre todos, las conclusiones son las siguientes:
- P2 es significativamente más lento que P1.
- P3 no es significativamente diferente de P1.
- P3 es significativamente más rápido que P2.
Dado que el efecto de la carga de trabajo no es significativo (p = 0.689), cualquier nivel de carga (L1, L2, L3) puede ser utilizado sin afectar significativamente el rendimiento del procesador
Recomendaciones:
Procesador recomendado: P3, ya que ofrece un rendimiento similar a P1 y significativamente mejor que P2. Nivel de carga recomendado: Cualquier nivel de carga (L1, L2, L3), ya que no hay diferencias significativas en los tiempos de procesamiento entre los diferentes niveles de carga.
Pregunta 2
a
Análisis de Correlación:
De esta matriz se observa que y tiene una correlación negativa muy baja con X1(−0.1308951) y una correlación positiva moderada con X2 (0.6516615)
Análisis de Regresión:
La pendiente no es significativamente diferente de cero (p-valor > 0.05), lo que indica que no hay evidencia suficiente para afirmar una relación lineal significativa entre y y X1
Regresión de y en función de X2 La pendiente de X2 tiene un p-valor cercano a 0.05, lo que sugiere una posible relación significativa
El R2ajustado indica que aproximadamente el 32.88% de la variabilidad en y puede ser explicada por X2
X2 tiene una relación más fuerte y posiblemente significativa con y,mientras que X1 no muestra una relación significativa con y
b
La correlacion de X1 y X2 es 0.6250541. Lo cual sugiere una correlación moderada a fuerte
Sabemos que VIF = 1/(1-R^2)
Podriamos asumir para este caso que el R^2 es 0.6250541 ^2 = 0.3906926
r2 = 0.6250541 ^2
VIF = 1/(1-r2);VIF## [1] 1.641208Dado que el VIF calculado es aproximadamente 1.641208, esto sugiere que no hay problemas graves de multicolinealidad entre X1 y X2
c
Hipotesis Nula: beta1 = beta2 = 0 Hipótesis alternativa : el modelo es significativo
El p-valor de la prueba F es 0.07999, que es mayor que 0.05. Esto significa que no podemos rechazar la hipótesis nula al nivel de significancia del 5%. Por lo tanto, no hay suficiente evidencia para afirmar que el modelo en su conjunto es significativo
En la parte a), se determinó que X1 no tiene una relación significativa con y, mientras que X2 muestra una relación significativa con Y
En la parte b), se determinó que no hay problemas graves de multicolinealidad
Los resultados de la prueba de significancia de la regresión no contradicen los hallazgos anteriores. Aunque X2 muestra una relación moderadamente significativa con y el modelo completo no es significativo a un nivel de alpha = 0.05. Esto puede deberse a la falta de poder estadístico dado el pequeño tamaño de la muestra (n = 8). Además, la falta de significancia del modelo completo se alinea con la observación de que X1 no es significativo # Pregunta 3
a
vacation = read.csv("vacation.csv")
head(vacation)## miles income age kids
## 1 902 41 26 0
## 2 491 31 38 3
## 3 1841 87 40 2
## 4 406 54 48 4
## 5 0 77 43 4
## 6 1899 70 55 2dim(vacation)## [1] 200 4# Cálculo de las millas recorridas
lm1 = lm(miles ~ income + age + kids,data=vacation)
summary(lm1)##
## Call:
## lm(formula = miles ~ income + age + kids, data = vacation)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1198.14 -295.31 17.98 287.54 1549.41
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -391.548 169.775 -2.306 0.0221 *
## income 14.201 1.800 7.889 2.10e-13 ***
## age 15.741 3.757 4.189 4.23e-05 ***
## kids -81.826 27.130 -3.016 0.0029 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 452.3 on 196 degrees of freedom
## Multiple R-squared: 0.3406, Adjusted R-squared: 0.3305
## F-statistic: 33.75 on 3 and 196 DF, p-value: < 2.2e-16confint(lm1,level = 0.95)["kids",]## 2.5 % 97.5 %
## -135.32981 -28.32302De acuerdo al summary podemos decir de la variable kids:
El coeficiente de la variable kids es -81.826, lo cual indica que, manteniendo constantes las variables income y age, cada niño adicional en el hogar está asociado con una disminución de aproximadamente 81.826 millas recorridas por año.
Intervalo de Confianza: El intervalo de confianza del 95% para este coeficiente es [-135.33, -28.32]. Esto significa que estamos un 95% seguros de que el verdadero efecto de un niño adicional en las millas recorridas por año se encuentra entre -135.33 y -28.32 millas, manteniendo constantes las otras dos variables (income y age).
Dado que el intervalo de confianza no incluye el valor 0, podemos concluir que el efecto de la variable kids es estadísticamente significativo al nivel del 5%. En otras palabras, hay evidencia suficiente para afirmar que la cantidad de niños tiene un efecto negativo y significativo en las millas recorridas por año por los hogares.
b
# Configurar la ventana gráfica para múltiples gráficos
par(mfrow = c(2,3))
# Gráfico de secuencia de los residuos
plot(lm1$residuals, main = "Run-Sequence Plot", ylab = "Residuos", xlab = "Índice")
# Gráfico de residuos vs valores ajustados
plot(lm1$fitted.values, lm1$residuals, main = "Residuals vs Fitted",
xlab = "Valores Ajustados", ylab = "Residuos")
abline(h = 0, col = "red")
# Extraer las variables del data frame 'vacation'
income <- vacation$income
age <- vacation$age
kids <- vacation$kids
# Gráfico de residuos vs INCOME
plot(income, lm1$residuals, main = bquote("Residuals vs"~income),
xlab = "INCOME", ylab = "Residuos")
abline(h = 0, col = "red")
# Gráfico de residuos vs AGE
plot(age, lm1$residuals, main = bquote("Residuals vs"~age),
xlab = "AGE", ylab = "Residuos")
abline(h = 0, col = "red")
# Gráfico de residuos vs KIDS
plot(kids, lm1$residuals, main = bquote("Residuals vs"~kids),
xlab = "KIDS", ylab = "Residuos")
abline(h = 0, col = "red")
i) Run-Sequence Plot No parece haber un patrón claro o tendencia en el
gráfico de secuencia de residuos, lo que sugiere que los residuos son
independientes a lo largo de las observaciones
Residuals vs Fitted No parece haber un patrón claro en los residuos en función de los valores ajustados. Sin embargo, hay algunos residuos que se desvían significativamente de 0, lo que podría sugerir la presencia de algunos outliers. La dispersión de los residuos parece ser relativamente constante a lo largo de los valores ajustados, lo que sugiere que no hay una fuerte evidencia de heterocedasticidad
Residuals vs income Los residuos no muestran un patrón claro en relación con income. La dispersión de los residuos parece ser relativamente constante a lo largo de los diferentes valores de income, lo que sugiere que no hay una fuerte evidencia de heterocedasticidad con respecto a esta variable
Residuals vs age Los residuos no muestran un patrón claro en relación con age. La dispersión de los residuos parece ser relativamente constante a lo largo de los diferentes valores de age, lo que sugiere que no hay una fuerte evidencia de heterocedasticidad con respecto a esta variable
Residuals vs kids
Los residuos no muestran un patrón claro en relación con kids. La dispersión de los residuos parece ser relativamente constante a lo largo de los diferentes valores de kids, lo que sugiere que no hay una fuerte evidencia de heterocedasticidad con respecto a esta variable
c
df_ordenado = vacation[order(vacation$income),]
df_p90_u90 = rbind(head(df_ordenado,90),tail(df_ordenado,90))
modelo_p90_u90 = lm(miles ~ income + age + kids, data = df_p90_u90)
summary(modelo_p90_u90)##
## Call:
## lm(formula = miles ~ income + age + kids, data = df_p90_u90)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1179.65 -303.73 24.95 280.24 1599.58
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -463.191 174.630 -2.652 0.00872 **
## income 13.978 1.811 7.719 8.43e-13 ***
## age 17.946 3.938 4.558 9.65e-06 ***
## kids -91.931 28.759 -3.197 0.00165 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 453.5 on 176 degrees of freedom
## Multiple R-squared: 0.3715, Adjusted R-squared: 0.3608
## F-statistic: 34.68 on 3 and 176 DF, p-value: < 2.2e-16GQ_t = lmtest::gqtest(modelo_p90_u90, alternative="two.sided", order.by = ~ df_p90_u90$income)
print(GQ_t)##
## Goldfeld-Quandt test
##
## data: modelo_p90_u90
## GQ = 3.1041, df1 = 86, df2 = 86, p-value = 3.28e-07
## alternative hypothesis: variance changes from segment 1 to 2d
hc3 = lmtest::coeftest(lm1,vcov. = sandwich::vcovHC(lm1,type="HC3"))# Extraer el coeficiente y el error estándar de la variable 'kids'
coef_kids <- hc3["kids", "Estimate"]
se_kids <- hc3["kids", "Std. Error"]
# Calcular el intervalo de confianza del 95%
alpha <- 0.05
z <- qnorm(1 - alpha/2)
ci_lower <- coef_kids - z * se_kids
ci_upper <- coef_kids + z * se_kids
# Mostrar el intervalo de confianza
cat("Intervalo de confianza del 95% para el efecto de 'kids' utilizando HC1: [", ci_lower, ", ", ci_upper, "]\n")## Intervalo de confianza del 95% para el efecto de 'kids' utilizando HC1: [ -139.9676 , -23.68529 ]e
# Crear un dataframe de los residuos y la variable INCOME
resid_data <- data.frame(log_e2 = log(lm1$residuals^2), INCOME = vacation$income)
# Ajustar un modelo de regresión para los residuos log-transformados
resid_mdl <- lm(log_e2 ~ INCOME, data = resid_data)
summary(resid_mdl)##
## Call:
## lm(formula = log_e2 ~ INCOME, data = resid_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.0059 -1.2400 0.5879 1.5703 3.2315
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.566282 0.580407 14.759 < 2e-16 ***
## INCOME 0.034446 0.008729 3.946 0.00011 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.258 on 198 degrees of freedom
## Multiple R-squared: 0.07291, Adjusted R-squared: 0.06823
## F-statistic: 15.57 on 1 and 198 DF, p-value: 0.0001103# Calcular la varianza estimada
h_est <- exp(resid_mdl$fitted.values)
# Calcular los pesos inversos basados en la estimación de la varianza
weights <- 1 / h_est
# Ajustar el modelo ponderado usando MCO con los pesos inversos
lm_wls <- lm(miles ~ income + age + kids, data = vacation, weights = weights)
summary(lm_wls)##
## Call:
## lm(formula = miles ~ income + age + kids, data = vacation, weights = weights)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -4.4099 -1.4571 0.0718 1.2464 5.0784
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -424.223 131.836 -3.218 0.001512 **
## income 13.733 1.657 8.289 1.79e-14 ***
## age 16.936 3.120 5.428 1.67e-07 ***
## kids -75.348 22.550 -3.341 0.000998 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.916 on 196 degrees of freedom
## Multiple R-squared: 0.4033, Adjusted R-squared: 0.3942
## F-statistic: 44.16 on 3 and 196 DF, p-value: < 2.2e-16# Calcular el intervalo de confianza del 95% para el coeficiente de KIDS
confint_wls <- confint(lm_wls, level = 0.95)["kids", ]
print(confint_wls)## 2.5 % 97.5 %
## -119.81908 -30.87671IC hallado en a: IC: 2.5 % 97.5 % -135.32981 -28.32302
IC hallado en d: [ -139.9676 , -23.68529 ]
IC hallado e: 2.5 % 97.5 % -119.81908 -30.87671
El intervalo más estrecho es el de MCG, seguido de MCO, y el más ancho es el de HCE3. Esto indica que el modelo ponderado por MCG proporciona la estimación más precisa (menor varianza de los errores).
Pregunta 4
a
traffic = read.csv("traffic.csv")
head(traffic)## CONGESTION TRAFFIC_VOLUME ROAD_CONDITION WEATHER
## 1 54.08679 2251.093 1 0
## 2 48.58948 2610.467 0 0
## 3 57.05467 2191.286 0 1
## 4 37.87363 2319.645 1 1
## 5 60.25308 2878.645 0 1
## 6 26.34060 531.758 1 1# Ajustar el modelo usando MCO
modelo_mco <- lm(CONGESTION ~ TRAFFIC_VOLUME + ROAD_CONDITION + WEATHER, data = traffic)
# Resumen del modelo para ver los coeficientes estimados y sus errores estándar
summary(modelo_mco)##
## Call:
## lm(formula = CONGESTION ~ TRAFFIC_VOLUME + ROAD_CONDITION + WEATHER,
## data = traffic)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.3787 -4.6428 -0.0118 5.0151 15.9809
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.923e+01 1.530e+00 32.175 < 2e-16 ***
## TRAFFIC_VOLUME 5.326e-03 6.866e-04 7.756 4.68e-13 ***
## ROAD_CONDITION -1.145e+01 1.007e+00 -11.371 < 2e-16 ***
## WEATHER -5.061e+00 1.003e+00 -5.045 1.03e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.091 on 196 degrees of freedom
## Multiple R-squared: 0.5351, Adjusted R-squared: 0.528
## F-statistic: 75.2 on 3 and 196 DF, p-value: < 2.2e-16TRAFFIC_VOLUME (beta1):Por cada incremento de mil vehículos en el volumen de tráfico diario, el índice de congestión diaria promedio aumenta en 0.005326 unidades, manteniendo constantes las otras variables. Este efecto es estadísticamente significativo
ROAD_CONDITION (beta2):Los días con buena condición de carretera (ROAD_CONDITION = 1) tienen, en promedio, un índice de congestión diaria 11.45 unidades menor que los días con mala condición de carretera (ROAD_CONDITION = 0), manteniendo constantes las otras variables. Este efecto es estadísticamente significativo
WEATHER (beta3):Los días soleados (WEATHER = 1) tienen, en promedio, un índice de congestión diaria 5.061 unidades menor que los días lluviosos (WEATHER = 0), manteniendo constantes las otras variables. Este efecto es estadísticamente significativo
b
# Obtener los residuos del modelo
residuos <- residuals(modelo_mco)
# Gráfico de residuos versus tiempo
plot(residuos, type = "o", main = "Residuos vs. Tiempo", ylab = "Residuos", xlab = "Tiempo")
abline(h = 0, col = "red")
# Gráfico de la función de autocorrelación de los residuos
acf(residuos, main = "Función de Autocorrelación de los Residuos")
Las barras de la ACF para los primeros lags serían significativamente
diferentes de cero Las primeras barras (lags) están significativamente
por encima de las líneas de significancia (líneas azules)
c
- Prueba de Durbin Watson: Hipótesis nula: No hay autocorrelación en los residuos Hipótesis alternativa: Hay autocorrelación en los residuos
dw_test <- dwtest(modelo_mco)
print(dw_test)##
## Durbin-Watson test
##
## data: modelo_mco
## DW = 0.79026, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0Interpretación: El p-valor es extremadamente pequeño (< 2.2e-16), lo que indica que podemos rechazar la hipótesis nula. Esto confirma la presencia de autocorrelación positiva en los residuos.
- Prueba de Ljung-Box Hipótesis nula: No hay autocorrelación en los residuos para cualquier lag Hipótesis alternativa: Hay autocorrelación en los residuos para al menos un lag
lb_test <- Box.test(residuos, lag = 6, type = "Ljung-Box")
print(lb_test)##
## Box-Ljung test
##
## data: residuos
## X-squared = 172.44, df = 6, p-value < 2.2e-16Interpretación: El p-valor es extremadamente pequeño (< 2.2e-16), lo que indica que podemos rechazar la hipótesis nula. Esto sugiere que hay autocorrelación en los residuos para al menos uno de los primeros 6 lags.
- Prueba de Breusch-Godfrey Hipótesis nula: No hay autocorrelación de orden en los residuos Hipótesis alternativa: Hay autocorrelación de orden en los residuos
bg_test <- bgtest(modelo_mco, order = 6)
print(bg_test)##
## Breusch-Godfrey test for serial correlation of order up to 6
##
## data: modelo_mco
## LM test = 87.446, df = 6, p-value < 2.2e-16Interpretación: El p-valor es extremadamente pequeño (< 2.2e-16), lo que indica que podemos rechazar la hipótesis nula. Esto sugiere la presencia de autocorrelación de orden hasta 6 en los residuos del modelo.
d
library(sandwich)
# Coeficientes y errores estándar HAC usando Newey-West
print(lmtest::coeftest(modelo_mco, sandwich::NeweyWest(modelo_mco, lag = 6))[, ])## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 49.231695624 2.0229545554 24.336531 3.834732e-61
## TRAFFIC_VOLUME 0.005325644 0.0007920864 6.723564 1.891536e-10
## ROAD_CONDITION -11.450884874 1.0847120351 -10.556613 6.488292e-21
## WEATHER -5.061374460 0.9702954152 -5.216323 4.621912e-07De lo hallado en a: Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.923e+01 1.530e+00 32.175 < 2e-16
TRAFFIC_VOLUME 5.326e-03 6.866e-04 7.756 4.68e-13
ROAD_CONDITION -1.145e+01 1.007e+00 -11.371 < 2e-16
WEATHER -5.061e+00 1.003e+00 -5.045 1.03e-06
Comentarios:
Coeficientes: Los coeficientes estimados son idénticos en ambos métodos porque el ajuste del modelo no cambia
Errores estándar: Los errores estándar ajustados con HAC (Newey-West) son generalmente mayores que los errores estándar originales, lo que refleja la corrección por heterocedasticidad y autocorrelación.
Significancia: Los p-valores son ligeramente mayores con HAC, pero todos los coeficientes siguen siendo altamente significativos
e
library(orcutt)# Aplicar el método de Cochrane-Orcutt
modelo_cochrane_orcutt <- cochrane.orcutt(modelo_mco)
summary(modelo_cochrane_orcutt)## Call:
## lm(formula = CONGESTION ~ TRAFFIC_VOLUME + ROAD_CONDITION + WEATHER,
## data = traffic)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.7531e+01 1.4516e+00 32.745 < 2.2e-16 ***
## TRAFFIC_VOLUME 5.8136e-03 4.8455e-04 11.998 < 2.2e-16 ***
## ROAD_CONDITION -1.0307e+01 6.9163e-01 -14.902 < 2.2e-16 ***
## WEATHER -4.9582e+00 6.8538e-01 -7.234 1.045e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.6304 on 195 degrees of freedom
## Multiple R-squared: 0.6999 , Adjusted R-squared: 0.6953
## F-statistic: 151.6 on 3 and 195 DF, p-value: < 9.976e-51
##
## Durbin-Watson statistic
## (original): 0.79026 , p-value: 3.591e-18
## (transformed): 2.30137 , p-value: 9.869e-01De lo hallado en a: Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.923e+01 1.530e+00 32.175 < 2e-16
TRAFFIC_VOLUME 5.326e-03 6.866e-04 7.756 4.68e-13
ROAD_CONDITION -1.145e+01 1.007e+00 -11.371 < 2e-16
WEATHER -5.061e+00 1.003e+00 -5.045 1.03e-06
De lo hallado en d: Estimate Std. Error t value Pr(>|t|) (Intercept) 49.231695624 2.0229545554 24.336531 3.834732e-61 TRAFFIC_VOLUME 0.005325644 0.0007920864 6.723564 1.891536e-10 ROAD_CONDITION -11.450884874 1.0847120351 -10.556613 6.488292e-21 WEATHER -5.061374460 0.9702954152 -5.216323 4.621912e-07
Comentarios: Los coeficientes estimados son muy similares entre los métodos, pero hay ligeras diferencias en los valores debido a la corrección por autocorrelación con Cochrane-Orcutt
Los errores estándar ajustados con HAC (Newey-West) son generalmente mayores que los errores estándar originales (MCO), lo que refleja la corrección por heterocedasticidad y autocorrelación.
Con Cochrane-Orcutt, los errores estándar son generalmente menores que los obtenidos con HAC, indicando que la corrección específica por autocorrelación de primer orden puede ser más efectiva en este caso
Los coeficientes siguen siendo altamente significativos en todos los métodos, pero los valores t son mayores con Cochrane-Orcutt, lo que sugiere que la corrección por autocorrelación de primer orden mejora la precisión de las estimaciones
Pregunta 5
a
library(tidyr)
library(leaps)library(ggplot2)
library(ggvis)##
## Attaching package: 'ggvis'## The following object is masked from 'package:ggplot2':
##
## resolutioncat_data = read.csv("cat_data.csv")
# Número de variables sin la predictora
nvars = dim(cat_data)[2]-1;nvars## [1] 17# Aplicando el método de mejores subconjuntos a los datos
regfit.full <- regsubsets(birds ~ .,data = cat_data, nvmax = nvars)
reg.summary <- summary(regfit.full)
reg.summary## Subset selection object
## Call: regsubsets.formula(birds ~ ., data = cat_data, nvmax = nvars)
## 17 Variables (and intercept)
## Forced in Forced out
## sex FALSE FALSE
## weight FALSE FALSE
## dryfood FALSE FALSE
## wetfood FALSE FALSE
## age FALSE FALSE
## owner.income FALSE FALSE
## daily.playtime FALSE FALSE
## fellow.cats FALSE FALSE
## owner.age FALSE FALSE
## house.area FALSE FALSE
## children.13 FALSE FALSE
## urban FALSE FALSE
## bell FALSE FALSE
## dogs FALSE FALSE
## daily.outdoortime FALSE FALSE
## daily.catnip FALSE FALSE
## neutered FALSE FALSE
## 1 subsets of each size up to 17
## Selection Algorithm: exhaustive
## sex weight dryfood wetfood age owner.income daily.playtime
## 1 ( 1 ) " " " " " " " " " " " " " "
## 2 ( 1 ) " " " " " " " " " " " " "*"
## 3 ( 1 ) " " " " " " " " " " " " "*"
## 4 ( 1 ) " " " " " " "*" " " " " "*"
## 5 ( 1 ) " " " " " " "*" " " " " "*"
## 6 ( 1 ) " " " " " " "*" " " " " "*"
## 7 ( 1 ) " " "*" " " "*" " " " " "*"
## 8 ( 1 ) " " "*" " " "*" " " " " "*"
## 9 ( 1 ) " " "*" " " "*" " " " " "*"
## 10 ( 1 ) " " "*" " " "*" " " " " "*"
## 11 ( 1 ) " " "*" "*" "*" " " " " "*"
## 12 ( 1 ) " " "*" "*" "*" "*" " " "*"
## 13 ( 1 ) "*" "*" "*" "*" "*" " " "*"
## 14 ( 1 ) "*" "*" "*" "*" "*" " " "*"
## 15 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## 16 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## 17 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## fellow.cats owner.age house.area children.13 urban bell dogs
## 1 ( 1 ) " " " " " " " " " " " " " "
## 2 ( 1 ) " " " " " " " " " " " " " "
## 3 ( 1 ) " " " " " " " " " " "*" " "
## 4 ( 1 ) " " " " " " " " " " "*" " "
## 5 ( 1 ) " " " " " " " " "*" "*" " "
## 6 ( 1 ) " " " " " " "*" "*" "*" " "
## 7 ( 1 ) " " " " " " "*" "*" "*" " "
## 8 ( 1 ) " " " " " " "*" "*" "*" " "
## 9 ( 1 ) " " " " "*" "*" "*" "*" " "
## 10 ( 1 ) "*" " " " " "*" "*" "*" "*"
## 11 ( 1 ) "*" " " " " "*" "*" "*" "*"
## 12 ( 1 ) "*" " " " " "*" "*" "*" "*"
## 13 ( 1 ) "*" " " " " "*" "*" "*" "*"
## 14 ( 1 ) "*" " " " " "*" "*" "*" "*"
## 15 ( 1 ) "*" " " " " "*" "*" "*" "*"
## 16 ( 1 ) "*" " " "*" "*" "*" "*" "*"
## 17 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## daily.outdoortime daily.catnip neutered
## 1 ( 1 ) "*" " " " "
## 2 ( 1 ) "*" " " " "
## 3 ( 1 ) "*" " " " "
## 4 ( 1 ) "*" " " " "
## 5 ( 1 ) "*" " " " "
## 6 ( 1 ) "*" " " " "
## 7 ( 1 ) "*" " " " "
## 8 ( 1 ) "*" " " "*"
## 9 ( 1 ) "*" " " "*"
## 10 ( 1 ) "*" " " "*"
## 11 ( 1 ) "*" " " "*"
## 12 ( 1 ) "*" " " "*"
## 13 ( 1 ) "*" " " "*"
## 14 ( 1 ) "*" "*" "*"
## 15 ( 1 ) "*" "*" "*"
## 16 ( 1 ) "*" "*" "*"
## 17 ( 1 ) "*" "*" "*"# Aplicando el método de mejores subconjuntos a los datos
regfit.full <- regsubsets(birds ~ .,data = cat_data, nvmax = nvars)
reg.summary <- summary(regfit.full)
reg.summary## Subset selection object
## Call: regsubsets.formula(birds ~ ., data = cat_data, nvmax = nvars)
## 17 Variables (and intercept)
## Forced in Forced out
## sex FALSE FALSE
## weight FALSE FALSE
## dryfood FALSE FALSE
## wetfood FALSE FALSE
## age FALSE FALSE
## owner.income FALSE FALSE
## daily.playtime FALSE FALSE
## fellow.cats FALSE FALSE
## owner.age FALSE FALSE
## house.area FALSE FALSE
## children.13 FALSE FALSE
## urban FALSE FALSE
## bell FALSE FALSE
## dogs FALSE FALSE
## daily.outdoortime FALSE FALSE
## daily.catnip FALSE FALSE
## neutered FALSE FALSE
## 1 subsets of each size up to 17
## Selection Algorithm: exhaustive
## sex weight dryfood wetfood age owner.income daily.playtime
## 1 ( 1 ) " " " " " " " " " " " " " "
## 2 ( 1 ) " " " " " " " " " " " " "*"
## 3 ( 1 ) " " " " " " " " " " " " "*"
## 4 ( 1 ) " " " " " " "*" " " " " "*"
## 5 ( 1 ) " " " " " " "*" " " " " "*"
## 6 ( 1 ) " " " " " " "*" " " " " "*"
## 7 ( 1 ) " " "*" " " "*" " " " " "*"
## 8 ( 1 ) " " "*" " " "*" " " " " "*"
## 9 ( 1 ) " " "*" " " "*" " " " " "*"
## 10 ( 1 ) " " "*" " " "*" " " " " "*"
## 11 ( 1 ) " " "*" "*" "*" " " " " "*"
## 12 ( 1 ) " " "*" "*" "*" "*" " " "*"
## 13 ( 1 ) "*" "*" "*" "*" "*" " " "*"
## 14 ( 1 ) "*" "*" "*" "*" "*" " " "*"
## 15 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## 16 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## 17 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## fellow.cats owner.age house.area children.13 urban bell dogs
## 1 ( 1 ) " " " " " " " " " " " " " "
## 2 ( 1 ) " " " " " " " " " " " " " "
## 3 ( 1 ) " " " " " " " " " " "*" " "
## 4 ( 1 ) " " " " " " " " " " "*" " "
## 5 ( 1 ) " " " " " " " " "*" "*" " "
## 6 ( 1 ) " " " " " " "*" "*" "*" " "
## 7 ( 1 ) " " " " " " "*" "*" "*" " "
## 8 ( 1 ) " " " " " " "*" "*" "*" " "
## 9 ( 1 ) " " " " "*" "*" "*" "*" " "
## 10 ( 1 ) "*" " " " " "*" "*" "*" "*"
## 11 ( 1 ) "*" " " " " "*" "*" "*" "*"
## 12 ( 1 ) "*" " " " " "*" "*" "*" "*"
## 13 ( 1 ) "*" " " " " "*" "*" "*" "*"
## 14 ( 1 ) "*" " " " " "*" "*" "*" "*"
## 15 ( 1 ) "*" " " " " "*" "*" "*" "*"
## 16 ( 1 ) "*" " " "*" "*" "*" "*" "*"
## 17 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## daily.outdoortime daily.catnip neutered
## 1 ( 1 ) "*" " " " "
## 2 ( 1 ) "*" " " " "
## 3 ( 1 ) "*" " " " "
## 4 ( 1 ) "*" " " " "
## 5 ( 1 ) "*" " " " "
## 6 ( 1 ) "*" " " " "
## 7 ( 1 ) "*" " " " "
## 8 ( 1 ) "*" " " "*"
## 9 ( 1 ) "*" " " "*"
## 10 ( 1 ) "*" " " "*"
## 11 ( 1 ) "*" " " "*"
## 12 ( 1 ) "*" " " "*"
## 13 ( 1 ) "*" " " "*"
## 14 ( 1 ) "*" "*" "*"
## 15 ( 1 ) "*" "*" "*"
## 16 ( 1 ) "*" "*" "*"
## 17 ( 1 ) "*" "*" "*"rsq <- as.data.frame(reg.summary$rsq)
names(rsq) <- "R2"
rsq %>%
ggvis(x=~ c(1:nrow(rsq)), y=~R2 ) %>%
layer_points(fill = ~ R2 ) %>%
add_axis("y", title = "R2") %>%
add_axis("x", title = "Número de Variables")# Grafica comparativa de la RSS, el R2 ajustado, Cp y BIC
par(mfrow = c(2,2))
plot(reg.summary$rss, xlab="Numero de Variables", ylab="RSS", type="l")
plot(reg.summary$adjr2, xlab="Numero de Variables", ylab="R2 Ajustado", type="l")
np.r2 <- which.max(reg.summary$adjr2) #valor optimo para este indicador
points(np.r2, reg.summary$adjr2[np.r2], col = "red", cex = 2, pch = 20)
plot(reg.summary$cp, xlab="Numero de Variables", ylab="Cp", type='l')
np.cp <- which.min(reg.summary$cp)
points(np.cp,reg.summary$cp[np.cp],col="red",cex=2,pch=20)
plot(reg.summary$bic, xlab="Numero de Variables", ylab = "BIC", type = 'l')
np.bic <- which.min(reg.summary$bic)
points(np.bic, reg.summary$bic[np.bic], col = "red", cex = 2, pch = 20)
# Las variables seleccionadas para el mejor modelo con un número de predictores
# R2
plot(regfit.full,scale="r2")
# R2 ajustado
plot(regfit.full,scale="adjr2")
# Cp de Mallow
plot(regfit.full,scale="Cp")
# BIC
plot(regfit.full,scale="bic")
# Coeficientes asociados al mejor modelo según el BIC
coefi <- coef(regfit.full, np.bic);coefi## (Intercept) weight wetfood daily.playtime
## 118.542381 -3.713249 -7.338219 -10.266126
## children.13 urban bell daily.outdoortime
## 5.247393 -12.287953 -51.449174 1.123477
## neutered
## -7.668978b
# Selección paso a paso hacia adelante
regfit.fwd = regsubsets(birds ~., data = cat_data, nvmax=nvars,method="forward")
summary(regfit.fwd)## Subset selection object
## Call: regsubsets.formula(birds ~ ., data = cat_data, nvmax = nvars,
## method = "forward")
## 17 Variables (and intercept)
## Forced in Forced out
## sex FALSE FALSE
## weight FALSE FALSE
## dryfood FALSE FALSE
## wetfood FALSE FALSE
## age FALSE FALSE
## owner.income FALSE FALSE
## daily.playtime FALSE FALSE
## fellow.cats FALSE FALSE
## owner.age FALSE FALSE
## house.area FALSE FALSE
## children.13 FALSE FALSE
## urban FALSE FALSE
## bell FALSE FALSE
## dogs FALSE FALSE
## daily.outdoortime FALSE FALSE
## daily.catnip FALSE FALSE
## neutered FALSE FALSE
## 1 subsets of each size up to 17
## Selection Algorithm: forward
## sex weight dryfood wetfood age owner.income daily.playtime
## 1 ( 1 ) " " " " " " " " " " " " " "
## 2 ( 1 ) " " " " " " " " " " " " "*"
## 3 ( 1 ) " " " " " " " " " " " " "*"
## 4 ( 1 ) " " " " " " "*" " " " " "*"
## 5 ( 1 ) " " " " " " "*" " " " " "*"
## 6 ( 1 ) " " " " " " "*" " " " " "*"
## 7 ( 1 ) " " "*" " " "*" " " " " "*"
## 8 ( 1 ) " " "*" " " "*" " " " " "*"
## 9 ( 1 ) " " "*" " " "*" " " " " "*"
## 10 ( 1 ) " " "*" "*" "*" " " " " "*"
## 11 ( 1 ) " " "*" "*" "*" "*" " " "*"
## 12 ( 1 ) " " "*" "*" "*" "*" " " "*"
## 13 ( 1 ) " " "*" "*" "*" "*" " " "*"
## 14 ( 1 ) "*" "*" "*" "*" "*" " " "*"
## 15 ( 1 ) "*" "*" "*" "*" "*" " " "*"
## 16 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## 17 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## fellow.cats owner.age house.area children.13 urban bell dogs
## 1 ( 1 ) " " " " " " " " " " " " " "
## 2 ( 1 ) " " " " " " " " " " " " " "
## 3 ( 1 ) " " " " " " " " " " "*" " "
## 4 ( 1 ) " " " " " " " " " " "*" " "
## 5 ( 1 ) " " " " " " " " "*" "*" " "
## 6 ( 1 ) " " " " " " "*" "*" "*" " "
## 7 ( 1 ) " " " " " " "*" "*" "*" " "
## 8 ( 1 ) " " " " " " "*" "*" "*" " "
## 9 ( 1 ) " " " " "*" "*" "*" "*" " "
## 10 ( 1 ) " " " " "*" "*" "*" "*" " "
## 11 ( 1 ) " " " " "*" "*" "*" "*" " "
## 12 ( 1 ) " " " " "*" "*" "*" "*" "*"
## 13 ( 1 ) "*" " " "*" "*" "*" "*" "*"
## 14 ( 1 ) "*" " " "*" "*" "*" "*" "*"
## 15 ( 1 ) "*" " " "*" "*" "*" "*" "*"
## 16 ( 1 ) "*" " " "*" "*" "*" "*" "*"
## 17 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## daily.outdoortime daily.catnip neutered
## 1 ( 1 ) "*" " " " "
## 2 ( 1 ) "*" " " " "
## 3 ( 1 ) "*" " " " "
## 4 ( 1 ) "*" " " " "
## 5 ( 1 ) "*" " " " "
## 6 ( 1 ) "*" " " " "
## 7 ( 1 ) "*" " " " "
## 8 ( 1 ) "*" " " "*"
## 9 ( 1 ) "*" " " "*"
## 10 ( 1 ) "*" " " "*"
## 11 ( 1 ) "*" " " "*"
## 12 ( 1 ) "*" " " "*"
## 13 ( 1 ) "*" " " "*"
## 14 ( 1 ) "*" " " "*"
## 15 ( 1 ) "*" "*" "*"
## 16 ( 1 ) "*" "*" "*"
## 17 ( 1 ) "*" "*" "*"# Mejor modelo según el Cp
plot(regfit.fwd,scale="Cp")
# Selección paso a paso hacia atrás
regfit.bwd = regsubsets(birds~.,data= cat_data, nvmax=nvars,method = "backward")
summary(regfit.bwd)## Subset selection object
## Call: regsubsets.formula(birds ~ ., data = cat_data, nvmax = nvars,
## method = "backward")
## 17 Variables (and intercept)
## Forced in Forced out
## sex FALSE FALSE
## weight FALSE FALSE
## dryfood FALSE FALSE
## wetfood FALSE FALSE
## age FALSE FALSE
## owner.income FALSE FALSE
## daily.playtime FALSE FALSE
## fellow.cats FALSE FALSE
## owner.age FALSE FALSE
## house.area FALSE FALSE
## children.13 FALSE FALSE
## urban FALSE FALSE
## bell FALSE FALSE
## dogs FALSE FALSE
## daily.outdoortime FALSE FALSE
## daily.catnip FALSE FALSE
## neutered FALSE FALSE
## 1 subsets of each size up to 17
## Selection Algorithm: backward
## sex weight dryfood wetfood age owner.income daily.playtime
## 1 ( 1 ) " " " " " " " " " " " " " "
## 2 ( 1 ) " " " " " " " " " " " " "*"
## 3 ( 1 ) " " " " " " " " " " " " "*"
## 4 ( 1 ) " " " " " " "*" " " " " "*"
## 5 ( 1 ) " " " " " " "*" " " " " "*"
## 6 ( 1 ) " " " " " " "*" " " " " "*"
## 7 ( 1 ) " " "*" " " "*" " " " " "*"
## 8 ( 1 ) " " "*" " " "*" " " " " "*"
## 9 ( 1 ) " " "*" " " "*" " " " " "*"
## 10 ( 1 ) " " "*" " " "*" " " " " "*"
## 11 ( 1 ) " " "*" "*" "*" " " " " "*"
## 12 ( 1 ) " " "*" "*" "*" "*" " " "*"
## 13 ( 1 ) "*" "*" "*" "*" "*" " " "*"
## 14 ( 1 ) "*" "*" "*" "*" "*" " " "*"
## 15 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## 16 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## 17 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## fellow.cats owner.age house.area children.13 urban bell dogs
## 1 ( 1 ) " " " " " " " " " " " " " "
## 2 ( 1 ) " " " " " " " " " " " " " "
## 3 ( 1 ) " " " " " " " " " " "*" " "
## 4 ( 1 ) " " " " " " " " " " "*" " "
## 5 ( 1 ) " " " " " " " " "*" "*" " "
## 6 ( 1 ) " " " " " " "*" "*" "*" " "
## 7 ( 1 ) " " " " " " "*" "*" "*" " "
## 8 ( 1 ) " " " " " " "*" "*" "*" " "
## 9 ( 1 ) " " " " " " "*" "*" "*" "*"
## 10 ( 1 ) "*" " " " " "*" "*" "*" "*"
## 11 ( 1 ) "*" " " " " "*" "*" "*" "*"
## 12 ( 1 ) "*" " " " " "*" "*" "*" "*"
## 13 ( 1 ) "*" " " " " "*" "*" "*" "*"
## 14 ( 1 ) "*" " " " " "*" "*" "*" "*"
## 15 ( 1 ) "*" " " " " "*" "*" "*" "*"
## 16 ( 1 ) "*" " " "*" "*" "*" "*" "*"
## 17 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## daily.outdoortime daily.catnip neutered
## 1 ( 1 ) "*" " " " "
## 2 ( 1 ) "*" " " " "
## 3 ( 1 ) "*" " " " "
## 4 ( 1 ) "*" " " " "
## 5 ( 1 ) "*" " " " "
## 6 ( 1 ) "*" " " " "
## 7 ( 1 ) "*" " " " "
## 8 ( 1 ) "*" " " "*"
## 9 ( 1 ) "*" " " "*"
## 10 ( 1 ) "*" " " "*"
## 11 ( 1 ) "*" " " "*"
## 12 ( 1 ) "*" " " "*"
## 13 ( 1 ) "*" " " "*"
## 14 ( 1 ) "*" "*" "*"
## 15 ( 1 ) "*" "*" "*"
## 16 ( 1 ) "*" "*" "*"
## 17 ( 1 ) "*" "*" "*"plot(regfit.bwd,scale="Cp")
Dado los resultados se llega a los mismo modelo elegido que en la parte a)