Reproducing Figure 2A-D

This notebook aims to reproduce Figures 2A-D from Shoaib M, Ramamohan V. Simulation modeling and analysis of primary health center operations. SIMULATION 98(3):183-208. (2022). https://doi.org/10.1177/00375497211030931.

Figure 2A

Figure 2B

Figure 2C

Figure 2D

Parameters

In these figures, we vary:

• Number of outpatients per day:
  • 170 (same as config 4)
  • 85
  • 65
• Average service time for outpatients - mean (SD):
  • 0.87 (0.21) (same as config 1)
  • 2.5 (0.5)
  • 5 (1)

To calculate inter-arrival times from those numbers per day, based on article description and the provided patient counts and equivalent IAT, understand the method for calculation to be round(60/(n/8.5)), where n is the number of arrivals per day. As such…

# Calculation of inter-arrival times
print(f'For 170 outpatients, use IAT (rounded to nearest int): {60/(170/8.5)}')
print(f'For 85 outpatients, use IAT (rounded to nearest int): {60/(85/8.5)}')
print(f'For 65 outpatients, use IAT (rounded to nearest int): {60/(65/8.5)}')

For 170 outpatients, use IAT (rounded to nearest int): 3.0
For 85 outpatients, use IAT (rounded to nearest int): 6.0
For 65 outpatients, use IAT (rounded to nearest int): 7.846153846153846

Set up

# To run model
import PHC

# To import results and produce figures
from reproduction_helpers import process_results
import pandas as pd
import os
import matplotlib.pyplot as plt
import numpy as np

# To speed up run time
from multiprocessing import Pool

# Additional package to record runtime of this notebook
import time
start = time.time()

# Paths to save image files to
output_folder = '../outputs'
fig2a_path = os.path.join(output_folder, 'fig2a.png')
fig2b_path = os.path.join(output_folder, 'fig2b.png')
fig2c_path = os.path.join(output_folder, 'fig2c.png')
fig2d_path = os.path.join(output_folder, 'fig2d.png')

Run model

As this is a variation on configuration 1 (which is the default parameters in PHC.py), we just need to input the varying number of outpatients and service time.

# Varying number of outpatients
arr_dict = [
    {
        'OPD_iat': 3,
        'rep_file': 'arr170'
    },
    {
        'OPD_iat': 6,
        'rep_file': 'arr85',
    },
    {
        'OPD_iat': 8,
        'rep_file': 'arr65',
    }
]

# Varying service time
serv_dict = [
    {
        'mean': 0.87,
        'sd': 0.21,
        'consult_boundary_1': 0.5,  # From PHC.py
        'consult_boundary_2': 0.3,  # From PHC.py
        'rep_file': 'serv087'
    },
    {
        'mean': 2.5,
        'sd': 0.5,
        'consult_boundary_1': 1,  # Guess
        'consult_boundary_2': 1,  # Guess
        'rep_file': 'serv25'
    },
    {
        'mean': 5,
        'sd': 1,
        'consult_boundary_1': 2,  # From config 4
        'consult_boundary_2': 2,  # From config 4
        'rep_file': 'serv5'
    }
]

Create each combination for the reproduction

dict_list = []
for arr in arr_dict:
    for serv in serv_dict:
        # Combine the dictionaries
        comb = {**arr, **serv}
        # Replace the file name
        comb['rep_file'] = f'''f2_{arr['rep_file']}_{serv['rep_file']}.xls'''
        # Save to list
        dict_list.append(comb)

len(dict_list)

9

# Append 's_' to all items
for i, d in enumerate(dict_list):
    dict_list[i] = {f's_{k}': v for k, v in d.items()}

# Preview example
dict_list[0]

{'s_OPD_iat': 3,
 's_rep_file': 'f2_arr170_serv087.xls',
 's_mean': 0.87,
 's_sd': 0.21,
 's_consult_boundary_1': 0.5,
 's_consult_boundary_2': 0.3}

Run the model (with parallel processing to reduce run time)

# Wrapper function to allow input of dictionary with pool
def wrapper(d):
    return PHC.main(**d)

