Classes Sociais, Renda e Clusters

Breve introdução

Se tem um negócio mal explicado é a questão da divisão em classes sociais no Brasil. A impressão que tenho é que cada um usa a metodologia que mais lhe interessa para explicar aquilo que quer.

O IBGE

O IBGE tem sua metodologia que divide as classes baseados no número de salários mínimos de um domicílio. Ok, bacana, se a renda é maior do que ‘X’ a família é classe C, se é maior do que ‘Y’ classe B e por aí vai. Até que é legal, mas utilizar a renda de uma família e não a per capita nos deixa com o problema Homer Simpson:

Problema Homer

ABEP

Esse aqui é bacana!!! Sabemos que a pesquisa do IBGE tem vários problemas, entre eles a imprecisão das suas medidas, então alguém teve a ideia de contar os bens de uma pessoa para definir sua classe social, a ideia até que é interessante, mas ao tentar resolver o problema criamos tantos outros que vai ser difícil descrever em um post só…

A tabela é mais ou menos assim:

ABEP

Terceirizando o problema

É óbvio que não temos uma solução para isso, então vamos tentar encarar esse problema de outra forma.

Imagine um alienígena que veio para nosso país, começou a observar a população brasileira e tomar notas para depois responder a seguinte pergunta: ‘Como subdividir a população brasileira em classes sociais’?

Em ‘machine learning’ de modo bem grosseiro podemos dividir os algoritmos em dois:

Supervisionados e não supervisionados.

No caso da tarefa que passamos para nosso colega de outro mundo, estamos falando de um problema não supervisionado, pois ninguém definiu para o alien as classes e como subdividi-las, nem o que é ou não mais importante nesse momento, simplesmente demos para o coitado o fardo de realizar tal classificação.

O modelo

Em poucas palavras, vamos tentar separar a população brasileira em grupos baseados simplesmente na distância entre essas pessoas. Então, imagine uma linha reta:

X——————————————–X

Agora, imagine que temos algumas pessoas distribuídas ao longo da mesma:

X—OOO——OOOO——OOO——X

Só de olhar para essa linha nós institivamente conseguimos separar em 3 grupos de pessoas. Agora imagine que ao invés de uma linha medida em metros ou quilômetros, estamos falando da Renda das famílias. Para famílias que ganham X, Y ou Z, a separação em grupos seria praticamente a mesma.

A diferença com o nosso modelo é que ao invés de calcularmos a distância em um espaço unidimensional (uma linha) ou bidimensional, calculamos a distância em um espaço n-dimensional.

Os dados

Para simular o problema do alienígena, vamos utilizar a PNAD de 2013, porém o triste é que com ela não é possível reproduzir o método da ABEP, já que não temos o número de vários itens, só uma variável binária dizendo se a família o possui ou não, então vamos jogar a PNAD praticamente inteira no modelo e criar nossa própria classificação híbrida, parte ABEP e parte IBGE.

Ao invés de contar pontos de item por item, criamos uma variável chamada BENS que soma 18 questões diferentes da PNAD para dar o número de bens das famílias, como TV, DVD, Fogão, Celular, 3/4G e etc. (A ABEP atribui pesos para cada item, nós não o faremos por uma questão de simplificação, mas um refinamento deste modelo poderá ter pesos diferentes.)

Além desta, utilizamos mais 8 variáveis da PNAD domicílios, e 03 variáveis construídas da PNAD pessoa que são ‘MAIS_INSTR’ que é igual ao número de anos de estudo da pessoa mais instruída do domicílio, número de filhos e número de empregados domésticos.

##  [1] "Total_de_moradores" "Comodos"            "Banheiros"         
##  [4] "n_TV_telafina"      "n_TV_Tubo"          "Automovel"         
##  [7] "Renda"              "N_FILHO"            "N_DOMESTICO"       
## [10] "INSTR"              "BENS"

Resultados – Iteração #1

Os resultados ficaram assim (PART é a participação da CLASSE no total):

## Source: local data frame [6 x 3]
## 
##   CLASSE PART      RENDA
## 1      A  0.2 84225.6966
## 2     B1  0.9 23496.5027
## 3     B2  3.0 12404.7848
## 4     C1  8.5  6536.8275
## 5     C2 21.2  3394.0439
## 6    D-E 66.1   899.4745

plot of chunk unnamed-chunk-4

Ficou bem diferente do que estamos acostumados a ver, acho que nosso amigo alienígena errou feio!

