Open Data: Berlin Kitas

Kitas per District

We start by looking into the first layer of geo-location granularity: district (there are 12 districts in Berlin).

fig, ax = plt.subplots()
kitas_df \
    .groupby(['district']) \
    .agg(n=('id', 'count')) \
    .reset_index(drop=False) \
    .sort_values('n', ascending=False) \
    .pipe((sns.barplot, 'data'), 
        x='n', 
        y='district',
        color=sns_c[0],
        edgecolor=sns_c[0],
        ax=ax
    )
ax.set(
    title='Number of Kitas per District in Berlin', 
    xlabel='number of kitas', 
    ylabel='district'
);

We generate the same plot but showing the shares instead.

fig, ax = plt.subplots()
kitas_df \
    .groupby(['district']) \
    .agg(n=('id', 'count')) \
    .assign(share= lambda x: x['n'] / x['n'].sum()) \
    .reset_index(drop=False) \
    .sort_values('n', ascending=False) \
    .pipe((sns.barplot, 'data'), 
        x='share', 
        y='district',
        color=sns_c[1],
        edgecolor=sns_c[1],
        ax=ax
    )
ax.xaxis.set_major_formatter(mtick.FuncFormatter(lambda y, _: '{:0.0%}'.format(y)))
ax.set(
    title='Share of Kitas per District in Berlin', 
    xlabel='share of kitas', 
    ylabel='district'
);

Pankow is a rather big district compared to Mitte and Friedrichshain-Kreuzberg area-wise. In order to run a fair comparison across districts, we need to control by size or population (which we will do bellow).

Kitas per PLZ

Now we dig deeper into plz level. A natural question is: How many Kitas each plz has? To answer it let us look into the distribution (again, in order to make a fair comparison we need to control by plz population or size).

fig, ax = plt.subplots()
kitas_df \
    .groupby(['district', 'plz']) \
    .agg(n=('name', 'count')) \
    .pipe((sns.histplot, 'data'),
        x='n', 
        bins=25, 
        kde=True,
        ax=ax
    )
ax.set(
    title='Number of Kitas per PLZ', 
     xlabel='number of kitas per plz',
     xlim=(0, None)
);

We see there is a long tail. Let us see the plz with more Kitas:

fig, ax = plt.subplots(figsize=(12, 6))
kitas_df \
    .groupby(['district', 'plz']) \
    .agg(n=('id', 'count')) \
    .reset_index(drop=False) \
    .assign(plz = lambda x: 'plz_' + x['plz']) \
    .sort_values('n', ascending=False) \
    .head(40) \
    .pipe((sns.barplot, 'data'),  
        x='plz',
        y='n',
        hue='district',
        dodge=False,
        ax=ax
    )
ax.tick_params(axis='x', labelrotation=90)
ax.set(
    title='Number of Kitas per PlZ in Berlin (Top 40)', 
     xlabel='plz', 
     ylabel='number of kitas', 
);

Let us locate plz = 10437 using Google maps.

This is a very trendy area in Berlin, specially for young families.

Next, let us plot the distribution of the number of Kitas per plz split by district.

g = kitas_df \
    .groupby(['district', 'plz']) \
    .agg(n=('name', 'count')) \
    .pipe((sns.displot, 'data'), 
        x='n',
        hue='district',
        kind='kde',
        rug=True, 
        height=5,
        aspect=1.3,
        palette='Paired'
    )
ax = g.axes.flatten()[0]
ax.set(
    title='Number of Kitas per PlZ', 
     xlabel='number of kitas per plz',
     xlim=(0, None)
);

fig, ax = plt.subplots(figsize=(12, 6))
kitas_df \
    .groupby(['district', 'plz']) \
    .agg(n=('id', 'count')) \
    .reset_index(drop=False) \
    .sort_values('district') \
    .pipe((sns.boxplot, 'data'), 
        x='district', 
        y='n',
        color=sns_c[0],
        saturation=1.0,
        ax=ax
    )
ax.tick_params(axis='x', labelrotation=90)
ax.set(
    title='Number of Kitas per PlZ in Berlin within each District', 
     ylabel='number of kitas per plz', 
);

It seems that, on average, Friedrichshain-Kreuzberg have more Kitas per plz and Pankow has the highest variance (again, we need to control by plz size to be more precise).

Spots per Kita

Let us now analyze the number of spots data.

# Some descriptive stats.
kitas_df['spots'].describe()

count 2736.000000 mean 67.617690 std 55.624663 min 0.000000 25% 25.000000 50% 45.000000 75% 95.000000 max 320.000000 Name: spots, dtype: float64

