Tutorial introductorio: Communities and Crime con tidymodels

Regresión lineal, modelos penalizados, Quarto y gráficas estáticas

Autor/a

Manuel Toral

1 Objetivo

Este tutorial muestra un flujo introductorio de modelado predictivo con tidymodels usando el dataset Communities and Crime del UCI Machine Learning Repository. El objetivo es predecir ViolentCrimesPerPop, una medida normalizada de delitos violentos per cápita, usando características socioeconómicas, demográficas y de vivienda de comunidades en Estados Unidos.

Trabajaremos únicamente con modelos lineales: 1. Regresión lineal ordinaria (lm). 2. Ridge regression (glmnet, mezcla = 0). 3. Lasso regression (glmnet, mezcla = 1). 4. Elastic net (glmnet, mezcla ajustable).

Nota ética y técnica

Este dataset incluye variables demográficas sensibles o potencialmente sensibles. El tutorial se centra en aprendizaje estadístico y flujo de trabajo reproducible; no debe interpretarse como una herramienta para tomar decisiones policiales, judiciales o de política pública sin una evaluación profunda de sesgos, contexto histórico, medición y efectos sociales.

2 Paquetes

Código

library(tidyverse)
library(tidymodels)
library(DT)
library(broom)

tidymodels_prefer()
theme_set(theme_minimal(base_size = 12))

3 1. Lectura de datos

Los datos originales usan ? para valores faltantes. Además, la documentación oficial separa el archivo .data de los nombres de columnas, así que definimos los nombres manualmente.

Buenas prácticas en lectura de datos

Al usar read_csv(), el argumento na = "?" es sumamente útil. Evita que R lea las columnas numéricas como texto (caracteres) al encontrarse con un símbolo desconocido. Esto nos ahorra pasos posteriores de coerción de tipos (parsing).

Código

col_names <- c(
  "state", "county", "community", "communityname", "fold",
  "population", "householdsize", "racepctblack", "racePctWhite",
  "racePctAsian", "racePctHisp", "agePct12t21", "agePct12t29",
  "agePct16t24", "agePct65up", "numbUrban", "pctUrban", "medIncome",
  "pctWWage", "pctWFarmSelf", "pctWInvInc", "pctWSocSec", "pctWPubAsst",
  "pctWRetire", "medFamInc", "perCapInc", "whitePerCap", "blackPerCap",
  "indianPerCap", "AsianPerCap", "OtherPerCap", "HispPerCap",
  "NumUnderPov", "PctPopUnderPov", "PctLess9thGrade", "PctNotHSGrad",
  "PctBSorMore", "PctUnemployed", "PctEmploy", "PctEmplManu",
  "PctEmplProfServ", "PctOccupManu", "PctOccupMgmtProf",
  "MalePctDivorce", "MalePctNevMarr", "FemalePctDiv", "TotalPctDiv",
  "PersPerFam", "PctFam2Par", "PctKids2Par", "PctYoungKids2Par",
  "PctTeen2Par", "PctWorkMomYoungKids", "PctWorkMom", "NumIlleg",
  "PctIlleg", "NumImmig", "PctImmigRecent", "PctImmigRec5",
  "PctImmigRec8", "PctImmigRec10", "PctRecentImmig", "PctRecImmig5",
  "PctRecImmig8", "PctRecImmig10", "PctSpeakEnglOnly",
  "PctNotSpeakEnglWell", "PctLargHouseFam", "PctLargHouseOccup",
  "PersPerOccupHous", "PersPerOwnOccHous", "PersPerRentOccHous",
  "PctPersOwnOccup", "PctPersDenseHous", "PctHousLess3BR", "MedNumBR",
  "HousVacant", "PctHousOccup", "PctHousOwnOcc", "PctVacantBoarded",
  "PctVacMore6Mos", "MedYrHousBuilt", "PctHousNoPhone", "PctWOFullPlumb",
  "OwnOccLowQuart", "OwnOccMedVal", "OwnOccHiQuart", "RentLowQ",
  "RentMedian", "RentHighQ", "MedRent", "MedRentPctHousInc",
  "MedOwnCostPctInc", "MedOwnCostPctIncNoMtg", "NumInShelters",
  "NumStreet", "PctForeignBorn", "PctBornSameState", "PctSameHouse85",
  "PctSameCity85", "PctSameState85", "LemasSwornFT", "LemasSwFTPerPop",
  "LemasSwFTFieldOps", "LemasSwFTFieldPerPop", "LemasTotalReq",
  "LemasTotReqPerPop", "PolicReqPerOffic", "PolicPerPop",
  "RacialMatchCommPol", "PctPolicWhite", "PctPolicBlack", "PctPolicHisp",
  "PctPolicAsian", "PctPolicMinor", "OfficAssgnDrugUnits",
  "NumKindsDrugsSeiz", "PolicAveOTWorked", "LandArea", "PopDens",
  "PctUsePubTrans", "PolicCars", "PolicOperBudg", "LemasPctPolicOnPatr",
  "LemasGangUnitDeploy", "LemasPctOfficDrugUn", "PolicBudgPerPop",
  "ViolentCrimesPerPop"
)

