Tarea 3 RLM
Angel Acevedo Pichardo
2026-09-16
#Primero importamos la data necesaria y verificamos la estructura
datos <- read.csv("data_rlm.csv")
str(datos)## 'data.frame': 78 obs. of 7 variables:
## $ municipio : chr "Adjuntas" "Aguada" "Aguadilla" "Aguas Buenas" ...
## $ ing_pc : int 10225 12755 7709 9413 17074 13032 17144 19549 13987 4670 ...
## $ bach_pct : num 15.6 21 21.3 20 28.4 ...
## $ desem_pct : num 9.89 9.25 17.13 9.98 13.04 ...
## $ banda_pct : num 25.7 34.5 28.4 29.7 35.8 ...
## $ pobreza_pct: num 28.9 24.1 30.3 31.5 10.1 ...
## $ tam_hogar : num 3.63 3.02 3.03 3.01 2.63 2.68 3.29 3.08 2.8 3.11 ...#Luego, exploramos los datos y observamos un mapa de calor de correlaciones
library(corrplot)
R <- cor(datos[, c("ing_pc", "bach_pct", "desem_pct", "banda_pct", "pobreza_pct", "tam_hogar")],
use = "pairwise.complete.obs")
corrplot(R,
method = "color",
type = "lower",
addCoef.col = "black",
tl.col = "black",
tl.srt = 45,
diag = FALSE)
#Ahora, creamos el modelo saturado con todas las variables disponibles y buscamos la mejor combinacion de variables.
saturado <- ing_pc ~ bach_pct + desem_pct + banda_pct + pobreza_pct + tam_hogar
reg <- lm(saturado, data = datos)
summary(reg)##
## Call:
## lm(formula = saturado, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4585.2 -1168.6 -125.4 1208.2 5224.9
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3757.23 3338.95 1.125 0.264211
## bach_pct 264.14 68.73 3.843 0.000259 ***
## desem_pct 20.05 74.85 0.268 0.789533
## banda_pct 125.32 45.30 2.766 0.007198 **
## pobreza_pct -210.64 47.68 -4.417 3.45e-05 ***
## tam_hogar 816.67 622.90 1.311 0.193997
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1982 on 72 degrees of freedom
## Multiple R-squared: 0.8457, Adjusted R-squared: 0.835
## F-statistic: 78.93 on 5 and 72 DF, p-value: < 2.2e-16library(leaps)
ajuste <- regsubsets(
saturado, data = datos,
nvmax = 5,
method = "exhaustive")
s <- summary(ajuste)
names(s)## [1] "which" "rsq" "rss" "adjr2" "cp" "bic" "outmat" "obj"#Entonces, pasamos a crear la formula del modelo, junto al modelo nulo y modelo saturado para encontrar la mejor combinacion de variables, comenzando por el forward stepwise selection.
form <- ing_pc ~ bach_pct + desem_pct + banda_pct + pobreza_pct + tam_hogar
m0 <- lm(ing_pc ~ 1, data = datos)
mS <-lm(saturado, data = datos)
m1_step <- step(m0,
scope = list(lower= ~1, upper = saturado),
direction = "forward",
k = log(nrow(datos)),
trace = TRUE)## Start: AIC=1328.25
## ing_pc ~ 1
##
## Df Sum of Sq RSS AIC
## + bach_pct 1 1438089133 395800604 1213.0
## + pobreza_pct 1 1331610382 502279355 1231.6
## + banda_pct 1 1121735582 712154155 1258.8
## <none> 1833889736 1328.2
## + desem_pct 1 8354792 1825534944 1332.2
## + tam_hogar 1 5821235 1828068502 1332.4
##
## Step: AIC=1213.01
## ing_pc ~ bach_pct
##
## Df Sum of Sq RSS AIC
## + pobreza_pct 1 77610119 318190485 1200.3
## + banda_pct 1 29535087 366265517 1211.3
## <none> 395800604 1213.0
## + tam_hogar 1 4923053 390877551 1216.4
## + desem_pct 1 2817657 392982947 1216.8
##
## Step: AIC=1200.34
## ing_pc ~ bach_pct + pobreza_pct
##
## Df Sum of Sq RSS AIC
## + banda_pct 1 28478762 289711723 1197.4
## <none> 318190485 1200.3
## + tam_hogar 1 4915101 313275383 1203.5
## + desem_pct 1 672748 317517737 1204.5
##
## Step: AIC=1197.39
## ing_pc ~ bach_pct + pobreza_pct + banda_pct
##
## Df Sum of Sq RSS AIC
## <none> 289711723 1197.4
## + tam_hogar 1 6489907 283221816 1200.0
## + desem_pct 1 17073 289694650 1201.7#Luego observamos backward stepwise selection
m2_step <- step(mS,
direction = "backward",
k = log(nrow(datos)),
trace = TRUE)## Start: AIC=1204.26
## ing_pc ~ bach_pct + desem_pct + banda_pct + pobreza_pct + tam_hogar
##
## Df Sum of Sq RSS AIC
## - desem_pct 1 282063 283221816 1200.0
## - tam_hogar 1 6754897 289694650 1201.7
## <none> 282939753 1204.3
## - banda_pct 1 30072895 313012647 1207.8
## - bach_pct 1 58048450 340988202 1214.5
## - pobreza_pct 1 76685336 359625088 1218.6
##
## Step: AIC=1199.98
## ing_pc ~ bach_pct + banda_pct + pobreza_pct + tam_hogar
##
## Df Sum of Sq RSS AIC
## - tam_hogar 1 6489907 289711723 1197.4
## <none> 283221816 1200.0
## - banda_pct 1 30053567 313275383 1203.5
## - bach_pct 1 61021640 344243456 1210.8
## - pobreza_pct 1 76514196 359736013 1214.3
##
## Step: AIC=1197.39
## ing_pc ~ bach_pct + banda_pct + pobreza_pct
##
## Df Sum of Sq RSS AIC
## <none> 289711723 1197.4
## - banda_pct 1 28478762 318190485 1200.3
## - bach_pct 1 62624091 352335814 1208.3
## - pobreza_pct 1 76553794 366265517 1211.3#Por ultimo, observamos el metodo híbrido
m3_step <- step(m0,
scope = list(lower= ~1, upper = saturado),
direction = "both",
k = log(nrow(datos)),
trace = TRUE)## Start: AIC=1328.25
## ing_pc ~ 1
##
## Df Sum of Sq RSS AIC
## + bach_pct 1 1438089133 395800604 1213.0
## + pobreza_pct 1 1331610382 502279355 1231.6
## + banda_pct 1 1121735582 712154155 1258.8
## <none> 1833889736 1328.2
## + desem_pct 1 8354792 1825534944 1332.2
## + tam_hogar 1 5821235 1828068502 1332.4
##
## Step: AIC=1213.01
## ing_pc ~ bach_pct
##
## Df Sum of Sq RSS AIC
## + pobreza_pct 1 77610119 318190485 1200.3
## + banda_pct 1 29535087 366265517 1211.3
## <none> 395800604 1213.0
## + tam_hogar 1 4923053 390877551 1216.4
## + desem_pct 1 2817657 392982947 1216.8
## - bach_pct 1 1438089133 1833889736 1328.2
##
## Step: AIC=1200.34
## ing_pc ~ bach_pct + pobreza_pct
##
## Df Sum of Sq RSS AIC
## + banda_pct 1 28478762 289711723 1197.4
## <none> 318190485 1200.3
## + tam_hogar 1 4915101 313275383 1203.5
## + desem_pct 1 672748 317517737 1204.5
## - pobreza_pct 1 77610119 395800604 1213.0
## - bach_pct 1 184088870 502279355 1231.6
##
## Step: AIC=1197.39
## ing_pc ~ bach_pct + pobreza_pct + banda_pct
##
## Df Sum of Sq RSS AIC
## <none> 289711723 1197.4
## + tam_hogar 1 6489907 283221816 1200.0
## - banda_pct 1 28478762 318190485 1200.3
## + desem_pct 1 17073 289694650 1201.7
## - bach_pct 1 62624091 352335814 1208.3
## - pobreza_pct 1 76553794 366265517 1211.3#Observamos los tres metodos
summary(m1_step)##
## Call:
## lm(formula = ing_pc ~ bach_pct + pobreza_pct + banda_pct, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4644.3 -995.6 -35.5 1360.0 5031.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6481.30 2522.94 2.569 0.012215 *
## bach_pct 270.45 67.62 3.999 0.000149 ***
## pobreza_pct -209.43 47.36 -4.422 3.3e-05 ***
## banda_pct 119.54 44.32 2.697 0.008657 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1979 on 74 degrees of freedom
## Multiple R-squared: 0.842, Adjusted R-squared: 0.8356
## F-statistic: 131.5 on 3 and 74 DF, p-value: < 2.2e-16summary(m2_step)##
## Call:
## lm(formula = ing_pc ~ bach_pct + banda_pct + pobreza_pct, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4644.3 -995.6 -35.5 1360.0 5031.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6481.30 2522.94 2.569 0.012215 *
## bach_pct 270.45 67.62 3.999 0.000149 ***
## banda_pct 119.54 44.32 2.697 0.008657 **
## pobreza_pct -209.43 47.36 -4.422 3.3e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1979 on 74 degrees of freedom
## Multiple R-squared: 0.842, Adjusted R-squared: 0.8356
## F-statistic: 131.5 on 3 and 74 DF, p-value: < 2.2e-16summary(m3_step)##
## Call:
## lm(formula = ing_pc ~ bach_pct + pobreza_pct + banda_pct, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4644.3 -995.6 -35.5 1360.0 5031.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6481.30 2522.94 2.569 0.012215 *
## bach_pct 270.45 67.62 3.999 0.000149 ***
## pobreza_pct -209.43 47.36 -4.422 3.3e-05 ***
## banda_pct 119.54 44.32 2.697 0.008657 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1979 on 74 degrees of freedom
## Multiple R-squared: 0.842, Adjusted R-squared: 0.8356
## F-statistic: 131.5 on 3 and 74 DF, p-value: < 2.2e-16#Ahora hacemos el proceso de seleccionar las variables por subconjuntos y seleccionamos el modelo óptimo que minimice BIC o maximice R^2-ajustado.
library(leaps)
reg <- lm(saturado, data = datos)
summary(reg)##
## Call:
## lm(formula = saturado, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4585.2 -1168.6 -125.4 1208.2 5224.9
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3757.23 3338.95 1.125 0.264211
## bach_pct 264.14 68.73 3.843 0.000259 ***
## desem_pct 20.05 74.85 0.268 0.789533
## banda_pct 125.32 45.30 2.766 0.007198 **
## pobreza_pct -210.64 47.68 -4.417 3.45e-05 ***
## tam_hogar 816.67 622.90 1.311 0.193997
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1982 on 72 degrees of freedom
## Multiple R-squared: 0.8457, Adjusted R-squared: 0.835
## F-statistic: 78.93 on 5 and 72 DF, p-value: < 2.2e-16s <- regsubsets(
saturado,
data = datos,
nvmax = 5,
method = "exhaustive"
)
suma <- summary(s)
suma## Subset selection object
## Call: regsubsets.formula(saturado, data = datos, nvmax = 5, method = "exhaustive")
## 5 Variables (and intercept)
## Forced in Forced out
## bach_pct FALSE FALSE
## desem_pct FALSE FALSE
## banda_pct FALSE FALSE
## pobreza_pct FALSE FALSE
## tam_hogar FALSE FALSE
## 1 subsets of each size up to 5
## Selection Algorithm: exhaustive
## bach_pct desem_pct banda_pct pobreza_pct tam_hogar
## 1 ( 1 ) "*" " " " " " " " "
## 2 ( 1 ) "*" " " " " "*" " "
## 3 ( 1 ) "*" " " "*" "*" " "
## 4 ( 1 ) "*" " " "*" "*" "*"
## 5 ( 1 ) "*" "*" "*" "*" "*"best_bic <- which.min(suma$bic)
best_r2 <- which.max(suma$adjr2)
best_bic## [1] 3best_r2## [1] 4suma$which## (Intercept) bach_pct desem_pct banda_pct pobreza_pct tam_hogar
## 1 TRUE TRUE FALSE FALSE FALSE FALSE
## 2 TRUE TRUE FALSE FALSE TRUE FALSE
## 3 TRUE TRUE FALSE TRUE TRUE FALSE
## 4 TRUE TRUE FALSE TRUE TRUE TRUE
## 5 TRUE TRUE TRUE TRUE TRUE TRUEpar(mfrow = c(1, 2))
plot(
suma$bic,
type = "b",
pch = 19,
xlab = "Número de predictores",
ylab = "BIC",
main = "Criterio BIC"
)
points(
best_bic,
suma$bic[best_bic],
pch = 19,
cex = 1.5
)
plot(
suma$adjr2,
type = "b",
pch = 19,
xlab = "Número de predictores",
ylab = "R²-ajustado",
main = "Criterio R²-ajustado"
)
points(
best_r2,
suma$adjr2[best_r2],
pch = 19,
cex = 1.5
)
suma$bic[best_bic]## [1] -126.5072suma$adjr2[best_r2]## [1] 0.8370999suma$which[best_bic, ]## (Intercept) bach_pct desem_pct banda_pct pobreza_pct tam_hogar
## TRUE TRUE FALSE TRUE TRUE FALSEsuma$which[best_r2, ]## (Intercept) bach_pct desem_pct banda_pct pobreza_pct tam_hogar
## TRUE TRUE FALSE TRUE TRUE TRUE#Ajustamos el mejor modelo BIC
coef_best <- coef(s,best_bic)
vars_best <- names(coef_best)[-1]
f_final <- as.formula(paste("ing_pc ~", paste(vars_best, collapse = "+")))
modelo_final <-lm(f_final, data = datos)
summary(modelo_final)##
## Call:
## lm(formula = f_final, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4644.3 -995.6 -35.5 1360.0 5031.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6481.30 2522.94 2.569 0.012215 *
## bach_pct 270.45 67.62 3.999 0.000149 ***
## banda_pct 119.54 44.32 2.697 0.008657 **
## pobreza_pct -209.43 47.36 -4.422 3.3e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1979 on 74 degrees of freedom
## Multiple R-squared: 0.842, Adjusted R-squared: 0.8356
## F-statistic: 131.5 on 3 and 74 DF, p-value: < 2.2e-16#Ajustamos el mejor modelo R^2
coef_best <- coef(s,best_r2)
vars_best <- names(coef_best)[-1]
f_final <- as.formula(paste("ing_pc ~", paste(vars_best, collapse = "+")))
modelo_final <-lm(f_final, data = datos)
summary(modelo_final)##
## Call:
## lm(formula = f_final, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4630.0 -1129.5 -83.4 1182.5 5207.2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4047.29 3138.40 1.290 0.201259
## bach_pct 267.16 67.36 3.966 0.000169 ***
## banda_pct 123.03 44.21 2.783 0.006848 **
## pobreza_pct -209.38 47.15 -4.441 3.12e-05 ***
## tam_hogar 790.86 611.48 1.293 0.199967
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1970 on 73 degrees of freedom
## Multiple R-squared: 0.8456, Adjusted R-squared: 0.8371
## F-statistic: 99.92 on 4 and 73 DF, p-value: < 2.2e-16#Durante el proceso de selección de variables, los tres métodos (forward, backward y híbrido) seleccionaron la misma combinación de variables independientes: el porcentaje de personas con bachillerato o más, el porcentaje de población bajo pobreza y el porcentaje de hogares con conexión de banda ancha.
#En la selección del modelo óptimo mediante los criterios BIC y R^2 ajustado, el BIC seleccionó tres predictores, en consistencia con los métodos de selección anteriores. Por otro lado, el R^2 ajustado seleccionó una cuarta variable, el tamaño promedio del hogar. Sin embargo, esta variable obtuvo un p-value de 0.199, por lo que no es estadísticamente significativa al nivel de 5 %. Por esta razón, se selecciona el modelo de tres predictores mediante el criterio BIC, ya que ofrece un modelo más parsimonioso y consistente con los resultados de los tres métodos de selección de variables.
#El modelo nos presenta un R^2 ajustado de 0.8356, lo que significa que aproximadamente el 83.56% de la variabilidad observada en ingreso per capita es explicada por las tres variables incluidas en el modelo
#Con esto, concluimos que el modelo final BIC es: #ing_pc = 6481.30 + 270.45(bach_pct) + 119.54(banda_pct) - 209.43(pobreza_pct)
```