Day 6 – Reproducing Shoaib and Ramamohan 2022

Note

Reproduced in-text results 3, 4 and 6. Progress with 5 and 7. Total time used: 25h 36m (64.0%)

Consensus

I emailed Figure 4 and in-text results 1 and 2 to Tom and Alison, who both confirmed that they agreed these to be successfully reproduced.

9.22-9.34: In-text result 2

On reflection, I realised that I had not looked at the result for staff nurse utilisation. Add this, and it was consistent with the text. Hence, successfully reproduced at 9.34.

import sys
sys.path.append('../')
from timings import calculate_times

# Minutes used prior to today
used_to_date = 1262

# Times from today
times = [
    ('9.22', '9.34')
]

calculate_times(used_to_date, times)

Time spent today: 12m, or 0h 12m
Total used to date: 1274m, or 21h 14m
Time remaining: 1126m, or 18h 46m
Used 53.1% of 40 hours max

9.37-10.32: In-text result 3

This result varies whether and for how long the doctor intervenes in delivery

50% delivery cases: no doctor
30% delivery cases: doctor but only one-third of typical intervention
20% delivery cases: doctor and full intervention

In PHC.py, there are two classes: Delivery_with_doctor() and Delivery_no_doc(). Ordinarily, these are used in the class Delivery(), arrivals during OPD hours will get delivery with doctor, and arrivals outside of hours will get without:

if 0 < (self.registration_time - N * 1440) < 480:
    Delivery_with_doctor(urgent=True)      # sets priority
else:
    Delivery_no_doc(urgent=True)

Hence, I think to have 50% with no doctor, we add a sample from a uniform distribution of 0 to 1, and then either side of 0.5 they will get a doctor or not, where currently it just allows Delivery_with_doctor(). This is based on the assumption that this scenario is just referring to delivery cases during OPD hours, when a doctor could have otherwise intervened. I set this up with a doctor_delivery_scenario parameter in main() and then changes to Delivery().

# If True for scenario where reduce doctor intervention...
if doctor_delivery_scenario:
    # Sample intervention probability (between 0 and 1)
    int_prob = sim.Uniform(0, 1)
    # Assign 50% to receive doctor, and 50% to not
    if int_prob >= 0.5:
        Delivery_with_doctor(urgent=True)
    else:
        Delivery_no_doc(urgent=True)
# If scenario not in place (i.e. base case, normal function),
# then just assign all to receive doctor...
else:
    Delivery_with_doctor(urgent=True)

Whilst doing this, changed comments by inputs to main() into a doc string for clarity as I am often adjusting these.

For the cases where a doctor is provided, look to then have 60% receive partial intervention and 40% receive full intervention. In Delivery_with_doctor(), doctor time is sample twice (depending on whether or not the model is in the warm-up period or not). Both as: doc_time = round(sim.Uniform(30, 60, 'minutes').sample(),2). This contributes to doctor utilisation as follows:

Patient.doc_service_time.append(doc_time)
doc_tot_time = sum(Patient.doc_service_time)
z = doc_tot_time/(420*days*doc_cap)
doc_occupancy.append(z) - and then doc_occupancy is used to calculate utilisation in the replication results spreadsheet

As such, I introduced a probability of receiving shorter or longer times around this Uniform sampling. I assumed 1/3 normal intervention to mean that the sampled time is simply divided by three.

# If we are in a scenario where we are reducing the time spent
# by doctors during delivery to one third of normal time, for 60%
# of delivery patients...
if doctor_delivery_scenario:
    # Sample to get probability of short intervention
    short_int_prop = sim.Uniform(0, 1)
    if short_int_prop >= 0.4:
        doc_time = doc_time / 3

However, this came back wrong, with very little change in utilisation:

import pandas as pd

pd.read_csv('intext3_wrong.csv')

	Output	Normal delivery	Reduced doctor intervention	Change
0	Doctor utilisation	1.142	1.131	-0.011
1	Staff nurse utilisation	0.323	0.321	-0.002

Reflections

As before, time spent working out how to implement scenarios from the article

10.41-11.16, 11.28-11.47: Problem-solving in-text result 3

Noticed I hadn’t add global doctor_delivery_scenario to Delivery_with_doctor, and that I hadn’t add the scenario to the warm-up section either. However, that seemed to have no impact -

import pandas as pd

pd.read_csv('intext3_wrong2.csv')

	Output	Normal delivery	Reduced doctor intervention	Change
0	Doctor utilisation	1.140	1.129	-0.011
1	Staff nurse utilisation	0.321	0.323	0.002

I tried running it with print statements after the 50% chance of doctor assignment, which confirmed that patients were being split between getting a doctor and not. I also ran with print statement around shortening of doctor time, which likewise confirmed that around 60% were having their times shortened.

