13 Segundo caso de negocio

Objetivos del capítulo

El lector, al finalizar este capítulo, estará en capacidad de:

  • Emplear las herramientas estudiadas en los capítulos anteriores para responder una pregunta de negocio que implique analítica predictiva.
  • Presentar los resultados de una regresión de manera gráfica empleando R con solución H.A.C..
  • Determinar cuál variable tiene más efecto sobre la variable explicativa empleando R en presencia de heteroscedasticidad.

13.1 Introducción

En los capítulos anteriores hemos estudiado las bases del modelo clásico de regresión múltiple, cómo encontrar el mejor modelo para hacer analítica diagnóstica o analítica predictiva. En este capítulo pondremos todos los elementos juntos para resolver un caso de negocio que implica analítica predictiva.

13.2 La pregunta de negocio

Una empresa de consultoría en seguridad en los Estados Unidos quiere contar con un modelo que le permita, al visitar a un alcalde de cualquier ciudad, poder tener un estimado de la tasa de crímenes violentos por cada 100 mil habitantes y como esta puede cambiar cuando se “intervengan” algunas variables. Es decir, esta empresa quiere tener una “fórmula” que le permita determinar cuál sería la tasa de crímenes violentos por cada 100 mil habitantes bajo diferentes escenarios. La empresa quiere ofrecer una nueva linea de servicios que le permita asesorar a los alcaldes en que variables intervenir para disminuir la tasa de crímenes violentos.

Para realizar esta tarea, contamos con una base de datos con información sociodemográfica y datos de crimen del FBI para diferentes comunidades de Estados Unidos, los datos se encuentran en el archivo DatosCaso2.csv133. Los datos son reales y fueron suministrada por Redmond & Baveja (2002) y tomados de la siguiente página https://archive.ics.uci.edu/ml/datasets/Communities+and+Crime.