Note that there is a Kita wih spots=0. I am not sure about the reason behind this.

kitas_df.query('spots == 0')

We plot the spots distribution:

spots_mean = kitas_df['spots'].describe()['mean']
spots_median = kitas_df['spots'].describe()['50%']

fig, ax = plt.subplots()
sns.histplot(x='spots', data=kitas_df, kde=True, color=sns_c[3], ax=ax)
ax.axvline(x=spots_mean, color=sns_c[1], linestyle='--', label=f'mean = {spots_mean: 0.0f}')
ax.axvline(x=spots_median, color=sns_c[4], linestyle='--', label=f'median = {spots_median: 0.0f}')
ax.legend(loc='upper right')
ax.set(
    title='Number of Spots per Kita', 
     xlabel='number of spots per kita',
     xlim=(0, None)
);

Here is a list of the top 20 kitas with more spots:

kitas_df \
    [['name', 'district', 'plz', 'spots']] \
    .sort_values('spots', ascending=False) \
    .head(20) 

Let us see the spots distributoin per district:

g = kitas_df \
    .pipe((sns.displot, 'data'), 
        x='spots',
        hue='district',
        kind='kde',
        rug=True, 
        height=5,
        aspect=1.3,
        palette='Paired'
    )
ax = g.axes.flatten()[0]
ax.set(
    title='Number of Spots per Kita', 
    xlabel='number of spots', 
    xlim=(0, None)
);

fig, ax = plt.subplots(figsize=(12, 6))
kitas_df \
    .sort_values('district') \
    .pipe((sns.boxplot, 'data'), 
        x='district', 
        y='spots', 
        color=sns_c[3],
        saturation=1.0,
        ax=ax
    )
ax.axhline(y=kitas_df['spots'].describe()['mean'], color=sns_c[1], linestyle='--', label='mean')
ax.axhline(y=kitas_df['spots'].describe()['50%'], color=sns_c[4], linestyle='--', label='median')
ax.tick_params(axis='x', labelrotation=90)
ax.legend(loc='upper left')
ax.set(
    title='Number of Spots per Kita', 
     ylabel='number of spots', 
);

Note that Pankow, Mitte and Friedrichshain-Kreuzberg have a lowe number of spots per Kita (these are the top 3 districts with mote Kitas). We want to see if there is a linear relation between the number of Kitas per plz and the number (median) of spots per Kita.

g = kitas_df \
    .groupby(['district', 'plz']) \
    .agg({'num_kitas_plz': np.mean, 'spots': np.median}) \
    .pipe((sns.jointplot, 'data'), 
        x='num_kitas_plz',  
        y='spots',
        kind='reg', 
        order=1,
        height=8,
        ratio=3, 
        marginal_kws=dict(color=sns_c[1])
    )
g.plot_joint(sns.kdeplot, zorder=0, levels=6, color=sns_c[3], alpha=0.7)
g.set_axis_labels(xlabel='number of kitas per plz', ylabel='median spots per plz')
g.fig.suptitle('Number of Kitas vs Spots per PLZ', y=1.03);

It seems there is a negative relation between them, i.e. plz with more Kitas seem to have less spots per Kita. Let us plot this per district:

g = kitas_df \
    .groupby(['district', 'plz'], as_index=False) \
    .agg({'num_kitas_plz': np.mean, 'spots': np.median}) \
    .pipe((sns.FacetGrid, 'data'),
        col='district', 
        col_wrap=4, 
        height=3,
    )
g.map(sns.regplot, 'num_kitas_plz', 'spots')
g.fig.suptitle('Number of Kitas vs Spots per PLZ', y=1.03);

Carriers

Now we want to study the carrier feature. First, let us count the number of kitas per carrier_type.

fig, ax = plt.subplots()
kitas_df \
    .groupby(['carrier_type']) \
    .agg(n=('id', 'count')) \
    .reset_index() \
    .sort_values('n', ascending=False) \
    .assign(share= lambda x: x['n'] / x['n'].sum()) \
    .pipe((sns.barplot, 'data'), 
        x='share', 
        y='carrier_type',
        color=sns_c[2],
        edgecolor=sns_c[2],
        ax=ax
    )
fmt = lambda y, _: f'{y :0.0%}'
ax.xaxis.set_major_formatter(mtick.FuncFormatter(fmt))
ax.set(
    title='Number of Kitas per Carrier Type', 
    xlabel='number of kitas', 
    ylabel='carrier type'
);

Here are some explanations on the carrier_type:

