Previo

Vamos a importar de nuevo de la ECOVID-ML, una nueva fuente desarrollada por INEGI. https://www.inegi.org.mx/investigacion/ecovidml/2020/

Vamos a llamar algunas librerías básicas, el tidyverse (que son muchas librerías) y sjlabelled que nos sirve para el manejo de etiquetas

if (!require("pacman")) install.packages("pacman") # instala pacman si se requiere
## Loading required package: pacman
pacman::p_load(tidyverse, readxl,haven, sjlabelled, foreign, janitor, srvyr) #carga los paquetes necesarios para esta práctica

Por si no tuviéramos cargada la base de datos, la volveremos a cargar

ecovid0420 <- read_dta("./datos/ecovid0420.dta")

También vamos usar la base de Índice de Competitividad Internacional ICI, desarrollado por el Instituto Mexicano de Competitividad. (véase http://imco.org.mx/indices/mexico-suenos-sin-oportunidad/)

ICI_2018 <- read_excel("./datos/ICI_2018.xlsx", sheet = "para_importar")
## New names:
## * `` -> ...128
## * `` -> ...129
## * `` -> ...132
## * `PIB (Paridad de Poder Adquisitivo)` -> `PIB (Paridad de Poder Adquisitivo)...135`
## * `PIB (Paridad de Poder Adquisitivo)` -> `PIB (Paridad de Poder Adquisitivo)...136`
## * ...

Continuación:variables cuantitativas

La media y la desviación estándar

Una de las medidas más comunes para establecer el centro de la distribución es el promedio o media aritmética. La suma de todos los valores de nuestra variable dividida entre el total de observaciones. La media tiene varias propiedades, como que si sumamos todas las desviaciones a este valor, la suma de ellas es cero.

\[ \bar{x}=\frac{\sum_{i=1}^{n} x_i}{n} \]

\[ \mu=\frac{\sum_{i=1}^{N} x_i}{N} \]

Para su cálculo, podemos hacerlo con la función “mean()” de base, pero podemos utilizar el comando “summarize” de dplyr() para obtenerlo como resultado de una operación después de otras a través de los “pipes”

ecovid0420 %>% 
  filter(clase2==1) %>% #filtro de casos 
    summarise(promedio=mean(pe10_1))
## # A tibble: 1 x 1
##   promedio
##      <dbl>
## 1       NA

Antes de calcular la desviación estándar, debemos calcular la varianza. Para ello, de nuevo necesitamos el concepto de desviación: desviación es la diferencia de un valor con respecto a una norma. Por lo general, asumimos esta norma como la media aritmética. Del mismo modo,la media aritmética de los cuadrados de las desviaciones de los valores de la variable con respecto a una constante cualquiera se hace mínima cuando dicha constante coincide con la media aritmética. De ahí que la varianza nos da una medida de distancia promedio, sin el problema que siempre dé cero, como pasaría si no la eleváramos a la cuadrado. A continuación presentamso sus fórmulas para la población y la muestra.

\[s^2=\frac{\sum_{i=1}^{n}(x_i-\bar{x})^2 }{n-1}\]

\[\sigma^2=\frac{\sum_{i=1}^{N}(x_i-\mu)^2 }{N}\]

La medida muestral es diferente a la poblacional en su denominador. Esto proviene de la corrección de Bessel, que corrige el sesgo estadístico en la estimación de la varianza poblacional.

La varianza es una medida muy importante pero díficil de interpretar. Debido que tenemos las unidades originales de nuestra variable: pesos al cuadrado, años al cuadrado, horas al cuadrado. De ahí, que sea importante sacarle raíz cuadrado:

\[s= \sqrt{s^2}=\sqrt{\frac{\sum_{i=1}^{n}(x_i-\bar{x})^2 }{n-1}}\]

La desviación estándar es entonces una medida de dispersión, que nos dice qué tan alejados están los datos de la media. Por lo que aporta mucha más información que la media sola. Por lo general las colocamos juntas:

ecovid0420 %>% 
  filter(clase2==1) %>% #filtro de casos 
    summarise(media=mean(pe10_1), # media
              sd=sd(pe10_1), # desviación estándar
              var=var(pe10_1)) #varianza
## # A tibble: 1 x 3
##   media    sd   var
##   <dbl> <dbl> <dbl>
## 1    NA    NA    NA

Estas funciones brindan las estimaciones muestrales. Si necesitas las estimaciones poblacionales, tenemos que hacer un artilugio de multiplicar por (N-1/N) para que se elimine el denominador (N-1) y quede multiplicado por N. Esto lo podemos hacer con la función “length()” para nuestro vector de análisis pe10_1

Asumiendo que tuviéramos una población y no una muestra

ecovid0420 %>% 
  filter(clase2==1) %>% #filtro de casos 
    summarise(media=mean(pe10_1), # media
              var.p=var(pe10_1)*(length(pe10_1)-1)/length(pe10_1), #varianza
              sd.p=sqrt(var(pe10_1)*(length(pe10_1)-1)/length(pe10_1))) # desviación estándar
## # A tibble: 1 x 3
##   media var.p  sd.p
##   <dbl> <dbl> <dbl>
## 1    NA    NA    NA

Las diferencias entre las estimaciones poblacionales y las muestrales son muy pocas porque nuestra muestra es grande y las diferencias en el denominador generan pocos cambios.

El resumen de cinco números y los gráficos de caja y brazos

La media es un medida muy popular, pero tiene un problema y es que está afectado por los valores atípicos.

Otra medida popular y más robusta a los “outliers” es la mediana. Ella representa el valor de la variable en posición central en un conjunto de datos ordenados. Es decir supera al 50% de los casos y su valor es superado el otro 50% restante.

ecovid0420 %>% 
  filter(clase2==1) %>% #filtro de casos 
    summarise(mediana=median(pe10_1))
## # A tibble: 1 x 1
##   mediana
##     <dbl>
## 1      NA

El que la media supere a la mediana, da información sobre el sesgo a la derecha que mantiene la distribución. Si los valores son iguales o muy cercanos, seguro estamos ante una distribución bastante simétrica; mientras que si la mediana supera a la media, ello da cuenta que existen valores a la izquierda de la distribución que la están sesgando, de ahí que podemos aducir que hay un sesgo negativo.

Cuando tenemos esta situación la media no es tan representativa y para comprender más nuestra distribución, necesitamos medidas que acompañen a una media de centro como la mediana. De ahí proviene la necesidad del resumen de cinco números:

ecovid0420 %>% 
  filter(clase2==1) %>% #filtro de casos 
    with(summary(pe10_1))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    1.00   16.00   35.00   32.65   48.00   99.00     431

Con el comando “summary()”, obtenemos estos seis números (se agrega la media), el resto es lo que conocemos como el resumen de cinco números. Incluye el mínimo y máximo en los extremos y otras dos medidas de posición: el cuartil 1 y el cuartil 3. El cuartil 1 es una medida de posición igual que la mediana que separa la población en un 25% inferior y un 75% superior; mientras que el cuartil 3 separa a la población en un 75% inferior y un 25% superior. Estas medidas nos dan un idea de cómo se distribuye nuestra variable, pero también son la base de unos de los gráficos más famosos: el gráfico de caja y brazos o “boxplot”.

El gráfico de caja y brazos (o caja y bigotes), también toma en cuenta el concenpto de rango intercuartílico (RIC), que es la diferencia entre el cuartil 1 y el cuartil 3, es decir, establece el rango donde se concentra el 50% de los datos.

\[ RIC= Q_3 - Q_1 \]

Otro concepto que utiliza el gráfico es el de atípicos, pero más allá de lo que ya habíamos hecho anteriormente (visualmente), propone unas medidas de límites inferior y superior:

\[ L_{inferior}= Q_1- 1.5 * RIC \]

\[ L_{superior}= Q_3+ 1.5 * RIC \]

Cualquier dato será atípico si es menor al \(L_{inferior}\) y mayor al \(L_{superior}\).

Para graficarlo en R tenemos:

ecovid0420 %>% 
  filter(clase2==1) %>% #filtro de casos 
    with(boxplot(pe10_1))

Tenemos tantos casos atípicos y un sesgo positivo tan grande que casi no podemos verlo.

Análisis bivariado: cuanti-cuali

Para hacer estadísticos para grupos, podemos agregar a nuestros códigos anteriores una línea en nuestros “pipes” que hemos utilizado anteriormente, y además de hacer esto revisáremos qué hacer cuando tenemos missings o valores perdidos:

ecovid0420 %>% 
  filter(clase2==1) %>% #vamos a quitar los que no reciben ingresos
  group_by(as_label(pb1))  %>%  # hace el agrupamiento para la variables categóricas
  summarise(media=mean(pe10_1, na.rm=T), # checa que ponemos que nos remueva los missings 
            sd=sd(pe10_1, na.rm=T),
            mediana=median(pe10_1, na.rm=T))
## `summarise()` ungrouping output (override with `.groups` argument)
## # A tibble: 2 x 4
##   `as_label(pb1)` media    sd mediana
##   <fct>           <dbl> <dbl>   <dbl>
## 1 Hombre           35.6  19.5      40
## 2 Mujer            28.6  19.0      28

Estadísticos con datos expandidos

# Conteo de casos 
ecovid0420 %>% 
  mutate(clase2=as_label(clase2)) %>%  # para usar etiquetas
         group_by(clase2) %>%  # agrupa la base
           tally() # da conteos para grupos
## # A tibble: 4 x 2
##   clase2     n
##    <dbl> <int>
## 1      1  2893
## 2      2   202
## 3      3  1080
## 4      4  1418

La ventaja de “tally()”, es que podemos ponerle un peso a su interior, y en lugar de contar casos puede sumar variables, tal como sucede con el factor de expansión:

# Conteo de factor de expansión
ecovid0420 %>% 
  mutate(clase2=as_label(clase2)) %>%  # para usar etiquetas
         group_by(clase2) %>%  # agrupa la base
           tally(fac_per) # suma el factor de expansión
## # A tibble: 4 x 2
##   clase2        n
##    <dbl>    <dbl>
## 1      1 32891935
## 2      2  2060938
## 3      3 13613258
## 4      4 19618299

Estos valores ya expandidos, coinciden con los tabulados de Inegi. Siempre es un buena práctica revisar nuestros resultados contra los tabulados publicados, para revisar si estamos calculando los elementos correctamente, o saber que nuestras diferencias se basan en algún cambio de criterio.

Algunas opciones de “janitor”, se pueden obtener para tener los totales como una nueva fila:

# Conteo de factor de expansión
ecovid0420 %>% 
  mutate(clase2=as_label(clase2)) %>%  # para usar etiquetas
         group_by(clase2) %>%  # agrupa la base
           tally(fac_per) %>%  # suma el factor de expansión
             adorn_totals("row") # agrega una fila con totales 
##  clase2        n
##       1 32891935
##       2  2060938
##       3 13613258
##       4 19618299
##   Total 68184430

Y finalmente, podemos calcular proporciones con “adorn_percentages”" y luego ponerle formato de “%”.

ecovid0420 %>% 
  mutate(clase2=as_label(clase2)) %>% 
  group_by(clase2) %>% 
  tally(fac_per) %>% 
   adorn_totals("row") %>% 
     adorn_percentages("all")  %>% 
      adorn_pct_formatting()
##  clase2      n
##       1  48.2%
##       2   3.0%
##       3  20.0%
##       4  28.8%
##   Total 100.0%

Para la media, las estimaciones de la muestra pueden diferir de acuerdo al factor de expansión que es parte del diseño muestral. En base tenemos la función que nos calcula la media expandida o ponderada:

ecovid0420 %>% 
  filter(clase2==1) %>% #sólo trabajadores
  group_by(as_label(pb1))  %>%  # hace el agrupamiento para la variables categóricas
  summarise(media=mean(pe10_1, na.rm=T),
            media_ponderada=
              weighted.mean(pe10_1, na.rm=T, w=fac_per)) # checa que ponemos que nos remueva los missings 
## `summarise()` ungrouping output (override with `.groups` argument)
## # A tibble: 2 x 3
##   `as_label(pb1)` media media_ponderada
##   <fct>           <dbl>           <dbl>
## 1 Hombre           35.6            35.8
## 2 Mujer            28.6            29.2

Vemos que la media ponderada es menor. Por lo que es importante cuando hacemos referencia a la población, retomar la estructura del diseño muestral

Para otras medidas, y también tiene capacidad de introducir el diseño muestral completo, podemos utilizar el paquete “srvyr”

ecovid0420 %>% 
    as_survey_design(weights = fac_per) %>% #establece los pesos
    mutate(pb1=as_label(pb1)) %>% #para mejor lectura de las etiquetas
      filter(clase2==1) %>% #filtro
        group_by(pb1)  %>%  # hace el agrupamiento para la variables categóricas
            summarise(media_ponderada=
                      survey_mean(pe10_1, na.rm=T))
## # A tibble: 2 x 3
##   pb1    media_ponderada media_ponderada_se
##   <fct>            <dbl>              <dbl>
## 1 Hombre            35.8              0.628
## 2 Mujer             29.2              0.743

Análisis bivariado: dos variables cuantitativas

Mide la fuerza de la relación, lineal si es de Pearson. Debemos recordar que la correlación puede tener un valor:

  • 1 es una correlación positiva perfecta
  • 0 es sin correlación (los valores no parecen vinculados en absoluto)
  • -1 es una correlación negativa perfecta)

    • Elementos importantes

    • El coeficiente de correlación, r, no hace ninguna distinción entre las variables explicativas y dependientes. No hace ninguna diferencia cuál variable por X y cuál se llama Y en el cálculo de la correlación.

    • r utiliza los valores estandarizados de las observaciones, r no cambia cuando cambiamos las unidades de medida de x, y, o ambos. La medición de la altura en pulgadas en lugar de centímetros y el peso en libras en lugar de kilogramos no cambia la correlación entre la altura y el peso. La correlación, r, en sí no tiene unidad de medida; sólo es un númer

    Aplicación con R

    Los datos de ICI Están un poco sucios. Primero los vamos a limpiar:

    ICI_2018<-clean_names(ICI_2018)

    Para observar la relación que pudiera existir entre la tasa de homicidios y la producción per cápita

    cor(ICI_2018$homicidios_dolosos, ICI_2018$pib_per_capita_139)
    ## [1] -0.4788776
    cor(ICI_2018$homicidios_dolosos, ICI_2018$pib_per_capita_139, use="pairwise")
    ## [1] -0.4788776

    ¡La correlación es negativa!

    Tipos de correlación

    Por default está la correlación de Pearson, pero en realidad podemos obtener otros tipos

    #Pearson - default