# Create a process pool that uses all CPUs
with Pool() as pool:
    # Run PHC.main() using each of inputs from config
    pool.map(wrapper, dict_list)

 No of replications done 0
 No of replications done 0
 No of replications done 0
 No of replications done 0
 No of replications done 0
 No of replications done 0
 No of replications done 0
 No of replications done 0
 No of replications done 0
 No of replications done 1
 No of replications done 1
 No of replications done 1
 No of replications done 1
 No of replications done 1
 No of replications done 2
 No of replications done 1
 No of replications done 2
 No of replications done 2
 No of replications done 3
 No of replications done 2
 No of replications done 1
 No of replications done 1
 No of replications done 2
 No of replications done 2
 No of replications done 1
 No of replications done 3
 No of replications done 4
 No of replications done 3
 No of replications done 3
 No of replications done 3
 No of replications done 5
 No of replications done 4
 No of replications done 3
 No of replications done 4
 No of replications done 2
 No of replications done 2
 No of replications done 4
 No of replications done 4
 No of replications done 6
 No of replications done 5
 No of replications done 5
 No of replications done 5
 No of replications done 2
 No of replications done 4
 No of replications done 7
 No of replications done 5
 No of replications done 6
 No of replications done 3
 No of replications done 6
 No of replications done 6
 No of replications done 8
 No of replications done 6
 No of replications done 5
 No of replications done 7
 No of replications done 7
 No of replications done 3
 No of replications done 9
 No of replications done 6
 No of replications done 7
 No of replications done 4
 No of replications done 8
 No of replications done 7
 No of replications done 8
 No of replications done 3
 No of replications done 7
 No of replications done 8
 No of replications done 9
 No of replications done 9
 No of replications done 4
 No of replications done 8
 No of replications done 5
 No of replications done 8
 No of replications done 9
 No of replications done 4
 No of replications done 9
 No of replications done 9
 No of replications done 5
 No of replications done 6
 No of replications done 5
 No of replications done 6
 No of replications done 7
 No of replications done 6
 No of replications done 7
 No of replications done 8
 No of replications done 7
 No of replications done 8
 No of replications done 9
 No of replications done 8
 No of replications done 9
 No of replications done 9

Create Figure 2A

# Import and process results
data_full = process_results([
    'f2_arr170_serv087', 'f2_arr85_serv087', 'f2_arr65_serv087',
    'f2_arr170_serv25', 'f2_arr85_serv25', 'f2_arr65_serv25',
    'f2_arr170_serv5', 'f2_arr85_serv5', 'f2_arr65_serv5'])

# Filter to doctor utilisation
data_2a = data_full.loc['doc occ']
data_2a

f2_arr170_serv087    0.311815
f2_arr85_serv087     0.225192
f2_arr65_serv087     0.204317
f2_arr170_serv25     0.640142
f2_arr85_serv25      0.389486
f2_arr65_serv25      0.327015
f2_arr170_serv5      1.141365
f2_arr85_serv5       0.640515
f2_arr65_serv5       0.517174
Name: doc occ, dtype: float64

def create_2a_2d(s, ylab, file, ylim=False):
    '''
    Creates Figure 2A or 2D (as both are very similar)

    Parameters:
    -----------
    s : pd.Series
        Series with mean result as values, and the model variant as index
    ylab : string
        Label for y axis
    file : string
        Path to save figure
    ylim : list, optional, default False
        If provided, gives the lower and upper limits for the Y axis

    Returns:
    --------
    matplotlib figure 
    '''
    # Reshape data so in appropriate format for plotting grouped bar chart
    names = [170, 85, 65]
    s5 = [s['f2_arr170_serv5'], s['f2_arr85_serv5'], s['f2_arr65_serv5']]
    s25 = [s['f2_arr170_serv25'], s['f2_arr85_serv25'], s['f2_arr65_serv25']]
    s87 = [s['f2_arr170_serv087'], s['f2_arr85_serv087'], s['f2_arr65_serv087']]

    data = pd.DataFrame(
        {'5 (1)': s5, '2.5 (0.5)': s25, '0.87 (0.21)': s87}, index=names)

    # Plot data
    ax = data.plot.bar(edgecolor='black', color='white', width=0.7)

    # Add patterns
    bars = ax.patches
    pattern = np.repeat(['..', '+++++', '\\\\\\\\'], 3)
    for bar, hatch in zip(bars, pattern):
        bar.set_hatch(hatch)
    ax.legend(title='Consultation time (min): mean (sd)')

    # Adjust figure
    plt.xlabel('Number of patients/day')
    plt.ylabel(ylab)
    if ylim:
        plt.ylim(ylim)
    plt.xticks(rotation=0)
    ax.grid(axis='y')
    ax.set_axisbelow(True)
    plt.savefig(file, bbox_inches='tight')
    plt.show()