url_data <- "https://archive.ics.uci.edu/ml/machine-learning-databases/communities/communities.data"

crime_raw <- read_csv(
  url_data,
  col_names = col_names,
  na = "?",
  col_types = cols(
    communityname = col_character(),
    .default = col_double()
  )
)

glimpse(crime_raw)

Rows: 1,994
Columns: 128
$ state                 <dbl> 8, 53, 24, 34, 42, 6, 44, 6, 21, 29, 6, 36, 25, …
$ county                <dbl> NA, NA, NA, 5, 95, NA, 7, NA, NA, NA, NA, NA, 21…
$ community             <dbl> NA, NA, NA, 81440, 6096, NA, 41500, NA, NA, NA, …
$ communityname         <chr> "Lakewoodcity", "Tukwilacity", "Aberdeentown", "…
$ fold                  <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ population            <dbl> 0.19, 0.00, 0.00, 0.04, 0.01, 0.02, 0.01, 0.01, …
$ householdsize         <dbl> 0.33, 0.16, 0.42, 0.77, 0.55, 0.28, 0.39, 0.74, …
$ racepctblack          <dbl> 0.02, 0.12, 0.49, 1.00, 0.02, 0.06, 0.00, 0.03, …
$ racePctWhite          <dbl> 0.90, 0.74, 0.56, 0.08, 0.95, 0.54, 0.98, 0.46, …
$ racePctAsian          <dbl> 0.12, 0.45, 0.17, 0.12, 0.09, 1.00, 0.06, 0.20, …
$ racePctHisp           <dbl> 0.17, 0.07, 0.04, 0.10, 0.05, 0.25, 0.02, 1.00, …
$ agePct12t21           <dbl> 0.34, 0.26, 0.39, 0.51, 0.38, 0.31, 0.30, 0.52, …
$ agePct12t29           <dbl> 0.47, 0.59, 0.47, 0.50, 0.38, 0.48, 0.37, 0.55, …
$ agePct16t24           <dbl> 0.29, 0.35, 0.28, 0.34, 0.23, 0.27, 0.23, 0.36, …
$ agePct65up            <dbl> 0.32, 0.27, 0.32, 0.21, 0.36, 0.37, 0.60, 0.35, …
$ numbUrban             <dbl> 0.20, 0.02, 0.00, 0.06, 0.02, 0.04, 0.02, 0.00, …
$ pctUrban              <dbl> 1.00, 1.00, 0.00, 1.00, 0.90, 1.00, 0.81, 0.00, …
$ medIncome             <dbl> 0.37, 0.31, 0.30, 0.58, 0.50, 0.52, 0.42, 0.16, …
$ pctWWage              <dbl> 0.72, 0.72, 0.58, 0.89, 0.72, 0.68, 0.50, 0.44, …
$ pctWFarmSelf          <dbl> 0.34, 0.11, 0.19, 0.21, 0.16, 0.20, 0.23, 1.00, …
$ pctWInvInc            <dbl> 0.60, 0.45, 0.39, 0.43, 0.68, 0.61, 0.68, 0.23, …
$ pctWSocSec            <dbl> 0.29, 0.25, 0.38, 0.36, 0.44, 0.28, 0.61, 0.53, …
$ pctWPubAsst           <dbl> 0.15, 0.29, 0.40, 0.20, 0.11, 0.15, 0.21, 0.97, …
$ pctWRetire            <dbl> 0.43, 0.39, 0.84, 0.82, 0.71, 0.25, 0.54, 0.41, …
$ medFamInc             <dbl> 0.39, 0.29, 0.28, 0.51, 0.46, 0.62, 0.43, 0.15, …
$ perCapInc             <dbl> 0.40, 0.37, 0.27, 0.36, 0.43, 0.72, 0.47, 0.10, …
$ whitePerCap           <dbl> 0.39, 0.38, 0.29, 0.40, 0.41, 0.76, 0.44, 0.12, …
$ blackPerCap           <dbl> 0.32, 0.33, 0.27, 0.39, 0.28, 0.77, 0.40, 0.08, …
$ indianPerCap          <dbl> 0.27, 0.16, 0.07, 0.16, 0.00, 0.28, 0.24, 0.17, …
$ AsianPerCap           <dbl> 0.27, 0.30, 0.29, 0.25, 0.74, 0.52, 0.86, 0.27, …
$ OtherPerCap           <dbl> 0.36, 0.22, 0.28, 0.36, 0.51, 0.48, 0.24, 0.18, …
$ HispPerCap            <dbl> 0.41, 0.35, 0.39, 0.44, 0.48, 0.60, 0.36, 0.21, …
$ NumUnderPov           <dbl> 0.08, 0.01, 0.01, 0.01, 0.00, 0.01, 0.01, 0.03, …
$ PctPopUnderPov        <dbl> 0.19, 0.24, 0.27, 0.10, 0.06, 0.12, 0.11, 0.64, …
$ PctLess9thGrade       <dbl> 0.10, 0.14, 0.27, 0.09, 0.25, 0.13, 0.29, 0.96, …
$ PctNotHSGrad          <dbl> 0.18, 0.24, 0.43, 0.25, 0.30, 0.12, 0.41, 0.82, …
$ PctBSorMore           <dbl> 0.48, 0.30, 0.19, 0.31, 0.33, 0.80, 0.36, 0.12, …
$ PctUnemployed         <dbl> 0.27, 0.27, 0.36, 0.33, 0.12, 0.10, 0.28, 1.00, …
$ PctEmploy             <dbl> 0.68, 0.73, 0.58, 0.71, 0.65, 0.65, 0.54, 0.26, …
$ PctEmplManu           <dbl> 0.23, 0.57, 0.32, 0.36, 0.67, 0.19, 0.44, 0.43, …
$ PctEmplProfServ       <dbl> 0.41, 0.15, 0.29, 0.45, 0.38, 0.77, 0.53, 0.34, …
$ PctOccupManu          <dbl> 0.25, 0.42, 0.49, 0.37, 0.42, 0.06, 0.33, 0.71, …
$ PctOccupMgmtProf      <dbl> 0.52, 0.36, 0.32, 0.39, 0.46, 0.91, 0.49, 0.18, …
$ MalePctDivorce        <dbl> 0.68, 1.00, 0.63, 0.34, 0.22, 0.49, 0.25, 0.38, …
$ MalePctNevMarr        <dbl> 0.40, 0.63, 0.41, 0.45, 0.27, 0.57, 0.34, 0.47, …
$ FemalePctDiv          <dbl> 0.75, 0.91, 0.71, 0.49, 0.20, 0.61, 0.28, 0.59, …
$ TotalPctDiv           <dbl> 0.75, 1.00, 0.70, 0.44, 0.21, 0.58, 0.28, 0.52, …
$ PersPerFam            <dbl> 0.35, 0.29, 0.45, 0.75, 0.51, 0.44, 0.42, 0.78, …
$ PctFam2Par            <dbl> 0.55, 0.43, 0.42, 0.65, 0.91, 0.62, 0.77, 0.45, …
$ PctKids2Par           <dbl> 0.59, 0.47, 0.44, 0.54, 0.91, 0.69, 0.81, 0.43, …
$ PctYoungKids2Par      <dbl> 0.61, 0.60, 0.43, 0.83, 0.89, 0.87, 0.79, 0.34, …
$ PctTeen2Par           <dbl> 0.56, 0.39, 0.43, 0.65, 0.85, 0.53, 0.74, 0.34, …
$ PctWorkMomYoungKids   <dbl> 0.74, 0.46, 0.71, 0.85, 0.40, 0.30, 0.57, 0.29, …
$ PctWorkMom            <dbl> 0.76, 0.53, 0.67, 0.86, 0.60, 0.43, 0.62, 0.27, …
$ NumIlleg              <dbl> 0.04, 0.00, 0.01, 0.03, 0.00, 0.00, 0.00, 0.02, …
$ PctIlleg              <dbl> 0.14, 0.24, 0.46, 0.33, 0.06, 0.11, 0.13, 0.50, …
$ NumImmig              <dbl> 0.03, 0.01, 0.00, 0.02, 0.00, 0.04, 0.01, 0.02, …
$ PctImmigRecent        <dbl> 0.24, 0.52, 0.07, 0.11, 0.03, 0.30, 0.00, 0.50, …
$ PctImmigRec5          <dbl> 0.27, 0.62, 0.06, 0.20, 0.07, 0.35, 0.02, 0.59, …
$ PctImmigRec8          <dbl> 0.37, 0.64, 0.15, 0.30, 0.20, 0.43, 0.02, 0.65, …
$ PctImmigRec10         <dbl> 0.39, 0.63, 0.19, 0.31, 0.27, 0.47, 0.10, 0.59, …
$ PctRecentImmig        <dbl> 0.07, 0.25, 0.02, 0.05, 0.01, 0.50, 0.00, 0.69, …
$ PctRecImmig5          <dbl> 0.07, 0.27, 0.02, 0.08, 0.02, 0.50, 0.01, 0.72, …
$ PctRecImmig8          <dbl> 0.08, 0.25, 0.04, 0.11, 0.04, 0.56, 0.01, 0.71, …
$ PctRecImmig10         <dbl> 0.08, 0.23, 0.05, 0.11, 0.05, 0.57, 0.03, 0.60, …
$ PctSpeakEnglOnly      <dbl> 0.89, 0.84, 0.88, 0.81, 0.88, 0.45, 0.73, 0.12, …
$ PctNotSpeakEnglWell   <dbl> 0.06, 0.10, 0.04, 0.08, 0.05, 0.28, 0.05, 0.93, …
$ PctLargHouseFam       <dbl> 0.14, 0.16, 0.20, 0.56, 0.16, 0.25, 0.12, 0.74, …
$ PctLargHouseOccup     <dbl> 0.13, 0.10, 0.20, 0.62, 0.19, 0.19, 0.13, 0.75, …
$ PersPerOccupHous      <dbl> 0.33, 0.17, 0.46, 0.85, 0.59, 0.29, 0.42, 0.80, …
$ PersPerOwnOccHous     <dbl> 0.39, 0.29, 0.52, 0.77, 0.60, 0.53, 0.54, 0.68, …
$ PersPerRentOccHous    <dbl> 0.28, 0.17, 0.43, 1.00, 0.37, 0.18, 0.24, 0.92, …
$ PctPersOwnOccup       <dbl> 0.55, 0.26, 0.42, 0.94, 0.89, 0.39, 0.65, 0.39, …
$ PctPersDenseHous      <dbl> 0.09, 0.20, 0.15, 0.12, 0.02, 0.26, 0.03, 0.89, …
$ PctHousLess3BR        <dbl> 0.51, 0.82, 0.51, 0.01, 0.19, 0.73, 0.46, 0.66, …
$ MedNumBR              <dbl> 0.5, 0.0, 0.5, 0.5, 0.5, 0.0, 0.5, 0.0, 0.0, 0.0…
$ HousVacant            <dbl> 0.21, 0.02, 0.01, 0.01, 0.01, 0.02, 0.01, 0.01, …
$ PctHousOccup          <dbl> 0.71, 0.79, 0.86, 0.97, 0.89, 0.84, 0.89, 0.91, …
$ PctHousOwnOcc         <dbl> 0.52, 0.24, 0.41, 0.96, 0.87, 0.30, 0.57, 0.46, …
$ PctVacantBoarded      <dbl> 0.05, 0.02, 0.29, 0.60, 0.04, 0.16, 0.09, 0.22, …
$ PctVacMore6Mos        <dbl> 0.26, 0.25, 0.30, 0.47, 0.55, 0.28, 0.49, 0.37, …
$ MedYrHousBuilt        <dbl> 0.65, 0.65, 0.52, 0.52, 0.73, 0.25, 0.38, 0.60, …
$ PctHousNoPhone        <dbl> 0.14, 0.16, 0.47, 0.11, 0.05, 0.02, 0.05, 0.28, …
$ PctWOFullPlumb        <dbl> 0.06, 0.00, 0.45, 0.11, 0.14, 0.05, 0.05, 0.23, …
$ OwnOccLowQuart        <dbl> 0.22, 0.21, 0.18, 0.24, 0.31, 0.94, 0.37, 0.15, …
$ OwnOccMedVal          <dbl> 0.19, 0.20, 0.17, 0.21, 0.31, 1.00, 0.38, 0.13, …
$ OwnOccHiQuart         <dbl> 0.18, 0.21, 0.16, 0.19, 0.30, 1.00, 0.39, 0.13, …
$ RentLowQ              <dbl> 0.36, 0.42, 0.27, 0.75, 0.40, 0.67, 0.26, 0.21, …
$ RentMedian            <dbl> 0.35, 0.38, 0.29, 0.70, 0.36, 0.63, 0.35, 0.24, …
$ RentHighQ             <dbl> 0.38, 0.40, 0.27, 0.77, 0.38, 0.68, 0.42, 0.25, …
$ MedRent               <dbl> 0.34, 0.37, 0.31, 0.89, 0.38, 0.62, 0.35, 0.24, …
$ MedRentPctHousInc     <dbl> 0.38, 0.29, 0.48, 0.63, 0.22, 0.47, 0.46, 0.64, …
$ MedOwnCostPctInc      <dbl> 0.46, 0.32, 0.39, 0.51, 0.51, 0.59, 0.44, 0.59, …
$ MedOwnCostPctIncNoMtg <dbl> 0.25, 0.18, 0.28, 0.47, 0.21, 0.11, 0.31, 0.28, …
$ NumInShelters         <dbl> 0.04, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, …
$ NumStreet             <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, …
$ PctForeignBorn        <dbl> 0.12, 0.21, 0.14, 0.19, 0.11, 0.70, 0.15, 0.59, …
$ PctBornSameState      <dbl> 0.42, 0.50, 0.49, 0.30, 0.72, 0.42, 0.81, 0.58, …
$ PctSameHouse85        <dbl> 0.50, 0.34, 0.54, 0.73, 0.64, 0.49, 0.77, 0.52, …
$ PctSameCity85         <dbl> 0.51, 0.60, 0.67, 0.64, 0.61, 0.73, 0.91, 0.79, …
$ PctSameState85        <dbl> 0.64, 0.52, 0.56, 0.65, 0.53, 0.64, 0.84, 0.78, …
$ LemasSwornFT          <dbl> 0.03, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ LemasSwFTPerPop       <dbl> 0.13, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ LemasSwFTFieldOps     <dbl> 0.96, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ LemasSwFTFieldPerPop  <dbl> 0.17, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ LemasTotalReq         <dbl> 0.06, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ LemasTotReqPerPop     <dbl> 0.18, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ PolicReqPerOffic      <dbl> 0.44, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ PolicPerPop           <dbl> 0.13, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ RacialMatchCommPol    <dbl> 0.94, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ PctPolicWhite         <dbl> 0.93, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ PctPolicBlack         <dbl> 0.03, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ PctPolicHisp          <dbl> 0.07, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ PctPolicAsian         <dbl> 0.10, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ PctPolicMinor         <dbl> 0.07, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ OfficAssgnDrugUnits   <dbl> 0.02, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ NumKindsDrugsSeiz     <dbl> 0.57, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ PolicAveOTWorked      <dbl> 0.29, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ LandArea              <dbl> 0.12, 0.02, 0.01, 0.02, 0.04, 0.01, 0.05, 0.01, …
$ PopDens               <dbl> 0.26, 0.12, 0.21, 0.39, 0.09, 0.58, 0.08, 0.33, …
$ PctUsePubTrans        <dbl> 0.20, 0.45, 0.02, 0.28, 0.02, 0.10, 0.06, 0.00, …
$ PolicCars             <dbl> 0.06, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ PolicOperBudg         <dbl> 0.04, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ LemasPctPolicOnPatr   <dbl> 0.90, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ LemasGangUnitDeploy   <dbl> 0.5, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ LemasPctOfficDrugUn   <dbl> 0.32, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, …
$ PolicBudgPerPop       <dbl> 0.14, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ ViolentCrimesPerPop   <dbl> 0.20, 0.67, 0.43, 0.12, 0.03, 0.14, 0.03, 0.55, …