• Arbeiterwohlfahrt: Workers welfare
• Betriebskita: Company daycare
• Deutscher Caritasverband: German Caritas Association
• Eigenbetrieb: Own operation
• EKT (Eltern-Initiativ-Kindertagesstätte): Parents initiative.
• Diakonisches Werk: Diaconal work
• private Kita ohne Abschluss der Rahmenvereinbarung: private daycare without signing the framework agreement
• Sonstiger freier Träger: Other private carrier

# Note that there is. a 1-1 correspondence between `carrier_number` and `carrier_name`.
kitas_df \
    .groupby('carrier_number')\
    .agg({'carrier_name': 'nunique'}) \
    .query('carrier_name > 1') \
    .shape

(0, 1)

Now we plot the carrier_type distriibution per district.

fig, ax = plt.subplots(figsize=(15, 6))
kitas_df \
    .groupby(['district', 'carrier_type']) \
    .agg(n=('id', 'count')) \
    .reset_index(drop=False) \
    .sort_values('n', ascending=False) \
    .pipe((sns.barplot, 'data'), 
        y='n', 
        x='district',
        hue='carrier_type',
        edgecolor='black',
        dodge=True,
        ax=ax
    )
ax.tick_params(axis='x', labelrotation=90)
ax.set(
    title='Number of Kitas per District in Berlin (per Carrier Type)', 
    xlabel='district', 
    ylabel='number of kitas'
);

We can also compute carrier_type share per district:

fig, ax = plt.subplots()
cmap = sns.light_palette(sns_c[0])
fmt = lambda y, _: f'{y :0.0%}'
kitas_df \
    .groupby(['district', 'carrier_type']) \
    .agg(n=('id', 'count')) \
    .reset_index(drop=False) \
    .assign(
        num_kitas_district = lambda x: x.groupby(['district'])['n'].transform(np.sum),
        share = lambda x: x['n'].div(x['num_kitas_district'])
    ) \
    .pivot(index='district', columns='carrier_type', values='share') \
    .fillna(0.0) \
    .pipe((sns.heatmap, 'data'), 
        vmin=0.0,
        vmax=1.0,
        cmap=cmap,
        linewidths=0.2, 
        linecolor='black',
        annot=True, 
        fmt='0.2%',
        cbar_kws={'format': mtick.FuncFormatter(fmt)},
        ax=ax
    )
ax.set(title='Carrier Type Share per District');

It is very interesting to see how high the share of EKT in Charlottenburg-Wilmersdorf and Friedrichshain-Kreuzberg.

Now we want to investigate the carriers coverage on the number of kitas. The questions we want to answer are:

1. Do carriers with more than one Kita dominate the Kita distribution?

2. Do we see a difference across districts?

To answer both question we need to compute how many Kitas each carrier has. In addition,we sort them out by the number of Kitas calculate the (cumulative) share on the number of Kitas. This will allow us to answer the first question.

carrier_df = kitas_df \
    .groupby(['carrier_name']) \
    .agg(num_kitas=('name', pd.Series.nunique)) \
    .sort_values('num_kitas', ascending=False) \
    .assign(
        single_kita = lambda x: x['num_kitas'] < 2,
        share = lambda x: x['num_kitas'] / x['num_kitas'].sum(),
        cumsum_share = lambda x: x['share'].cumsum(),
        idx= lambda x: range(x.shape[0])
    )

carrier_df.head()

Let us compute the distribution of carriers with a single Kita vs carriers with more than one Kita.

carrier_df \
    .groupby('single_kita')\
    .agg(n=('idx', 'count')) \
    .assign(share=lambda x: x['n'].div(x['n'].sum()))

Next we plot the cumulative share data on the number of Kitas:

k = 200 
k_cumsum_share = carrier_df.query(f'idx == {k - 1}')['cumsum_share'].values[0]
idx_single_kita = carrier_df.query('single_kita').iloc[0]['idx']
cumsum_share_single_kita = carrier_df.query('single_kita').iloc[0]['cumsum_share']

fig, ax = plt.subplots()
sns.lineplot(x='idx', y='cumsum_share', data=carrier_df, color=sns_c[0], label='cumulative sum share', ax=ax)
ax.axvline(x=idx_single_kita, color=sns_c[4], linestyle='--', label=f'{idx_single_kita} carriers have more than one Kita ({cumsum_share_single_kita: 0.2%})')
ax.axhline(y=cumsum_share_single_kita, color=sns_c[4], linestyle='--')
ax.axvline(x=k, color=sns_c[3], linestyle='--', label=f'{k} carriers account for {k_cumsum_share: 0.2%} of the Kitas')
ax.axhline(y=k_cumsum_share , color=sns_c[3], linestyle='--')
ax.legend(loc='lower right')
ax.set(
    title='Cumulative Share of Top Kita Carriers', 
    xlabel='carrier rank (based on number of Kitas)',
    ylabel='cumsum share'
);

