Laboratorio I: Conceptos de estadiśtica descriptiva

Introducción

Esta práctica inicial de laboratorio tiene como objetivo realizar una primera aproximación al Lenguaje R, utilizando el enfoque de análisis exploratorio de datos sobre un dataset, a efectos de repasar conceptos fundamentales de estadística descriptiva.

Estadistica Descriptiva

A partir del dataset MPI_subnational.csv (Multidimensional Poverty Measures), se solicita trabajar sobre las siguientes consignas:

Fuente de datos

1. Exploración de datos

Explore y explique en que consiste el dataset utilizando herramientas de exploración de datos.

1.a. Releve las características de los atributos.

Para comenzar cargamos las librerías necesarias para esta práctica:

library(pacman)
pacman::p_load(tidyverse, modeest, WVPlots, DT, plotly, GGally)
options(warn=-1)

Luego cargamos el dataset a analizar:

read_csv <- function(name) {
  read.csv(
    paste('https://raw.githubusercontent.com/adrianmarino/dm/master/datasets/', name, '.csv', sep=''), 
    header = TRUE, 
    sep = ','
  )
}

dataset = read_csv('MPI_subnational')

Renombramos las columnas del dataset a nombres mas cortos y reemplazamos valores nulos por ceros:

dataset = dataset %>%
  rename(
    country_code   = ISO.country.code,
    country        = Country,
    sub_nat_region = Sub.national.region,
    world_region   = World.region,
    mpi_national   = MPI.National,
    mpi_regional   = MPI.Regional,
    hc_regional    = Headcount.Ratio.Regional,
    iod_regional   = Intensity.of.deprivation.Regional
  ) %>%
  replace_na(list(
    iod_regional = 0,
    mpi_national = 0,
    hc_regional  = 0,
    mpi_regional = 0
  ))

Creamos una nueva tabla con las columnas numéricas normalizadas, que usaremos para comparar distribuciones:

numeric_columns = c("mpi_national", "mpi_regional", "hc_regional", "iod_regional")

normalize <- function(x) { (x - min(x)) / ( max(x) - min(x) ) }

dataset.n <- dataset %>% mutate_at(numeric_columns, normalize)

summary(dataset.n)

##  country_code         country          sub_nat_region     world_region      
##  Length:984         Length:984         Length:984         Length:984        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##   mpi_national     mpi_regional      hc_regional      iod_regional   
##  Min.   :0.0000   Min.   :0.00000   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:0.1002   1st Qu.:0.07124   1st Qu.:0.1260   1st Qu.:0.5455  
##  Median :0.2805   Median :0.20833   Median :0.3429   Median :0.6008  
##  Mean   :0.3307   Mean   :0.28405   Mean   :0.4059   Mean   :0.6210  
##  3rd Qu.:0.4958   3rd Qu.:0.45901   3rd Qu.:0.6740   3rd Qu.:0.6838  
##  Max.   :1.0000   Max.   :1.00000   Max.   :1.0000   Max.   :1.0000

Defininos las funciones que vamos a luego utilizar para graficar:

pie_plot <- function(data, title) {
  values = table(data)
  labels = paste(names(values), " (", values, ")", sep="")

  pie(values, labels = labels, main=title)  
}

ggpie_plot <- function(df, seg_label="", sum_label="") {
  p <- ggplot(df, aes(x="", y=Frequency, fill=Value)) +
    geom_bar(stat="identity", width=1, color="white") +
    coord_polar("y", start=0) +
    geom_text(
      aes(label = Frequency), 
      position = position_stack(vjust = 0.5), 
      color = "white"
    )+
    labs(
      x = NULL, 
      y = NULL, 
      fill = seg_label, 
      title = paste(sum_label, "por", seg_label, sep=" ")
    )  +
    theme_void()
  p
}

