Taller 1: Modelos de Regresión Lineal (Econometría I)
[Juan Esteban Reyes Rodríguez]
Octubre 2025
1 1. Planteamiento y Carga de Datos
Se utilizará la base de datos crime1 del paquete
wooldridge. La variable dependiente seleccionada es
narr86 (número de arrestos en 1986) y se buscará explicarla
en función de variables de disuasión (p.ej., pcnv,
avgsen) y variables socioeconómicas (p.ej.,
qemp86, inc86).
1.1 1.1. Carga de Paquetes y Datos
# Fijar semilla para reproducibilidad
set.seed(123)
# Carga de paquetes
library(wooldridge)
library(tidyverse)
library(broom)
library(lmtest)
library(sandwich)
library(car)
library(ggplot2)
library(tseries)
# Cargar datos
data("crime1")
df <- as_tibble(crime1)1.2 1.2. Limpieza y Preparación de Variables
# Filtrar NAs en variables clave y crear transformaciones
df <- df %>%
filter(!is.na(narr86), !is.na(pcnv), !is.na(avgsen), !is.na(inc86)) %>%
mutate(
# Transformación de ingresos para interpretación ($100s)
inc86_100 = inc86 / 100,
# Transformación logarítmica de Y (con ajuste de +0.1 para ceros)
ln_narr86 = log(narr86 + 0.1)
)
glimpse(df)## Rows: 2,725
## Columns: 18
## $ narr86 <int> 0, 2, 1, 2, 1, 0, 2, 5, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, …
## $ nfarr86 <int> 0, 2, 1, 2, 1, 0, 2, 3, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, …
## $ nparr86 <int> 0, 0, 0, 1, 0, 0, 1, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pcnv <dbl> 0.38, 0.44, 0.33, 0.25, 0.00, 1.00, 0.44, 0.75, 0.33, 0.23, …
## $ avgsen <dbl> 17.6, 0.0, 22.8, 0.0, 0.0, 0.0, 0.0, 0.0, 10.9, 0.0, 0.0, 31…
## $ tottime <dbl> 35.2, 0.0, 22.8, 0.0, 0.0, 0.0, 0.0, 0.0, 21.8, 0.0, 0.0, 63…
## $ ptime86 <int> 12, 0, 0, 5, 0, 0, 0, 0, 9, 0, 0, 12, 0, 0, 0, 3, 0, 0, 0, 0…
## $ qemp86 <dbl> 0.0, 1.0, 0.0, 2.0, 2.0, 4.0, 0.0, 0.0, 0.0, 3.0, 4.0, 0.0, …
## $ inc86 <dbl> 0.0, 0.8, 0.0, 8.8, 8.1, 97.6, 0.0, 0.0, 0.0, 16.7, 162.5, 0…
## $ durat <dbl> 0.0, 0.0, 11.0, 0.0, 1.0, 0.0, 1.0, 3.0, 19.3, 0.0, 0.0, 0.0…
## $ black <int> 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, …
## $ hispan <int> 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ born60 <int> 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, …
## $ pcnvsq <dbl> 0.1444, 0.1936, 0.1089, 0.0625, 0.0000, 1.0000, 0.1936, 0.56…
## $ pt86sq <int> 144, 0, 0, 25, 0, 0, 0, 0, 81, 0, 0, 144, 0, 0, 0, 9, 0, 0, …
## $ inc86sq <dbl> 0.00000, 0.64000, 0.00000, 77.44000, 65.61001, 9525.75977, 0…
## $ inc86_100 <dbl> 0.000, 0.008, 0.000, 0.088, 0.081, 0.976, 0.000, 0.000, 0.00…
## $ ln_narr86 <dbl> -2.30258509, 0.74193734, 0.09531018, 0.74193734, 0.09531018,…2 2. Exploración y Modelo Base
2.1 2.1. Estadísticos Descriptivos
df %>%
summarize(
Observaciones = n(),
Media_Arrestos = round(mean(narr86), 3),
DE_Arrestos = round(sd(narr86), 3),
Media_ProbCondena = round(mean(pcnv), 3),
Media_Sentencia = round(mean(avgsen), 3)
)## # A tibble: 1 × 5
## Observaciones Media_Arrestos DE_Arrestos Media_ProbCondena Media_Sentencia
## <int> <dbl> <dbl> <dbl> <dbl>
## 1 2725 0.404 0.859 0.358 0.6322.2 2.2. Estimación Modelo Base (M0)
[cite_start]Se estima el modelo lineal por MCO[cite: 25].
m0 <- lm(narr86 ~ pcnv + avgsen + tottime + qemp86 + inc86_100 +
durat + black + hispan + born60, data = df)
# [cite_start]Reporte de coeficientes, SE, R2, etc. [cite: 26]
summary(m0)##
## Call:
## lm(formula = narr86 ~ pcnv + avgsen + tottime + qemp86 + inc86_100 +
## durat + black + hispan + born60, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.8859 -0.4569 -0.2549 0.2895 11.5419
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.5469505 0.0434979 12.574 < 2e-16 ***
## pcnv -0.1416106 0.0405248 -3.494 0.000483 ***
## avgsen -0.0072566 0.0122603 -0.592 0.553980
## tottime 0.0048661 0.0093500 0.520 0.602796
## qemp86 -0.0382997 0.0152748 -2.507 0.012221 *
## inc86_100 -0.1544451 0.0343909 -4.491 7.39e-06 ***
## durat -0.0006318 0.0039640 -0.159 0.873373
## black 0.3204411 0.0455846 7.030 2.61e-12 ***
## hispan 0.1819218 0.0398049 4.570 5.09e-06 ***
## born60 -0.0216075 0.0334261 -0.646 0.518058
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.832 on 2715 degrees of freedom
## Multiple R-squared: 0.06514, Adjusted R-squared: 0.06204
## F-statistic: 21.02 on 9 and 2715 DF, p-value: < 2.2e-163 3. Validación de Supuestos (Modelo Base)
[cite_start]Se realizan los diagnósticos clásicos sobre el modelo
m0[cite: 6].
3.1 3.1. Linealidad / Especificación (Prueba RESET)
# Prueba RESET [cite: 28, 44]
resettest(m0, power = 2:3, type = "fitted")##
## RESET test
##
## data: m0
## RESET = 5.9206, df1 = 2, df2 = 2713, p-value = 0.0027183.2 3.2. Homocedasticidad (Prueba BP y Gráfico)
# Prueba de Breusch-Pagan [cite: 30, 45]
bptest(m0)##
## studentized Breusch-Pagan test
##
## data: m0
## BP = 53.42, df = 9, p-value = 2.433e-08# Gráfico Residuos vs Ajustados [cite: 30, 45]
ggplot(data = NULL, aes(x = fitted(m0), y = resid(m0))) +
geom_point(alpha = 0.5) +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
labs(title = "Residuos vs Valores Ajustados - Modelo Base",
x = "Valores Ajustados", y = "Residuos") +
theme_minimal()
3.3 3.3. Normalidad de Residuos (Jarque-Bera y QQ-Plot)
# Prueba Jarque-Bera [cite: 31, 46]
tseries::jarque.bera.test(resid(m0))##
## Jarque Bera Test
##
## data: resid(m0)
## X-squared = 111059, df = 2, p-value < 2.2e-16# Gráfico QQ-plot [cite: 31, 46]
qqnorm(resid(m0), main = "QQ-Plot de Residuos - Modelo Base")
qqline(resid(m0), col = "red")
3.4 3.4. Multicolinealidad (VIF)
# Factores de Inflación de Varianza (VIF) [cite: 32, 47]
vif(m0)## pcnv avgsen tottime qemp86 inc86_100 durat black hispan
## 1.009288 7.279196 7.301600 2.381162 2.066087 1.312411 1.105507 1.061937
## born60
## 1.0165124 4. Soluciones Econométricas
Los diagnósticos sugieren problemas de especificación (RESET) y heterocedasticidad (Breusch-Pagan). [cite_start]Se probarán re-especificaciones[cite: 34].
4.1 4.1. Modelos Alternativos
[cite_start]Se centran variables para la interacción y se prueba un modelo log-lineal (dado que Y >= 0)[cite: 37, 56].
# Centrar variables para interacción
df <- df %>%
mutate(
pcnv_c = pcnv - mean(pcnv),
avgsen_c = avgsen - mean(avgsen)
)
# Modelo con interacción
m1 <- lm(narr86 ~ pcnv_c + avgsen_c + pcnv_c:avgsen_c +
qemp86 + inc86_100 + durat + black + hispan + born60, data = df)
# Modelo log-lineal (M2)
m2 <- lm(ln_narr86 ~ pcnv_c + avgsen_c + pcnv_c:avgsen_c +
qemp86 + inc86_100 + durat + black + hispan + born60, data = df)4.2 4.2. Comparación y Selección de Modelo
comparacion <- tibble(
Modelo = c("Base (m0)", "Interacción (m1)", "Log-lineal (m2)"),
R2_ajustado = c(summary(m0)$adj.r.squared,
summary(m1)$adj.r.squared,
summary(m2)$adj.r.squared),
RESET_p = c(resettest(m0)$p.value,
resettest(m1)$p.value,
resettest(m2, power=2:3, type="fitted")$p.value)
)
# El modelo M2 (log-lineal) muestra el mejor p-valor para RESET
comparacion## # A tibble: 3 × 3
## Modelo R2_ajustado RESET_p
## <chr> <dbl> <dbl>
## 1 Base (m0) 0.0620 0.00272
## 2 Interacción (m1) 0.0623 0.00275
## 3 Log-lineal (m2) 0.0764 0.002265 5. Modelo Final y Corrección de Heterocedasticidad
Se selecciona m2 como modelo final. Se vuelve a revisar
la heterocedasticidad.
final_model <- m2
bptest(final_model)##
## studentized Breusch-Pagan test
##
## data: final_model
## BP = 222.34, df = 9, p-value < 2.2e-16El p-valor de Breusch-Pagan sigue siendo bajo, indicando que la heterocedasticidad persiste. [cite_start]Se aplicará la corrección de Errores Estándar Robustos (Eicker-Huber-White) para una inferencia válida[cite: 35, 48].
# Coeficientes con errores clásicos (para comparar)
tidy(final_model)## # A tibble: 10 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -1.58 0.0560 -28.2 3.22e-153
## 2 pcnv_c -0.368 0.0566 -6.51 8.89e- 11
## 3 avgsen_c -0.00266 0.00637 -0.417 6.76e- 1
## 4 qemp86 0.00774 0.0211 0.367 7.13e- 1
## 5 inc86_100 -0.308 0.0475 -6.47 1.14e- 10
## 6 durat 0.0122 0.00548 2.23 2.56e- 2
## 7 black 0.449 0.0631 7.12 1.39e- 12
## 8 hispan 0.252 0.0550 4.59 4.71e- 6
## 9 born60 -0.00336 0.0462 -0.0726 9.42e- 1
## 10 pcnv_c:avgsen_c 0.0142 0.0238 0.596 5.51e- 1# Coeficientes con Errores Estándar Robustos (HC1)
coef_robustos_test <- coeftest(final_model, vcov = vcovHC(final_model, type = "HC1"))
# Tabla final con IC del 95% robustos
tabla_final_robusta <- tidy(coef_robustos_test, conf.int = TRUE)
tabla_final_robusta## # A tibble: 10 × 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -1.58 0.0584 -27.0 1.85e-142 -1.69 -1.46
## 2 pcnv_c -0.368 0.0517 -7.12 1.37e- 12 -0.470 -0.267
## 3 avgsen_c -0.00266 0.00705 -0.378 7.06e- 1 -0.0165 0.0112
## 4 qemp86 0.00774 0.0209 0.370 7.12e- 1 -0.0333 0.0488
## 5 inc86_100 -0.308 0.0406 -7.57 4.95e- 14 -0.387 -0.228
## 6 durat 0.0122 0.00600 2.04 4.16e- 2 0.000466 0.0240
## 7 black 0.449 0.0708 6.34 2.61e- 10 0.310 0.588
## 8 hispan 0.252 0.0572 4.41 1.06e- 5 0.140 0.365
## 9 born60 -0.00336 0.0458 -0.0732 9.42e- 1 -0.0933 0.0865
## 10 pcnv_c:avgsen_c 0.0142 0.0315 0.450 6.53e- 1 -0.0477 0.07606 6. Interpretación Económica (Modelo Final)
6.1 6.1. Efectos Marginales de la Interacción
Se calcula el efecto marginal de pcnv (probabilidad de
condena) sobre narr86 para diferentes niveles de
avgsen (sentencia).
b <- coef(final_model)
mu_avgsen <- mean(df$avgsen) # Media de la sentencia original
# Calcular efecto para distintos niveles de sentencia (avgsen)
efectos <- tibble(
Sentencia_meses = c(0, round(mu_avgsen, 1), 12),
# Efecto marginal sobre ln(Y)
Efecto_marginal = b["pcnv_c"] + b["pcnv_c:avgsen_c"] * (Sentencia_meses - mu_avgsen),
# Cambio porcentual en Y por 10 puntos (0.1) de aumento en pcnv
Cambio_porcentual_10pts = round((exp(Efecto_marginal * 0.1) - 1) * 100, 2)
)
efectos## # A tibble: 3 × 3
## Sentencia_meses Efecto_marginal Cambio_porcentual_10pts
## <dbl> <dbl> <dbl>
## 1 0 -0.377 -3.7
## 2 0.6 -0.369 -3.62
## 3 12 -0.207 -2.056.2 6.2. Gráfico de Interacción
Visualización del efecto de pcnv sobre los arrestos
esperados (narr86) en diferentes niveles de sentencia
(avgsen).
# Preparar datos para predicción
niveles_avgsen <- c(0, 6, 12)
pcnv_range <- seq(min(df$pcnv), max(df$pcnv), length.out = 50)
# Crear una grilla con valores promedio para otras variables
grid_data <- expand.grid(
pcnv_c = pcnv_range - mean(df$pcnv),
avgsen_c = niveles_avgsen - mean(df$avgsen),
qemp86 = mean(df$qemp86),
inc86_100 = mean(df$inc86_100),
durat = mean(df$durat),
black = 1, # Fijo para un grupo
hispan = 0,
born60 = 1
) %>%
mutate(
pcnv = pcnv_c + mean(df$pcnv),
avgsen = avgsen_c + mean(df$avgsen),
Grupo = paste("Sentencia =", avgsen, "meses")
)
# Predecir ln(narr86) y convertir a narr86
pred <- predict(final_model, newdata = grid_data)
grid_data$narr86_pred <- exp(pred) - 0.1
# Graficar
ggplot(grid_data, aes(x = pcnv, y = narr86_pred, color = Grupo)) +
geom_line(linewidth = 1) +
labs(title = "Interacción: Probabilidad de Condena vs Arrestos",
subtitle = "Por nivel de severidad de sentencia (Modelo Log-Lineal)",
x = "Probabilidad de Condena (pcnv)",
y = "Número Esperado de Arrestos (narr86)",
color = "Sentencia") +
theme_minimal() +
theme(legend.position = "bottom")
7 7. Exportación de Resultados
Se exportan las tablas principales requeridas para los entregables[cite: 16].
# Tabla de coeficientes robustos
write.csv(tabla_final_robusta, "tabla_coeficientes_final.csv", row.names = FALSE)
# Tabla de efectos marginales
write.csv(efectos, "tabla_efectos_marginales.csv", row.names = FALSE)