5. Demographic data analysis

Here, we do an exploratory data analysis of the Brazilian PNADC - Pesquisa Nacional por Amostra de Domicílios Contínua.

# Download the PNADC dataset using the pynad package
# !pip install pynad
# !pynad pnadc microdados
# !pynad pnadc metadados
## (then unzip everything or the desired files)

import numpy as np
import scipy.stats as stats
import pandas as pd
import scipy as sp
import matplotlib.pyplot as plt
import pickle
import torch
import torch.nn as nn
import torch.nn.functional as F
import torch.optim as optim
import json
from collections import Counter

%matplotlib inline

# nan corrected Counter
def ncCounter(x):
    isnan = np.isnan(x)
    res = Counter(x[np.logical_not(isnan)])
    res[np.nan] = np.sum(isnan)
    return res

# weighted cumulative empirical distribution plot

def ecdf_plot(x, weights, ax, *args, **kwargs):
    x = np.array(x)
    w = np.array(weights)
    sort_indx = np.argsort(x)
    xc = np.concatenate(([min(x)], x[sort_indx], [max(x)*1.1]))
    yc = np.cumsum(w[sort_indx])/w.sum()
    yc = np.concatenate(([0, 0], yc))
    ax.step(xc, yc, *args, **kwargs)

# weighted boxplot function

def wboxplot(x, weights, ax):
    x = np.array(x)
    w = np.array(weights)

    sort_indx = np.argsort(x)
    x = x[sort_indx]
    w = w[sort_indx]

    wcs = np.cumsum(w)
    wsum = w.sum()
    median = x[np.searchsorted(wcs, 0.50 * wsum)]
    q1 = x[np.searchsorted(wcs, 0.25 * wsum)]
    q3 = x[np.searchsorted(wcs, 0.75 * wsum)]

    iqr = q3 - q1
    whishi = q3 + 1.5 * iqr
    whishi = max(whishi, x.max())
    whislo = q1 - 1.5 * iqr
    whislo = min(whislo, x.min())

    bxpstats = [{
        'med': median,
        'q1': q1,
        'q3': q3,
        'iqr': iqr,
        'whishi': whishi,
        'whislo': whislo,
    }]

    ax.bxp(bxpstats, showfliers=False)

## load the variables descriptions

with open("variaveis.pnadc.trimestral.json", 'r') as f:
    variables = json.load(f)

print(list(variables))
variables['v2009']

# You can translate a variable description to English with:
# !pip install translate
# from translate import Translator
# translator= Translator(to_lang="en", from_lang='pt')
# translator.translate(str(variables['v2009']))