Por outro lado, uma informação bem interessante que é possível derivar deste gráfico é que o brasil é um país extremamente desigual e a renda média da classe ‘A1’ acabou ficando em torno de $84k por mês!

We are the 99% – Iteração #2

A desigualdade no Brasil é tão grande, mas tão grande que enviesou completamente o modelo, porém até o documento da ABEP alerta sobre essa desigualdade indicando que para casos onde a renda é superior a R$30.000 o modelo de classes não ajuda muito.

Para resolver isso, vamos adotar o parâmetro ‘occupy Wall Street’ no nosso modelo e simplesmente ignorar o 1% mais rico da população, afinal, eles já venceram a corrida dos ratos.

Essas famílias são automaticamente classe ‘A’, da mesma forma que o pessoal de renda zero é classe ‘E’, porém como a PNAD contém várias omissões e erros, vamos ignorar o pessoal de renda Zero também.

## Source: local data frame [6 x 3]
## 
##   CLASSE  PART     RENDA
## 1      A  4.14 12657.533
## 2     B1  5.09  6973.986
## 3     B2  8.94  4720.220
## 4     C1 12.61  3274.278
## 5     C2 19.16  2245.416
## 6    D-E 50.06  1060.438

plot of chunk unnamed-chunk-6

Melhorou bastante! Agora está até lembrando as classificações mais tradicionais, apesar que ainda nosso modelo coloca mais pessoas nas classes D e E do que o que estamos acostumados, isso pode ser um problema da amostra ou político, afinal é bem conveniente classificar a maioria da população como vivendo no ‘meio do bolo’ para dar uma perpétua sensação de mobilidade social, porém os dados pintam um quadro diferente.

Ficando mais técnico – Controle da variância

Sabemos que as variáveis utilizadas no modelo podem ter impactos diferentes na classe social, por isso na hora de montar o modelo de clusterização (k-means) decidimos deixar as variáveis no estado mais puro o possível (direto da base), porém algumas variáveis podem roubar o show em um modelo de clusterização.

Vamos pegar um exemplo de livro texto.

plot of chunk unnamed-chunk-7

Lindo, certo? Em um espaço bidimensional, a relação entre as variáveis X e Y praticamente implora por uma clusterização, note que o modelo até comete alguns erros classificando os ‘O’s e os ‘+’s mas ele acabou fazendo um ótimo trabalho na clusterização no final das contas.

Agora vamos olhar outro exemplo

plot of chunk unnamed-chunk-8

Ao olhar para este segundo exemplo conseguimos até separar em dois clusters, mas essa separação visual é feita praticamente só com a informação do eixo ‘Y’, isso faz com que Y enviese o resultado por uma condição matemática e não real.

Apesar de controlarmos por várias variáveis no nosso modelo, isso pode acontecer, por exemplo com o número de empregados domésticos que vai de 0 a 2 no máximo, versus a renda que vai de 0 a mais de 150.000.

Pré-processamento

Para tentar tratar esse problema vamos utilizar técnicas básicas de pré-processamento centralizando e colocando as variáveis na mesma escala, dessa forma todas as variáveis do nosso modelo ficarão com média ‘0’ e variância ‘1’.

##                   VAR MEDIA VARIANCIA
## 1  Total_de_moradores     0         1
## 2             Comodos     0         1
## 3           Banheiros     0         1
## 4       n_TV_telafina     0         1
## 5           n_TV_Tubo     0         1
## 6           Automovel     0         1
## 7               Renda     0         1
## 8             N_FILHO     0         1
## 9         N_DOMESTICO     0         1
## 10              INSTR     0         1
## 11               BENS     0         1

Variáveis ajustadas – Iteração #3

## Source: local data frame [6 x 3]
## 
##   CLASSE  PART    RENDA
## 1      A  6.71 7656.665
## 2     B1 10.29 3998.186
## 3     B2 18.88 3000.812
## 4     C1 12.63 1976.285
## 5     C2 19.08 1849.001
## 6    D-E 32.40 1483.811

plot of chunk unnamed-chunk-11

Esse novo modelo ficou com uma distribuição melhor ainda que o anterior, mas o gráfico de Renda x Classe acabou ficando com vários outliers, indicando que agora com as variáveis ajustadas a renda não tem uma importância tão grande quanto tinha antes. É bem possível que a influência das outras variáveis na classe social seja maior. Ótimo, era mais ou menos essa a intenção quando misturamos as metodologias do IBGE e ABEP.

Análise dos outliers