cor(ICI_2018$homicidios_dolosos, ICI_2018$pib_per_capita_139, 
    use = "pairwise", method = "pearson")
    ## [1] -0.4788776
    #Tau-Kendall
cor(ICI_2018$homicidios_dolosos, ICI_2018$pib_per_capita_139, 
    use = "pairwise", method = "kendall")
    ## [1] -0.3754153
    #Rho-Spearman
cor(ICI_2018$homicidios_dolosos, ICI_2018$pib_per_capita_139, 
    use = "pairwise", method = "spearman")
    ## [1] -0.5203866

    Con dplyr:

    ICI_2018 %>% 
  with(cor(homicidios_dolosos,pib_per_capita_139))
    ## [1] -0.4788776

    Etiquetado de variables

    Adenda de lo visto en el curso

    #Revisamos los valores de la variable

table(ecovid0420$pos_ocu)
    ## 
##    0    1    2    3 
## 2700 2050  796   47
    # De acuerdo a la documentación de la base de datos tenemos:

ecovid0420$pos_ocu<-set_labels(ecovid0420$pos_ocu,
                               labels=c("No Aplica", # para 0
                                        "Subordinadxs", # para el 1
                                        "Cuenta Propia", # para el 2
                                        "Sin pago")) # Para el 3

    Así establecemos estiquetas de valores.

    Para etiquetar variable:

    set_label(ecovid0420$pos_ocu)<-"Posición en la ocupación"

    Comando “paste”, pega el valor de algún objeto. Así verificamos que tengamos esta información

    paste(get_label(ecovid0420$pos_ocu))
    ## [1] "Posición en la ocupación"

    Ejercicio

    • Realice un análisis descriptivo de dos variables de alguna de las bases de datos (ICI_2018 o ecovid0420)

    • PRIMERO: realice el análisis descriptivo de manera UNIVARIADA de cada una de las variables que escoja.

    • SEGUNDO: realice el análisis descriptivo BIVARIADO.

    Mande el Envíe el script utilizando la liga: https://forms.gle/NECSbcCRjnTBRTMP6