Nota técnica: Datos normalizados

En este dataset, la mayoría de variables numéricas ya vienen normalizadas al intervalo 0–1. Aun así, en el flujo de tidymodels usaremos step_normalize() después de remover variables de varianza cero. Esto ayuda a que los coeficientes penalizados de glmnet sean estrictamente comparables y estables.

4 2. Auditoría de valores faltantes

Código

missing_tbl <- crime_raw |>
  summarise(across(everything(), ~ sum(is.na(.)))) |>
  pivot_longer(everything(), names_to = "variable", values_to = "n_na") |>
  mutate(prop_na = n_na / nrow(crime_raw)) |>
  filter(n_na > 0) |>
  arrange(desc(prop_na))

DT::datatable(
  missing_tbl,
  rownames = FALSE,
  options = list(pageLength = 10, scrollX = TRUE),
  caption = "Variables con valores faltantes"
)

Código

p_missing <- missing_tbl |>
  slice_max(prop_na, n = 25) |>
  mutate(variable = fct_reorder(variable, prop_na)) |>
  ggplot(aes(x = prop_na, y = variable)) +
  geom_col() +
  scale_x_continuous(labels = scales::percent) +
  labs(
    title = "Variables con más valores faltantes",
    x = "Proporción de valores faltantes",
    y = NULL
  )