Vamos olhar para a renda para tentar entender o comportamento desses outliers.

plot of chunk unnamed-chunk-12

Até que no final das contas eles estão bem-comportados.

Próximos passos

Até que para uma análise de fim de semana o modelo ficou bem interessante. Alguns caminhos que podem ser seguidos daqui para frente é tentar analisar com esse modelo a evolução do conceito de classe social ao longo dos anos e tentar derivar algum conceito de mobilidade social.

Um exemplo disso é que de acordo com a metodologia da ABEP, se uma família comprar um carro ela marca mais pontos e pode ‘subir de classe’, porém se todo mundo fizer isso as famílias continuam semelhantes entre elas, o que não justifica uma mudança de classe social.

Uma forma interessante de analisar a ‘Tragédia dos comuns’…

Advertisements

A PNAD, os DVDS e as Florestas (Aleatórias)

Motivação

Em jun/15, em uma conversa de ‘boteco’ alguém sugeriu que eu comprasse um DVD da Galinha pintadinha para minha filha recém-nascida e a primeira coisa que pensei era: – ‘Ótimo, aí eu penduro o DVD no móbile dela’, pois não tenho onde reproduzir esse troço…

Isso me lembrou a época de ouro das locadoras de vídeo, onde era praticamente um passeio de família ir a locadora pegar um ou dois filmes para assistir no fim de semana e gigantes como a Blockbuster quando chegaram ao Brasil mudaram a cara deste mercado para sempre, da locadora de porta de garagem para algo grandioso, cheio de opções e principalmente sem filas de espera…

Porém no final de 2013 foi oficialmente declarada a morte da Blockbuster link, uma morte mais do que esperada, mas ainda assim o fim de uma era…

…Com isso, vamos fazer uma viagem no tempo, pegar os dados da PNAD 2013 e tentar responder a seguinte pergunta:

“Com tantas maneiras de assistir filmes, e as ruínas das locadoras espalhadas pelo país, o que leva uma família a ter um aparelho de DVD (em 2013)?”

Vamos tentar montar um modelo para explicar isso…

A PNAD

Para quem gosta de trabalhar com dados, o IBGE é um prato cheio. A PNAD é uma pesquisa divertidíssima onde temos uma bela amostra de dados (longe de ser perfeita, mas boa para algumas análises) e perguntas o suficiente para a realização das mais variadas análises.

Estávamos passando pelos dicionários da pesquisa de 2013 recentemente e encontramos a seguinte questão:

  • ‘Tem aparelho de DVD’?

Os dados

Não vamos nos preocupar em explicar como extraímos os dados da PNAD, as instruções estão no site do IBGE neste link e como isso por si só já dá um belo trabalho, talvez um dia no futuro escreveremos um post sobre como fazer isso.

A última PNAD anual foi realizada em 2013, então é essa pesquisa que vamos utilizar. Atualmente existe a PNAD contínua, mas para o efeito que queremos capturar nesse exercício a de 2013 está mais do que ótima.

Factor engineering

Essa parte é relativamente trabalhosa, pois temos que:

  1. Transformar variáveis numéricas em fatores. Olhar para o dicionário da PNAD é um exercício de paciência. Por exemplo:

‘Tem televisão em cores?’

2 = SIM
4 = NÃO

‘Tem aparelho de DVD?’

1 = SIM
3 = NÃO

PORQUEEEEEEEEEEEE?????

É muito difícil normalizar isso? Que tal 1=sim e 0=não? Ou, para facilitar, 2 = sim, 4 = não….sempre!

  1. Criar algumas variáveis que vão ajudar na nossa modelagem, como por exemplo ‘Acesso_INTERNET’ que combina as questões de banda larga e discada em uma variável
  2. Vamos usar os dados do questionário dos domicílios, mas no dicionário das pessoas vamos criar variáveis como EH_FILHO que é igual ao número de filhos de uma família, EH_FILHO_6 que é igual a quantidade de filhos menores do que 6 anos e FAIXA_IDADE_FILHO que é a média de idade dos filhos de um domicílio

Os datasets

A primeira coisa a fazer é criar os datasets de treino e teste, só que temos um problema, considerando que a amostra da PNAD domicílios é de mais de 148.000 observações onde 116.000 são válidas e cada um tem seu peso, na teoria teríamos que repetir cada domicílio um número de vezes igual ao seu peso, ou seja, se o peso médio é 560, temos mais de 64.000.000 domicílios.

O dataset do nosso modelo vai consistir em 100.000 observações selecionadas aleatoriamente do universo de 64.000.000

