Análisis de temas: Estimación

Ahora que tenemos nuestros datos en una matriz de términos, podemos empezar a hacer nuestros análisis.

library(topicmodels)
library(tm)
dtm <- readRDS("./dta/dtm.RDS")

El paquete topicmodels ofrece una función LDA para la estimación de Latent Dirichlet Allocation models. Esta función toma tres argumentos principales: 1. una matriz de términos como la que acabamos de construir, 2. el número de temas que creemos que existen en nuestra base de datos, y 3. el método de estimación. Con ello, nuestro modelo tendría el siguiente aspecto.

lda_model <- LDA(dtm,
              5,
              method="Gibbs")
lda_model

La estimación mediante Gibbs sampler es una forma común de estimar este tipo de modelos, aunque los fundamentos de estos métodos de simulación nos llevaría por derroteros muy alejados de los contenidos del curso. Nos conformaremos por ahora con fijar valores que son razonables pero que seguramente requerirán ajuste para situaciones reales. Lo importante ahora es entender intuitivamente qué hace cada uno de los parámetros y cómo manipularlos si nuestra simulación no funciona correctamente.

Una versión más completa que incorpora detalles de cómo queremos que se ejecute la simulación es la siguiente:

k <- 5
burnin <- 4000
iter <- 2000
thin <- 500
seed <- 2003
nstart <- 1
best <- TRUE
ctrl <- list(nstart=nstart,
            seed=seed,
            best=best,
            burnin=burnin,
            iter=iter,
            thin=thin)

lda_model <- LDA(dtm,
                k,
                method="Gibbs",
                control=ctrl)

Aquí estamos fijando el número de iteraciones que queremos descartar de las estimación final porque forman parte de la transición a la distribución estacionaria (burnin), el número de iteraciones que queremos escoger después de que el modelo haya terminado el burn-in period (iter) y el offset entre iteraciones (thin). El parámetro seed es la semilla del generador de números pseudo-aleatorios.

Ahora que hemos ajustado el modelo, podemos ver la distribución de términos en cada tema

terms(lda_model, 5)

O podemos ver también qué tema se asocia con cada uno de los documentos, si asumimos que el tema de mayor probabilidad es al mismo tiempo el tema común de un documento.

head(topics(lda_model))

Una visión más completa la podemos ver en la matriz gamma que muestra la probabilidad de cada tema en cada documento.

as.data.frame(lda_model@gamma)

Por supuesto, el problema de la estrategia que hemos utilizado es que hemos fijado un número de temas que puede no ser el correcto ¿Cómo escoger cuál es el número de temas si no tenemos información ex ante? La mejor estrategia, al igual que vimos en el caso de los modelos de agrupamiento, es estimar el modelo para un gran número de temas y comprobar qué tan bien cada supuesto ajusta los datos.

El límite en este caso es que cada uno de los modelos consume mucho tiempo. Para esquivar ese problema, lo que haremos será estimar en paralelo un modelo con un número diferente de temas. Dependiendo del número de núcleos que tenga nuestra máquina, esto puede rebajar muy considerablemente el tiempo que tenemos que esperar.

Empezaremos por cargar en nuestra sesión los paquetes de R para realizar computaciones en paralelo.

library(parallel)
library(doMC); library(doParallel)
cl <- makeCluster(detectCores() - 1)
registerDoParallel(cl)

En este caso, hemos definido un cluster de tantos núcleos como tenga la máquina en la que estamos trabajando, reservando uno para el sistema operativo.

Ahora podemos estimar nuestro modelo para una serie de modelos que asumen un número de temas entre 2 y 29.

mlist <- foreach(i=seq(2, 29, 3),
                 .packages="topicmodels",
                 .export=c("ctrl", "dtm")) %dopar% {
    out <- LDA(dtm,
               i,
               method="Gibbs",
               control=ctrl)
    return(out)
}

Aquí estamos usando la función foreach que es similar a for pero intenta enviar cada iteración a un proceso diferente. El número de temas es el índice de la iteración y solo tenemos que tener en cuenta que tenemos que exportar explícitamente los objetos que el cuerpo de la iteración necesita para poder trabajar. Pensad en ello como lo que una nueva sesión de R necesitaría para poder hacer el análisis que hemos solicitado.

Lo que tenemos es una lista mlist que contiene, en cada elemento, el resultado de estimar un modelo LDA con un número diferente de temas ¿Cómo podeoms escoger entre estas alternativas? Las formas más comunes son fijándose en la verosimilitud de cada modelo o en su perplejidad, que es una medida de la capacidad del modelo para capturar la estructura de nuevos datos. Idealmente, haríamos esta estimación usando validación cruzada, que veremos en la segunda mitad del curso, pero por ahora nos limitaremos a ver cómo recuperar estas dos cantidades.

mlogLik <- as.data.frame(as.matrix(lapply(mlist, logLik)))
## mperp <- as.data.frame(as.matrix(lapply(mlist, perplexity)))

Y podemos crear gráficos para ver con más claridad qué tan bien funciona cada modelo:

plot(seq(2, 29, 3),
     unlist(mlogLik),
     xlab="Número de temas",
     ylab="Log-Verosimilitud")