La base de datos contiene 1994 observaciones y las siguientes 101 variables:

  • population: población para la comunidad: (numérico - decimal)
  • householdsize: media de personas por hogar (numérico - decimal)
  • racepctblack: porcentaje de población afroamericana (numérico - decimal)
  • racePctWhite: porcentaje de población caucásica (numérico - decimal)
  • racePctAsian: porcentaje de población de origen asiático (numérico - decimal)
  • racePctHisp: porcentaje de población de origen hispano (numérico - decimal)
  • agePct12t21: porcentaje de población de 12-21 años (numérico - decimal)
  • agePct12t29: porcentaje de población de 21-29 años (numérico - decimal)
  • agePct16t24: porcentaje de población de 16-24 años (numérico - decimal)
  • agePct65up: porcentaje de población de 65 años o más (numérico - decimal)
  • numbUrban: número de personas que viven en zonas clasificadas como urbanas (numérico - decimal)
  • pctUrban: porcentaje de personas que viven en zonas clasificadas como urbanas (numérico - decimal)
  • medIncome: ingreso familiar medio (numérico - decimal)
  • pctWWage: porcentaje de hogares con ingresos salariales en 1989 (numérico - decimal)
  • pctWFarmSelf: porcentaje de hogares con ingresos agrícolas o por cuenta propia en 1989 (numérico - decimal)
  • pctWInvInc: porcentaje de hogares con ingresos por inversiones/alquileres en 1989 (numérico - decimal)
  • pctWSocSec: porcentaje de hogares con ingresos de la seguridad social en 1989 (numérico - decimal)
  • pctWPubAsst: porcentaje de hogares con ingresos de la asistencia pública en 1989 (numérico - decimal)
  • pctWRetire: porcentaje de hogares con ingresos por jubilación en 1989 (numérico - decimal)
  • medFamInc: renta familiar mediana (difiere de la renta de los hogares no familiares) (numérico - decimal)
  • perCapInc: renta per cápita (numérico - decimal)
  • whitePerCap: renta per cápita de los caucásicos (numérico - decimal)
  • blackPerCap: renta per cápita de los afroamericanos (numérico - decimal)
  • IndianPerCap: renta per cápita de los nativos americanos (numérico - decimal)
  • AsianPerCap: renta per cápita de las personas de origen asiático (numérico - decimal)
  • OtherPerCap: renta per cápita de las personas con “otra” herencia (numérico - decimal)
  • HispPerCap: renta per cápita de las personas de origen hispano (numérico - decimal)
  • NumUnderPov: número de personas por debajo del nivel de pobreza (numérico - decimal)
  • PctPopUnderPov: porcentaje de personas bajo el nivel de pobreza (numérico - decimal)
  • PctLess9thGrade: porcentaje de personas de 25 años o más con menos de 9 grados de educación (numérico - decimal)
  • PctNotHSGrad: porcentaje de personas de 25 años o más que no tienen estudios secundarios (numérico - decimal)
  • PctBSorMore: porcentaje de personas de 25 años o más con una licenciatura o educación superior (numérico - decimal)
  • PctUnemployed: porcentaje de personas de 16 años o más, que forman parte de la población activa y están desempleadas (numérico - decimal)
  • PctEmploy: porcentaje de personas de 16 años o más que están empleadas (numérico - decimal)
  • PctEmplManu: porcentaje de personas de 16 años o más que trabajan en la industria manufacturera (numérico - decimal)
  • PctEmplProfServ: porcentaje de personas de 16 años o más que trabajan en servicios profesionales (numérico - decimal)
  • PctOccupManu: porcentaje de personas de 16 años o más que trabajan en la industria manufacturera (numérico - decimal)
  • PctOccupMgmtProf: porcentaje de personas de 16 años o más empleadas en ocupaciones de gestión o profesionales (numérico - decimal)
  • MalePctDivorce: porcentaje de hombres divorciados (numérico - decimal)
  • MalePctNevMarr: porcentaje de hombres que nunca se han casado (numérico - decimal)
  • FemalePctDiv: porcentaje de mujeres divorciadas (numérico - decimal)
  • TotalPctDiv: porcentaje de población divorciada (numérico - decimal)
  • PersPerFam: número medio de personas por familia (numérico - decimal)
  • PctFam2Par: porcentaje de familias (con hijos) encabezadas por dos padres (numérico - decimal)
  • PctKids2Par: porcentaje de niños en viviendas familiares con dos padres (numérico - decimal)
  • PctYoungKids2Par: porcentaje de niños de 4 años o menos en hogares con dos padres (numérico - decimal)
  • PctTeen2Par: porcentaje de niños de 12 a 17 años en hogares con dos padres (numérico - decimal)
  • PctWorkMomYoungKids: porcentaje de madres de niños de 6 años o menos que trabajan (numérico - decimal)
  • PctWorkMom: porcentaje de madres de niños menores de 18 años que trabajan (numérico - decimal)
  • NumIlleg: número de niños nacidos de madres nunca casadas (numérico - decimal)
  • PctIlleg: porcentaje de hijos nacidos de personas nunca casadas (numérico - decimal)
  • NumImmig: número total de personas que se sabe que han nacido en el extranjero (numérico - decimal)
  • PctImmigRecent: porcentaje de inmigrantes que han inmigrado en los últimos 3 años (numérico - decimal)
  • PctImmigRec5: porcentaje de inmigrantes que han inmigrado en los últimos 5 años (numérico - decimal)
  • PctImmigRec8: porcentaje de inmigrantes que han inmigrado en los últimos 8 años (numérico - decimal)
  • PctImmigRec10: porcentaje de inmigrantes que han inmigrado en los últimos 10 años (numérico - decimal)
  • PctRecentImmig: porcentaje de la población que ha inmigrado en los últimos 3 años (numérico - decimal)
  • PctRecImmig5: porcentaje de población que ha inmigrado en los últimos 5 años (numérico - decimal)
  • PctRecImmig8: porcentaje de población que ha inmigrado en los últimos 8 años (numérico - decimal)
  • PctRecImmig10: porcentaje de población que ha inmigrado en los últimos 10 años (numérico - decimal)
  • PctSpeakEnglOnly: porcentaje de personas que sólo hablan inglés (numérico - decimal)
  • PctNotSpeakEnglWell: porcentaje de personas que no hablan bien el inglés (numérico - decimal)
  • PctLargHouseFam: porcentaje de hogares familiares que son grandes (6 o más) (numérico - decimal)
  • PctLargHouseOccup: porcentaje de hogares ocupados que son grandes (6 o más personas) (numérico - decimal)
  • PersPerOccupHous: media de personas por hogar (numérico - decimal)
  • PersPerOwnOccHous: media de personas por hogar ocupado por el propietario (numérico - decimal)
  • PersPerRentOccHous: media de personas por hogar de alquiler (numérico - decimal)
  • PctPersOwnOccup: porcentaje de personas en hogares ocupados por el propietario (numérico - decimal)
  • PctPersDenseHous: porcentaje de personas en viviendas densas (más de 1 persona por habitación) (numérico - decimal)
  • PctHousLess3BR: porcentaje de viviendas con menos de 3 dormitorios (numérico - decimal)
  • MedNumBR: número medio de dormitorios (numérico - decimal)
  • HousVacant: número de viviendas vacías (numérico - decimal)
  • PctHousOccup: porcentaje de viviendas ocupadas (numérico - decimal)
  • PctHousOwnOcc: porcentaje de hogares ocupados por el propietario (numérico - decimal)
  • PctVacantBoarded: porcentaje de viviendas vacías que están tapiadas (numérico - decimal)
  • PctVacMore6Mos: porcentaje de viviendas desocupadas que han estado desocupadas más de 6 meses (numérico - decimal)
  • MedYrHousBuilt: año medio de construcción de las viviendas (numérico - decimal)
  • PctHousNoPhone: porcentaje de viviendas ocupadas sin teléfono (en 1990, esto era raro) (numérico - decimal)
  • PctWOFullPlumb: porcentaje de viviendas sin instalaciones completas de acueducto (numérico - decimal)
  • OwnOccLowQuart: viviendas ocupadas por sus propietarios - valor del cuartil inferior (numérico - decimal)
  • OwnOccMedVal: viviendas ocupadas por sus propietarios - valor medio (numérico - decimal)
  • OwnOccHiQuart: vivienda ocupada por el propietario - valor del cuartil superior (numérico - decimal)
  • RentLowQ: vivienda de alquiler - renta del cuartil inferior (numérico - decimal)
  • RentMedian: vivienda de alquiler - renta mediana (variable censal H32B del fichero STF1A) (numérico - decimal)
  • RentHighQ: vivienda de alquiler - renta del cuartil superior (numérico - decimal)
  • MedRent: alquiler bruto medio (variable censal H43A del fichero STF3A - incluye los servicios públicos) (numérico - decimal)
  • MedRentPctHousInc: alquiler bruto medio como porcentaje de los ingresos del hogar (numérico - decimal)
  • MedOwnCostPctInc: costo medio de los propietarios como porcentaje de los ingresos del hogar - para propietarios con hipoteca (numérico - decimal)
  • MedOwnCostPctIncNoMtg: costo medio de los propietarios como porcentaje de los ingresos del hogar - para los propietarios sin hipoteca (numérico - decimal)
  • NumInShelters: número de personas en refugios para personas sin hogar (numérico - decimal)
  • NumStreet: número de personas sin hogar contabilizadas en la calle (numérico - decimal)
  • PctForeignBorn: porcentaje de personas nacidas en el extranjero (numérico - decimal)
  • PctBornSameState: porcentaje de personas nacidas en el mismo estado en el que viven actualmente (numérico - decimal)
  • PctSameHouse85: porcentaje de personas que viven en la misma casa que en 1985 (5 años antes) (numérico - decimal)
  • PctSameCity85: porcentaje de personas que viven en la misma ciudad que en 1985 (5 años antes) (numérico - decimal)
  • PctSameState85: porcentaje de personas que viven en el mismo estado que en 1985 (5 años antes) (numérico - decimal)
  • LandArea: superficie terrestre en millas cuadradas (numérico - decimal)
  • PopDens: densidad de población en personas por milla cuadrada (numérico - decimal)
  • PctUsePubTrans: porcentaje de personas que utilizan el transporte público para desplazarse (numérico - decimal)
  • LemasPctOfficDrugUn: porcentaje de agentes asignados a unidades de drogas (numérico - decimal)
  • ViolentCrimesPerPop: número total de delitos violentos por cada 100 mil habitantes (numérico - decimal) (Variable dependiente)