p_missing

4.1 Limpieza inicial

Quitamos identificadores no predictivos y variables policiales/LEMAS con muchísimos valores faltantes.

Decisiones de modelado: ¿Imputar o eliminar?

Esta decisión replica una estrategia común para este dataset: retirar columnas que tienen alrededor de 80% o más de ausencias. Cuando una variable tiene tantos datos faltantes, imputarla (ej. con la media, KNN o Random Forest) puede introducir más ruido que señal. Es más seguro descartarlas o, en un escenario avanzado, crear una variable binaria (is_missing) si creemos que la ausencia de información no es aleatoria.

Código

crime <- crime_raw |>
  select(
    -state, -county, -community, -communityname, -fold,
    -starts_with("Lemas"), -starts_with("Polic"),
    -contains("PctPolic"), -contains("OtherPerCap"),
    -PolicBudgPerPop, -OfficAssgnDrugUnits,
    -NumKindsDrugsSeiz, -PolicAveOTWorked,
    -RacialMatchCommPol
  )

crime |>
  summarise(across(everything(), ~ sum(is.na(.)))) |>
  pivot_longer(everything()) |>
  filter(value > 0)

# A tibble: 0 × 2
# ℹ 2 variables: name <chr>, value <int>

5 3. Exploración visual

5.1 Variable objetivo