gplot_hist <- function(
  values,
  ylab = "Frecuencia",
  name = "",
  line_size=1.05,
  truncated_mean_value=0.05,
  binwidth=1,
  linetype="solid"
) {
  df = as.data.frame(values)

  p <- ggplot(df, aes(x=values)) +
    geom_histogram(aes(y=..density..), color="darkblue", fill="lightblue", binwidth=binwidth) +
    geom_density(alpha=0.2, size=line_size) +

    # Plot measures of central tendency...
    geom_vline(aes(xintercept = mean(values), color='Media'), linetype = linetype, size=line_size) +
    geom_vline(aes(xintercept = mean(values, truncated_mean_value), color='Media Truncada'), linetype=linetype, size=line_size) +
    geom_vline(aes(xintercept = median(values), color='Mediana'), linetype=linetype, size=line_size) +
    geom_vline(aes(xintercept = max(values), color='Máximo'), linetype=linetype, size=line_size) +
    geom_vline(aes(xintercept = min(values), color='Mínimo'), linetype=linetype, size=line_size)

  for(var in mfv(df$values)) {
      p <- p + geom_vline(aes(xintercept = var, color='Moda'), linetype=linetype, size=line_size)
  }
  
  p <- p + scale_color_manual(
    name = "Medidas de tendencia central", 
    values = c(
      'Media'   = "black", 
      'Media Truncada' = 'wheat3',
      'Mediana' = 'red',
      'Máximo'  = 'darkolivegreen4',
      'Mínimo'  = 'darkgoldenrod1',
      'Moda' = 'blue'
    )
  )
  
  p <- p + labs(x=name, y = ylab, title = paste("Histograma", name, sep=" - "))

  p
}

plot_hist <- function(
  values, 
  name, 
  xlab, 
  ylab = "Frecuencia",
  lwd = 3,
  legend_ = TRUE,
  missing_value = 0,
  truncated_mean_value=0.05
) {
  # Fill missing values...
  values[is.na(values)] <- missing_value
  
  # Plot histogram..
  hist(
    values, 
    col  = "deepskyblue", 
    main = sprintf("Distribución - %s", name), 
    xlab = name,
    ylab = ylab,
    freq = FALSE
  )
  
  # Plot measures of central tendency...
  lines(density(values),                         col = "chocolate3",      lwd = lwd)
  abline(v = mean(values),                       col = "black",           lwd = lwd)
  abline(v = mean(values, truncated_mean_value), col = "wheat3",          lwd = lwd)
  abline(v = median(values),                     col = "red",             lwd = lwd)
  abline(v = mfv(values),                        col = "blue",            lwd = lwd)
  abline(v = max(values),                        col = "darkolivegreen4", lwd = lwd)
  abline(v = min(values),                        col = "darkgoldenrod1",  lwd = lwd)
  
  # Plot legend...
  if (isTRUE(legend_)) {
    legend(
      x = "topright",
      c(
        "Densidad", "Media","Media Truncada", 
        "Mediana", "Moda", "Máximo", "Mínimo"
      ),
      col = c(
        "chocolate3", "black", "wheat3", "red", 
        "blue", "darkolivegreen4", "darkgoldenrod1"
      ),
      lwd = c(lwd, lwd, lwd, lwd, lwd, lwd, lwd),
      cex = 1
    )
  }
}

box_plot <- function(data, horizontal = TRUE, xlab="", ylab="") {
  boxplot(
    data,
    xlab=xlab, 
    ylab=ylab,
    horizontal = horizontal,
    las=1,
    cex.lab=0.8, 
    cex.axis=0.6,
    pars=list(boxlwd = 2, boxwex=.8),
    col=colors()
  )
}

plot_heatmap <- function(data) {
  plot_ly(
    z = data, 
    y = colnames(data),
    x = rownames(data),
    colors = colorRamp(c("white", "red")),
    type = "heatmap"
  )
}

show_table <- function(table, page_size = 6, filter = 'top') {
  datatable(
    table, 
    rownames = FALSE, 
    filter=filter, 
    options = list(page_size = page_size, scrollX=T)
  )
}

A continuación se carga el dataset en una tabla para tener una primera vista de sus datos:

show_table(dataset[,2:8])

Veamos con que variables contamos y cual es su tipo. Para esto mostramos un resumen de la estructura del dataset:

str(dataset)