Responder esta pregunta de negocio implicaba realizar analítica predictiva al tener que encontrar una “fórmula” que permita determinar cuál sería la tasa de crímenes violentos por cada 100 mil habitantes para diferentes valores de las variables explicativas (escenarios).

13.3 El plan

Recordemos que nuestra primera tarea siempre es trazar una ruta analítica para responder la pregunta de negocio. Para este momento, ya deben estar intuyendo que es muy probable que exista un problema de heteroscedasticidad dado que la base de datos con que trabajaremos corresponde a datos de corte transversal. Esto implicará que nuestra ruta tendrá los siguientes pasos:

  • Encontrar diferentes modelos candidatos a ser el mejor modelo.
  • Determinar si existe un problema de heteroscedasticidad en los modelos candidatos.
  • Limpiar los modelos candidatos de variables no significativas.
  • Comparar los modelos candidatos para seleccionar el mejor modelo para hacer analítica descriptiva empleando Validación cruzada de \(k\) iteraciones.
  • Determinar si existe un problema de multicolinealidad.
  • Generar escenarios para mostrar como funcionaría la herramienta.

Empecemos a ejecutar esa ruta analítica para resolver la pregunta de negocio.

13.4 Detección de posibles modelos

Lean los datos empleando la función read.csv() , guárdenlos en el objeto datos.caso2 y eliminen la primera variable que no es relevante (se carga una variable “X” con el número de la observación. Esto puede cambiar de computador a computador.)