['ano', 'trimestre', 'uf', 'capital', 'rm_ride', 'upa', 'estrato', 'v1008', 'v1014', 'v1016', 'v1022', 'v1023', 'v1027', 'v1028', 'v1029', 'posest', 'v2001', 'v2003', 'v2005', 'v2007', 'v2008', 'v20081', 'v20082', 'v2009', 'v2010', 'v3001', 'v3002', 'v3002a', 'v3003', 'v3003a', 'v3004', 'v3005', 'v3005a', 'v3006', 'v3006a', 'v3007', 'v3008', 'v3009', 'v3009a', 'v3010', 'v3011', 'v3011a', 'v3012', 'v3013', 'v3013a', 'v3013b', 'v3014', 'v4001', 'v4002', 'v4003', 'v4004', 'v4005', 'v4006', 'v4006a', 'v4007', 'v4008', 'v40081', 'v40082', 'v40083', 'v4009', 'v4010', 'v4012', 'v40121', 'v4013', 'v40132', 'v40132a', 'v4014', 'v4015', 'v40151', 'v401511', 'v401512', 'v4016', 'v40161', 'v40162', 'v40163', 'v4017', 'v40171', 'v401711', 'v4018', 'v40181', 'v40182', 'v40183', 'v4019', 'v4020', 'v4021', 'v4022', 'v4024', 'v4025', 'v4026', 'v4027', 'v4028', 'v4029', 'v4032', 'v4033', 'v40331', 'v403311', 'v403312', 'v40332', 'v403321', 'v403322', 'v40333', 'v403331', 'v4034', 'v40341', 'v403411', 'v403412', 'v40342', 'v403421', 'v403422', 'v4039', 'v4039c', 'v4040', 'v40401', 'v40402', 'v40403', 'v4041', 'v4043', 'v40431', 'v4044', 'v4045', 'v4046', 'v4047', 'v4048', 'v4049', 'v4050', 'v40501', 'v405011', 'v405012', 'v40502', 'v405021', 'v405022', 'v40503', 'v405031', 'v4051', 'v40511', 'v405111', 'v405112', 'v40512', 'v405121', 'v405122', 'v4056', 'v4056c', 'v4057', 'v4058', 'v40581', 'v405811', 'v405812', 'v40582', 'v405821', 'v405822', 'v40583', 'v405831', 'v40584', 'v4059', 'v40591', 'v405911', 'v405912', 'v40592', 'v405921', 'v405922', 'v4062', 'v4062c', 'v4063', 'v4063a', 'v4064', 'v4064a', 'v4071', 'v4072', 'v4072a', 'v4073', 'v4074', 'v4074a', 'v4075a', 'v4075a1', 'v4076', 'v40761', 'v40762', 'v40763', 'v4077', 'v4078', 'v4078a', 'v4082', 'vd2002', 'vd2003', 'vd2004', 'vd3004', 'vd3005', 'vd3006', 'vd4001', 'vd4002', 'vd4003', 'vd4004', 'vd4004a', 'vd4005', 'vd4007', 'vd4008', 'vd4009', 'vd4010', 'vd4011', 'vd4012', 'vd4013', 'vd4014', 'vd4015', 'vd4016', 'vd4017', 'vd4018', 'vd4019', 'vd4020', 'vd4023', 'vd4030', 'vd4031', 'vd4032', 'vd4033', 'vd4034', 'vd4035', 'vd4036', 'vd4037']

{'parte': 'parte 2 - características gerais dos moradores',
 'colunas': [92, 94, 3],
 'bytes': 2,
 'quesito': 9,
 'desc': 'idade do morador na data de referência',
 'periodo': '1º tri/2012 - atual',
 'valores': {'0 a 130': 'idade (em anos)'}}

# text and v1XXX -> "Parte 1 - Identificação e controle"
# v2XXX -> "Parte 2 - Características gerais dos moradores"
# v3XXX -> "Parte 3 - Características de educação para os moradores de 5 anos ou mais de idade"
# v4XXX -> "Parte 4 - Características de trabalho das pessoas de 14 anos ou mais de idade"
# v4009 <= x <= v4064a -> "2 - Pessoas ocupadas"
# v4071 <= x <= v4082 -> "3 - Pessoas não ocupadas"
# vdXXXX -> "Variáveis Derivadas"

# Some interesting ones:

# v1008 numero selecao municipio
# v1014 painel - grupo da amostra
# v1016 Número da entrevista no domicílio (1 a 5)

# v2001 Número de pessoas no domicílio (1 a 30)
# v2003 Número de ordem da pessoa no domicilio (1 a 30)
# v2005 Condição da pessoa no domicílio

# v1028 peso
# v1029 projeção populacional

# v2009 Analfabetismo

# v4009 numero de empregos
# v403312 renda emprego principal

# Load datasets

df211 = pd.read_csv("microdados.pnadc.trimestral.2021.1.csv")
df194 = pd.read_csv("microdados.pnadc.trimestral.2019.4.csv")

# number of jobs per person
(
 ncCounter(df211.v4009[np.logical_and(df211.v2009>=18, df211.v2009<60)]),
 ncCounter(df194.v4009[np.logical_and(df194.v2009>=18, df194.v2009<60)])
)

