Merge pull request #10 from pythonhealthdatascience/dev

Dev

Author: amyheather

README.md

@@ -3,39 +3,113 @@
 # Stroke capacity planning model: python DES RAP
 
 [![python](https://img.shields.io/badge/-Python_3.13.1-blue?logo=python&logoColor=white)](https://www.python.org/)
-![licence](https://img.shields.io/badge/Licence-MIT-green.svg?labelColor=gray)
+![licence](https://img.shields.io/badge/🛡️_Licence-MIT-green.svg?labelColor=gray)
 [![Tests](https://github.com/pythonhealthdatascience/stroke_rap_python/actions/workflows/tests.yaml/badge.svg)](https://github.com/pythonhealthdatascience/stroke_rap_python/actions/workflows/tests.yaml)
 [![Linting](https://github.com/pythonhealthdatascience/stroke_rap_python/actions/workflows/lint.yaml/badge.svg)](https://github.com/pythonhealthdatascience/stroke_rap_python/actions/workflows/lint.yaml)
-[![ORCID: Heather](https://img.shields.io/badge/ORCID_Amy_Heather-0000--0002--6596--3479-brightgreen)](https://orcid.org/0000-0002-6596-3479)
+[![ORCID](https://img.shields.io/badge/ORCID_Amy_Heather-0000--0002--6596--3479-A6CE39?&logo=orcid&logoColor=white)](https://orcid.org/0000-0002-6596-3479)
 
 </div>
 
 This repository applies the [Python DES RAP Template](https://github.com/pythonhealthdatascience/rap_template_python_des) to a real-life example:
 
 > Monks T, Worthington D, Allen M, Pitt M, Stein K, James MA. A modelling tool for capacity planning in acute and community stroke services. BMC Health Serv Res. 2016 Sep 29;16(1):530. doi: [10.1186/s12913-016-1789-4](https://doi.org/10.1186/s12913-016-1789-4). PMID: 27688152; PMCID: PMC5043535.
 
+Model diagram:
+
 ![](images/stroke_rehab_design.png)
 
 <br>
 
 ## Installation
 
-TBC
+Clone the repository locally:
+
+```
+git clone https://github.com/pythonhealthdatascience/stroke_rap_python.git
+cd stroke_rap_python
+```
+
+Use the provided `environment.yaml` file to set up a Python environment with `conda`:
 
-<!-- TODO: Provide instructions for installing dependencies and setting up the environment. -->
+```
+conda env create --file environment.yaml
+conda activate
+```
+
+The provided `environment.yaml` file is a snapshot of the environment used when creating the repository, including specific package versions. You can update this file if necessary, but be sure to test that everything continues to work as expected after any updates. Also note that some dependencies are not required for modelling, but instead served other purposes, like running `.ipynb` files and linting.
+
+As an alternative, a `requirements.txt` file is provided which can be used to set up the environment with `virtualenv`. This is used by GitHub actions, which run much faster with a virtual environment than a conda environment. However, we recommend locally installing the environment using conda, as it will also manage the Python version for you. If using `virtualenv`, it won't fetch a specific version of Python - so please note the version listed in `environment.yaml`.
 
 <br>
 
 ## How to run
 
-TBC
+The simulation code is provided as a **package** within `simulation/`. There are notebooks executing the model and analysing the results in `notebooks/`.
+
+To run the model with base parameters once or with replications:
+
+```
+from simulation.parameters import Param
+from simulation.runner import Runner
+
+param = Param()
+runner = Runner(param=param)
+
+single_result = runner.run_single(run=0)
+rep_results = runner.run_reps()
+```
+
+Example altering the model parameters:
+
+```
+from simulation.parameters import Param, ASUArrivals, RehabRouting
+from simulation.runner import Runner
+
+# Modified one of the arrival rates, some routing probabilities, and the
+# number of replications
+param = Param(
+    asu_arrivals=ASUArrivals(tia=10),
+    rehab_routing=RehabRouting(neuro_esd=0.2, neuro_other=0.8),
+    number_of_runs=10
+)
+runner = Runner(param=param)
+rep_results = runner.run_reps()
+```
+
+### Generating the results from the article
+
+The original study used Simul8. Each of the outputs from that article have been replicated in this repository using Python:
+
+* Figure 1. Simulation probability density function for occupancy of an acute stroke unit.
+* Figure 3. Simulated trade-off between the probability that a patient is delayed and the no. of acute beds available.
+* Table 2. Likelihood of delay. Current admissions versus 5% more admissions.
+* Table 3. Results of pooling of acute and rehab beds.
+* Supplementary Table 1. Likelihood of delay. Current admissions versus No Complex neurological patients.
+* Supplementary Table 3. Likelihood of delay. Current admissions versus ring fenced acute stroke beds.
+
+To generate these, simply execute `notebooks/analysis.ipynb`.
+
+#### Examples
+
+**Figure 1**
+
+Original:
+
+![](docs/article/fig1.png)
+
+From this repository:
+
+![](outputs/occupancy_freq_asu.png)
+
+**Figure 3**
+
+Original:
 
-<!-- Provide step-by-step instructions and examples.
+![](docs/article/fig3.png)
 
-Clearly indicate which files will create each figure in the paper. Hypothetical example:
+From this repository:
 
-* To generate **Figures 1 and 2**, execute `notebooks/base_case.ipynb`
-* To generate **Table 1** and **Figures 3 to 5**, execute `notebooks/scenario_analysis.ipynb` -->
+![](outputs/delay_prob_asu.png)
 
 <br>
 

docs/article_monks_etal_2016.pdf renamed to docs/article/article_monks_etal_2016.pdf

@@ -805,4 +805,52 @@ I then add new tests based on the python template and on the tests ran in https:
 
 I add methods to check parameter validity in `Param`. This flagged that rehab other routing probability don't sum to 100% (88% and 13%) - but this is as described in the paper, and presumed to be due to rounding, so altered the validation test to allow.
 
-> 💡 When explain tests, could do all in one section, like Tests > Back tests, Tests > Functional tests, Tests > Unit tests - and then on each of those pages, it's like, if you have parameter validation... if you have warm-up... etc. etc. suggesting tests could include.
+> 💡 When explain tests, could do all in one section, like Tests > Back tests, Tests > Functional tests, Tests > Unit tests - and then on each of those pages, it's like, if you have parameter validation... if you have warm-up... etc. etc. suggesting tests could include.
+
+## Scenario logic
+
+Having successfully implemented the base model generating Figure 1 and 3 (as in https://github.com/pythonhealthdatascience/llm_simpy/), I then moved on to the scenarios from Monks et al. 2016. These were:
+
+0. **Current admissions** Current admission levels; beds are reserved for either acute or rehab patients
+1. **5% more admissions** A 5% increase in admissions across all patient subgroups.
+2. **Pooling of acute and rehab beds** The acute and rehab wards are co-located at same site. Beds are pooled and can be used by either acute or rehabilitation patients. Pooling of the total bed stock of 22 is compared to the pooling of an increased bed stock of 26.
+3. **Partial pooling of acute and rehab beds** The acute and rehab wards are co-located at same site. A subset of the 26 beds are pooled and can be used by either acute or rehab patients.
+4. **No complex-neurological cases** Complex neurological patients are excluded from the pathway in order to assess their impact on bed requirements
+
+### Scenario 1 and 4
+
+As from the supplementary:
+
+> "Scenarios investigating increased demand multiply the mean arrival rates (supplied in main text) by the appropriate factor. To exclude a particular patient group, the mean inter-arrival time for that group is multiplied by a large number such that no arrivals will occur in the modelled time horizon."
+
+Hence, it is understood that:
+
+* **Scenario 1** can be achieved by multiplying all patient IAT by 1.05 (see below: 0.95).
+* **Scenario 4** can be achieved by multiplying IAT for complex neurological patients by a very high number (e.g. 10,000,000) - and can add a test which checks no patients are complex neurological.
+
+## Using multiple replications
+
+Altered `Runner` to output summary tables from across replications, and switched to using these for Figures 1 and 3.
+
+## Scenario 1 + Table 2
+
+Ran scenario 1 in `analysis.ipynb` and created Table 2. Noticed some differences. 
+
+I find scenario with same bed number, probability of delay drops. They find that it goes up.
+
+Thinking through the logic, scenario has more arrivals -> wards more full -> expect delays for lower max bed numbers earlier -> expect higher probability of delay for lower bed numbers.
+
+I then realised my mistake! I had actually lower admissions, as I'd multiplied IAT by 1.05, when I should've multiplied by 0.95.
+
+## Scenario 4 + supplementary table 1
+
+Ran scenario, adjusted function from table 2 so it could be used to make this table too.
+
+## Scenario logic
+
+We now have two remaining scenarios:
+
+* Scenario 2: **Pooling of acute and rehab beds** The acute and rehab wards are co-located at same site. Beds are pooled and can be used by either acute or rehabilitation patients. Pooling of the total bed stock of 22 is compared to the pooling of an increased bed stock of 26.
+* Scenario 3: **Partial pooling of acute and rehab beds** The acute and rehab wards are co-located at same site. A subset of the 26 beds are pooled and can be used by either acute or rehab patients.
+
+It took quite a while to understand the formula and how to implement them.

docs/article/fig1.png

@@ -0,0 +1,9 @@
+# Data dictionary for `parameters.csv`
+
+| Column | Data type | Description | Possible values |
+| - | - | - | - |
+| unit | str | Hospital unit | `asu`: Acute Stroke Unit<br>`rehab`: Rehabilitation Unit (post-acute recovery care) |
+| parameter | str | Type of operational metric | `iat`: Inter-arrival time (time between patient admissions)<br>`los`: Length of stay (duration from admission to discharge)<br>`routing`: Transition probability between care pathways  |
+| type | str | Patient classification or care transition path | **For `iat`/`los`:**<br>`stroke`: Stroke patients<br>`stroke_esd`: Stroke patients transferred to Early Supported Discharge<br>`stroke_no_esd`: Stroke patients not transferred to Early Supported Discharge<br>`tia`: Transient Ischemic Attack patients<br>`neuro`: Complex neurological patients<br>`other`: Other patient types<br><br>**For `routing`:**<br>`[diagnosis]_rehab`: Probability of transferring to rehabilitation unit<br>`[diagnosis]_esd`: Probability of Early Supported Discharge<br>`[diagnosis]_other`: Probability of other discharge pathways |
+| mean | float | **For `iat`:** Mean days between admissions<br>**For `los`:** Mean days in unit<br>**For `routing`:** Probability (0-1 scale) | - |
+| sd | float | Standard deviation of the mean | - |

docs/article/fig2.png

@@ -0,0 +1,38 @@
+unit,parameter,type,mean,sd
+asu,iat,stroke,1.2,NA
+asu,iat,tia,9.3,NA
+asu,iat,neuro,3.6,NA
+asu,iat,other,3.2,NA
+rehab,iat,stroke,21.8,NA
+rehab,iat,neuro,31.7,NA
+rehab,iat,other,28.6,NA
+asu,los,stroke_no_esd,7.4,8.61
+asu,los,stroke_esd,4.6,4.8
+asu,los,tia,1.8,2.3
+asu,los,neuro,4,5
+asu,los,other,3.8,5.2
+rehab,los,stroke_no_esd,28.4,27.2
+rehab,los,stroke_esd,30.3,2un3.1
+rehab,los,tia,18.7,23.5
+rehab,los,neuro,27.6,28.4
+rehab,los,other,16.1,14.1
+asu,routing,stroke_rehab,0.24,NA
+asu,routing,stroke_esd,0.13,NA
+asu,routing,stroke_other,0.63,NA
+asu,routing,tia_rehab,0.01,NA
+asu,routing,tia_esd,0.01,NA
+asu,routing,tia_other,0.98,NA
+asu,routing,neuro_rehab,0.11,NA
+asu,routing,neuro_esd,0.05,NA
+asu,routing,neuro_other,0.84,NA
+asu,routing,other_rehab,0.05,NA
+asu,routing,other_esd,0.1,NA
+asu,routing,other_other,0.85,NA
+rehab,routing,stroke_esd,0.4,NA
+rehab,routing,stroke_other,0.6,NA
+rehab,routing,tia_esd,0,NA
+rehab,routing,tia_other,1,NA
+rehab,routing,neuro_esd,0.09,NA
+rehab,routing,neuro_other,0.91,NA
+rehab,routing,other_esd,0.13,NA
+rehab,routing,other_other,0.88,NA