datos.caso2 <- read.csv("./Data/DatosCaso2.csv", sep=",")

datos.caso2 <- datos.caso2[,-1]

El siguiente paso es encontrar los mejores modelos empleando las estrategias de regresión paso a paso forward, backward y combinada con el AIC, con el valor p y el \(R^2\) ajustado. Como lo discutimos en el Capítulo 10, dado que es posible que exista heteroscedasticidad (más adelante lo demostraremos) es mejor no emplear el criterio del valor p.

Entonces tendremos 6 opciones de algoritmos y criterios de selección que se presentan en el Cuadro 13.1. Los modelos que obtenemos los guardaremos con los nombres que se presentan en el Cuadro 13.1.

Cuadro 13.1: Modelos a estimar con los diferentes algoritmos
Nombre del objeto Algoritmo Criterio
modelo1.C2 Forward \(R^2\) ajustado
modelo3.C2 Forward AIC
modelo4.C2 Backward \(R^2\) ajustado
modelo6.C2 Backward AIC
modelo7.C2 Both \(R^2\) ajustado
modelo9.C2 Both AIC
Fuente: elaboración propia.

Partamos de estimar los modelos lineales con todas las variables potenciales (max.model).

# modelo con todas las variables
max.model <- lm(ViolentCrimesPerPop ~.,data = datos.caso2)

Ahora procedamos a encontrar modelos candidatos para ser los mejores modelos. Empecemos con la estrategia stepwise Forward.

13.4.1 Stepwise forward

Empecemos por el modelo 1: algoritmo Forward y criterio de \(R^2\) ajustado. Estímenlo, guarda lo en el objeto modelo1.C2 y ahora constatemos si el modelo 1 tiene o no heteroscedasticidad. Parte del código no se presenta de manera intencional; ¡intenta reproducir todos los resultados!

Para realizar la prueba de Breusch-Pagan, debemos constatar si el supuesto de normalidad se cumple. Para esto extraigamos los residuales con la función residuals() y efectúemos las pruebas de normalidad con la función ks.test() del paquete base (Prueba de Kolmogorov-Smirnov) y la función jarque.bera.test() del paquete tseries (Trapletti & Hornik, 2019). (Prueba de Jarque-Bera) 

# se extraen los residuales
res.modelo1.C2 <- residuals(modelo1.C2)
# pruebas de normalidad Kolmogorov-Smirnov
ks.test(res.modelo1.C2, y= pnorm)
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  res.modelo1.C2
## D = 0.38847, p-value < 2.2e-16
## alternative hypothesis: two-sided
# pruebas de normalidad Kolmogorov-Smirnov
library(tseries)
jarque.bera.test(res.modelo1.C2)
## 
##  Jarque Bera Test
## 
## data:  res.modelo1.C2
## X-squared = 1185.3, df = 2, p-value < 2.2e-16

Los residuales de este modelo no siguen una distribución normal. Por eso no son confiables los resultados de la prueba de Breusch-Pagan tradicional. Deberíamos entonces emplear la versión studentizada de la prueba propuesta por Koenker (1981).

# se carga la librería
library(lmtest)
# prueba de Breusch-Pagan studentizada
bptest(modelo1.C2, studentize = TRUE)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo1.C2
## BP = 307.37, df = 65, p-value < 2.2e-16

Con un 99% de confianza podemos rechazar la hipótesis de homoscedasticidad. Y para la prueba de White obtendremos el mismo resultado (¡Inténtelo!).

Por tanto, al eliminar las variables no significativas podemos emplear la función construida en el Capítulo10: remueve.no.sinifica.HC3() .

El modelo 1 después de la eliminación de las variables explicativas no significativas (con la corrección HC3) se reporta en el Cuadro 13.2.

#remueve las variables no significativas HC3
modelo1.C2 <- remueve.no.sinifica.HC3(modelo1.C2, 0.05)

Ahora sigamos con el modelo 2: algoritmo Forward y criterio AIC. Pueden encontrar que este modelo también tiene heteroscedasticidad (¡Háganlo!). Ahora realicemos la corrección HC3 y eliminemos las variables no significativas. El resultado se reporta en el Cuadro 13.2.

Cuadro 13.2: Modelo seleccionado el algoritmo stepwise forward con corrección HC3
  Modelo 1 (HC3) (\(R^2\) aj) Modelo 2 (HC3) (AIC)
(Intercept) 0.65*** 0.75***
  (0.11) (0.10)
racepctblack 0.20*** 0.17***
  (0.04) (0.05)
agePct12t29 -0.26*** -0.39***
  (0.08) (0.11)
