As the primary function of the model is to determine the number of children in each health status at a given point in time, the model’s core consists of a basic Markov model. This category of model has been used for various types of probabilistic modeling over the past 20 years, including epidemiological modeling, especially when the model has a defined number of outcomes, or states.11 The Markov model used for the treatment coverage model for Yemen was composed of four different health statuses (Healthy, Sick, In-treatment, and Deceased) for children to move between based on a set of probabilistic transition rates as depicted in Figure 1. This means that by just knowing these probabilities the model can produce a breakdown of the number of children who are estimated to become sick or seek and receive treatment. The model rests on the following assumptions:
- Any sick person who decides to seek treatment will contact a healthcare provider in community, in a mobile team or at a health facility, and will be treated;
- To reach the In-treatment health status one must have previously resided in the Sick health status; and
- The Deceased health status is an absorbing state, meaning once someone enters it, they cannot transition out to another health status again.
These assumptions mean that children will not be turned away for treatment, they must be sick before they can be treated, and if a child dies, they remain deceased for the duration of the simulation. Furthermore, as all children that reach care are assumed to be treated, this means that the coverage rate is really the percent of sick children who are able to reach points of care, as all will be treated.
Table1: Transition Rate Probability Table
Health Status Transition
|
Transition Probability
|
Healthy to Healthy (pHH)
|
.7601
|
Healthy to Sick (pHS)
|
.2388
|
Healthy to In-treatment (pHT)
|
0
|
Healthy to Deceased (pHD)
|
.0011
|
Sick to Healthy (pSH)
|
.6968
|
Sick to Sick (pSS)
|
(1-(Movement function+6.2712)/9)
|
Sick to In-treatment (pST)
|
(Movement function-.01296834)/9
|
Sick to Deceased (pSD)
|
.001441
|
In-treatment to Healthy (pTH)
|
.9902
|
In-treatment to Sick (pTS)
|
.00961
|
In-treatment to Treated (pTT)
|
0
|
In-treatment to Deceased (pTD)
|
.00019
|
Deceased to Healthy (pDH)
|
0
|
Deceased to Sick (pDS)
|
0
|
Deceased to In-treatment (pDT)
|
0
|
Deceased to Deceased (pDD)
|
1
|
List of the probabilities of a child transitioning from one state of health to another in the model.
Determining the base transition rates
The health facility data collected through UNICEF programs in-country provided insight into the number of cases of diarrheal disease among children under five years old that were seen and treated (pST), the mortality rate among children treated with Zn/ORS solutions, the efficacy of Zn/ORS solutions in treating the disease and the percentage of the population the monitored facilities are estimated to have served. As the UNICEF data can only track consultations, our model has been designed to also allow for potential double counting by tracking the number of cases seen and not the number of individual children who fall sick. This information was filtered to examine only diarrheal diseases and to eliminate any incomplete reports where the total number of reported children treated for diarrheal diseases did not equal the sum of all the individual children reported to have diarrheal diseases. Only completed data was used to not skew the value of the resulting transition probabilities and to match the assumption that all who reach a health facility will be treated. Coupled with the most recent overall mortality rates in the literature from the WHO and an estimated incidence rate for Yemen from El Bcheraoui et. al, the transition rates were calculated to allow the model’s output to match the previously mentioned program data.2,12 The final transition probabilities were as listed in Table 1. An overview of the terms used, as well as their values and meaning are provided in Table 2. These are described in more detail later in the manuscript.
Relationship between incidence rate, Rotavirus and Healthy to Sick transition rates
An incidence rate of seven episodes of diarrheal disease among children under five per person-year estimated by El Bcheraoui et. al was used as the initial assumed metric for the healthy to sick transition rates to match.2 This was then modified to account for general sources of diarrheal diseases and in particular Rotavirus, the most common source which is noted as the cause of between 35-60 % of enteric diseases in children under five years old.13 Rotavirus infections confer natural immunity, protecting against 87% of subsequent severe diarrhea cases which increases with each subsequent infection.13 Thus, for n rotavirus infection after the first infection from rotavirus the likelihood that a child is not protected against severe diarrhea can be roughly estimated at a probability of 0.13n.13 With these factors implemented the general transition rate from healthy to sick with Rotavirus or non-Rotavirus diarrheal disease was adjusted to attempt to match previously estimated incidence rates as discussed further in the validation section.
Variables and Parameters used throughout the model: weather, conflict and their combined effects on transportation infrastructure
The health information system in Yemen was weak even prior to the current crisis, with a reliance on paper-based reporting and challenges with completeness and timeliness of data received at governorate and national levels. With the onset of the conflict and humanitarian crisis which has now brought the health system to the brink of collapse, these gaps in the available health system data are even more notable. This issue affects the data on childhood diarrheal disease as with all other routine administrative data tracked through the health information system. To mitigate this concern about quality of routine administrative data on childhood diarrheal diseases, other types of data outside of traditional health data were incorporated, including the occurrence of conflict events, weather patterns, and the effects of both of these on the traversable nature of transportation infrastructure (roads and bridges). The factors influence pST, the probability determining the likelihood that someone will be willing and able to seek out treatment and thus move from the Sick health status to the In-treatment health status, and conversely, the probability that a sick person will remain sick via pss since a larger proportion of sick people seeking treatment should result in a smaller proportion not seeking treatment. The extent of conflict events (represented in the model as the state of the infrastructure) and the estimated weather conditions for that time period allows for the creation of an equation incorporating the time varying nature of these data types into the probabilistic decision process of the model. Thus, the model is not only probabilistic but also time varying.
The state of the weather and the quality of the roads and few bridges in Yemen, which are already affected by and dependent on the state of the conflict and so serve as a proxy for the severity of the conflict, have been shown to be related functions. These factors were selected due to the importance that climate change (in the form of precipitation and storms) and conflict (in the form of destroyed transportation infrastructure and unwillingness or inability to cross lines of conflict) have to the average Yemeni. The weather component is assumed to be a sinusoidal function with a period of one year. This weather function was created as a normalized function of best fit from monthly precipitation data for Yemen from 1901 to 2016 published by the World Bank.9 Clear weather will not alter movement but as rainfall increases travel becomes more difficult and therefore less likely. For our model this means the range of unweighted probabilities could vary between 0, which in the context of the model would indicate a perfect storm during the heaviest rainfall of the year preventing anyone from traveling or reaching their destination, to 1, which would result in anyone being able to access roads leading to health services during the driest week of the year. A constant threat of airstrikes and fighting will further reduce this function. Based on data published by the Yemen Data Project Organization indicating that there are 21.6 airstrikes for every 1000 people, it is further multiplied by a set constant of 97.84%.14 This probability represents an individual’s likelihood of not being in an airstrike as there will be .0216 airstrikes per person. If viewed as a percent, then on average a person has a 2.16% chance of being affected or a 97.84% chance of not being affected. Like the rains, this can deter people from traveling and further destroy infrastructure, albeit much more directly.
It is estimated that Yemen has about 50,000 km of roads.15 As the country does not have a functioning rail system to aid in transportation, these roads are the primary means of transportation for people within the country.15,16 Despite the importance of these roads, reports from the World Bank estimated that only 28% were all-weather paved before the conflict’s start, with only 11% of rural roads being paved.15 These nonpaved road segments have been reported by USAID to be damaged both by airstrikes and seasonal rains which has led to flooding, hindering travel during the rainy season.17 The model was constructed on the assumption that the rains predominantly affect the coastal regions to the south and west of the country and so only 50% of the unpaved roads potentially experience wash out from the worst of the rains. The bridges have also been affected by both the weather and the war: air strikes and fighting have destroyed some bridges, while many others have their access restricted, according to maps published by the World Food Program and the Yemen Logistics Cluster.10 These restricted bridges were noted to have detours that “may not be accessible during the rainy season”.10 Thus, when the rains begin the weather function discussed previously will start to decrease. When the weather function drops below a tunable threshold called “bridgecuttoff” in the model, a percentage of the roads corresponding to passable bridge detours that were previously opened will close. This mimics the idea that once the rain intensity increases to this specific threshold, routes that were previously accessible due to dry season detours around destroyed bridges will no longer be traversable, decreasing the number of open routes. The percentage of bridges that would close was designed to scale with the weather and vary from its maximum and minimum values due to the conflict, providing an alternative to “all or nothing” values for the bridges.
Different sections of the same continuous road may have different conditions, with some sections of the same road being in better or worse conditions than others. To improve specificity, the roads shown on the maps published by the Logistic Cluster and World Food Program were subdivided into individual road segments. A road segment was defined as a unique stretch of road that connects between marked communities or another road segment. Each segment could be either open (traversable) or closed (non-traversable). The percentage of open road segments is the output of the Road function. This function is composed of the first five numbered functions below.
Equation 1 determines the maximum value the road function can produce given the conflict and weather. The maximum value will occur when the weather is clear and there are no restrictions due to bridges or the associated detours being closed from rain. However, while its maximum value occurred when there were no weather based restrictions on travel, the conflict is still ongoing. As this is a significant factor in the state of travel infrastructure, the max road value was ultimately determined from the extent of the conflict. This Roadmax value serves as a proxy for the severity of the fighting as these values were based on how extensive the conflict was during the period examined; the weather provided more regular variation determining how useable these open road segments were. (see Equation 1 in the Supplementary Files)
With the base road segment conditions from the fighting set in equation 1, the minimum road segment value was calculated according to equation 2. This value was the percent of road segments left open from the conflict (Roadmax) being directly degraded further by weather conditions, specifically heavy, seasonal rainfall. This was calculated as the percent of road segments left functional after the assumed max wash-out of the previously open road segments, minus the percent of road segments that have a dry season bridge detour as the heaviest rains would make these detours unusable. (see Equation 2 in the Supplementary Files)
As previously noted, the state of all the bridges are not simply purely opened or purely closed. Based on the severity of weather, the percentage of these bridges that are inaccessible will change. As weather worsens beyond a point, more bridges and detours will close and vice versa. This is only true for severe weather events. Bridges will retain their normal degree of functionality up to a certain degree of weather severity, the tunable “bridgecutoff” value discussed previously. This leads the function depicting the number of bridges down to have different behaviors depending on weather conditions as shown in Equation 3. (see Supplementary Files)
With these components all determined, the road function could be fully built and is described in Equation 4. The minimum function (min) in Equation 4 ensures the Road function never produces a value beyond its maximum allowed value as determined by Equation 1. All previous equations are contained within this function including the best fit weather function Y(t). When combined with a scaling factor this function allows for user input, allowing it to be more finely tuned beyond the national average if needed. (see Equation 4 in the Supplementary Files)
While the roads are a major factor in a person’s willingness to travel to seek treatment, the weather too may play a role in this decision. As such, the overall transition probability movement function YT(t) was then created as a weighted average of the weather function, Y(t), and the infrastructure and conflict function, Road(t), as seen in Equation 5. As YT(t) is a component of a probability, the weighting was necessary to ensure the sum of the weather and road functions never resulted in a case of greater than 100% probability. (see Equation 5 in the Supplementary Files)
Combining transition rates and variable parameters to estimate treatment coverage for children with diarrhea
With the weather and conflict (i.e. road) components combined to form the transition rate function for the probability that a person will travel to seek treatment and thus change from the Sick health status to the In-treatment health status, it became possible to estimate the coverage of treatment services for children under five with diarrheal diseases. From the assumptions, each sick person that reached a treatment center (pST) was treated with Zn /ORS. As this transition occurred only after a person became sick, the estimation of those seeking treatment (the coverage rate) was calculated from Equation 6. (see Supplementary Files)
The model is designed that it takes one iteration (i.e. one simulated week) for a sick child to move into the In-treatment health status. A subset of the Sick population at week t would move to In-treatment at week (t+1). Thus, to calculate the treatment coverage of sick children, the In-treatment population at week (t+1) must be divided by the Sick population from week t.
Determining the model’s population size
To ensure the model output accurately matched the UNICEF provided data, the model’s population size was next determined. The under-five population was estimated at 5.16 million children.18 Half of this population was unable to reach a fixed health facility before the war started.2 With this in mind, these pre-conflict health facilities were assumed to be all functional as there was no conflict to close them. To determine what proportion of this 50% of the population was reaching care during the conflict we next examined the health facilities that were open during the conflict. According to the HeRAMS data, 45% of health facilities are fully functional and 38% are partially functional. For the purposes of this model, it was assumed that a partially operating health facility was roughly equivalent to half of a fully functioning one, and thus 38% partially functioning facilities is equivalent to 19% fully functioning. This leads to a final estimation of 64% functionality of health facilities. 53% of the facilities in the data set on facilities treating children provided by UNICEF had complete data on the number of diarrheal disease cases they saw for children under five years old. This data set consisted of about 41% of all the primary care facilities monitored by UNICEF. This information was combined to estimate the population size that has been able to reach health facilities on the ground and produce the primary data. In order to finalize transition rates in the model and still match the data outputs, the model would also need to use this same sample size which was estimated at 5,156,951 x 50% x 64.26% x 53% x 40.82% = 358,498 children.
Table 2: Table of Terms
Object Notation
|
Summary
|
Value
|
Susceptible Population
|
Total under 5 child population
|
358,498
|
Cycles
|
Number of weeks the simulation runs for
|
52
|
Maxroad
|
Max value of roads
|
variable
|
Roadmin
|
Minimum value of roads due to washout
|
equation
|
Bridge
|
Percent of open road segments with open bridges
|
variable
|
Nonbridgemin
|
Minimum road value with operating bridges
|
equation
|
Bridgecuttoff
|
Weather threshold where lower quality roads may begin to experience washout
|
0.5
|
Scalemax
|
Parameter determining how long infrastructure stays at max value
|
variable
|
Weatherweight
|
Weighting of weather on transition rate
|
0.3
|
Roadweight
|
Weighting of infrastructure on transition rate
|
0.7
|
Y(t)
|
Weather function
|
Equation
|
Road(t)
|
Road and bridge infrastructure function
|
Equation
|
YT(t)
|
Transition probability equation from Sick to Treated
|
Equation
|
Overview of the variables and parameters used throughout the model.
Deterministic Sensitivity Analysis
A 1-way simple, deterministic sensitivity analysis, one of the most widely used methods to study model uncertainty, was then performed for “maxroad”, “Roadmin” “Bridge”, and two parameters from the function Y(t), namely the conflict multiplier and the weather amplitude. 19,20 This analysis is useful as by examining the stability in the model and its parameters, we can further improve confidence in the model.21 For this 1-way deterministic analysis, each parameter was modulated between 10% and 200%
of its base value, shown in Table 3 except for the weather amplitude which could only be increased to 190% due to probability constraints in the model. The model estimation of the number of children seeking care was repeated the same way as described above for each parameter throughout the sweep, identifying which pieces of data had the greatest impact on both the maximum weekly coverage and minimum weekly coverage.
Table 3: Sensitivity Analysis Parameters
Parameter
|
Base Value
|
Values Swept
|
Maxroad
|
.4430
|
.0443-.886
|
Roadmin
|
.2479
|
.0248-.4958
|
Bridge
|
.0804
|
.00804-.1608
|
Air Strike Proportion
|
.0216
|
.00216-.0432
|
Weather Amplitude
|
7.174
|
.7174-13.6306
|
Summary of parameter baselines and ranges swept for a 1-way deterministic sensitivity analysis