Código

p_y <- crime |>
  ggplot(aes(x = ViolentCrimesPerPop)) +
  geom_histogram(bins = 40, fill = "steelblue", color = "white") +
  labs(
    title = "Distribución de delitos violentos per cápita",
    x = "ViolentCrimesPerPop",
    y = "Número de comunidades"
  )

p_y

5.2 Correlaciones con la respuesta

Código

cor_target <- crime |>
  summarise(across(-ViolentCrimesPerPop, ~ cor(.x, ViolentCrimesPerPop))) |>
  pivot_longer(everything(), names_to = "variable", values_to = "correlation") |>
  mutate(abs_correlation = abs(correlation)) |>
  arrange(desc(abs_correlation))

p_cor <- cor_target |>
  slice_max(abs_correlation, n = 25) |>
  mutate(variable = fct_reorder(variable, correlation)) |>
  ggplot(aes(x = correlation, y = variable, fill = correlation > 0)) +
  geom_col(show.legend = FALSE) +
  scale_fill_manual(values = c("tomato", "steelblue")) +
  labs(
    title = "Predictores con mayor correlación con la respuesta",
    x = "Correlación de Pearson",
    y = NULL
  )

p_cor

5.3 Relaciones bivariadas estáticas

Código

top_vars <- cor_target |>
  slice_max(abs_correlation, n = 9) |>
  pull(variable)