pctUrban 0.04*** 0.04***
  (0.01) (0.01)
pctWWage -0.26*** -0.21***
  (0.06) (0.05)
pctWFarmSelf 0.05** 0.03
  (0.02) (0.02)
pctWInvInc -0.16** -0.19***
  (0.05) (0.06)
pctWRetire -0.06* -0.10***
  (0.03) (0.03)
medFamInc 0.18*  
  (0.09)  
whitePerCap -0.19** -0.18***
  (0.07) (0.05)
indianPerCap -0.04* -0.04*
  (0.02) (0.02)
OtherPerCap 0.05** 0.05**
  (0.02) (0.02)
PctPopUnderPov -0.09* -0.19***
  (0.04) (0.05)
PctEmploy 0.20*** 0.20***
  (0.06) (0.05)
MalePctNevMarr 0.23*** 0.14*
  (0.06) (0.06)
FemalePctDiv -0.29** 0.31
  (0.11) (0.29)
TotalPctDiv 0.34** -0.80
  (0.11) (0.50)
PctKids2Par -0.33*** -0.34***
  (0.08) (0.08)
PctWorkMom -0.13*** -0.15***
  (0.03) (0.03)
PctIlleg 0.14** 0.15**
  (0.05) (0.05)
NumImmig -0.17*  
  (0.07)  
PctNotSpeakEnglWell -0.15**  
  (0.05)  
PctLargHouseFam -0.10*  
  (0.04)  
PersPerOccupHous 0.27***  
  (0.07)  
PersPerRentOccHous -0.22***  
  (0.05)  
PctPersOwnOccup -0.63**  
  (0.20)  
PctPersDenseHous 0.28*** 0.17***
  (0.07) (0.04)
HousVacant 0.17*** 0.27***
  (0.04) (0.06)
PctHousOccup -0.06*  
  (0.02)  
PctHousOwnOcc 0.57**  
  (0.19)  
PctVacMore6Mos -0.05* -0.06**
  (0.02) (0.02)
OwnOccLowQuart -0.44**  
  (0.15)  
OwnOccMedVal 0.32*  
  (0.15)  
RentLowQ -0.23*** -0.22***
  (0.04) (0.04)
MedRent 0.27*** 0.25***
  (0.05) (0.05)
MedOwnCostPctIncNoMtg -0.10*** -0.08***
  (0.02) (0.02)
NumStreet 0.20** 0.19**
  (0.06) (0.06)
PctForeignBorn 0.14**  
  (0.04)  
racePctWhite   -0.06
    (0.05)
numbUrban   -0.21*
    (0.08)
MalePctDivorce   0.58*
    (0.24)
MedRentPctHousInc   0.06*
    (0.03)
PctEmplManu   -0.04*
    (0.02)
AsianPerCap   0.03
    (0.02)
agePct16t24   0.14
    (0.09)
PctVacantBoarded   0.04
    (0.03)
PctBSorMore   0.09*
    (0.04)
LemasPctOfficDrugUn   0.03
    (0.02)
MedOwnCostPctInc   -0.05
    (0.03)
PctUsePubTrans   -0.03
    (0.02)
R2 0.69 0.68
Adj. R2 0.68 0.68
Num. obs. 1993 1993
***p < 0.001; **p < 0.01; *p < 0.05

13.4.2 Stepwise backward

De manera similar en el Cuadro 13.3 se presentan los resultados de emplear el algoritmo stepwise backward y tras limpiar las variables no significativas con la corrección HC3 (con un 95% de confianza). Previamente pueden mostrar que estos dos modelos tienen un problema de heteroscedasticidad.

Cuadro 13.3: Modelos seleccionados con el algoritmo stepwise backward con corrección HC3
  Modelo 3 (HC3) (\(R^2\) aj) Modelo 4 (HC3) (AIC)
(Intercept) 0.64*** 0.75***
  (0.10) (0.05)
racepctblack 0.18***  
  (0.04)  
agePct12t29 -0.22**  
  (0.08)  
pctUrban 0.04***  
  (0.01)  
pctWWage -0.26***  
  (0.06)  
pctWFarmSelf 0.05**  
  (0.02)  
pctWInvInc -0.17**  
  (0.06)  
pctWRetire -0.07**  
  (0.03)  
whitePerCap -0.11*  
  (0.05)  
indianPerCap -0.03*  
  (0.02)  
OtherPerCap 0.05**  
  (0.02)  
PctPopUnderPov -0.14**  
  (0.05)  
PctEmploy 0.16**  
  (0.06)  
PctOccupMgmtProf 0.10*  
  (0.04)  
MalePctDivorce 0.32***  
  (0.10)  
MalePctNevMarr 0.20***  
  (0.06)  
