Cómo analizar la ENVIPE en R

De acuerdo con el INEGI, a nivel nacional 73% percibe a su estado como “inseguro”, pero si uno analiza la base de datos, la cifra es 68%.¿Por qué existen esas diferencias?

[También revisa Cómo analizar la ENVIPE en Stata]
La Encuesta Nacional de Victimización y Percepción sobre Seguridad Pública (ENVIPE) del INEGI es una fuente de datos muy valiosa para los expertos en seguridad pública. Sin embargo, es común que al analizar la base de datos de la encuesta, algunos no puedan replicar los resultados que proporciona el INEGI. Por poner un ejemplo, de acuerdo a los tabulados del INEGI, a nivel nacional 24.5% de los individuos percibe a su estado como “seguro” y 73.2% como “inseguro”; pero si uno realiza un tabulado directamente de la base de datos, los porcentajes son 29.3% y 68%, respectivamente.¿Por qué existe esa discrepancia?

La razón es que la ENVIPE es una encuesta recolectada a través de un diseño muestral complejo. A diferencia de lo que asumen las herramientas estadísticas más comunes, las observaciones de la encuesta no tienen la misma probabilidad de selección ni son independientes. Omitir el diseño de la muestra puede arrojar estimaciones sesgadas y subestimar la varianza de los estimadores de interés.

El INEGI advierte que el análisis de la ENVIPE “es responsabilidad exclusiva del usuario”, pero la gran mayoría de las personas con sólidas habilidades para análisis de datos no están necesariamente familiarizados con estimaciones para muestras complejas.

Mi propósito en este post es presentar algunos fragmentos de código en R para replicar dos tablas de los Tabulados Básicos del INEGI. El código puede servir de referencia para otros análisis. El post presenta tres ejemplos: uno sobre percepción de seguridad, otro sobre incidencia, prevalencia y tazas del delito de secuestro, y el último sobre estimación de homicidios. El primer ejemplo hace estimaciones de proporciones a nivel individual para subpoblaciones (entidades federativas). Los otros dos tratan sobre estimaciones a nivel del hogar para totales, y razones. En mi opinión, estos ejemplos cubren los aspectos más generales para poder analizar prácticamente cualquier pregunta de la encuesta.

Como en cualquier análisis, es conveniente que el lector tenga nociones sobre el análisis de muestras complejas. Algunos textos introductorios son:

Ejemplo 1: Percepción de inseguridad en los estados

La tabla de abajo muestra las estimaciones a nivel estatal del número (“absolutos”) y porcentaje (“relativos”) de individuos que perciben a su estado como seguro o inseguro. Como se observa en la primera fila de la tabla, a nivel nacional 24.5% percibe a su estado como seguro y 73.2% como inseguro.

V_percepcion_seguridad_2015 5.10

Tabla 1. Percepción de inseguridad por entidad federativa

Ahora analicemos la base de datos disponible en este link. Primero carguemos los datos en R:

rm(list = ls(all = TRUE)) #clear workspace
library(foreign)
data = read.dbf('TPer_Vic1.DBF', as.is=TRUE)

Ahora tabulemos las respuestas de la pregunta sobre percepción de seguridad (AP4_3_3; véase el documento sobre la estructura de la base de datos). Como se observa, los porcentajes son distintos a los de la Tabla 1: 29.3% se siente seguro y 68% inseguro.

> round(prop.table(table(data$AP4_3_3))*100,1)
   1    2    9 
29.3 68.0  2.7 

La diferencia radica en que los tabulados del INEGI incorporan las características del diseño muestral de la encuesta. Para fines de estimación, las más importantes son:

  • Conglomerado: La variable UPM_DIS identifica el conglomerado o unidad primaria de muestreo.
  • Estrato: La variable EST_DIS identifica el estrato al que pertenece la unidad primaria de muestreo.
  • Ponderador a nivel individual: La variable FAC_ELE contiene el factor de expansión (ponderador) que incorpora la probabilidad de selección del individuo y ajuste por no-respuesta, omisión de delitos y proyección de la población. De acuerdo con la metodología de la encuesta, este ponderador es “requerido para estimar resultados de las preguntas de percepción de la seguridad pública y la victimización de la población de 18 años y más.”