scatter_tbl <- crime |>
  select(ViolentCrimesPerPop, all_of(top_vars)) |>
  pivot_longer(-ViolentCrimesPerPop, names_to = "variable", values_to = "value")

p_scatter <- scatter_tbl |>
  ggplot(aes(
    x = value,
    y = ViolentCrimesPerPop
  )) +
  geom_point(alpha = 0.35) +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  facet_wrap(~ variable, scales = "free_x") +
  labs(
    title = "Relaciones lineales bivariadas: top predictores",
    x = "Valor normalizado del predictor",
    y = "ViolentCrimesPerPop"
  )

p_scatter

6 4. Separación train/test

Código

set.seed(42)
split <- initial_split(crime, prop = 0.80, strata = ViolentCrimesPerPop)
train <- training(split)
test  <- testing(split)

Nota técnica: Estratificación en regresión

Aunque ViolentCrimesPerPop es continua, initial_split(..., strata = ViolentCrimesPerPop) agrupa los datos en cuantiles (bins) de forma interna para asegurar que tanto el conjunto de entrenamiento como el de prueba mantengan la misma distribución de la variable objetivo. Esto reduce el riesgo de tener un modelo entrenado sin conocer los valores más extremos.

7 5. Receta de preprocesamiento

Usaremos una receta común para todos los modelos.