TotalPctDiv -0.24*  
  (0.11)  
PctKids2Par -0.36***  
  (0.08)  
PctWorkMom -0.12***  
  (0.03)  
PctIlleg 0.14**  
  (0.05)  
NumImmig -0.17*  
  (0.07)  
PctNotSpeakEnglWell -0.13*  
  (0.05)  
PctLargHouseOccup -0.14**  
  (0.04)  
PersPerOccupHous 0.36***  
  (0.07)  
PersPerRentOccHous -0.22***  
  (0.05)  
PctPersOwnOccup -0.59**  
  (0.20)  
PctPersDenseHous 0.27***  
  (0.07)  
HousVacant 0.16***  
  (0.04)  
PctHousOccup -0.05*  
  (0.02)  
PctHousOwnOcc 0.55**  
  (0.19)  
PctVacMore6Mos -0.04*  
  (0.02)  
OwnOccLowQuart -0.09*  
  (0.04)  
RentLowQ -0.22***  
  (0.04)  
MedRent 0.27***  
  (0.05)  
MedOwnCostPctIncNoMtg -0.09***  
  (0.02)  
NumStreet 0.20**  
  (0.06)  
PctForeignBorn 0.13**  
  (0.04)  
PctFam2Par   -0.58***
    (0.04)
racePctAsian   -0.07**
    (0.02)
PctSameState85   -0.05*
    (0.02)
agePct12t21   -0.10***
    (0.03)
MedYrHousBuilt   0.07***
    (0.02)
racePctWhite   -0.30***
    (0.03)
numbUrban   0.24***
    (0.04)
PersPerFam   0.10**
    (0.03)
blackPerCap   -0.05*
    (0.02)
PctSameCity85   0.07*
    (0.03)
RentHighQ   0.08***
    (0.02)
PctWorkMomYoungKids   -0.08***
    (0.02)
PctHousLess3BR   0.06*
    (0.03)
R2 0.68 0.63
Adj. R2 0.68 0.63
Num. obs. 1993 1993
***p < 0.001; **p < 0.01; *p < 0.05

13.4.3 Combinando forward y backward

Y finalmente, el Cuadro 13.4 se presentan los resultados de emplear el algoritmo combinado y tras limpiar las variables no significativas y la corrección HC3 (con un 95% de confianza). (Recuerden hacer las pruebas de heteroscedasticidad para cada modelo)

Cuadro 13.4: Modelo seleccionado el algoritmo stepwise forward y backward con corrección HC3
  Modelo 5 (HC3) (\(R^2\) aj) Modelo 6 (HC3) (AIC)
(Intercept) 0.53*** 0.63***
  (0.10) (0.10)
racepctblack 0.17*** 0.19***
  (0.03) (0.03)
agePct12t29 -0.20** -0.20**
  (0.07) (0.07)
pctUrban 0.04*** 0.04***
  (0.01) (0.01)
pctWWage -0.25*** -0.26***
  (0.05) (0.05)
pctWFarmSelf 0.05** 0.04**
  (0.02) (0.02)
pctWInvInc -0.14** -0.18***
  (0.05) (0.05)
pctWRetire -0.07* -0.08**
  (0.03) (0.03)
indianPerCap -0.03* -0.03*
  (0.02) (0.02)
OtherPerCap 0.05** 0.05**
  (0.02) (0.02)
PctPopUnderPov -0.11* -0.10*
  (0.05) (0.04)
PctEmploy 0.14* 0.14**
  (0.06) (0.05)
MalePctDivorce 0.12** 0.37***
  (0.04) (0.09)
MalePctNevMarr 0.21*** 0.19***
  (0.06) (0.06)
PctKids2Par -0.29*** -0.38***
  (0.08) (0.08)
PctWorkMom -0.11*** -0.12***
  (0.03) (0.02)
PctIlleg 0.15** 0.16**
  (0.05) (0.05)
NumImmig -0.16*  
  (0.07)  
PctNotSpeakEnglWell -0.15**  
  (0.05)  
PctLargHouseFam -0.11** -0.14***
  (0.04) (0.04)
PersPerOccupHous 0.34*** 0.15**
  (0.07) (0.05)
PersPerRentOccHous -0.24***  
  (0.05)  
PctPersOwnOccup -0.61**  
  (0.20)  
PctPersDenseHous 0.26*** 0.27***
  (0.07) (0.04)
HousVacant 0.15*** 0.23***
  (0.04) (0.05)
PctHousOccup -0.05*  
  (0.02)  
PctHousOwnOcc 0.55**  
  (0.19)  
OwnOccLowQuart -0.12***  
  (0.03)  
