Forecasting the medical workforce: a stochastic agent-based simulation approach

Abstract

Starting in the 50s, healthcare workforce planning became a major concern for researchers and policy makers, since an imbalance of health professionals may create a serious insufficiency in the health system, and eventually lead to avoidable patient deaths. As such, methodologies and techniques have evolved significantly throughout the years, and simulation, in particular system dynamics, has been used broadly. However, tools such as stochastic agent-based simulation offer additional advantages for conducting forecasts, making it straightforward to incorporate microeconomic foundations and behavior rules into the agents. Surprisingly, we found no application of agent-based simulation to healthcare workforce planning above the hospital level. In this paper we develop a stochastic agent-based simulation model to forecast the supply of physicians and apply it to the Portuguese physician workforce. Moreover, we study the effect of variability in key input parameters using Monte Carlo simulation, concluding that small deviations in emigration or dropout rates may originate disparate forecasts. We also present different scenarios reflecting opposing policy directions and quantify their effect using the model. Finally, we perform an analysis of the impact of existing demographic projections on the demand for healthcare services. Results suggest that despite a declining population there may not be enough physicians to deliver all the care an ageing population may require. Such conclusion challenges anecdotal evidence of a surplus of physicians, supported mainly by the observation that Portugal has more physicians than the EU average.

Appendix

1.1 A Mathematical description of the model

The ABM proposed can be loosely approximated by a set of differential equations. Keep in mind, however, that the differential equations that define the stock and flow equations emerge out of the agents’ actions, and are thus not hardcoded in the model. In fact, they provide an aggregate overview of the number of physicians at each state, and not a precise account for each individual physician. Under this formulation it is not possible, for instance, to know the age of an individual physician or if that particular physician is retired, only the aggregate number of physicians retired by year and speciality. It is also not possible, or extremely laborious, to incorporate microeconomic foundations. Notwithstanding, this alternative formulation may be of interest to those well versed in math or with experience in SD simulation, for whom it will be easier to follow through.

Let t be the time period ($t \in T=\{0,\dots ,T\}$), s the state the physician is in ($s \in \mathcal {S} =\{1,2,3,\dots \}$) from the list of states described in Fig. 2 and i the physician specialty ($i \in \mathcal {I} = \{0,1,2,\dots ,52\}$). The stock S of physicians by specialty i at time t and state s can be computed by integrating the following equation:

$$ S_{is}(t) = {{\int}_{0}^{T}} [\operatorname{Inflows}_{is}(t) - \operatorname{Outflows}_{is}(t)]\operatorname{dt} + S_{is}(0) $$

(4)

The Inflows function represents all positive flows of physicians with specialty i accruing at each period t in state s, while Outflows represents negative flows to the stock i at period t and state s. S _{i s} (0) represents the initial stock of physicians in state s with specialty i. Which inflows and outflows affect the stock S will depend on the state s the physician is in. For instance, dropping out of school only applies to agents still in training (see Fig. 2). The stock equations for each state s will then be an adaptation of Eq. 4.

For instance, let us consider the case of students undergoing training. These students are not specialized yet, and so we assume i=0. The growth rate in the stock of students in medical school is given by:

$$ \frac{S_{1}(t)}{\operatorname{dt}} = M(t) + (x - d - m_{a})\cdot S_{1}(t) - \delta\cdot S_{1}(t), $$

(5)

where M(t) is the fixed number of fresh students joining medical school, x is the net migration rate (i.e. immigrations - emigrations), d is the dropout rate, m is the mortality rate at age a and δ⋅S ₁(t) is the share of already enrolled students moving to the next state (S ₂). The mortality rate is obtained from a life table, more specifically q _a , the probability that someone aged exactly a will die before reaching (a+1).

The stock equations for the remaining states are similar, except for one detail. Since “Med School” is the initial state, there are no agents transitioning from a previous state, only an initial stock and a constant positive flow M of fresh students arriving at each period t. However, subsequent states will have previous states. For instance, medical students transition from state “Med School” to “IAC”. Hence, we should also include a term representing all the positive flows arising from the previous state as follows:

$$ \frac{S_{s}(t)}{\operatorname{dt}} = (x - d - m_{a})\cdot S_{s}(t) - \delta\cdot S_{s}(t) + \omega\cdot S_{s-1}(t), $$

(6)

where ω⋅S _s−1(t) is the share of physicians moving from the previous state to current state s and δ⋅S _s (t) represents the share of physicians moving onwards to the next state.

To make it clear, let us instantiate a particular case. Consider state S ₂, corresponding to “IAC”, the practice year between medical school S ₁ and residency training S ₃. We assume students are uniformly distributed throughout college years, and so the probability of completing medical school S ₁ and moving to next state S ₂ can be approximated by the probability of being in the last year of studies. However, the time required to complete studies is not deterministic. Although medical school takes no less than 6 years to complete, some students take longer. We model the probability of shifting from state by first considering that the number of years λ necessary to fulfil the requirements follows a probabilistic distribution. In this case, we assume the exponential distribution, as it exhibits the best fit to the data:

$$ \lambda \sim \operatorname{Exponential}. $$

(7)

Then, assuming all years are equally laborious, the probability of transitioning to the next state is equal to the probability of being in the final year, which is equal to:

$$ \omega = p_{1,2} = \frac{1}{\lambda}, $$

(8)

and the share of physicians moving from the previous state S ₁ to the current state S ₂ is then approximated by:

$$ \omega S_{1}(t) = \frac{S_{1}(t)}{\lambda}. $$

(9)

We also need to compute the probability of moving to state S ₃ to obtain the complete equation. Transition to S ₃ implies obtaining a vacancy at a residency program. There are several types of residencies corresponding to each medical specialty, and the probability of moving to any residency program can be approximated by the ratio between the number of residency vacancies available for year t, n(t), and the number of students wishing to apply in that year, ω⋅S ₂(t):

$$ p_{2,3}= \left\{ \begin{array}{ll} n(t)/\omega S_{2}(t) & \quad \text{if } n(t) < \omega S_{2}(t) \\ 1 & \quad \text{if } n(t) \geq \omega S_{2}(t) \end{array} \right. $$

(10)

Consequently, the number of physicians moving to residency, δ⋅S ₃, can be approximated by p _2,3⋅S ₃(t).

When the physician is already employed the equations change slightly. To calculate the number of physicians with specialty i in the public sector (s=5) at time t we need to integrate the following equation:

$$ \frac{S_{i5}(t)}{\operatorname{dt}} = (x - y - m_{a} - r)\cdot S_{5i}(t) + \omega\cdot S_{4i}(t), $$

(11)

where y is the net migration rate from the public to the private sector and r is the rate of physicians retiring at year t. This is required since, at any time, the physician may move back and forth between the public and the private sector. When calibrating the model we assume that the net transfer rate is negative with regards to the public sector. After working for some years the physician retires. The retirement age is drawn from a normal distribution:

$$ R\ \sim\ \mathcal{N}(\mu,\,\sigma^{2}), $$

(12)

and calibrated according to the distributions specified in Table 4.

The remaining differential equations are trivial to build following the same procedure.