Código

crime_recipe <- recipe(ViolentCrimesPerPop ~ ., data = train) |>
  step_zv(all_predictors()) |>
  step_normalize(all_predictors())

crime_recipe

Concepto Clave: Fuga de Información (Data Leakage)

¿Por qué normalizar dentro de la receta y no antes de separar los datos? Si aplicáramos la normalización a todo el dataset original, la media y desviación estándar utilizadas para escalar los datos incluirían información del conjunto de validación/prueba. Esto causaría un sesgo optimista en las métricas. Al usar step_normalize() en la receta, estas transformaciones se estiman estrictamente con el set de entrenamiento.

8 6. Remuestreo y métricas

La validación cruzada estima el desempeño del modelo de forma más estable que una sola partición. En cada fold, el modelo se entrena con una parte de los datos y se evalúa en la parte restante.

Código

set.seed(42)
folds <- vfold_cv(train, v = 10, strata = ViolentCrimesPerPop)

reg_metrics <- metric_set(rmse, rsq, mae)

9 7. Modelo 1: regresión lineal ordinaria

Empezamos con la regresión lineal ordinaria porque funciona como modelo base. Antes de introducir penalizaciones, conviene saber qué tan bien predice una combinación lineal simple.

\[\widehat{y} = \beta_0 + \beta_1x_1 + \beta_2x_2 + \cdots + \beta_px_p\]

Código

lm_spec <- linear_reg() |>
  set_engine("lm")

lm_wf <- workflow() |>
  add_recipe(crime_recipe) |>
  add_model(lm_spec)

set.seed(42)
lm_res <- fit_resamples(
  lm_wf,
  resamples = folds,
  metrics = reg_metrics,
  control = control_resamples(save_pred = TRUE)
)

collect_metrics(lm_res)

# A tibble: 3 × 6
  .metric .estimator   mean     n std_err .config             
  <chr>   <chr>       <dbl> <int>   <dbl> <chr>               
1 mae     standard   0.0960    10 0.00190 Preprocessor1_Model1
2 rmse    standard   0.136     10 0.00267 Preprocessor1_Model1
3 rsq     standard   0.660     10 0.0206  Preprocessor1_Model1

10 8. Modelos penalizados

Después del modelo base añadimos ridge, lasso y elastic net. Estos modelos siguen siendo lineales, pero agregan una penalización para controlar la magnitud de los coeficientes, lo cual es vital con datos altamente colineales.

El Espectro de la Regularización: Ridge vs. Lasso

Ridge (L2 - mixture = 0): Penaliza la suma de los coeficientes al cuadrado. Reduce todos los coeficientes hacia cero de forma proporcional, pero nunca los elimina por completo. Ideal para variables muy correlacionadas.
Lasso (L1 - mixture = 1): Penaliza la suma del valor absoluto de los coeficientes. Puede llevar los coeficientes exactamente a cero, actuando como un método automático de selección de características.
Elastic net: Una combinación lineal de ambas penalizaciones gobernada por el parámetro mixture.

10.1 Grid de penalización

Código

set.seed(42)
penalty_grid <- grid_regular(
  penalty(range = c(-4, 0)),
  levels = 30
)

10.2 Ridge

Código

ridge_spec <- linear_reg(penalty = tune(), mixture = 0) |>
  set_engine("glmnet")

ridge_wf <- workflow() |>
  add_recipe(crime_recipe) |>
  add_model(ridge_spec)

set.seed(42)
ridge_res <- tune_grid(
  ridge_wf,
  resamples = folds,
  grid = penalty_grid,
  metrics = reg_metrics
)

10.3 Lasso

Código

lasso_spec <- linear_reg(penalty = tune(), mixture = 1) |>
  set_engine("glmnet")

lasso_wf <- workflow() |>
  add_recipe(crime_recipe) |>
  add_model(lasso_spec)

set.seed(42)
lasso_res <- tune_grid(
  lasso_wf,
  resamples = folds,
  grid = penalty_grid,
  metrics = reg_metrics
)

10.4 Elastic net

Código

elastic_spec <- linear_reg(penalty = tune(), mixture = tune()) |>
  set_engine("glmnet")