I looked across all the outputs and could spot only minor differences between the normal and reduced delivery results.

I then tried setting the implementation up in a different way - instead of assigning people to the Delivery_no_doc() class, setting the intervention time to 0 in the Delivery_with_doctor() class. However, this made no difference:

pd.read_csv('intext3_wrong3.csv')

	Output	Normal delivery	Reduced doctor intervention	Change
0	Doctor utilisation	1.144	1.129	-0.015
1	Staff nurse utilisation	0.326	0.320	-0.006

On reflection, I wondered if perhaps I had misinterpreted this scenario - and that actually, this scenario was on top of the change in admin, and that the language was ambiguous (that actually, it could be describing doing this seperately or in addition). This then resolved the issue, with results matching up to the paper. Hence, feel it is successfully reproduced at 11.47.

Reflections

Potential misinterpretation of scenarios from a text description which can be ambiguous.

# Minutes used prior to this item
used_to_date = 1274

# Times from today
times = [
    ('9.37', '10.32'),
    ('10.41', '11.16'),
    ('11.28', '11.47')
]

calculate_times(used_to_date, times)

Time spent today: 109m, or 1h 49m
Total used to date: 1383m, or 23h 3m
Time remaining: 1017m, or 16h 57m
Used 57.6% of 40 hours max

11.55-12.03, 12.06-12.11: Reproducing in-text result 4

It is again unclear whether the intervention in this case - adding a doctor - is in addition to or seperate to the prior interventions. As such, I will try running both options (just on top of base case, or in addition to admin and delivery changes).

Either result matches up with the description, so present both as possible answers, but consder to be succesfully reproduced at 12.11

# Minutes used prior to this item
used_to_date = 1383

# Times from today
times = [
    ('11.55', '12.03'),
    ('12.06', '12.11')
]

calculate_times(used_to_date, times)

Time spent today: 13m, or 0h 13m
Total used to date: 1396m, or 23h 16m
Time remaining: 1004m, or 16h 44m
Used 58.2% of 40 hours max

Times below: Start on in-text result 5 (and amend result 3’s table)

Times:

12.14-12.18
12.21-12.37
12.41-12.56
13.38-13.53
15.10-15.23

Re-ran with less delivery as an option alone for in-text result 3 (as this was missing from the table but would provide helpful context). Whilst that ran, started on in-text result 5.

This result appears to build on configuration 1 but with:

High demand conditions (2 inpatient and childbirth cases per day - as in Figure 3A-D with (2/2/2))
Varying the number of inpatient beds (6 to 4)

However, changing the bed number with inpatient_bed_n had no impact on results, with utilisation at 18% for 4 and 6 beds.

This was based on ipd bed occ - although I noticed that ipd occ had results of 0.196 and then 0.295, which - if that is referring to utilisation - is much closer to the paper result of 33%. I explored how each of these were calculated:

ipd bed occ:

This was what I used for inpatient bed utilisation in Table 6 and Figure 3C
ipd bed occ is equal to bed_util
bed_util is based on bed_time - bed_util.append(bed_time/(days*1440*6))
bed_time is used in Delivery_no_doc(), Delivery_with_doctor(), IPD_no_doc(), and IPD_with_doc(). In each case, it is bed_time += y, where y is sampling from a distribution (e.g. y = sim.Uniform(240, 1440, 'minutes').sample(), y = sim.Triangular(60, 360, 180, 'minutes').bounded_sample(0))

ipd occ:

Created from bed_occupancy_list
That is based on bed.occupancy: bed_occupancy_list.append(round(bed.occupancy.mean(), 2))
bed is made from sim.Resource("Bed", capacity=inpatient_bed_n). It is used in Delivery_no_doc(), Delivery_with_doctor(), IPD_no_doc(), and IPD_with_doc()
In Salabim, a Resource occupancy is described as “occupancy (=claimed_quantity / capacity)”, which fits with being related to utilisation.

Hence, both methods should be providing utilisation of the beds (since both are based on time the beds are in use over the total possible time) - yet they have differing results - particularly when set to 4 beds.

pd.read_csv('intext5_wrong.csv', index_col=0)

	in5_6bed	in5_4bed
ipd bed occ	0.185185	0.185522
ipd occ	0.196000	0.295000

ipd occ is based on the bed SimPy resource, which is what we altered to change the number of beds. Hence, it seems the bed capacity has been altered for that, but not altered somewhere for ipd bed occ.

I looked back at the calculation of ipd bed occ and noticed that bed utilisation has 6 beds hard coded into it - bed_util.append(bed_time/(days*1440*6)). When I amended the script to adjust the bed capacity via a parameter input to main(), I had not noticed this, and had only changed it in the simpy Resource. I amended this accordingly and re-ran. Unfortunately, this did not fix it.