DICA DO DIA: Nunca tente criar um vetor de 64.000.000 de observações no R e selecionar 100.000 linhas aleatoriamente!!!

Depois de travar o PC mudamos de estratégia…

Como?

Vamos assumir um universo de 3 domicílios com pesos 1.000.000, 50.000.000 e 13.000.000, criamos um vetor de pesos acumulados, 1.000.000, 51.000.000 e 64.000.000, fazemos 100.000 sorteios aleatórios de 1 a 64.000.000 e simplesmente contamos o número de itens no vetor de pesos acumulados menores do que o sorteio mais 1, por exemplo:

Sorteio: 10; número < 0; posição = 1
Sorteio: 1.100.000; numero < 1; posição = 2
Sorteio: 64.000.000; numero < 2; posição = 3

Isso faz com que a probabilidade se selecionarmos um domicílio seja igual seu peso no total.

Treino e teste

Agora que temos as 100.000 observações aleatórias, vamos dividir em dois datasets, um para treinar o nosso modelo e outro para testar seus resultados, faremos isso com a razão de 70% treino e 30% teste

 

O modelo

Para o modelo temos um problema de classificação, ou seja, para um certo grupo de inputs queremos uma resposta 0 ou 1 para classificar um domicílio como ‘Tem DVD’ ou ‘Não tem DVD’, esses inputs são algumas variáveis socioeconômicas dos domicílios e dos seus habitantes.

library(randomForest)
modelFit <- randomForest(Tem_DVD ~.,
                         data=train,
                         do.trace=T, 
                         ntree=200,
                         importance=TRUE,
                         mtry=10)
plot(modelFit,main='randomForest error rate')

plot of chunk unnamed-chunk-4

O pessoal costuma ir para o Overkill rodando milhares de árvores, mas a verdade é que não precisamos de tantas para estabilizar o erro, além disso demora uma eternidade! É possível ver pelo gráfico acima que depois de 100 árvores o erro já está mais do que controlado e o ganho é marginal.

Mostrando o pau

Agora que matamos a cobra, vamos aplicar nosso modelo no dataset de testes e ver se o mesmo ficou aceitável ou se somente capturou o ruído do dataset de treino.

pred <- predict(modelFit,newdata=test)
CM <- confusionMatrix(test$Tem_DVD, pred)
CM
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction     0     1
##          0  5573  2835
##          1  1578 20013
##                                           
##                Accuracy : 0.8529          
##                  95% CI : (0.8488, 0.8569)
##     No Information Rate : 0.7616          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.6179          
##  Mcnemar's Test P-Value : < 2.2e-16       
##                                           
##             Sensitivity : 0.7793          
##             Specificity : 0.8759          
##          Pos Pred Value : 0.6628          
##          Neg Pred Value : 0.9269          
##              Prevalence : 0.2384          
##          Detection Rate : 0.1858          
##    Detection Prevalence : 0.2803          
##       Balanced Accuracy : 0.8276          
##                                           
##        'Positive' Class : 0               
## 

TMI [Too much information] – Mais para frente faremos um post explicando detalhes do que é cada um desses resultados, por enquanto vamos focar em duas coisas:

Accuracy : 0.8529

O modelo tem uma precisão de 85%, ou seja, com base em variáveis socioeconômicas é possível prever com 85% de precisão se um domicílio tem ou não DVD. O intervalo de confiança é bem pequeno (95% CI : (0.8488, 0.8569)), o que é bom.

CM$table
##           Reference
## Prediction     0     1
##          0  5573  2835
##          1  1578 20013

Temos dois tipos de erro, quando a referência diz que um domicilio tem um DVD e o modelo diz que não e quando a referência = não e o modelo = sim. Como diria um colega, se o modelo diz algo e a realidade diz outro é a realidade que tem um problema…

Importância das variáveis

plot of chunk unnamed-chunk-7

Um negócio bacana das florestas é poder observar a importância das variáveis e para esse modelo Renda vence de longe como sendo a mais importante, depois a UF (afinal estamos no Brasil e Renda+UF explicam praticamente tudo!) e só depois a média de idade dos filhos.

Uma coisa interessante é que as variáveis de filhos que criamos para este modelo (EH_FILHO = número de filhos no domicílio e EH_FILHO_6 = Número de filhos com 6 anos ou menos) não fazem muita diferença no modelo…

Pelo visto ninguém vai comprar um aparelho de DVD só pela Galinha Pintadinha…