We have two main takeaways from this plot:

1. 356 carriers, which are ~ 30% of the total number of carriers, cover ~ 67% of the number of Kitas.

2. The top 200 carriers, which are ~ 16% of the total number of carriers, cover more than half (~ 55%) of the number of Kitas.

Hence, the carriers. with more than one Kita have a very high coverage.

# We merge the carriers data to. the Kitas data.
kitas_df = pd.merge(
    left=kitas_df, 
    right=carrier_df['single_kita'].reset_index(), 
    on='carrier_name', 
    how='left'
)

Now let us answer the second question by looking into the districts.

fig, ax = plt.subplots()
kitas_df \
    .groupby(['district', 'single_kita']) \
    .agg(n=('id', 'count')) \
    .reset_index(drop=False) \
    .sort_values('n', ascending=False) \
    .pipe((sns.barplot, 'data'), 
        x='n', 
        y='district',
        hue='single_kita',
        edgecolor='black',
        palette=[sns_c[3], sns_c[0]],
        ax=ax
    )
ax.set(
    title='Number of Kitas per District in Berlin (Carrier with/without Single Kita)', 
    xlabel='number of kitas', 
    ylabel='district'
);

We now plot the relative share over districts of number of Kitas ans spots (which are of course correlated):

carrier_kitas_df = kitas_df \
    .groupby(['district', 'single_kita']) \
    .agg(n=('id', 'count'), spots=('spots', np.sum)) \
    .reset_index(drop=False) \
    .sort_values('district') \
    .assign(
        total_n = lambda x: x.groupby('district')['n'].transform(np.sum),
        total_spots = lambda x: x.groupby('district')['spots'].transform(np.sum),
        share_n = lambda x: x['n'].div(x['total_n']),
        share_spots = lambda x: x['spots'].div(x['total_spots']),
    )

fig, ax = plt.subplots(1, 2, constrained_layout=True)
cmap0 = sns.light_palette(sns_c[1])
cmap1 = sns.light_palette(sns_c[4])
fmt = lambda y, _: f'{y :0.0%}'
carrier_kitas_df \
    .pivot(index='district', columns='single_kita', values='share_n') \
    .pipe((sns.heatmap, 'data'), 
        vmin=0.0,
        vmax=1.0,
        cmap=cmap0,
        linewidths=0.2, 
        linecolor='black',
        annot=True, 
        fmt='0.2%',
        cbar_kws={'format': mtick.FuncFormatter(fmt)},
        ax=ax[0]
    )

carrier_kitas_df \
    .pivot(index='district', columns='single_kita', values='share_spots') \
    .pipe((sns.heatmap, 'data'), 
        vmin=0.0,
        vmax=1.0,
        cmap=cmap1,
        linewidths=0.2, 
        linecolor='black',
        annot=True, 
        fmt='0.2%',
        cbar_kws={'format': mtick.FuncFormatter(fmt)},
        ax=ax[1]
    )

ax[0].set(title='Number of Kitas Share', xlabel='carrier single kita')
ax[1].set(title='Number of Spots Share', xlabel='carrier single kita', ylabel='');

• In Charlottenburg-Wilmersdorf the share of carriers with/without single Kitas is quite balanced whereas in Marzahn-Hellersdorf and Treptow-Köpenick is quite unbalanced.

This ends the an initial exploratory data analysis of the Kitas data. There are still many interesting things to look, but we will leave them for a second iteration.

Linear Regression Model

A natural question is: How do the relations above depend on the district? We of course expect different demographic distributions over different areas of the city. In this first iteration we run a simple linear regression model to estimate the variance explained of num_kitas_plz as a linear function of einwohner and district.

import statsmodels.formula.api as smf