Estas variables son incorporadas al análisis con la función svydesign:

library(survey)
design = svydesign(id=~UPM_DIS,strata=~EST_DIS, weights=~data$FAC_ELE, data=data)

Luego, se pueden obtener los estimadores con la función svymean. Como se observa abajo, los porcentajes son iguales a los de la Tabla 1.

> svymean(~AP4_3_3, design)
             mean     SE
AP4_3_31 0.244594 0.0022
AP4_3_32 0.732087 0.0023
AP4_3_39 0.023319 0.0008

También podemos obtener el número “absoluto” de individuos que consideran seguro o inseguro a su estado con el comando svytotal. De nuevo, como se observa en la tabla de abajo, las cifras coinciden con las de la Tabla 1:

> svytotal(~AP4_3_3, design)
            total     SE
AP4_3_31 19790834 190069
AP4_3_32 59235506 322330
AP4_3_39  1886823  61281

Para obtener el porcentaje y número absoluto de individuos por estado, es necesario hacer inferencias para subpoblaciones con la función svyby:

ENT = substr(data$UPM, 1, 2) #Entidad
total.edo = svyby(~AP4_3_3, by=ENT, design=design, svytotal)
prop.edo = svyby(~AP4_3_3, by=ENT, design=design, svymean)

Así, por ejemplo, el objeto prop.edo contiene –por estado– la proporción de individuos que se sienten seguros, inseguros, o rechazaron contestar la pregunta. El objeto también incluye los errores estándar de cada proporción. Como se observa abajo, los porcentajes coinciden con los de la Tabla 1:

> head(prop.edo)
   by  AP4_3_31  AP4_3_32    AP4_3_39 se.AP4_3_31 se.AP4_3_32 se.AP4_3_39
01 01 0.5464177 0.4317832 0.021799133  0.01338317  0.01346271 0.004186313
02 02 0.4314004 0.5324062 0.036193386  0.01241809  0.01220218 0.003954120
03 03 0.3675080 0.6182497 0.014242278  0.01074400  0.01181397 0.003078726
04 04 0.4223809 0.5367231 0.040895998  0.01348049  0.01307536 0.005051368
05 05 0.2182016 0.7492093 0.032589059  0.01145402  0.01214261 0.004358251
06 06 0.4248653 0.5653547 0.009779912  0.01163603  0.01146453 0.002362011

Ahora repliquemos la Tabla 1. Los únicos datos que nos faltan incluir es el total de la población mayor de 18 años por estado. Para hacerlo, sumemos el ponderador individual por estado:

pob.edo = aggregate(data$FAC_ELE, by=list(ENT), sum)

Juntando toda la información, podemos elaborar una tabla que es idéntica a la del INEGI:

> cbind(Poblacion.18.y.mas = pob.edo[,2], 
+       Seguro.Absolutos = total.edo[,2], 
+       Seguro.Relativos = round(prop.edo[,2]*100, 1), 
+       Inseguro.Absolutos = total.edo[,3], 
+       Inseguro.Relativos = round(prop.edo[,3]*100, 1))
      Poblacion.18.y.mas Seguro.Absolutos Seguro.Relativos Inseguro.Absolutos Inseguro.Relativos
 [1,]             826776           451765             54.6             356988               43.2
 [2,]            2324817          1002927             43.1            1237747               53.2
 [3,]             513401           188679             36.8             317410               61.8
 [4,]             609644           257502             42.2             327210               53.7
 [5,]            1931814           421525             21.8            1447333               74.9
 [6,]             495710           210610             42.5             280252               56.5
 [7,]            3171336          1244180             39.2            1732141               54.6
 [8,]            2452925           614762             25.1            1805410               73.6
 [9,]            6885392          1406594             20.4            5403354               78.5