## 'data.frame':    984 obs. of  8 variables:
##  $ country_code  : chr  "AFG" "AFG" "AFG" "AFG" ...
##  $ country       : chr  "Afghanistan" "Afghanistan" "Afghanistan" "Afghanistan" ...
##  $ sub_nat_region: chr  "Badakhshan" "Badghis" "Baghlan" "Balkh" ...
##  $ world_region  : chr  "South Asia" "South Asia" "South Asia" "South Asia" ...
##  $ mpi_national  : num  0.295 0.295 0.295 0.295 0.295 0.295 0.295 0.295 0.295 0.295 ...
##  $ mpi_regional  : num  0.387 0.466 0.3 0.301 0.325 0.313 0.319 0.25 0.245 0.384 ...
##  $ hc_regional   : num  67.5 79.3 59.7 55.7 61 65.1 61.4 49.4 47.4 74.6 ...
##  $ iod_regional  : num  57.3 58.8 50.3 54.1 53.3 48.1 52 50.6 51.6 51.5 ...

Descripción de variables/columnas

• country_code: Código único de país.
• country: Nombre del país.
• sub_national_region: Regiones dentro de cada país.
• world_region: Regionales globales,
• mpi_national: Tasa de pobreza nacional.
• mpi_regional: Tasa de pobreza regional.
• hc_regional: Porcentaje de población pobre por región.
• iod_regional: Promedio de la distancia por debajo de la línea de pobreza por región.

Tipo de variables/columnas

• Categóricos o Nominales:
  • country_code
  • country
  • sub_national_region
  • world_region
• Escalas de razón o ratio:
  • mpi_national
  • mpi_regional
  • cp_regional
  • iod_regional

1.b. Represente gráficamente la cantidad de ciudades agrupados por región y pais.

cities.count.by.region <- dataset %>%
  group_by(world_region) %>%
  summarise(count = n_distinct(sub_nat_region)) %>%
  arrange(world_region)

cities.count.by.region.country <- dataset %>%
  group_by(world_region, country) %>%
  summarise(count = n_distinct(sub_nat_region), .groups = "keep") %>%
  arrange(world_region, country)

parent_regions  <- cities.count.by.region$world_region
regions         <- cities.count.by.region.country$world_region
countries       <- cities.count.by.region.country$country

cities.count         <- sum(cities.count.by.region.country$count)
country.cities.count <- cities.count.by.region.country$count
region.cities.count  <- cities.count.by.region$count

labels  <- c(c(parent_regions), c(countries))
parents <- c(rep(" ", length(parent_regions)), regions)
values  <- c(c(region.cities.count), c(country.cities.count))

# print(paste(length(labels), length(parents), length(values), sep=" - "))
plot_ly(labels = labels, parents = parents, values = values, type = 'sunburst', maxdepth=3)

2. Medidas de posición

Calcule las medidas de posición para los atributos numéricos y agrupe los cálculos de acuerdo a la Región.

2.a. Ordene los resultados del MPI resultante y concluya al respecto. Help(order).

 

Medida de posición para el indice de pobreza regional

show_table(
  dataset %>% 
    group_by(world_region) %>% 
    summarise(
      Media   = mean  (mpi_regional),
      Mediana = median(mpi_regional),
      Moda    = paste ("[", paste(mfv(mpi_regional), collapse = ', '), "]"),
      Min     = min   (mpi_regional),
      Max     = max   (mpi_regional) 
    ) %>%
    arrange(desc(Media))
)

Observaciones

• Utilizando la medida mas confiable de la tabla (Mediana), se evidencia que la región Africana es la más pobre, seguida de Asia (por una diferencia de 0.1) y por debajo las demás regiones, las cuales parecen estar en niveles mas cercanos de pobreza entre si.
• También podemos observar se manteiene la misma tendencia con el indicador Máximo. Es decir, que la sub región nacional con mas pobreza en el mundo esta en el continente Africano, seguido por el sudoeste de Asia.
• Los mínimos registran la misma tendencia que los máximos.
• Para alguna regiones se puede apreciar que hay varias modas. Esto indica que hay variaciones en las sub regiones con los mismos niveles de pobreza. Tenemos varios sub niveles de pobreza dentro del mismo país, a diferencia de otros países donde todos son pobres por igual (Nivel).

 

Medida de posición para el indice de pobreza nacional

show_table(
  dataset %>% 
    group_by(world_region) %>% 
    summarise(
      Media   = mean  (mpi_national),
      Mediana = median(mpi_national),
      Moda    = paste ("[", paste(mfv(mpi_national), collapse = ', '), "]"),
      Min     = min   (mpi_national),
      Max     = max   (mpi_national) 
    ) %>%
    arrange(desc(Media))
)

Observaciones

