Usando clustering para identificar cursos no Prouni

O que é clustering?

Clustering é uma classe de algoritimos não-supervisionados para classificação de observações. Existem vários tipos, cores e tamanho de técnicas de clustering, mas essa bonita variedade vai ficar para outro dia porque o foco de hoje é a abordagem de distância centrada.

“Agrupamento de observações”

A visualização é razoavelmente clara, clusters são literalmente agrupamentos. Com base em alguns critérios dependentes do algoritimo a ser utilizado, você classifica uma observação em um ou outro agrupamento (exceto nos modelos fuzzy, mas isso fica para outro dia).

Clustering k-means como um problema de otimização

Um problema de otimização irrestrita tem, a grosso modo, duas features. A função objetivo a ser maximizada ou minimizada e o instrumento com o qual atingir tal objetivo. Aqueles familiarizados com o canônico método de estimação por Mínimos Quadrados Ordinários vão reconhecer alguma semelheança.

K-means, ao invés de minimizar quadrado dos resíduos, minimiza a soma do quadrado da distância dentro do cluster (WCSS, em inglês). Nossos instrumentos são \(k\), o número de agrupamentos e \(S_i\), os conjuntos que dão qual elemento está em qual agrupamento. Podem parecer instrumentos redundantes à primeira vista. Pense que para um mesmo número de agrupamentos, é possível ter combinações de conjuntos com WCSSs diferentes.

Algumas definições antes. \(k\) é o número de clusters, \(S_i\) é conjunto de elementos do i-ésimo cluster, \(\mu_i\) é a média do i-ésimo cluster.

\[ \text{arg min}_S \sum_{i=1}^k \sum_{x_j \in S_i} || x_j - \mu_i ||^2 \]

O leitor atento percebeu que \(k\) não aparece aqui como um instrumento do problema, mas sim como um parâmetro dado. Bem, aí está uma das peculiaridades de k-means, nós escolhemos o k. É uma tarefa que tem um pouco de ciência e muita arte, vou me aprofundar um pouco nela mais à frente.

Os dados

A amostra que temos é do ProUni de 2017 e conta com algo em torno de 32 mil observações. Já tive o trabalho de limpar a base para vossa apreciação e vou deixa-la disponível aqui e o código que contém tudo aqui. Vamos primeiro explorar nossa amostra com a ajuda do ggplot2.

final %>%
  ggplot(aes(x = mensalidade, y = nota, 
             color = medicina, show.legend = FALSE)) +
  geom_point()+
  xlab("Mensalidade do curso no ProUni")+
  ylab("Nota de Corte do curso no ProUni")+
  labs(col="Medicina")

final %>%
ggplot(aes(x = mensalidade)) + 
  xlim(0,2500) +
  geom_histogram(aes(y=..density..), binwidth = 50, fill = "#3385ff") +
  xlab("Mensalidade do curso no ProUni") + 
  ylab("") +
  scale_y_continuous(labels = percent) +
  geom_vline(aes(xintercept = mean(mensalidade, na.rm=T)),   
             color="black", linetype="dashed", size=1)

final %>%
  ggplot(aes(x = mensalidade, fill = medicina)) + 
  xlab("Mensalidade do curso no ProUni") + 
  ylab("") +
  geom_histogram(aes(y=..density..), binwidth = 300) +
  scale_y_continuous(labels = percent) 

Aqui observamos três coisas muito interessantes. A primeira é que notas de corte seguem muito bem uma distribuição normal exceto pela regra que impõe nota de corte mínima de 450 no ProUni e Sisu. É o tipo de coisa em que seria legal aplicar um Teste de Densidade de McCrary (2006). Depois, que Medicina tem um padrão de distribuição de notas bem diferente do resto.

final %>%
  ggplot(aes(x = nota)) +
  geom_histogram(aes(y = ..density..), binwidth = 10, fill = "#3385ff") +
  scale_y_continuous(labels = percent) + 
  xlab("Nota de corte do curso no ProUni")

final %>%
  ggplot(aes(x = nota, fill = medicina)) +
  geom_histogram(aes(y = ..density..), binwidth = 10) +
  scale_y_continuous(labels = percent) + 
  xlab("Nota de corte do curso no ProUni")

Como escolher \(k\)?

Essa é a pergunta de um milhão de dólares, honestamente. Eu encontrei dois principais métodos, um é computacionalmente exigente e preciso, o outro é computacionalmente simples e depende mais de interpretação.

A primeira e mais complicada é a Gap Statistic (Tbishirani, Walther e Hastie, 2001). O método envolve algumas computações com bootstrap, então exige uma máquina preparada. Só consegui rodar usando um servidor, então evite esse método se tiver um computador normal (ou até mesmo um pessoal de alta qualidade). Em qualquer caso, a implementação desse método é a função cluster::clusGap.

O segundo método é o do “Cotovelo”. Não é muito sofisticado, mas é potente. Plotamos o WCSS como uma função de \(k\) e procuramos por uma inflexão na curva. Onde ela tiver um “cotovelo”, é provavelmente o \(k\) mais adequado. A função a seguir implementa o gráfico:

wssplot <- function(data, nc=15, seed=1234){
  wss <- (nrow(data)-1)*sum(apply(data,2,var))
  for (i in 2:nc){
    set.seed(seed)
    wss[i] <- sum(kmeans(data, centers=i)$withinss)}
  
  plot(1:nc, wss, type="b", xlab="Number of Clusters",
       ylab="Within groups sum of squares")}

final$completo = complete.cases(final)

final.cotovelo = data.frame(nota = final$nota, 
                            mensalidade = final$mensalidade,
                            completo = final$completo)

final.cotovelo = final.cotovelo[final$completo == TRUE,]
final.cotovelo$completo = NULL

wssplot(final.cotovelo, 
          nc = 6) 

Agora que escolhemos o número 3, podemos finalmente ver se o modelo classifica bem cursos de medicina.

Os finalmentes, rodando o modelo e resultados

Tendo 3 como o número mágico, podemos finalmente rodar o modelo. kmeans é um comando nativo do R. Ele retorna um objeto de classe kmeans e os agrupamentos em específico estão no vetor cluster. É interessante visualizar o resultado através do primeiro gráfico - agora com color coding baseado nos clusters, mas uma visualização específica pode ser melhor - afinal, é interessante ver os agrupamentos. Para isso, vamos usar cluster::clusplot:

analise_kmeans <- kmeans(final.cotovelo, 
                          centers = 3)

##### Agora visualize os resultados

clusplot(final, analise_kmeans$cluster,
                        main='Procurando por 3 agrupamentos no ProUni',
                            color = TRUE,
                              shade = TRUE,
                                lines = 0)

final %>%
  ggplot(aes(x = mensalidade, y = nota,
             colour = factor(analise_kmeans$cluster), 
             show.legend = FALSE)) +
  geom_point() +
  xlab("Mensalidade do curso no ProUni") +
  ylab("Nota de Corte do curso no ProUni") +
  labs(col = "Agrupamento")

Se o leitor fizer o exercício de replicar esse post, vai poder ver que o algoritimo identificou todos os 122 cursos de medicina da amostra e inseriu por engano no mesmo cluster 44 cursos que não são de medicina.

A performance é aceitável? Deixamos a resposta para outro post.