[10,]            1145726           342569             29.9             778912               68.0
[11,]            3683034          1145236             31.1            2386772               64.8
[12,]            2218531           254308             11.5            1924735               86.8
[13,]            1927642           670860             34.8            1196991               62.1
[14,]            5261794          1520348             28.9            3664604               69.6
[15,]           11458236           910688              7.9           10382800               90.6
[16,]            3015726           532797             17.7            2406763               79.8
[17,]            1319461           161878             12.3            1138051               86.3
[18,]             815389           326954             40.1             452675               55.5
[19,]            3519863           963442             27.4            2489875               70.7
[20,]            2586646           511555             19.8            2010110               77.7
[21,]            3943603          1125434             28.5            2657074               67.4
[22,]            1335421           613569             45.9             676465               50.7
[23,]            1034056           361165             34.9             630525               61.0
[24,]            1774437           495586             27.9            1227394               69.2
[25,]            2000620           536895             26.8            1409542               70.5
[26,]            1973052           722915             36.6            1232477               62.5
[27,]            1558732           155998             10.0            1385646               88.9
[28,]            2405118           264643             11.0            2090075               86.9
[29,]             834657           330639             39.6             494021               59.2
[30,]            5472670           963321             17.6            4407725               80.5
[31,]            1424466           927914             65.1             481647               33.8
[32,]             992468           153576             15.5             802782               80.9

También podemos graficar los resultados con la siguiente instrucción:

library(arm)
abbrev = c('Ags','BC','BCS','Camp','Chis','Chih','Coah','Col','DF','Dgo','Gto',
            'Gro','Hgo','Jal','Méx','Mich','Mor','Nay','NL','Oax','Pue','Qro',
            'QR','SLP','Sin','Son','Tab','Tams','Tlax','Ver.','Yuc','Zac')
coefplot(rev(prop.edo[,3]*100), rev(prop.edo[,5]*100), varnames=rev(abbrev), main='% Inseguro')
abline(v=svymean(~AP4_3_3, design)[2]*100, lty='dashed')

Los puntos de la gráfica de abajo muestran el % de individuos que perciben a su estado como inseguro, mientras que las barras horizontales representan intervalos de confianza al 95%. La línea vertical punteada indica el porcentaje a nivel nacional.

inseguro_estatal

Ejemplo 2: Estimación de secuestro

El segundo ejemplo trata sobre la estimación de la incidencia y prevalencia del delito de secuestro. La Tabla de abajo muestra los resultados del INEGI:

secuestro_2015

Tabla 2. Hogares, victimas y número de secuestros (ENVIPE, 2015)

Comencemos por cargar la base de datos en R:

rm(list = ls(all = TRUE)) #clear workspace
library(foreign)
data = read.dbf('TPer_Vic2.DBF', as.is=TRUE)

La estimación del delito de secuestro es un poco compleja. Primero, el encuestador pregunta al individuo seleccionado si algún miembro del hogar fue secuestrado antes de 2014.  Luego le pregunta si un miembro del hogar fue secuestrado durante 2014. Naturalmente, esta última pregunta es la de interés; el propósito de la pregunta previa es reducir la tendencia del informante a reportar delitos que ocurrieron antes del periodo de referencia (telescopeo). Luego, el encuestador verifica con una serie de preguntas si las víctimas de secuestro efectivamente vivían y compartían los alimentos en el hogar donde se realizó la entrevista. Podemos usar esta serie de preguntas para construir nuestra medición validada del número de delitos de secuestro, el número de víctimas de secuestro, y el número de hogares donde hubo al menos una víctima de secuestro:

#6.12 ¿Cuántas veces sufrió el secuestro la (NÚMERO ORDINAL) víctima del hogar?
foo = data.matrix(data[,c('AP6_12_1','AP6_12_2','AP6_12_3','AP6_12_4',
                           'AP6_12_5','AP6_12_6','AP6_12_7')])
data$hogares = as.numeric(rowSums(foo, na.rm=TRUE) > 0)
data$victimas = rowSums(foo > 0, na.rm=TRUE)
data$secuestros = rowSums(foo, na.rm=TRUE)

