Tijdens de coronaperiode gaf ik weinig aandacht aan dit blog. Het lijkt me leuk om met open data een voorspelmodel te maken over de toekomstige ontwikkeling van het aantal bijstandsuitkeringen. Daarvoor gaan we eerst data ophalen bij het CBS. Het werkloosheidspercentage voegt mogelijk kwaliteit toe aan de forecast. Dit artikel beschrijft het downloaden van CBS data. Daarbij willen we inzicht in het percentage van de variantie van het aantal bijstandsuitkeringen verklaart kan worden via het werkloosheidspercentage. Hiervoor berekenen we de R Squared van de correlatie coefficient.
Python packages in dit artikel
import pandas as pd # data wrangling en analyse
from pandas.tseries.offsets import MonthEnd # bepalen maandeinde
from scipy import stats # correlatie bereken
import matplotlib.pyplot as plt # data visualisatie
import matplotlib.ticker as mtick # format % y-as
import seaborn as sns # data visualisatie
import cbsodata # ophalen data van CBS
# Instellingen
pd.set_option('display.max_colwidth', None) # volledige kolombreedte
pd.set_option('display.max_columns', None) # toon alle kolommen
sns.set_theme() # seaborn standaard thema voor grafieken
sns.set(font_scale=1.2) # groter lettertype
Open data downloaden met cbsodata
Voor Python gebruikers is het package cbsodata beschikbaar via pip. Zelf zoek ik eerst online naar voorbeelden van beschikbare datasets. Het is ook mogelijk om een lijst te downloaden met get_table_list(). Elke dataset wordt geïdentificeerd door een eigen unieke sleutel: table_id. Zo kan je ‘Bijstandsuitkeringen per regio’ downloaden met 82015ned en ‘arbeidsdeelname en werkloosheid’ met 82015ned.
cbs_df = pd.DataFrame(cbsodata.get_table_list())
Het downloaden van een dataset doe je met get_data()
bijstandscijfers_df = pd.DataFrame(
cbsodata.get_data(
table_id='82015ned'))
arbeidsdeelname_landelijk_df = pd.DataFrame(
cbsodata.get_data(
table_id='80590ned'))
Metadata beschrijft de data eigenschappen
Het fijne van CBS open data is dat je een beschrijving van de data kan opvragen. Hieronder gaan we de meta data ophalen over de dataset Arbeidsdeelname en werkloosheid per maand.
meta = pd.DataFrame(cbsodata.get_meta('80590ned', 'DataProperties'))
Verschillende onderwerpen
Tabel ‘80590ned’ bevat cijfers over de beroepsbevolking, werkzame beroepsbevolking, werkloze beroepsbevolking en werkloosheidspercentage. Het filteren van de onderwerpen uit de metadata doe je zo:
meta[meta.Type=='TopicGroup'][['Title','Description']].head(4)
meta[meta.Type=='TopicGroup'][['Title','Description']].head(4)
Meetwaarden selecteren
Een onderwerp wordt gevolgd door de bijbehorende meetwaarden. Je kan dus slicen vanaf de index van het onderwerp om de gewenste meetwaarden te vinden.
meta[['Type','Key','Title','Description','Datatype',
'Unit','Decimals','Default']].iloc[12:,:].head(3)
De naamgeving van de meetwaarden bevat het soort cijfer, maar niet het onderwerp. Tijdens het laatste stukje data wrangling vervangen we de kolomnamen met de functionele naam. Deze staat in de kolom ‘Title’ van de metadata.
Dimensies
De meetwaarden worden door het CBS berekent vanuit verschillende invalshoeken. Per dimensie wordt een nieuw record toegevoegd aan de dataset met de meetwaarden. Elke dimensie heeft een eigen kolom, die je kan gebruiken om te filteren. Bijvoorbeeld om de arbeidsdeelname cijfers te selecteren voor een bepaalde leeftijdsgroep. De dataset bevat ook een record met de cijfers over de volledige populatie.
meta[meta.Type.isin(
['Dimension','TimeDimension'])
][['Key','Title']]
Landelijke maandcijfers
De dataset bevat cijfers over de volgende tijdsvakken: jaar, maand en kwartaal. Het is een beetje onhandig dat de datasets geen periode einddatum bevatten. Daarom converteren we de maandelijkse tijdsvakken, ‘2003 april’ naar de maand einddatum. De maand einddatum slaan we op de kolom ‘periode einddatum’ om haar vervolgens op te nemen in visualisaties en filters.
arbeidsdeelname_landelijk_df[
'Perioden'].unique()[15:25]
maanden = {'januari':'01','februari':'02','maart':'03','april':'04','mei':'05','juni':'06','juli':'07'
,'augustus':'08','september':'09','oktober':'10','november':'11','december':'12'}
def cbs_perioden_converteren(df):
# Periode type, let op! gaat ervan uit dat het altijd maand,kwartaal of jaar is
df['periode_type'] = 'maand'
df.loc[df['Perioden'].str.len()==4,'periode_type'] = 'jaar'
df.loc[df['Perioden'].str.contains('kwartaal',case=False),'periode_type'] = 'kwartaal'
# Maandeinddatum bepalen en datum dimensies toevoegen
# Let op! nu alleen nog voor periode type maand
df['periode_einddatum'] = df['Perioden'].str.slice(0,4) + '-' + df.Perioden.str.slice(5,99).replace(maanden) + '-' + '01'
df.loc[df['periode_type'].isin(['jaar','kwartaal']),'periode_einddatum'] =''
df['periode_einddatum'] = pd.to_datetime(df['periode_einddatum'], format='%Y-%m-%d', errors='ignore') + MonthEnd(0)
return df
bijstandscijfers_df = cbs_perioden_converteren(bijstandscijfers_df)
arbeidsdeelname_landelijk_df = cbs_perioden_converteren(arbeidsdeelname_landelijk_df)
Vervolgens kunnen we de landelijkecijfers per maand selecteren door te filteren op maandelijkse perioden en de juiste dimensie waarden. In de code hieronder zie je de filters, waarmee je de landelijke cijfers selecteert.
bijstandscijfers_df.query("RegioS=='Nederland' & UitkeringenInRelatieMetDeAOW=='Aan personen tot de AOW-leeftijd' & periode_type=='maand'", inplace=True)
arbeidsdeelname_landelijk_df.query("Geslacht=='Totaal' & Leeftijd == '15 tot 75 jaar' & periode_type == 'maand'", inplace=True)
Laatste stukje data wrangling
Voor we verder gaan met data analyse en visualisatie vereenvoudigen we deze stappen door de datasets te combineren. Het laatste stukje data wrangling heeft de volgende taken:
- Variabelen selecteren
- Technische namen hernoemen naar logische namen met .rename()
- Samenvoegen van de datasets met .merge()
- Sorteren op datum met .sort_values()
# Aantal bijstandsuitkeringen
bijstandscijfers_df = bijstandscijfers_df[['periode_einddatum','AantalBijstandsuitkeringen_1']].copy()
bijstandscijfers_df.rename(columns={'AantalBijstandsuitkeringen_1':'Aantal bijstandsuitkeringen'},inplace=True)
# Werkloosheidspercentage
# Geen zin om alle kolomnamen over te tikken, dan slicen :-)
arbeidsdeelname_meetwaarden = list(arbeidsdeelname_landelijk_df.columns)[10:12]
c = arbeidsdeelname_meetwaarden + ['periode_einddatum']
arbeidsdeelname_landelijk_df = arbeidsdeelname_landelijk_df[c]
arbeidsdeelname_landelijk_df.rename(
columns={'NietSeizoengecorrigeerd_7':'Werkloosheidspercentage',
'Seizoengecorrigeerd_8':'Werkloosheidspercentage seizoengecorrigeerd'},
inplace=True)
# Samenvoegen
kerncijfers_df = bijstandscijfers_df.merge(arbeidsdeelname_landelijk_df, on='periode_einddatum', how='left')
kerncijfers_df.sort_values(by='periode_einddatum', inplace=True)
kerncijfers_df.set_index('periode_einddatum', inplace=True)
Daarna kijken of het is gelukt 🙂
kerncijfers_df.tail(4)
Samenhang tussen tijdsreeksen verkennen met lijndiagram
Met beschrijvende statistieken en methoden kan je relatief snel veel gegevens tonen over jouw data. Er is echter niets beters dan een data visualisatie om je verhaal te vertellen of om feedback te vragen van domeinkenners. Vaak zie je al snel opvallende patronen of samenhang met bepaalde events. Bijvoorbeeld de corona hobbel vanaf begin 2020. Uit onderstaande grafiek kan je aflezen dat er voor corona een vertraging zit tussen werkgelegenheid en het aantal bijstandsuitkeringen. Hoe lang denk jij dat de vertraging duurt?
# Canvas
fig_dims = (16, 6)
fig, ax = plt.subplots(figsize=fig_dims)
fig.canvas.draw()
# Grafiek
ax1 = sns.lineplot(data=kerncijfers_df, x='periode_einddatum', y='Aantal bijstandsuitkeringen',
color='tab:orange', linewidth=3, ci=None,
label='Lopende bijstandsuitkeringen')
ax1.set_ylabel('Bijstandsuitkeringen', fontsize=14)
ax1.set_xlabel('')
ax1.yaxis.grid(False)
ax2 = ax.twinx()
ax2 = sns.lineplot(data=kerncijfers_df, x='periode_einddatum', y='Werkloosheidspercentage',
color='tab:gray', linewidth=3, ci=None, linestyle=(0,(1,1)),
label='Werkloosheid')
ax2.set_ylabel('Werkloosheid', fontsize=14)
ax2.yaxis.grid(False)
ax2.yaxis.set_major_formatter(mtick.PercentFormatter(decimals=0,
symbol=' %',
is_latex=False))
# Legenda
ax1.legend(title='', bbox_to_anchor=(0.01, 1.07, 0, 0), ncol=1,
loc='upper left', borderpad=0, borderaxespad=0, frameon=False,
fontsize=12)
ax2.legend(title='', bbox_to_anchor=(0.25, 1.07, 0, 0), ncol=1,
loc='upper left', borderpad=0, borderaxespad=0, frameon=False,
fontsize=12)
# Titel
transform = ax.transAxes
title_strings = ["VERTRAAGDE SAMENHANG ",
" BIJSTANDSUITKERINGEN ",
"EN ",
" WERKLOOSHEIDSPERCENTAGE ",
"TOT COVID-19"]
colors = ['black',
'tab:orange',
'black',
'tab:gray',
'black']
x_pos = 0.0
for string, col in zip(title_strings, colors):
t = ax.text(x_pos, 1.1, string, transform=transform, color=col, fontsize=14, fontweight='bold')
bbox = t.get_window_extent().transformed(transform.inverted())
x_pos = bbox.x1
x_pos = 0.0
Correlatie en R² (R squared) berekenen met stats uit SciPy
Als we naar de grafiek kijken dan lijkt er een vertraagde samenhang te zijn tussen het aantal bijstandsuitkeringen en het werkloosheidspercentage. Tijdens de coronacrisis lijkt deze samenhang te veranderen. Met SciPy kan je in Python correlatie berekenen volgens verschillende methodes Pearson, Spearman en Kendall’s tau. Hieronder maken we gebruik van pearsonr om de correlatie te bereken. Dit geeft ons de volgende inzichten:
- Voor corona sterke correlatie (0,87) met het werkloosheidspercentage van 24 maanden terug. Deze correlatie is erg laag tijdens corona en waarschijnlijk toeval (hoge p-waarde)
- Tijdens corona is correlatie met het actuele werkloosheidspercentage erg sterk. De samenhang is dus inderdaad (tijdelijk) veranderd
# Correlatie berekenen met stats uit SciPy
# Periode voor corona
kerncijfers_df['lag'] = kerncijfers_df.Werkloosheidspercentage.shift(periods=24)
stats.pearsonr(kerncijfers_df.dropna().loc['2011-01-01':'2020-02-01','Aantal bijstandsuitkeringen'],
kerncijfers_df.dropna().loc['2011-01-01':'2020-02-01','lag'])
# Vanaf corona
kerncijfers_df['lag'] = kerncijfers_df.Werkloosheidspercentage.shift(periods=24)
stats.pearsonr(kerncijfers_df.dropna().loc['2020-02-01':'2099-01-01','Aantal bijstandsuitkeringen'],
kerncijfers_df.dropna().loc['2020-02-01':'2099-01-01','lag'])
stats.pearsonr(kerncijfers_df.dropna().loc['2020-02-01':'2099-01-01','Aantal bijstandsuitkeringen'],
kerncijfers_df.dropna().loc['2020-02-01':'2099-01-01','Werkloosheidspercentage'])
Correlatiematrix met Pandas is eenvoudig
Het telkens opnieuw tikken van dezelfde code is omslachtig. Na twee jaar thuiswerken in coronatijd toont mijn toetsenbord samenhang met de grafiek. Je kan ook Pandas gebruiken om een correlatiematrix op te stellen. Het bepalen van voorliggende waarden (lag) in een ‘for’ loop maakt de code nog een stapje efficienter.
Als een toetsenbord kan glimlachen, dan zie je nu een smiley.
# Correlatie berekenen met Pandas
correlatie_matrix = kerncijfers_df.loc['2011-01-01':'2020-02-01',
['Aantal bijstandsuitkeringen','Werkloosheidspercentage']]
lags=[3, 6,12,24,36]
for l in lags:
c = 'ww_lag_' + str(l)
correlatie_matrix[c] = correlatie_matrix.Werkloosheidspercentage.shift(periods=l)
correlatie_matrix = correlatie_matrix.dropna().corr(method='pearson')[
['Aantal bijstandsuitkeringen','Werkloosheidspercentage']]
correlatie_matrix.round(2)
Mooi maken kan met een heatmap via seaborn.
cmap = sns.cubehelix_palette(start=.5, rot=-.5, as_cmap=True)
sns.heatmap(correlatie_matrix.iloc[1:10,0:1], annot=True,
cmap=cmap, square=True, xticklabels='')
plt.show()
Nu nog de R² (R squared) bepalen. Dit percentage wordt gebruikt als indicatie hoeveel de variatie X, de variatie van Y verklaart. Uit de berekening van R² volgt dat het werkloosheidspercentage van 24 maanden terug, voor 76% samenvalt met het aantal bijstandsuitkeringen gekeken naar de periode voorafgaande aan de coronacrisis.
# Correlatie
c = stats.pearsonr(kerncijfers_df.dropna().loc['2011-01-01':'2020-02-01','Aantal bijstandsuitkeringen'],
kerncijfers_df.dropna().loc['2011-01-01':'2020-02-01','lag'])
# Verklaarde variantie berekenen
r_squared = c[0] * c[0]
round(r_squared,2)
What’s next?
Nu je weet wat de correlatie is komt de vraag op: wat doe je ermee? Hoewel een verband in dit voorbeeld zeer aannemelijk is, toont een sterke correlatie geen oorzaak gevolg aan. Het is een indicatie van de sterkte van de onderlinge samenhang. Ik kies ervoor om het werkloosheidspercentage, met een lag van 24 maanden, op te nemen in een voorspelmodel voor toekomstige bijstandscijfers. Het ‘verklaart’ 76% van de variantie en is door de lag beschikbaar voor datapunten, die in de toekomst liggen (horizon van de voorspelling). In het volgende artikel over forecasting beschrijf ik hoe dat uitpakt.