• Se registra un comportamiento muy similar al indicador anterior.
• En Asia las regiones menos pobres son mas pobre que las menos pobres en Africa. Es decir, la base de la pobreza es mucho mas baja que Africa a nivel país. Por supuesto que si comparamos con el mismo indice a nivel regional, vamos a encontrar que hay regiones en Africa que inicialmente son mas pobres, pero hay que destacar que a nivel regional, Asia inicia con una base mas pobre que Africa. Los demás países menos pobres inician con un mínimo muy bajo, lo cual es de esperar.

 

Medida de posición para el porcentaje de población pobre por región (Head count)

show_table(
  dataset %>% 
    group_by(world_region) %>% 
    summarise(
      Media   = mean  (hc_regional),
      Mediana = median(hc_regional),
      Moda    = paste ("[", paste(mfv(hc_regional), collapse = ', '), "]"),
      Min     = min   (hc_regional),
      Max     = max   (hc_regional) 
    ) %>%
    arrange(desc(Media))
)

Observaciones

• Esta variable se comporta muy similar a las dos anteriores. Las regiones mas pobres tiene un head count mas alto.

 

Medida de posición para el promedio de la distancia por debajo de la línea de pobreza por región

show_table(
  dataset %>% 
    group_by(world_region) %>% 
    summarise(
      Media   = mean  (iod_regional),
      Mediana = median(iod_regional),
      Moda    = paste ("[", paste(mfv(iod_regional), collapse = ', '), "]"),
      Min     = min   (iod_regional),
      Max     = max   (iod_regional) 
    ) %>%
    arrange(desc(Media))
)

Observaciones

• Sigo viendo la misma tendencia de crecimiento hacia la regiones mas pobres.

 

2.b. Grafique las variables y observe su comportamiento (graph: barplot, pie & hist).

p <- gplot_hist(dataset$mpi_national, name="Indice de pobreza nacional", binwidth=0.06)
ggplotly(p)

Observaciones

• Es una distribución bimodal.
• La distribución parece tener un sesgo positivo.
• Pareciera indicar que hay mas regiones menos pobre y un menor número de regiones que son mas pobres.

p <- gplot_hist(dataset$mpi_regional, name="Indice de pobreza regional", binwidth=0.06)
ggplotly(p)

Observaciones

• La distribución parece tener un sesgo positivo.
• Pareciera indicar que hay mas regiones menos pobre y un menor número de regiones que son mas pobres.

p <- gplot_hist(dataset$hc_regional, name="% regional de población pobre", binwidth=4)
ggplotly(p)

Observaciones

• La distribución parece tener un sesgo positivo.
• Según este gráfico, a medida que aumenta en nivel de pobreza va disminuyendo su población.

p <- gplot_hist(dataset$iod_regional, name="Intensidad de pobreza regional", binwidth=4)
ggplotly(p)

Observaciones

• Es una distribución bimodal.
• La distribución tiene muchos valores faltantes (entre 0 y 30).
• La distribución parece tener un sesgo positivo muy leve.
• Al ser una distribución parecida a la normal, entiendo que los niveles de personas por debajo de la linea de pobreza se mantienen en valores medios.

3. Medidas de dispersión. Calcular el desvío estándar, la varianza y el rango para cada una de las variables.

show_table(
  dataset %>% 
    group_by(world_region) %>% 
    summarise(
      Varianza=var(mpi_regional),
      Dispersion=sd(mpi_regional),
      Rango=max(mpi_regional) - min(mpi_regional)
    ) %>%
    arrange(desc(Varianza))
)

3.a. Realice diagramas de cajas y scatterplot’s. Documente las conclusiones.

 

Boxplot de variables numericas

p <- dataset.n %>% select(mpi_national, mpi_regional, hc_regional, iod_regional) %>%
  pivot_longer(., cols = c(mpi_national, mpi_regional, hc_regional, iod_regional), names_to = "Variables", values_to = "Frecuencia") %>%
  ggplot(aes(x = Variables, y = Frecuencia, fill = Variables)) +
  geom_boxplot()

ggplotly(p)

Boxplot para cada variable numerica segmentado por world region

ggplot(data = dataset.n, aes(x = mpi_national, y=world_region, fill=world_region)) +
  geom_boxplot(alpha=0.4) +
  ggtitle("MPI National por World Region")