RentLowQ -0.21*** -0.22***
  (0.04) (0.04)
MedRent 0.28*** 0.20***
  (0.05) (0.04)
MedOwnCostPctIncNoMtg -0.10*** -0.08***
  (0.02) (0.02)
NumStreet 0.21*** 0.18**
  (0.06) (0.05)
PctForeignBorn 0.14**  
  (0.04)  
TotalPctDiv   -0.32**
    (0.10)
PctVacMore6Mos   -0.04*
    (0.02)
PctLess9thGrade   -0.09**
    (0.03)
NumIlleg   -0.18*
    (0.07)
R2 0.68 0.68
Adj. R2 0.68 0.68
Num. obs. 1993 1993
***p < 0.001; **p < 0.01; *p < 0.05

13.5 Selección del mejor modelo para predecir

En resumen, contamos con 6 modelos (Ver Cuadros 13.2, 13.3 y 13.4). Todos los modelos se encuentran anidados en el modelo 1 y cada modelo es diferente.

Ahora comparemos el comportamiento de los modelos para predecir fuera de muestra. Emplearemos el método de \(k\) iteraciones para hacer la validación cruzada. Empleemos, el paquete caret (Kuhn, 2020) y las funciones trainControl() y train() (tal como lo discutimos en el Capítulo 12). Además empleemos 5 iteraciones (\(k = 5\)).

# se fija el método de k iteraciones
train.control <- trainControl(method = "cv", number = 5)

# se realiza la evaluación fuera de muestra
fold.modelo1= train(formula(modelo1.C2), data = datos.caso2, 
                    method = "lm", trControl = train.control)
fold.modelo2= train(formula(modelo2.C2), data = datos.caso2, 
                    method = "lm", trControl = train.control)
fold.modelo3= train(formula(modelo3.C2), data = datos.caso2, 
                    method = "lm", trControl = train.control)
fold.modelo4= train(formula(modelo4.C2), data = datos.caso2, 
                    method = "lm", trControl = train.control)
fold.modelo5= train(formula(modelo5.C2), data = datos.caso2, 
                    method = "lm", trControl = train.control)
fold.modelo6= train(formula(modelo6.C2), data = datos.caso2, 
                    method = "lm", trControl = train.control)

Resumamos los resultados usando el \(RMSE\) y el \(MAE\) de cada modelo como se presenta en el Cuadro 13.5.

Cuadro 13.5: Metricas de precisión de los 6 modelos con 5 iteraciones
RMSE MAE
Modelo 1 0.133 0.093
Modelo 2 0.134 0.094
Modelo 3 0.133 0.093
Modelo 4 0.142 0.100
Modelo 5 0.134 0.094
Modelo 6 0.134 0.094
Fuente: elaboración propia.

Aunque las métricas no varían mucho entre modelos, de acuerdo con la \(RMSE\) el mejor modelo es el 1 y de acuerdo al MAE el mejor modelo es el 3. Noten que en este caso puede tener sentido penalizar errores de predicción grandes y favorecer los errores pequeños. Así, emplearemos el criterio del \(RMSE\) que sugiere el modelo 1. Este modelo se construyó a partir del algoritmo Forward y empleando el criterio de \(R^2\) ajustado.

Ahora, antes de pasar a emplear el modelo, procedamos a determinar si el mejor modelo tiene o no multicolinealidad. Noten que el \(VIF\) depende de la matriz de varianzas y covarianzas de los estimadores MCO y en presencia de heteroscedasticidad, éstos no son confiables. Así no podemos emplearlos en esta situación. Así que nos concentraremos en la prueba de Belsley et al. (1980) también conocida como la prueba Kappa.

El siguiente código permite hacer la prueba.

XTX <- model.matrix(modelo1.C2)
e <- eigen(t(XTX) %*% XTX)

lambda.1 <- max(e$val)
lambda.k <- min(e$val)
kappa <- sqrt(lambda.1/lambda.k)
kappa
## [1] 270.928

Este estadístico es muy grande (\(\kappa\) = 270.9280298). Esta prueba sugiere la existencia de un problema serio de multicolinealidad. Solucionar este problema no será necesario, pues la multicolinealidad alta nos está ayudando a mejorar el \(R^2\) del modelo y dado que nuestro objetivo es hacer predicciones, no será necesario solucionar el problema.

13.6 Predicciones (escenarios)

El modelo construido nos permite generar los escenarios que estaba buscando la empresa consultora. Recuerden que la pregunta de negocio implicaba tener una “fórmula” que le permita determinar ¿cuál sería la tasa de crímenes violentos por cada 100 mil habitantes bajo diferentes escenarios?.

