Análisis de conglomerados

Intentemos capturar ahora datos relacionados con Seguridad y Defensa para cada municipio en Colombia. Como vamos a obtener datos de endpoints de diferentes bases de datos, lo que haremos será usar un pequeño archivo de texto que contiene las URL de las que queremos extraer información. También crearemos una pequeña base de datos que nos sirva de referencia.

library(readr)
library(httr)
library(jsonlite)
library(plyr); library(dplyr)

rl <- readLines("endpoints.txt")
dbase <- gsub(".*Seguridad-y-Defensa/(.*)-2017/(.*)", "\\1", rl)
endpoint <- gsub(".*Seguridad-y-Defensa/(.*)-2017/(.*)", "\\2", rl)

endpoints <- data.frame("endpoint"=endpoint, "database"=dbase)
head(endpoints)

Ahora recogeremos los datos de cada uno de esos endpoints. Para ello usaremos un par de funciones que nos ayudarán en el proceso.

creds <- readLines("./credentials-datos-abiertos.txt")
baseurl <- "https://www.datos.gov.co/resource/%s.json"

pasar_pagina <- function(i, block=1000) {
    return(list("$limit"=block,
                "$offset"=block*(i - 1))) 
}

collect_data <- function(url, maxpag=10) {
    status <- TRUE
    i <- 1
    res <- list()
    while (status) {
        print(sprintf("Recogiendo datos de pagina %s", i))
        r <- GET(url,
                 query=pasar_pagina(i),
                 add_headers("X-App-Token"=creds))
        status <- status_code(r) == 200
        if (status) {
            print("...Guardando datos")
            if (length(content(r)) > 0) {
                res[[i]] <- fromJSON(content(r, "text"))
            } else {
                print("I'm stopping now!")
                status <- FALSE
            }
        }
        if (i >= maxpag | !status) { ## Para no eternizarnos
            status <- FALSE
        }
        ## Sys.sleep(2)
        i <- i + 1
    }
    return(do.call(rbind.fill, res))
}

Lo que haremos será llamar cada uno de esos endpoints, recoger hasta 5 páginas de información y agrupar la información por departamento y municipio en diferentes elementos de una lista:

out <- list()
for (i in 1:nrow(endpoints)) {
    URL <- sprintf(baseurl, endpoints$endpoint[i])
    out[[i]] <- as.data.frame(collect_data(URL, 5))
    out[[i]] <- out[[i]] %>%
        group_by(departamento, municipio) %>%
        summarize(n=n())
    out[[i]]$var <- endpoints$database[i]
}
str(out)

Ahora juntamos todas las bases en un único data.frame

out <- do.call(rbind, out)

y lo reorganizamos de tal manera que cada variable sea una nueva columna

dt <- reshape2::dcast(out, departamento ~ var, value.var="n", fun.aggregate=sum)
departamento <- dt$departamento
dt$departamento <- NULL

Lo que hemos conseguido después de estos pasos es una base datos en los que cada fila es una fila y cada columna representa un tipo de incidente.

Como hemos visto, los algoritmos de conglomerado, en esencia, lo que hacen es estructurar una matriz de distancias entre observaciones, poniendo juntos casos que están cerca y separando casos que están lejos. Para ver la intuición detrás de este proceso, podemos empezar por calcular las distancias entre tres observaciones de la base de datos.

dist(dt[1:3, ])

Como vemos el par 1/3 está “cerca”, pero el par 1/2 y el par 2/3 está “lejos”. Esto parece consistente con los datos en bruto. Si calculamos la distancia para cada una de las variables vemos que efectivamente hay poca discrepancia entre 1 y 3 en casi todas las dimensiones:

dt[1, ] - dt[3, ]

Pero que sin embargo, los otros dos pares muestran diferencias mucho más elevadas.

dt[1, ] - dt[2, ]
dt[2, ] - dt[3, ]

Una observación que no debería pasar desapercibida es que algunos de estos delitos son mucho más frecuentes que otros. Por ejemplo, hurto de motocicletas es mucho más frecuente que delitos de terrorismo. Por eso, la distancia entre dos casos sigue implícitamente lo que ocurra en esa variable. Sin embargo, eso no tiene en cuenta que una diferencia de 5 delitos de terrorismo es mucho más significativa que una diferencia en 5 hurtos de motocicleta, simplemente porque un delito de terrorismo ocurre con muy escasa frecuencia.

