<jupyter_start><jupyter_text>Calculate doctor capacity<jupyter_code>import pandas as pd
import numpy as np
import scipy.sparse
from os.path import join
from datetime import datetime
import matplotlib.pyplot as plt
data_path = 'data'  
file_name = 'matched_and_imputed_doctors.csv'
# load the matched and filtered doctor data (only doctors with > 1000 patients / year)
opening_hours = pd.read_csv(join(data_path, file_name))

# calculate quarterly patient numbers (as in Ruggiero's PNAS paper)
opening_hours['N_total_quarterly_patients'] = \
            opening_hours['N_total_yearly_patients'] / 4
opening_hours['N_total_quarterly_patients'] = \
            opening_hours['N_total_quarterly_patients'].round().astype(int)
opening_hours['N_unique_quarterly_patients'] = \
            opening_hours['N_unique_yearly_patients'] / 4
opening_hours['N_unique_quarterly_patients'] = \
            opening_hours['N_unique_quarterly_patients'].round().astype(int)

specs = opening_hours['specialization'].unique()<jupyter_output><empty_output><jupyter_text>Model 1: flat capacity estimate<jupyter_code># capacity increase factor: all capacities are increased by 10%
C = 0.1
opening_hours['capacity_total_flat'] = \
            opening_hours['N_total_quarterly_patients'] * (1 + C)
opening_hours['capacity_total_flat'] = \
            opening_hours['capacity_total_flat'].round().astype(int)
opening_hours['capacity_unique_flat'] = \
            opening_hours['N_unique_quarterly_patients'] * (1 + C)
opening_hours['capacity_unique_flat'] = \
            opening_hours['capacity_unique_flat'].round().astype(int)<jupyter_output><empty_output><jupyter_text>Model 2: opening hour based capacity For opening hour intervals of 5 opening hours: calculate the median patient capacity (patients per quarter) of the upper 10% of doctors. Assign this median capacity to all doctors that do not belong to the doctors within the upper 10% of capacity.<jupyter_code>opening_hours['capacity_total_hour-based'] = opening_hours['N_total_quarterly_patients']
opening_hours['capacity_unique_hour-based'] = opening_hours['N_unique_quarterly_patients']

# threshold for the upper percentile: the "maximum capacity" will be calculated
# as the median capacity of doctors in the upper (1-threshold) % of 
# doctors in a given speciality and hour bracket. The capacities of the doctors
# in the lower threshold % will be set to this maximum capacity
threshold = 0.9   ###0.9 original main setting

# hour bracket for which upper capacities will be calculated
hour_range_interval = 5
hour_ranges = np.arange(0, opening_hours['total_weekly_hours'].max(), hour_range_interval)

# iterate through all specializations
for spec in specs:
    # extract data from doctors from the current specialization
    spec_data = opening_hours[opening_hours['specialization'] == spec]
    # dictionary of the upper percentile medians in each hour-bracket for
    # the given speciality
    upper_percentile_medians = {i:np.nan for i in hour_ranges}
    
    # iterate through all hour brackets
    for iv_start in hour_ranges:
        # extract the data from doctors whose opening hours fall in the given
        # hour bracket
        interval_data = spec_data[ \
            (spec_data['total_weekly_hours'] >= iv_start) & \
            (spec_data['total_weekly_hours'] < iv_start + hour_range_interval)]
        
        # do the same calculation for both "unique" and "total" patient counts
        for pat_type in ['total', 'unique']:
            quantile = interval_data['N_{}_quarterly_patients'\
                                        .format(pat_type)].quantile(threshold)

            # calculate the median capacity of the upper (1-threshold) % of docs
            upper_median = interval_data[interval_data['N_{}_quarterly_patients'\
                    .format(pat_type)] >= quantile]['N_{}_quarterly_patients'\
                    .format(pat_type)].median()
            # find the indices of all doctors of a given speciality and hour
            # bracket that are not in the uper (1-threshold) % of capacities
            lower_index = interval_data[interval_data['N_{}_quarterly_patients'\
                                    .format(pat_type)] < quantile].index
            
            # if the number of doctors in the lower threshold % of capacities is
            # non-zero: set their capacities to the maximum capacity
            if len(lower_index) > 0:
                opening_hours.loc[lower_index, 'capacity_{}_hour-based'\
                                .format(pat_type)] = upper_median

            
