Erlang loss refusal probability, annual refusals and achieved occupancy of a six-bed unit

Signature

a = lambda * L_bar; P_block = (a^6 / 720) / (1 + a + a^2 / 2 + a^3 / 6 + a^4 / 24 + a^5 / 120 + a^6 / 720); R = 365 * lambda * P_block; o_ach = a * (1 - P_block) / 6

Function

Erlang loss sizing function for hospital beds with refused admissions

Computational function

• Computational function for Erlang B refusals, beds for a refusal target and cost-minimising beds

  Takes a ward's arrival rate per day, its mean stay, a refusal target and the yearly cost of a staffed bed and of a refused admission, and returns the Erlang loss probability, refusals a year and achieved occupancy for every number of beds up to a ceiling, the smallest number of beds that meets the refusal target and the number of beds that minimises the total cost of beds and refusals. It computes the loss probability of HE-FM-CAPP-001 for any number of beds by the Erlang B recursion B_0 = 1, B_k = a B_(k-1) / (k + a B_(k-1)), which avoids the large powers and factorials of the direct formula, and its target and cost inputs differ from the formula's symbols. The total cost of N beds is c_b N plus c_r times the refusals at N beds, and one more bed is worth adding while the refusals it avoids, times c_r, exceed c_b. Because each extra bed avoids fewer refusals than the one before, the first bed that fails this test marks the cost minimum (derived; it follows from the strict convexity of the loss function in the number of beds that Tuenter reports from Jagers and van Doorn, and a grid check of 75 cases found no exception). The calculator line writes the recursion out for six beds and returns the net saving of the sixth bed; the searches over N are in the Excel, R and Python code.

  Inputs and outputs: lambda is the arrival rate in patients per day, refused patients included; L_bar is the mean stay in days; a is the offered load in beds; B_1 to B_6 are the loss probabilities with one to six beds; R_5 and R_6 are the refusals a year with five and six beds; c_b is the yearly cost of one staffed bed and c_r the cost of one refused admission, in GBP; S_6 is the net saving of the sixth bed in GBP a year, positive when the bed is worth adding. In the code target is the refusal target as a proportion, Nmax the largest number of beds searched, and the outputs are the table of loss probability, refusals and occupancy by number of beds, the beds for the target and the cost-minimising beds.

  Assumption: The assumptions of HE-FM-CAPP-001 hold for every number of beds (HE-AS-CAPP-001 to HE-AS-CAPP-004). Each staffed bed costs c_b a year whether used or not and each refusal costs c_r, both constant, the article's illustrative convention; Gorunescu and colleagues balance the cost of empty beds against the cost of delayed patients instead; that convention would change the cost line of the code here but not the recursion (derived).

  Worked example (De Bruin and colleagues' coronary care unit at GBP 10,000 a refusal): With the unit's 3.52 arrivals a day and 1.694-day stay, illustrative costs of GBP 150,000 a bed and GBP 10,000 a refusal, the sixth bed cuts refusals by 122.7 a year, to 337.0, worth GBP 1.23 million against its GBP 150,000, so it pays; the code puts the cost minimum at 11 beds (computed here). lambda = 3.52; L_bar = 1.694; c_b = 150000; c_r = 10000; a = 5.96288; B_1 = 0.8564; B_2 = 0.7186; B_3 = 0.5882; B_4 = 0.4672; B_5 = 0.3578; B_6 = 0.2623; R_5 = 459.70; R_6 = 337.02; S_6 = 1076836

  Worked example (Same unit at GBP 1,000 a refusal): At an illustrative GBP 1,000 a refusal the sixth bed saves GBP 122,684 of refusal costs, less than its GBP 150,000, and the code puts the cost minimum at 4 beds (computed here). lambda = 3.52; L_bar = 1.694; c_b = 150000; c_r = 1000; a = 5.96288; B_1 = 0.8564; B_2 = 0.7186; B_3 = 0.5882; B_4 = 0.4672; B_5 = 0.3578; B_6 = 0.2623; R_5 = 459.70; R_6 = 337.02; S_6 = -27316

  Worked example (Article's acute medical unit, for reading against the code): With 20 arrivals a day and a 4.0-day stay, the offered load is 80 beds and 95 beds refuse 1.14% of arrivals, 83.0 patients a year, at 83.3% occupancy. Refusals are 99.7, 83.0 and 68.5 a year at 94, 95 and 96 beds, so at GBP 150,000 a bed and GBP 10,000 a refusal the 95th bed saves GBP 167,000 and the 96th GBP 145,000, and the cost minimum is 95 beds; at GBP 5,000 it is 89 and at GBP 20,000 it is 99. The unit one-tenth the size, 2 arrivals a day, refuses 12.2% with 10 beds, 88.8 patients a year, and needs 15 beds and 52.8% occupancy to hold refusals below 1%, all as in the article. This example is for reading against the R and Python code; the scalar lines above run on the site.

  Excel: With named cells ArrivalsPerDayCF, MeanStayCF, BlockTarget, BedCost and RefusalCost, enter the bed numbers 0 to 199 in A2:A201 and name the range BedCount. Enter 1 in B2 and =ArrivalsPerDayCF*MeanStayCF*B2/(A3+ArrivalsPerDayCF*MeanStayCF*B2) in B3, fill it down to B201 and name B2:B201 ErlangBColumn. Enter =365*ArrivalsPerDayCF*B2 in C2 and =BedCost*A2+RefusalCost*C2 in D2, fill both down and name them RefusalsColumn and TotalCostColumn. A cell named BedsForTarget holds =INDEX(BedCount,MATCH(TRUE,INDEX(ErlangBColumn<=BlockTarget,0),0)) and a cell named CostMinBeds holds =INDEX(BedCount,MATCH(MIN(TotalCostColumn),TotalCostColumn,0)). With a number of beds in BedsInPlan, =ArrivalsPerDayCF*MeanStayCF*(1-INDEX(ErlangBColumn,BedsInPlan+1))/BedsInPlan returns its achieved occupancy. The INDEX form with a column argument of 0 makes MATCH read the comparison as an array without array entry, so the formulas work in every Excel version; wards above 199 beds need longer columns.

  R: erlang_ward <- function(lambda, L, Nmax = 200, target = 0.05, cb = 150000, cr = 10000) { a <- lambda*L; B <- numeric(Nmax+1); B[1] <- 1; for (k in 1:Nmax) B[k+1] <- a*B[k]/(k+a*B[k]); N <- 0:Nmax; R <- 365*lambda*B; occ <- ifelse(N > 0, a*(1-B)/N, NA); TC <- cb*N+cr*R; list(table = data.frame(N = N, B = B, R = R, occ = occ), beds_target = N[which(B <= target)[1]], beds_costmin = N[which.min(TC)]) } Called as w <- erlang_ward(20, 4, target = 0.01) it gives w$table$B[96] of 0.01137, w$table$R[95:97] of 99.69, 82.99 and 68.51, w$beds_costmin of 95 and w$beds_target of 96, and erlang_ward(2, 4, target = 0.01)$beds_target is 15.

  Python: def erlang_ward(lam, L, nmax=200, target=0.05, cb=150000, cr=10000): a = lam*L; B = [1.0]; [B.append(a*B[-1]/(k+a*B[-1])) for k in range(1, nmax+1)]; R = [365*lam*b for b in B]; tc = [cb*n+cr*r for n, r in enumerate(R)]; return B, R, next(n for n, b in enumerate(B) if b <= target), min(range(nmax+1), key=tc.__getitem__) Standard library only; erlang_ward(20, 4, target=0.01) returns B[95] of 0.01137, beds for the target of 96 and cost-minimising beds of 95, and erlang_ward(20, 4, cr=5000)[3] is 89.

  Test (Article's acute medical unit at 95 beds): With ArrivalsPerDayCF 20 and MeanStayCF 4, the 96th cell of ErlangBColumn, the loss probability with 95 beds, holds 0.01137 within 0.00005. Expected result TRUE. FALSE shows the Poisson overflow at a mean of 80, 0.0446, in place of the recursion, a column that starts from 0 in B2, which returns 0 in every row, and a doubled (0.0227) or halved (0.0057) cell. Excel check =ABS(INDEX(ErlangBColumn,96)-0.01137)<0.00005

  Test (Article's cost-minimising beds): With the same inputs, BedCost 150,000 and RefusalCost 10,000, CostMinBeds is 95, as in the article. Expected result TRUE. FALSE shows a refusal cost entered as 5,000, which gives 89, and a doubled (190) or halved (47.5) cell. Excel check =CostMinBeds=95

  Test (De Bruin and colleagues' merged cardiac units): With ArrivalsPerDayCF 8.37 (3.52, 3.06 and 1.79 a day added), MeanStayCF 1.96, the weighted mean stay in their Table 7, and BlockTarget 0.05, BedsForTarget is 22, the beds they report for 5% refused admissions after merging, against 10, 12 and 7 (29 in all) for the three units apart. Expected result TRUE. FALSE shows a target typed as 5 in place of 0.05, which returns 0 beds, and a doubled (44) or halved (11) cell. Excel check =BedsForTarget=22

  Common error (Factorials in one cell for a large ward): Writing the loss formula with FACT and powers in one cell fails above 170 beds, because FACT returns #NUM! for arguments above 170, and powers of the load overflow for large wards in the same way. The recursion keeps every term at or below 1.

  Common error (One occupancy target for every unit size): Sizing the article's small unit at 85% occupancy gives 10 beds and 12.2% refused, about one patient in eight, while the same rule gives 1.14% for the 95-bed unit; the beds for a stated refusal target come from BedsForTarget, not from the occupancy rule.

  Source: de Bruin AM, Bekker R, van Zanten L, Koole GM. Dimensioning hospital wards using the Erlang loss model. Annals of Operations Research. 2010;178(1):23-43. Section 6.1 and Table 6, the number of required beds per ward for 2%, 5% and 10% blocking (10 beds at 5% for the coronary care unit), and section 6.3 and Table 7, 29 beds for 5% refused admissions at the three units apart and 22 after merging. Tuenter HJH. Erlang's loss function and a property of the exponential function. arXiv preprint math/0612052, 2006. Preprint section 2, equation (6), the recursion in reciprocal form with initial value B(0, lambda) = 1, and the statement that Jagers and van Doorn showed the continued loss function to be strictly convex in the number of circuits. Gorunescu F, McClean SI, Millard PH. A queueing model for bed-occupancy management and planning of hospitals. Journal of the Operational Research Society. 2002;53(1):19-24. Abstract, on optimising the number of beds to keep the delay probability low and optimising the average cost per day by balancing costs of empty beds against costs of delayed patients.

  a = lambda * L_bar; B_1 = a / (1 + a); B_2 = a * B_1 / (2 + a * B_1); B_3 = a * B_2 / (3 + a * B_2); B_4 = a * B_3 / (4 + a * B_3); B_5 = a * B_4 / (5 + a * B_4); B_6 = a * B_5 / (6 + a * B_5); R_5 = 365 * lambda * B_5; R_6 = 365 * lambda * B_6; S_6 = c_r * (R_5 - R_6) - c_b

Try this function

Implementations

• Excel

  Six-bed Erlang loss probability, refusals and occupancy from named cells

  With the arrival rate per day in a cell named ArrivalsPerDay and the mean stay in days in MeanStay, the four formulas return the offered load in a cell named OfferedLoad, the refusal probability in BlockProbSixBeds, the refusals a year in RefusalsPerYear and the achieved occupancy in AchievedOccupancy. For de Bruin and colleagues' coronary care unit, with 3.52 and 1.694, they return 5.96288, 0.2623, 337.02 and 0.7331.

  =ArrivalsPerDay*MeanStay; =(OfferedLoad^6/720)/(1+OfferedLoad+OfferedLoad^2/2+OfferedLoad^3/6+OfferedLoad^4/24+OfferedLoad^5/120+OfferedLoad^6/720); =365*ArrivalsPerDay*BlockProbSixBeds; =OfferedLoad*(1-BlockProbSixBeds)/6

Assumptions

• Poisson arrivals at a constant rate in the six-bed loss model

  Patients arrive one at a time and independently at a constant mean rate lambda, as de Bruin and colleagues assume, and the formula describes the steady state reached when that rate has held long enough. Green reports that the Poisson process represents unscheduled arrivals to intensive care units, obstetrics units and emergency departments well; booked elective admissions and weekday patterns break the assumption and call for a time-dependent model or simulation. A loss system reaches a steady state at any load, unlike a queue, because refused patients do not accumulate (derived).

• Refused patients leave the six-bed unit and do not wait

  A patient who finds all six beds occupied is refused and leaves the unit, which in practice means diversion to another hospital or admission to another ward, and does not wait or come back later. If patients who find the unit full wait for a bed instead, the measure that matters is the probability of delay of the M/M/s model in Green's equation 11-1, not P_block.

• Any length of stay distribution with the same mean in the six-bed loss model

  De Bruin and colleagues note that the loss model is insensitive to the length of stay distribution and valid for general service times, so only the mean stay L_bar enters. Stays are independent of each other and of the arrival process.

• Six interchangeable staffed beds and arrivals counted as demand

  All six beds can take any patient of the unit and no extra bed is opened when the unit is full, and lambda counts every patient who needs a bed, refused patients included. Hospital systems usually record only admissions, so de Bruin and colleagues estimate lambda from the average number of occupied beds, which the hospital system does record, by solving their expression (iii) numerically.

Worked examples

• De Bruin and colleagues' six-bed coronary care unit in 2006

  De Bruin and colleagues' coronary care unit had 6 operational beds in 2006, an estimated 3.52 arrivals a day and a mean stay of 1.694 days, so the offered load is 5.96 beds. The formula refuses 26.2% of arrivals and the unit runs at 73.3% occupancy, the refused fraction and occupancy rate in their Table 7; that is about 337 refused patients a year (computed here).

  lambda = 3.52; L_bar = 1.694; a = 5.96288; P_block = 0.2623; R = 337.02; o_ach = 0.7331

• Six-bed unit held to 5% refused admissions

  Illustrative inputs of 1.48 arrivals a day and a mean stay of 2 days give an offered load of 2.96 beds, at which the six-bed unit refuses 5.0% of arrivals, about 27 a year, and runs at 46.9% occupancy, below the 50% that de Bruin and colleagues give for a small unit held to 5% refused admissions.

  lambda = 1.48; L_bar = 2; a = 2.96; P_block = 0.05; R = 27; o_ach = 0.4687

• Overloaded six-bed unit with five times the load it has beds for

  Illustrative inputs of 6 arrivals a day and a mean stay of 5 days give an offered load of 30 beds. The unit is almost always full, with achieved occupancy of 96.2%, and refuses 80.8% of arrivals, about 1,769 patients a year. As the load grows the achieved occupancy approaches 1 and P_block approaches 1 minus 6 / a (derived), a limiting case in which a full unit says nothing about how many patients it turns away.

  lambda = 6; L_bar = 5; a = 30; P_block = 0.8076; R = 1768.62; o_ach = 0.9621

Common errors

• Recorded admissions entered as the arrival rate of the six-bed loss model

  Hospital systems record admissions, not refused patients, so admissions understate demand wherever beds run short. For de Bruin and colleagues' coronary care unit, entering the roughly 2.60 patients a day actually admitted (computed here from 3.52 arrivals with 26.2% refused) in place of 3.52 gives a refusal probability of 0.147 in place of 0.262 and about 139 refusals a year in place of 337. De Bruin and colleagues estimate the arrival rate from the average number of occupied beds instead.

• Poisson overflow used in place of the Erlang loss probability

  The probability that a Poisson count with mean a exceeds six, the complement of the Poisson sum up to the capacity that HE-CF-ARRP-001 uses for hourly arrivals, counts refused patients as if they still occupied beds. For the coronary care unit it gives 0.388 against the loss probability of 0.262, about 48% too high (computed here).

• Offered load over beds read as the achieved occupancy of a loss unit

  Dividing the offered load by the beds ignores the bed-days refused patients never use. For the coronary care unit a / 6 is 0.994, while the achieved occupancy is 0.733 (computed here), so a plan that reads the first figure as occupancy hides a quarter of its patients being turned away.

• Refusal probability read as the chance of waiting for a bed

  P_block counts patients turned away. If patients who find the unit full wait instead, the M/M/s probability of delay in Green's equation 11-1, which assumes exponential stays, applies and is higher; for the illustrative six-bed unit with a load of 2.96 it is 0.094 against a refusal probability of 0.050 (computed here). The article's 95-bed unit shows the same gap, a 6.8% chance of waiting against 1.14% refused.

Sources

• De Bruin and colleagues' Erlang loss formula and occupancy rate for hospital wards

  de Bruin AM, Bekker R, van Zanten L, Koole GM. Dimensioning hospital wards using the Erlang loss model. Annals of Operations Research. 2010;178(1):23-43. Section 4, which assumes that an arriving patient who finds all beds occupied is blocked and leaves the system, and section 4.1, The M/G/c/c queuing model, with expression (i) for the fraction of patients blocked, P_c = ((lambda mu)^c / c!) / sum over k from 0 to c of (lambda mu)^k / k!, where mu is the expected length of stay and c the operational beds, the note that the model is insensitive to the length of stay distribution, expression (ii) for the occupancy rate, one minus P_c times lambda mu, divided by c, and the name offered load for lambda mu. Section 4.2 estimates lambda from the average number of occupied beds with expression (iii). Table 6 gives the coronary care unit's 6 operational beds and 3.52 daily arrivals, and Table 7 its mean stay of 1.694 days, 73.3% occupancy and 26.2% refused admissions.

  View source →

• Tuenter on Erlang's loss function as the probability that all circuits are occupied

  Tuenter HJH. Erlang's loss function and a property of the exponential function. arXiv preprint math/0612052, version 1, 2 December 2006. Preprint section 2, equation (4), which defines Erlang's loss function B(n, lambda) as lambda^n / n! divided by the sum over j from 0 to n of lambda^j / j!, the probability that calls arriving in a Poisson stream of rate lambda to a link of n circuits, each call requesting unit service time, find all circuits occupied and are lost.

  View source →

• Green on pure loss models and the Erlang delay alternative

  Green L. Queueing analysis in healthcare. In: Hall RW, editor. Patient Flow: Reducing Delay in Healthcare Delivery. New York: Springer; 2006. p. 281-307. Section 2.2.3, Some useful extensions of the M/M/s model, which calls the finite capacity case with no waiting room for arrivals who find all servers busy a pure loss model; section 2.2.2, The M/M/s model, which gives the probability of delay in equation 11-1 for systems in which patients wait; and section 2.2.1, The Poisson process, which reports the Poisson process as a good representation of unscheduled arrivals to ICUs, obstetrics units and EDs.

  View source →

Canonical Identity