Luego, declaramos el diseño muestral con la función svydesign como en la sección anterior. Sin embargo, utilizaremos un ponderador distinto: FAC_HOG. De acuerdo con el documento metodológico de la encuesta, éste es el “ponderador que se utiliza para estimar resultados de las preguntas que se refieren a los hogares y la población en general“.

#Ponderador y diseno
library(survey)
design = svydesign(id=~UPM_DIS,strata=~EST_DIS, weights=~data$FAC_HOG, data=data)

#Estimacion
point.hogares = svytotal(~hogares, design) #Estimacion: 88,523
point.victimas =svytotal(~victimas, design) #Estimacion: 99,747
point.secuestros =svytotal(~secuestros, design) #Estimacion: 102,883
tasa.hogares = svyratio(~hogares, rep(1, nrow(data)), design)

Luego, obtenemos nuestros estimadores de interés. Nótese que son idénticos a las cifras de la Tabla 2 (en amarillo):

> print(point.hogares)
        total   SE
hogares 88523 8741
> print(point.victimas)
         total    SE
victimas 99747 10070
> print(point.secuestros)
            total    SE
secuestros 102883 10199
> print(tasa.hogares)
Ratio estimator: svyratio.survey.design2(~hogares, rep(1, nrow(data)), design)
Ratios=
                   
hogares 0.002723332
SEs=
                [,1]
hogares 0.0002688849

La Tabla 2 también incluye intervalos de confianza, margen de error, y coeficiente de variación. Podemos replicar el contenido de la tabla con el siguiente código:

#Replicar Tablas
#---------------
#Intervalos de confianza
ci.hogares = confint(point.hogares, level=0.9)
ci.victimas = confint(point.victimas, level=0.9)
ci.secuestros =  confint(point.secuestros, level=0.9)

#Coeficientes de variacion
cv.hogares = cv(point.hogares)
cv.victimas = cv(point.victimas)
cv.secuestros = cv(point.secuestros)

tabla1 = rbind(
  c(point.hogares[1], point.victimas[1], point.secuestros[1]),
  c(paste(round(ci.hogares), collapse='-'), paste(round(ci.victimas), collapse='-'), paste(round(ci.secuestros), collapse='-')),
  round(c((point.hogares[1]-ci.hogares[1])/point.hogares[1],
          (point.victimas[1]-ci.victimas[1])/point.victimas[1],
          (point.secuestros[1]-ci.secuestros[1])/point.secuestros[1])*100),
  round(c(cv.hogares, cv.victimas, cv.secuestros)*100)
)
rownames(tabla1) = c('Estimacion','Intervalo','Margen','CV')

El resultado es el siguiente:

> print(tabla1)
           hogares        victimas       secuestros    
Estimacion "88523"        "99747"        "102883"      
Intervalo  "74145-102901" "83183-116311" "86107-119659"
Margen     "16"           "17"           "16"          
CV         "10"           "10"           "10"

Podemos hacer un ejercicio similar para estimar la tasa de hogares víctimas de secuestro (columna azul en la Tabla 2):

cv.tasa.hogares = cv(tasa.hogares)
ci.tasa.hogares = confint(tasa.hogares, level=0.9)
tabla2 = c(round(tasa.hogares$ratio*100000), 
          paste(round(ci.tasa.hogares*100000), collapse='-'),
          round(((tasa.hogares$ratio-ci.tasa.hogares[1])/tasa.hogares$ratio)*100),
          round(cv.tasa.hogares*100)
)
names(tabla2) = c('Estimacion','Intervalo','Margen','CV')

La tabla que resulta es la siguiente:

> print(as.matrix(tabla2))
           [,1]     
Estimacion "272"    
Intervalo  "228-317"
Margen     "16"     
CV         "10"     

Ejemplo 3: Estimación de Homicidios

En 2013 el INEGI reportó la estimación del número de homicidios a partir de los datos de la ENVIPE, pero ha dejado de hacerlo para la ENVIPE 2014 y 2015. Sin embargo, podemos hacerlo de manera similar a la estimación de secuestro con las siguientes líneas de código:

rm(list = ls(all = TRUE)) #clear workspace

#Leer datos
library(foreign)
data = read.dbf('TPer_Vic2.DBF', as.is=TRUE)