pd.read_csv('intext5_wrong2.csv', index_col=0)

	in5_6bed	in5_4bed
ipd bed occ	1.122112	1.115264
ipd occ	0.198000	0.295000

I then searched for all occurrences of the number 6, but identified nothing further that related to the number of beds.

Struggling to identify a solution, I moved onto the next item in the scope, and plan to return to this later.

Reflections

Confusion from duplicate outputs.

Hard-coding of parameters (which can then be tricky to spot afterwards to change).

15.24-15.55: Reproducing in-text result 6

This result is a variant of Figure 2B IAT 3, where we assign administrative work to the staff nurse to decrease the utilisation of the NCD nurse.

In the normal model, the admin work is assigned to the doctor and the NCD nurse. We previously add a various to move the work from the doctor to the staff nurse. Similar to last time, we add another variant, this time moving work from the NCD nurse to the staff nurse. I amended the name of the boolean controlling it for the in-text results 3/4 (for clarity, distinguishing between that scenario and this one), and re-ran that notebook.

The resultant NCD utilisation was 100%, as in the paper, so consider to be successfully reproduced at 15.55.

# Minutes used prior to this item
used_to_date = 1396

# Times from today
times = [
    ('12.14', '12.18'),
    ('12.21', '12.37'),
    ('12.41', '12.56'),
    ('13.38', '13.53'),
    ('15.10', '15.23'),
    ('15.24', '15.55')
]

calculate_times(used_to_date, times)

Time spent today: 94m, or 1h 34m
Total used to date: 1490m, or 24h 50m
Time remaining: 910m, or 15h 10m
Used 62.1% of 40 hours max

15.58-16.44: Starting on in-text result 7

This scenario appears to build on-top of the change in admin work (“in addition to the administrative work”). In it, 10% of the NCD cases are assigned to the staff nurse instead of the NCD nurse.

In PHC.py, I can see that the NCD nurses do admin and do OPD appointments. The appointments are in the class Patient(), and if OPD_PatientGenerator.OPD_List[OPD_PatientGenerator.patient_count]["Age"] >= 30, then it requests an NCD nurse (yield self.request((NCD_Nurse, 1))). Their time spent is then saved to:

Patient.NCD_Nurse_time_list.append(x)
Patient.NCD_Nusre_1_time += x
ncd_time += x

This code is repeated before and after warm-up (with the difference being in whether results are saved). As such, to implement this scenario, it appears I need to change this section to instead:

Request a staff nurse - from other sections, can see this is yield self.request(staff_nurse)
Save the time to the staff nurse time - from other sections, can see this is:

Main.NT_list.append(temp)
Main.SN_time += z
Main.staff_nurse_IPD += temp

The NT_list is what I modified for in-text result 2.
SN_time appears to be legacy code, as it is not used in the results section.
staff_nurse_IPD is used to store staff nurse time specifically spent on IPD (as opposed to delivery or ANC). It is just used in the large results spreadsheet. However, this is specific to IPD appointments rather than OPD.

I sampled from a uniform distribution to control assignment of 10% of cases to the staff nurse.

Upon running the model, I found these modifications reduced the NCD nurse utilisation markedly (99% to 84%), although not as much as in the paper (100% to 71%).

pd.read_csv('intext7_close.csv')

	Output	Normal	Admin from NCD to staff nurse	Admin and 10% OPD appointments from NCD to staff nurse
0	NCD nurse utilisation	1.236	0.982	0.835

I tried adding a scenario with just appointments (and not the admin change), just in case this could relate to misinterpretation of the next (although wouldn’t expect it to go lower, so wasn’t optimistic). As expected, this was not the solution, and the results were higher.

pd.read_csv('intext7_wrong.csv')

	Output	Normal	Admin from NCD to staff nurse	10% OPD appointments from NCD to staff nurse	Admin and 10% OPD appointments from NCD to staff nurse
0	NCD nurse utilisation	1.237	0.989	1.071	0.834

Timings

# Minutes used prior to today
used_to_date = 1262

# Times from today
times = [
    ('9.22', '9.34'),
    ('9.37', '10.32'),
    ('10.41', '11.16'),
    ('11.28', '11.47'),
    ('11.55', '12.03'),
    ('12.06', '12.11'),
    ('12.14', '12.18'),
    ('12.21', '12.37'),
    ('12.41', '12.56'),
    ('13.38', '13.53'),
    ('15.10', '15.23'),
    ('15.24', '15.55'),
    ('15.58', '16.44')
]

calculate_times(used_to_date, times)

Time spent today: 274m, or 4h 34m
Total used to date: 1536m, or 25h 36m
Time remaining: 864m, or 14h 24m
Used 64.0% of 40 hours max