opening_hours['capacity_total_hour-based'] = \
        opening_hours['capacity_total_hour-based'].round().astype(int)
opening_hours['capacity_unique_hour-based'] = \
        opening_hours['capacity_unique_hour-based'].round().astype(int)
# visulaize capacities
fig, axes = plt.subplots(5, 3, figsize=(20, 15))
for spec, ax in zip(specs, axes.flatten()):
    spec_data = opening_hours[opening_hours['specialization'] == spec]
    median = spec_data['capacity_total_hour-based'].median()
    ax.hist([spec_data['capacity_total_hour-based'], spec_data['N_total_quarterly_patients'],\
             spec_data['capacity_total_hour-based'] - spec_data['N_total_quarterly_patients']],
             rwidth = 0.8, bins=range(int(spec_data['capacity_total_hour-based'].min()/1000)*1000, 10000, 500),
             label=['total capacity', 'quarterly patients', 'free capacity'])
    ax.plot([median, median], ax.get_ylim(), label='median capacity')
    ax.set_title(spec)
    ax.legend()
    
axes[4, 1].remove()
axes[4, 2].remove()
fig.tight_layout()<jupyter_output><empty_output><jupyter_text>Model 3: maximum capacity For doctors that are not yet in the top 25 percentile of patient capacity, assign the median value of the top 25 percentile as capacity (regardless of opening hours).<jupyter_code>def get_upper_median(data, spec, pat_type):
    spec_data = data[data['specialization'] == spec]
    quant_75 = spec_data['N_{}_quarterly_patients'.format(pat_type)].quantile(0.75)
    quant_median = spec_data[spec_data['N_{}_quarterly_patients'.format(pat_type)] >=\
                             quant_75]['N_{}_quarterly_patients'.format(pat_type)].median()
    return quant_median, quant_75
opening_hours['capacity_total_max'] = opening_hours['N_total_quarterly_patients']
opening_hours['capacity_unique_max'] = opening_hours['N_unique_quarterly_patients']

for spec in specs:
    for pat_type in ['total', 'unique']:
        upper_median, quant_75 = get_upper_median(opening_hours, spec, pat_type)
        print('[{}, {}]: upper median {}, 75th percentile {}'.format(spec, pat_type, upper_median, quant_75))
        lower_index = opening_hours[(opening_hours['specialization'] == spec) & \
                                    (opening_hours['N_{}_quarterly_patients'.format(pat_type)] < quant_75)].index
        opening_hours.loc[lower_index, 'capacity_{}_max'.format(pat_type)] = upper_median
        
opening_hours['capacity_total_max'] = opening_hours['capacity_total_max'].round().astype(int)
opening_hours['capacity_unique_max'] = opening_hours['capacity_unique_max'].round().astype(int)
fig, axes = plt.subplots(5, 3, figsize=(20, 15))
for spec, ax in zip(specs, axes.flatten()):
    spec_data = opening_hours[opening_hours['specialization'] == spec]
    median = spec_data['capacity_total_max'].median()
    ax.hist([spec_data['capacity_total_max'], spec_data['N_total_quarterly_patients'],
             spec_data['capacity_total_max'] - spec_data['N_total_quarterly_patients']],
             rwidth = 0.8, bins=range(int(spec_data['capacity_total_max'].min()/1000)*1000, 10000, 500),
             label=['total capacity', 'quarterly patients', 'free capacity'])
    ax.plot([median, median], ax.get_ylim(), label='median')
    ax.set_title(spec)
    ax.legend()
    
axes[4, 1].remove()
axes[4, 2].remove()
fig.tight_layout()<jupyter_output><empty_output><jupyter_text>Save the data<jupyter_code>today = datetime.now().date()
opening_hours['adj_index'] = opening_hours.index
opening_hours.to_csv('data/matched_and_imputed_doctors_with_capacity_threshold{}.csv'.format(threshold),
                     index=False)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11131364685867055, 'fraction_numerical': 0.01300279499331632, 'mean_word_length': 4.1373283395755305, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 26, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '544165', 'license': 'MIT', 'n_tokens_mistral': 2669, 'n_tokens_neox': 2383, 'n_words': 549, 'path': 'Code/calculate_capacity_threshold.ipynb', 'repo_name': 'MichiKal/healthcare-resilience', 'size': '8229'}