#Validacion de homicidios
foo = data.matrix(data[,c('AP6_21_1','AP6_21_2','AP6_21_3','AP6_21_4',
                           'AP6_21_5','AP6_21_6','AP6_21_7')])

#Calcular número de homicidios y hogares victima
data$homicidios = rowSums(foo == 1, na.rm = TRUE)
data$hogares = data$homicidios > 0

#Disenio muestral
library(survey)
design = svydesign(id=~UPM_DIS,strata=~EST_DIS, weights=~as.numeric(data$FAC_HOG), data=data)

Como se observa, el número estimado de homicidios es 24,917 su intervalo de confianza (95%) es 14,476-35,358. La estimación del número de hogares con al menos una víctima de homicidio es 22,824 (14,568-31,080).

> #Estimacion
> svytotal(~homicidios, design) #Estimacion: 24,917
           total     SE
homicidios 24917 5326.9
> confint(svytotal(~homicidios, design))
              2.5 %   97.5 %
homicidios 14476.44 35357.56
> 
> svytotal(~hogares, design) #Estimacion: 22,824
                total       SE
hogaresFALSE 32482575 105568.8
hogaresTRUE     22824   4212.2
> confint(svytotal(~hogares, design))
                   2.5 %      97.5 %
hogaresFALSE 32275663.89 32689486.11
hogaresTRUE     14568.27    31079.73

4 thoughts on “Cómo analizar la ENVIPE en R

  1. Pingback: Cómo analizar la ENVIPE en Stata | Políticamente Correcto

  2. Muchas gracias por el excelente articulo.
    Estoy intentado estimar los Intervalos de confianza
    para el Módulo de delitos (“TMod_Vic.dbf”). Pero solo tengo casi lo mismos resultado como el INEGI. Por ejemplo por mi CI es 9114345.06- 10148596.94, de INEGI es: 9 110 178-10 152 764. Tengo el mismo problema por 2014. Mi codigo es: setwd(“D:/Bjoern_backup/ENVIPE/ENVIPE2015”)
    require(foreign)
    require(survey)

    vic <- read.dbf("TMod_Vic.dbf", as.is=TRUE)

    vic <- subset(vic, BPCOD != "03") # quitar pintas, etc.

    design <- svydesign(id = ~UPM_DIS,
    weights = ~vic$FAC_DEL,
    strata= ~EST_DIS,
    data = vic)

    options(survey.lonely.psu="adjust")

    delitos=svytotal(~ BPCOD, design)
    confint(svytotal(~ BPCOD, design), level = 0.90, df = degf(design))
    sum(delitos, na.rm = FALSE)

    • Hola!
      En lugar de usar subset, sustituye por NAs:

      library(foreign)
      data = read.dbf('TMod_Vic.DBF', as.is=TRUE)
      table(as.numeric(data$BPCOD))
      data$BPCOD[data$BPCOD=='03'] = NA
      library(survey)
      design = svydesign(id=~UPM_DIS,strata=~EST_DIS, weights=~data$FAC_DEL, data=data)
      options(survey.lonely.psu='adjust')
      confint(svytotal(~BPCOD, design, na.rm=TRUE), level=0.90)
      

      Como verás, los resultados que arroja son idénticos a los del INEGI:

      > confint(svytotal(~BPCOD, design, na.rm=TRUE), level=0.90)
                     5 %        95 %
      BPCOD01  505026.79   585077.21
      BPCOD02 3279706.34  3537427.66
      BPCOD04 1970689.40  2129486.60
      BPCOD05 9110177.77 10152764.23
      BPCOD06 1105537.39  1280200.61
      BPCOD07 1609691.48  1879058.52
      BPCOD08 1480878.20  1915885.80
      BPCOD09 7608697.20  8330472.80
      BPCOD10 2988501.33  3661304.67
      BPCOD11  975208.23  1343279.77
      BPCOD12   59432.25    99093.75
      BPCOD13  563718.58   926621.42
      BPCOD14   45646.39   135327.61
      BPCOD15   47572.58    81579.42
      

Leave a Reply