elastic_wf <- workflow() |>
  add_recipe(crime_recipe) |>
  add_model(elastic_spec)

set.seed(42)
elastic_grid <- grid_regular(
  penalty(range = c(-4, 0)),
  mixture(range = c(0.05, 1)),
  levels = c(20, 6)
)

set.seed(42)
elastic_res <- tune_grid(
  elastic_wf,
  resamples = folds,
  grid = elastic_grid,
  metrics = reg_metrics
)

11 9. Comparación de modelos

Código

model_results <- list(
  lm = lm_res,
  ridge = ridge_res,
  lasso = lasso_res,
  elastic_net = elastic_res
)

comparison_tbl <- map_dfr(
  model_results,
  collect_metrics,
  .id = "modelo"
) |>
  select(modelo, .metric, mean, std_err, n, .config) |>
  arrange(.metric, mean)

p_compare <- comparison_tbl |>
  filter(.metric %in% c("rmse", "mae")) |>
  group_by(modelo, .metric) |>
  slice_min(mean, n = 1, with_ties = FALSE) |>
  ungroup() |>
  ggplot(aes(x = modelo, y = mean, fill = modelo)) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~ .metric, scales = "free_y") +
  labs(
    title = "Mejor resultado de cada modelo según RMSE y MAE",
    x = NULL,
    y = "Valor promedio en validación cruzada"
  )

p_compare

12 10. Seleccionar mejor modelo y evaluar en test

Elegimos el modelo con menor RMSE promedio en validación cruzada.

Código

best_elastic <- select_best(elastic_res, metric = "rmse")

final_elastic_wf <- finalize_workflow(elastic_wf, best_elastic)

final_fit <- last_fit(
  final_elastic_wf,
  split = split,
  metrics = reg_metrics
)

test_preds <- collect_predictions(final_fit)

Código

p_pred <- test_preds |>
  ggplot(aes(
    x = ViolentCrimesPerPop,
    y = .pred
  )) +
  geom_abline(linetype = 2, color = "red") +
  geom_point(alpha = 0.55) +
  coord_equal() +
  labs(
    title = "Predicción vs. valor observado en test",
    x = "Valor observado",
    y = "Valor predicho"
  )

p_pred

13 11. Interpretación de coeficientes

Para modelos lineales penalizados con glmnet, los coeficientes dependen del valor de penalty. Extraemos los coeficientes del workflow final.

Código

final_fit_obj <- extract_workflow(final_fit)

coef_tbl <- final_fit_obj |>
  extract_fit_parsnip() |>
  tidy() |>
  filter(term != "(Intercept)") |>
  mutate(abs_estimate = abs(estimate)) |>
  arrange(desc(abs_estimate))

p_coef <- coef_tbl |>
  slice_max(abs_estimate, n = 25) |>
  mutate(term = fct_reorder(term, estimate)) |>
  ggplot(aes(
    x = estimate,
    y = term,
    fill = estimate > 0
  )) +
  geom_col(show.legend = FALSE) +
  scale_fill_manual(values = c("steelblue", "tomato")) +
  labs(
    title = "Coeficientes más grandes en valor absoluto",
    x = "Coeficiente estimado",
    y = NULL
  )

p_coef

Cuidado al interpretar: Correlación vs. Causalidad

Un coeficiente positivo o negativo fuerte no significa causalidad. En este flujo, el modelo aprende asociaciones lineales condicionadas al conjunto de predictores disponibles. Las variables están altamente entrelazadas e interrelacionadas, representando dimensiones sociales y económicas muy complejas. Para hacer inferencia causal real, se requerirían diseños de estudio, control estricto de variables confusoras (confounders) y pruebas estadísticas mucho más profundas.

14 12. Resumen

En este tutorial construimos un flujo completo con Quarto y tidymodels: 1. Lectura y limpieza de datos UCI. 2. Auditoría de faltantes y manejo inicial. 3. Visualización estática de la variable objetivo y predictores. 4. Separación train/test estratificada para mantener la representatividad. 5. Receta reproducible de preprocesamiento, cuidando de evitar fugas de información. 6. Validación cruzada mediante K-Folds. 7. Comparación de modelos lineales ordinarios frente a penalizados (Ridge, Lasso, Elastic Net). 8. Evaluación final en el entorno de Test.

El resultado no debe verse sólo como “un modelo que predice crimen”, sino como una práctica de modelado responsable: entender los datos, sus límites, su contexto y la incertidumbre detrás de cada predicción.