Parece más razonable estandarizar los datos de tal manera que todas las variables estén en la misma escala. R provee la función scale, pero podemos hacerlo manualmente para que sea más claro:

sdt <- apply(dt, 2, function(x) (x-mean(x))/sd(x))
head(sdt)

Ahora podemos recalcular las mismas distancias que antes

dist(sdt[1:3, ])

Vemos que la misma jerarquía se mantiene pero ahora las distancias se acortan mucho.

Análisis de conglomerados

Métodos aglomerativos y divisivos

Ahora podemos pasar al análisis de nuestros datos

library(cluster)

Empezaremos por construir nuestros aglomerados usando un método que empieza cada observación como unidad independiente y va agrupando casos en función de su similitud:

agmodel <- agnes(sdt)
agmodel

Vemos la ordenación de los objetos y un resume de la altura, pero quizás lo más práctico sea mirar el dendrograma resultante:

plot(agmodel, ask=FALSE, which.plot=2)

Vemso que dos grupos aparecen con mucha claridad. Por una parte, los departamentos 2, 30 y 14, que, sin ser del todo parecidos entre ellos, son muy diferentes a los demás. Por otra, todas las demás observaciones.

departamento[c(2, 30, 14)]

Podemos comprobar manualmente la lógica de esta diferencia. He añadido dos observacioens más para hacer más clara la diferencia:

dist(sdt[c(2, 30, 14, 1, 3), ])

Estamos por tanto ante tres departamentos con relativamente pocos delitos.

El mismo proceso lo podemos hacer a la inversa, empezando con un único grupo y dividiendolo hasta llegar a observaciones individuales

dimodel <- diana(sdt)
dimodel
plot(dimodel, ask=FALSE, 3)

Como vemos no hace ninguna diferencia práctica, como suele ser el caso.

Método de k-medias

En el modelo anterior, hemos obtenido el número de grupos inspeccionando el dendrograma después de calcular distancias entre cada observación. Otra alternativa es fijar de antemano el número de conglomerados y decidir cómo se dividirían las observaciones en ese caso.

La función kmeans nso permite pasar una matriz de datos y un número de grupos (o un vector de centros) y calcula la división óptima de los casos asumiendo ese número como correcto.

Asumamos que existan 4 grupos en la población.

kmeans(sdt, 4)

La función nos ofrece, en primer lugar, el centro de cada uno de los grupos (que se leería como el vector fila). También nos dice a qué grupo pertenece cada observación y datos sobre la suma de cuadrados explicada por esta asignación.

¿Cómo decidir si nuestro supuesto es correcto? El mejor modo es hacer diferentes análisis bajo diferentes supuestos y comprobar los resultados. Por ejemplo, podemos clasificar los datos en una secuencia de 1 a 10 grupos:

K <- 10
m <- vector("list", K)
for (k in 1:K) {
    m[[k]] <- kmeans(sdt, centers=k)
}

Ahora podemos ver cómo evoluciona la descomposición de la suma de cuadrados:

bw <- ldply(m, function(x) data.frame("b"=x$betweenss, "w"=x$tot.withinss))
plot(1:K, bw$b, "l", col="blue")
lines(1:K, bw$w, col="red")

Vemos que a partir de dos grupos, el incremento marginal en la coherencia de los grupos y la distancia entre los mismos es mucho menor. Es decir, a partir de 2 grupos, no merece la pena seguir dividiendo.

Hay otros criterios para tomar esta decisión.

library(fpc)
model <- kmeansruns(sdt, krange=1:10, criterion="ch")
model$crit
model$bestk

o

model <- kmeansruns(sdt, krange=1:10, criterion="asw")
model$crit
model$bestk

Lo que vemos es que da igual qué criterio escojamos, dos grupos son suficientes para representar la estructura de los datos. Además, si usamos este modelo, con sus centros, para predecir qué observación cae en cada grupo vemos que

m <- kmeans(dt, 2)
which(m$cluster == 2)

que son las mismas observaciones que vimos antes.

saveRDS(dt, "dta/datos-criminalidad.RDS")