ggplot(data = dataset.n, aes(x = mpi_regional, y=world_region, fill=world_region)) +
  geom_boxplot(alpha=0.4) +
  ggtitle("MPI Regional por World Region")

ggplot(data = dataset.n, aes(x = hc_regional, y=world_region, fill=world_region)) +
  geom_boxplot(alpha=0.4) +
  ggtitle("Head Count Regional por World Region")

ggplot(data = dataset.n, aes(x=iod_regional, y=world_region, fill=world_region)) +
  geom_boxplot(alpha=0.4) +
  ggtitle("IOD regional por World Region")

Diagramas de dispersion para todas las medidas numericas por world region

plot <- PairPlot(
  dataset.n[,4:8], 
  colnames(dataset.n)[5:8],
  " ",
  group_var = "world_region", 
  palette=NULL
) +
ggplot2::scale_color_manual(values=unique(as.factor(dataset.n$world_region)))

ggplotly(plot)

A continuación tomaremos las gráficas donde se encuentre una relación medianamente clara entre variables:

Relación entre hc_regional y mpi_regional

plot3 <- ggplot(dataset.n, aes(x = hc_regional, y = mpi_regional)) +
  geom_point(aes(shape = world_region, color = world_region)) +
  geom_smooth(method = 'loess', formula = 'y ~ x')
ggplotly(plot3)

Observaciones

• Al parecer hay una relación positiva, es decir, a medida que aumenta hc_regional también aumenta mpi_regional.
• Por otro lado, la curva tiene pocos valores dispersos, lo que indica que la relación entre variables es fuerte.

 

Relación entre iod_regional y mpi_regional

plot1 <- ggplot(dataset.n, aes(x = iod_regional, y = mpi_regional)) +
  geom_point(aes(shape = world_region, color = world_region)) +
  geom_smooth(method = 'loess', formula = 'y ~ x')
ggplotly(plot1)

Observaciones

• Es la relación positiva.
• La curva tiene algunos pocos valores dispersos, pero igualmente la relación parece ser fuerte.

 

Relación entre iod_regional y hc_regional

plot2 <- ggplot(dataset.n, aes(x = iod_regional, y = hc_regional)) +
  geom_point(aes(shape = world_region, color = world_region)) +
  geom_smooth(method = 'loess', formula = 'y ~ x')
ggplotly(plot2)

Observaciones

• Es la relación positiva.
• La curva tiene varios valores dispersos, pero igualmente conserva su forma.
• La relación parece ser medianamente fuerte.
• Al parecer, en valores de hc_regional mayores a 90 la relación se vuelve negativa.

3.b. ¿Qué variable es la que presenta mayor dispersión? Tenga en cuenta que cada variable puede estar expresada en diferentes unidades y magnitudes.

La relaciones mas débiles que tiene mayor dispersión son las siguientes:

• mpi_region vs. mpi_nacional
• hc_regional vs. mpi_nacional
• iod_regional vs. mpi_nacional

4. Medidas de asociación.

Calcular el coeficiente de correlación de todas las variables y explique el resultado. ¿Qué tipo de gráficos describen mejor esta relación entre las variables?

 

Matríz de correlación de pearson

cor_matrix = cor(dataset[5:8])
cor_matrix[upper.tri(cor_matrix)] <- NA
cor_matrix

##              mpi_national mpi_regional hc_regional iod_regional
## mpi_national    1.0000000           NA          NA           NA
## mpi_regional    0.8591325    1.0000000          NA           NA
## hc_regional     0.8555896    0.9839779   1.0000000           NA
## iod_regional    0.8052553    0.9347065   0.8946479            1

plot_heatmap(cor_matrix)

Nota

• Rojo alta correlación.
• Blanco correlación media.
• Verde baja correlación.

 

Observaciones

• Al parecer hay una alta correlación entre las variables mpi_regional y hc_regional. Ya que ambas variable se comportan de la la misma manera, se puede descartar a una de las dos. Particularmente me quedaría con mpi_region, ya que es una medida mas directa y representativa del dataset.
• Hay una correlación no tan fuerte pero marcada entre las variables mpi_regional y iod_regional. En este caso también descartaría iod_regional ya que se comporta prácticamente igual que la variable mpi_regional.

 

---

Realizado por Adrian Marino

adrianmarino@gmail.com