Queuing network models for panel sizing in oncology

Abstract

Motivated by practices and issues at the British Columbia Cancer Agency (BCCA), we develop queuing network models to determine the appropriate number of patients to be managed by a single physician. This is often referred to as a physician’s panel size. The key features that distinguish our study of oncology practices from other panel size models are high patient turnover rates, multiple patient and appointment types, and follow-up care. The paper develops stationary and non-stationary queuing network models corresponding to stabilized and developing practices, respectively. These models are used to determine new patient arrival rates that ensure practices operate within certain performance thresholds. Data from the BCCA are used to calibrate and illustrate the implications of these models.

Appendix: Aggregation for the multi-class open queueing network

Here we describe the technique used to aggregate the parameters in the multi-class open queueing network in Sect. 5. See [3] for more details on the method.

Let $J'=\{\hbox {np},\hbox {ap},\hbox {ip}\}$ be the set of queues in the network and let $j\in J'$. We start by computing the first and second moments of the aggregate service time ($D_\mathrm{ap}^+$) for patients in the active patient queue. This amounts to a weighted average of the service time of each patient type as follows:

$$\begin{aligned} \mathbb {E}\left[ D_\mathrm{ap}^+\right]= & {} \frac{1}{\sum \nolimits _{i=1}^c\varLambda _{\mathrm{ap},i}}\displaystyle \sum \limits _{i=1}^c \varLambda _{\mathrm{ap},i}\mathbb {E}[D_{\mathrm{ap},i}], \\ V\left[ \left( D_\mathrm{ap}^+\right) \right]= & {} \frac{1}{\sum \nolimits _{i=1}^c\varLambda _{\mathrm{ap},i}}\displaystyle \sum \limits _{i=1}^c \varLambda _{\mathrm{ap},i} V[(D_{\mathrm{ap},i})], \end{aligned}$$

where V[X] is the variance of random variable X and $\varLambda _{j,i}$ is the arrival rate of patient type i to queue $j\in J'$.

From these aggregate values, the squared coefficient of variance for the service time in the active patient queue ($\hbox {SCV}_{\mathrm{ap,s}}$, note that the subscript ap indicates the active patient queue and the subscript s indicates service) can be obtained as follows:

$$\begin{aligned} \hbox {SCV}_\mathrm{ap,s}= & {} \frac{1}{\varLambda _\mathrm{ap}^+ \left( \mathbb {E}[D_\mathrm{ap}^+]\right) ^2} \nonumber \\&\times \displaystyle \sum _{i=1}^c \left( \varLambda _{\mathrm{ap},i} \mathbb {E}[D_{\mathrm{ap},i}]^2 \left( \left( \frac{\sqrt{V[D_{\mathrm{ap},i}]}}{(\mathbb {E}[D_{\mathrm{ap},i}]) }\right) ^2 +1 \right) -1\right) . \end{aligned}$$

(12)

The aggregate mean arrival rate and the aggregate routing probabilities are, respectively, $\varLambda _j^+=\sum _{i=1}^c\varLambda _{j,i}$, $j\in J'$, and $r_\mathrm{ap,np}^+=1$, $r_\mathrm{ap,ip}^+=1/\varLambda _\mathrm{ap}^+\sum _{i=1}^c\varLambda _{\mathrm{ap},i} r_{(\mathrm{ap},i),\mathrm{ip}}$, $r_\mathrm{ip,ap}^+=1/\varLambda _\mathrm{ip}^+\sum _{i=1}^c \sum _{k=1}^c \varLambda _{\mathrm{ip},i} r_{\mathrm{ip},(\mathrm{ap},k)}$. Since all new patients go to queue np, the external (new) patient arrival rate $\varLambda _0^+=\varLambda$ and $r_{0,j}^+=\mathbf{1}\{j=\hbox {np}\}$.

At this point, the c patient classes are aggregated into a single class and we now consider the network to be a single class open queueing network with the aggregate parameters described above. To analyze the single class open queueing network, we next determine the SCV for the interarrival times to each queue ($\hbox {SCV}_{j,a}$) as follows:

$$\begin{aligned} \hbox {SCV}_{\mathrm{np},a}= & {} \alpha _\mathrm{np}, \nonumber \\ \hbox {SCV}_{\mathrm{ap},a}= & {} \alpha _\mathrm{ap}+\hbox {SCV}_{\mathrm{np},a}\beta _\mathrm{np,ap}+ \hbox {SCV}_{\mathrm{ip},a}\beta _\mathrm{ip,ap}, \nonumber \\ \hbox {SCV}_{\mathrm{ip},a}= & {} \alpha _\mathrm{ip}+\hbox {SCV}_{\mathrm{ap},a}\beta _\mathrm{ap,ip}, \end{aligned}$$

(13)

where $\alpha _j$ and $\beta _{i,j}$ are constants depending on the input data:

$$\begin{aligned} \alpha _j= & {} 1+w_j\left( I_j\left( \frac{\sqrt{V[\varLambda ]}}{\varLambda } \right) ^2-1\right) \\&+ w_j\left( \displaystyle \sum _{k\in J'}\frac{\varLambda _k^+ r_{k,j}^+}{\varLambda _j^+}\left( \left( 1-r_{k,j}^+\right) +r_{k,j}^+\rho _k^2x_k \right) \right) , \\ \beta _{k,j}= & {} w_j r_{k,j}^+\frac{\varLambda _k^+ r_{k,j}^+}{\varLambda _j^+}\left( 1-\rho ^2_k\right) ,\quad k\in J', \end{aligned}$$

where $I_\mathrm{np}=1, I_\mathrm{ap}=I_\mathrm{ip}=0$ and

$$\begin{aligned} w_{j}= & {} \left( 1+4(1-\rho _{j})^2(v_{j}-1)\right) ^{-1}, \\ v_{j}= & {} \left( \displaystyle \sum _{k\in J' \cup \{0\}} \left( \frac{\varLambda _k^+ r_{k,j}^+}{\varLambda _j^+} \right) ^2 \right) ^{-1}, \\ \rho _{j}= & {} \frac{\varLambda _{j}^+\mathbb {E}[D_{j}^+]}{m_{j}}, \\ x_{j}= & {} 1+m_{j}^{-0.5}(\max [\hbox {SCV}_{j,s},0.2]-1),\quad j\in J'. \end{aligned}$$

Since $m_\mathrm{np}$ and $m_\mathrm{ip}$ are large, then $x_\mathrm{np}\approx 1,x_\mathrm{ip}\approx 1,$ and $\rho _\mathrm{np}<1, \rho _\mathrm{ip}<1$.

Using the SCV for the arrival process (13) and the SCV for the service time (12), the expected waiting time in the active patient queue is approximated as follows (see [16]):

$$\begin{aligned} \mathbb {E}[W]\approx \frac{\hbox {SCV}_{\mathrm{ap},a}+\hbox {SCV}_{\mathrm{ap},s}}{2} \frac{\rho _\mathrm{ap}^{\left( \sqrt{2(m_\mathrm{ap}+1)}-1 \right) }}{m_\mathrm{ap}(1-\rho _\mathrm{ap})}\mathbb {E}\left[ D_\mathrm{ap}^+\right] . \end{aligned}$$