model = smf.ols(formula='num_kitas_plz ~ einwohner + C(district)', data=plz_df)
result = model.fit()
print(result.summary())

    OLS Regression Results                            
    ==============================================================================
    Dep. Variable:          num_kitas_plz   R-squared:                       0.455
    Model:                            OLS   Adj. R-squared:                  0.423
    Method:                 Least Squares   F-statistic:                     14.00
    Date:                Sat, 19 Sep 2020   Prob (F-statistic):           6.33e-21
    Time:                        10:44:02   Log-Likelihood:                -696.16
    No. Observations:                 214   AIC:                             1418.
    Df Residuals:                     201   BIC:                             1462.
    Df Model:                          12                                         
    Covariance Type:            nonrobust                                         
    ===========================================================================================================
                                                  coef    std err          t      P>|t|      [0.025      0.975]
    -----------------------------------------------------------------------------------------------------------
    Intercept                                   0.8849      1.529      0.579      0.563      -2.130       3.900
    C(district)[T.Friedrichshain-Kreuzberg]     6.0192      2.399      2.509      0.013       1.289      10.750
    C(district)[T.Lichtenberg]                 -3.8870      2.365     -1.644      0.102      -8.550       0.776
    C(district)[T.Marzahn-Hellersdorf]         -6.0663      2.304     -2.633      0.009     -10.610      -1.523
    C(district)[T.Mitte]                        2.2961      1.903      1.206      0.229      -1.457       6.049
    C(district)[T.Neukölln]                    -1.1441      2.067     -0.554      0.581      -5.219       2.931
    C(district)[T.Pankow]                       1.1004      2.062      0.534      0.594      -2.966       5.167
    C(district)[T.Reinickendorf]               -3.6003      2.117     -1.700      0.091      -7.776       0.575
    C(district)[T.Spandau]                     -3.0675      2.293     -1.338      0.183      -7.589       1.454
    C(district)[T.Steglitz-Zehlendorf]         -3.5661      1.932     -1.846      0.066      -7.376       0.244
    C(district)[T.Tempelhof-Schöneberg]        -2.4554      1.751     -1.402      0.162      -5.908       0.997
    C(district)[T.Treptow-Köpenick]            -1.0303      2.120     -0.486      0.627      -5.210       3.149
    einwohner                                   0.0008   8.15e-05      9.586      0.000       0.001       0.001
    ==============================================================================
    Omnibus:                       13.969   Durbin-Watson:                   1.951
    Prob(Omnibus):                  0.001   Jarque-Bera (JB):               30.795
    Skew:                           0.237   Prob(JB):                     2.06e-07
    Kurtosis:                       4.797   Cond. No.                     2.02e+05
    ==============================================================================
    
    Warnings:
    [1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
    [2] The condition number is large, 2.02e+05. This might indicate that there are
    strong multicollinearity or other numerical problems.

We now look into the model predictions.

# Collect predictions.
plz_model_df = pd.concat([plz_df, result.get_prediction(plz_df).summary_frame()], axis=1)
# Compute errors.
plz_model_df = plz_model_df.assign(error = lambda x: x['mean'] - x['num_kitas_plz'])

fig, ax = plt.subplots(figsize=(10, 8))
sns.scatterplot(
    x='num_kitas_plz',
    y='mean',
    data=plz_model_df,
    hue='district',
    palette='Paired',
    edgecolor='black',
    linewidth=0.5,
    size='spots_plz',
    ax=ax
)
ax.legend(**{'loc': 'center left', 'bbox_to_anchor': (1, 0.5)})
ax.set(
    title='Number of Kitas Linear Model',
    xlabel='number of kitas',
    ylabel='prediction'
);

Finally, we plot the errors distribution:

error_mean = plz_model_df['error'].mean()
error_std = plz_model_df['error'].std()

fig, ax = plt.subplots(nrows=2, ncols=1, figsize=(12, 7), constrained_layout=True)
sns.histplot(x='error', data=plz_model_df, color=sns_c[9], kde=True, ax=ax[0])
ax[0].axvline(x=error_mean, color='black', linestyle='--', label=f'$\mu$')
ax[0].axvline(x=error_mean + 2*error_std, color='gray', linestyle='--', label=f'$\mu\pm2\sigma$')
ax[0].axvline(x=error_mean - 2*error_std, color='gray', linestyle='--')
sns.scatterplot(x='einwohner', y='error', data=plz_model_df, hue='district', palette='Paired', edgecolor='black', ax=ax[1])
ax[1].axhline(y=error_mean, color='black', linestyle='--', label=f'$\mu$')
ax[1].axhline(y=error_mean + 2*error_std, color='gray', linestyle='--', label=f'$\mu\pm2\sigma$')
ax[1].axhline(y=error_mean - 2*error_std, color='gray', linestyle='--')
ax[0].legend()
ax[1].legend(**{'loc': 'center left', 'bbox_to_anchor': (1, 0.5)})
ax[0].set(title='Model Erros Distribution');

The model is not perfect, but is a good place to start for a deeper analysis: ideally more data about the district demographics will be added and Bayesian methods like the ones presented in Statistical Rethinking by Richard McElreath.