Por ejemplo, supongamos que queremos tener un escenario en el que todas las variables se encuentran en sus medias.

Recuerden que tenemos heteroscedasticidad y el modelo no cumple el supuesto de normalidad. Por eso deberemos emplear la técnica de bootstrapping para crear los intervalos de confianza para las predicciones. Esto lo podemos hacer empleando la función ic.boot.predic() que construimos en la sección 12.7.

 data.temp = as.data.frame(t(apply(datos.caso2, 2, mean)))

set.seed(123445)
ic.boot.predic(modelo1.C2, data.temp, R= 1000, alpha = 0.05)
##       fit.1    lwr.2.5%   upr.97.5% 
##  0.23798294 -0.04723444  0.59553341

Ahora la empresa de consultoría puede jugar con este modelo mostrando diferentes escenarios. Por ejemplo, supongamos que visitamos el municipio de la fila 120 de la base de datos. Este municipio tiene las siguientes características (\(\bf{x}_{120}\)).

modelo1.C2$model[120,]
##     ViolentCrimesPerPop racepctblack agePct12t29 pctUrban pctWWage pctWFarmSelf
## 120                0.68         0.09        0.34        1     0.49          0.2
##     pctWInvInc pctWRetire medFamInc whitePerCap indianPerCap OtherPerCap
## 120       0.56       0.26      0.55        0.88         0.24        0.36
##     PctPopUnderPov PctEmploy MalePctNevMarr FemalePctDiv TotalPctDiv
## 120           0.24      0.59           0.73         0.81        0.82
##     PctKids2Par PctWorkMom PctIlleg NumImmig PctNotSpeakEnglWell
## 120        0.59       0.43     0.29     0.15                 0.3
##     PctLargHouseFam PersPerOccupHous PersPerRentOccHous PctPersOwnOccup
## 120             0.2             0.05               0.07            0.12
##     PctPersDenseHous HousVacant PctHousOccup PctHousOwnOcc PctVacMore6Mos
## 120             0.25       0.15         0.77          0.08            0.5
##     OwnOccLowQuart OwnOccMedVal RentLowQ MedRent MedOwnCostPctIncNoMtg
## 120              1            1     0.43    0.42                  0.17
##     NumStreet PctForeignBorn
## 120      0.84           0.75

Y la tasa de crímenes violentos por cada 100 mil habitantes de ese municipio es de:

datos.caso2[120,101]
## [1] 0.68

Ahora, se le podría proponer al alcalde una política pública para aumentar el porcentaje de personas de 16 años o más que están empleadas (variable PctEmploy) en 5 puntos porcentuales, de tal manera que pase de 59% a 64%. Entonces, en ese caso se esperaría que la tasa de crímenes violentos por cada 100 mil habitantes de ese municipio cambiaría a:

# se extraen los datos para  el municipio 10 
data.escenario <- modelo1.C2$model[120,]
# se modifica los datos del municipio 10 para ajustarse al escenario
data.escenario["PctEmploy"] <- data.escenario["PctEmploy"] +0.05

# predicción

ic.boot.predic(modelo1.C2, data.escenario, R= 1000, alpha = 0.05)
##   fit.120  lwr.2.5% upr.97.5% 
## 0.5044719 0.1684555 0.8725683

Noten que esto implicaría que el intervalo de confianza de la predicción sigue conteniendo el valor original de la tasa de crímenes violentos por cada 100 mil habitantes, por tanto ese escenario no cambiaría la tasa esperada. Ahora los consultores pueden jugar con esta fórmula para proponer diferentes escenarios a cada alcalde para disminuir la tasa de crímenes violentos o si es del caso saber ¿cuáles municipios no tiene sentido visitar?, pues no se puede hacer mucho para disminuir dicha tasa. De aquí en adelante, la imaginación es el límite.

Referencias

Belsley, D. A., Kuh, E., & Welsch, R. E. (1980). Regression Diagnostics: Identifying Influential Data and Sources of Collinearity. John Wiley & Sons, Inc. https://doi.org/10.1002/0471725153
Koenker, R. (1981). A note on studentizing a test for heteroscedasticity. Journal of Econometrics, 17(1), 107–112.
Kuhn, M. (2020). Caret: Classification and regression training. https://CRAN.R-project.org/package=caret
Redmond, M., & Baveja, A. (2002). A data-driven software tool for enabling cooperative information sharing among police departments. European Journal of Operational Research, 141(3), 660–678.
Trapletti, A., & Hornik, K. (2019). Tseries: Time series analysis and computational finance. https://CRAN.R-project.org/package=tseries

  1. Los datos se pueden descargar de la página web del libro: https://www.icesi.edu.co/editorial/modelo-clasico/↩︎