(Counter({1.0: 108488, 2.0: 3176, 3.0: 181, nan: 74750}),
 Counter({1.0: 200165, 2.0: 7031, 3.0: 447, nan: 110179}))

# income distribution plots for main job
dff = df211[np.logical_and(df211.v2009>=18, df211.v2009<60)].copy()

dff.loc[dff.v403312.isna(), 'v403312'] = 0
dff.loc[dff.v403312 < 89, 'v403312'] = 89 # assume everyone gets at least 89 BRL

# some other possible filters
# dff = dff.loc[dff.v403312 >= 5000]
# dff = dff.loc[dff.v403312 <= 8000]

print(np.average(dff.v403312, weights=dff.v1028))
print(np.sqrt(np.average(dff.v403312**2, weights=dff.v1028) - np.average(dff.v403312, weights=dff.v1028)**2))

fig, axes = plt.subplots(1,2,figsize=(14,5))
ecdf_plot(dff.v403312, dff.v1028, axes[0])
axes[0].set_xscale('symlog')
wboxplot(dff.v403312, dff.v1028, axes[1])
axes[1].set_yscale('log')
axes[1].set_ylim(10**-2, 10**6)
axes[1].set_xticks([])

1502.2946012826508
3060.08505198942

[]

# illiteracy rate for people of legal age

adultos = df211[df211.v2009>=18]
assert adultos.v3001.isna().sum()==0
np.average(adultos.v3001==2, weights=adultos.v1028)

0.05986003807566674

# person-dwelling-interview unique id (id unico da pessoa, domicilio, entrevista):
# ['upa', 'v1008', 'v1014', 'v1016', 'v2003']

# dwelling-interview unique id
# id unico do domicilio, entrevista
# ['upa', 'v1008', 'v1014', 'v1016']

print(len(df211))
len(np.unique(df211[['upa', 'v1008', 'v1014', 'v1016', 'v2003']].to_numpy(),axis=0))

319898

319898

# Counter of "situation" the person hold in the house
# e.g.:
# '1': 'person responsible for the household',
# '2': 'spouse or partner of different sex',
# '3': 'spouse or partner of the same sex',
# '4': 'child of both',
# '5': 'child of the spouse only',
# '6': 'stepson'
print(variables['v2005'])
dict(sorted(Counter(df211.v2005).items()))

{'parte': 'parte 2 - características gerais dos moradores', 'colunas': [81, 82, 2], 'bytes': 1, 'quesito': 5, 'desc': 'condição no domicílio', 'periodo': '1º tri/2012 - atual', 'valores': {'1': 'pessoa responsável pelo domicílio', '2': 'cônjuge ou companheiro(a) de sexo diferente', '3': 'cônjuge ou companheiro(a) do mesmo sexo', '4': 'filho(a) do responsável e do cônjuge', '5': 'filho(a) somente do responsável', '6': 'enteado(a)', '7': 'genro ou nora', '8': 'pai, mãe, padrasto ou madrasta', '9': 'sogro(a)', '10': 'neto(a)', '11': 'bisneto(a)', '12': 'irmão ou irmã', '13': 'avô ou avó', '14': 'outro parente', '15': 'agregado(a) - não parente que não compartilha despesas', '16': 'convivente - não parente que compartilha despesas', '17': 'pensionista', '18': 'empregado(a) doméstico(a)', '19': 'parente do(a) empregado(a) doméstico(a)'}}

{1: 110286,
 2: 67112,
 3: 261,
 4: 67069,
 5: 35531,
 6: 3408,
 7: 2836,
 8: 6969,
 9: 969,
 10: 14345,
 11: 423,
 12: 4933,
 13: 326,
 14: 4192,
 15: 580,
 16: 563,
 17: 9,
 18: 84,
 19: 2}

# average number of "children of either" per house (but not weight-adjusted!)
np.mean(df211.groupby(['upa', 'v1008', 'v1014', 'v1016'])[['v2005']].agg(lambda x: sum(x==4)+sum(x==5)+sum(x==6)))

v2005    0.96121
dtype: float64