create_2a_2d(data_2a, ylab='''Doctor's utilisation''', file=fig2a_path)

Create Figure 2B

Import and process data - note: used 0.87 since that is PHC1 and this is described as variants on that but, as paper notes and as I’ve observed (not shown here), “the doctor’s consultation time does not impact this outcome”, so it ultimately doesn’t matter which is chosen.

# Import and process results
data_2b = process_results([
    'f2_arr170_serv087', 'f2_arr85_serv087', 'f2_arr65_serv087'])

# Rename columns
data_2b.columns = (
    data_2b.columns.str.removeprefix('f2_arr').str.removesuffix('_serv087'))

# Preview data
data_2b.head()

	170	85	65
OPD patients	44022.000000	22156.400000	16744.600000
IPD patients	182.600000	184.000000	183.400000
ANC patients	357.700000	354.700000	368.400000
Del patients	359.300000	371.500000	373.900000
OPD Q wt	0.012598	0.006986	0.005684

Produce figure

ax = data_2b.loc['NCD occ'].plot.bar(
    edgecolor='black', color='white', hatch='..')
plt.xlabel('Number of patients/day')
plt.ylabel('''NCD nurse's utilisation''')
plt.ylim(0, 1.4)
plt.xticks(rotation=0)
ax.grid(axis='y')
ax.set_axisbelow(True)
plt.savefig(fig2b_path, bbox_inches='tight')
plt.show()

Create Figure 2C

Import and process data

# Import and process results
data_2c = process_results([
    'f2_arr170_serv5', 'f2_arr85_serv5', 'f2_arr65_serv5'])

# Rename columns
data_2c.columns = (
    data_2c.columns.str.removeprefix('f2_arr').str.removesuffix('_serv5'))

# Preview data
data_2c.head()

	170	85	65
OPD patients	43965.100000	22121.700000	16680.200000
IPD patients	178.600000	177.200000	185.500000
ANC patients	360.900000	356.900000	358.800000
Del patients	373.300000	NaN	370.900000
OPD Q wt	6.936889	0.540399	0.275889

Produce figure

# Create figure
ax = data_2c.loc['OPD Q wt'].plot.bar(
    edgecolor='black', color='white', hatch='..')
plt.xlabel('Number of patients/day')
plt.ylabel('''OPD average waiting time (minutes)''')
plt.xticks(rotation=0)
plt.ylim(0, 8)
ax.grid(axis='y')
ax.set_axisbelow(True)
plt.savefig(fig2c_path, bbox_inches='tight')
plt.show()

Create Figure 2D

# Filter to pharmacy waiting time
data_2d = data_full.loc['Pharmacy Q wt']
data_2d

f2_arr170_serv087    2.435524
f2_arr85_serv087     0.463927
f2_arr65_serv087     0.293510
f2_arr170_serv25     2.497280
f2_arr85_serv25      0.456076
f2_arr65_serv25      0.286186
f2_arr170_serv5      1.282328
f2_arr85_serv5       0.341688
f2_arr65_serv5       0.215585
Name: Pharmacy Q wt, dtype: float64

create_2a_2d(data_2d, ylab='Pharmacy average waiting time (minutes)',
             ylim=[0, 3], file=fig2d_path)

Run time

# Find run time in seconds
end = time.time()
runtime = round(end-start)

# Display converted to minutes and seconds
print(f'Notebook run time: {runtime//60}m {runtime%60}s')

Notebook run time: 3m 27s