Gravity accessibility with a population demand adjustment for two practices and three neighbourhoods

Signature

V_X = P_1 * d_1X^(-beta) + P_2 * d_2X^(-beta) + P_3 * d_3X^(-beta); V_Y = P_1 * d_1Y^(-beta) + P_2 * d_2Y^(-beta) + P_3 * d_3Y^(-beta); A_1 = S_X * d_1X^(-beta) / V_X + S_Y * d_1Y^(-beta) / V_Y; A_2 = S_X * d_2X^(-beta) / V_X + S_Y * d_2Y^(-beta) / V_Y; A_3 = S_X * d_3X^(-beta) / V_X + S_Y * d_3Y^(-beta) / V_Y

Function

Spatial accessibility index function from provider capacity, population and travel times

Computational function

• Computational function for gravity accessibility scores, mean and spread for any number of sites

  Takes the capacities of any number of provider locations, the populations of any number of residential locations, the matrix of travel times between them and a friction coefficient, and returns each provider's demand potential, each location's gravity score and the population-weighted mean and variance of the scores, scaled to a stated number of people. It applies HE-FM-ACCS-002 to every provider and location at once, so its inputs are a matrix and two lists rather than the formula's named scalars; the mean and variance are as in HE-CF-ACCS-001. The calculator line repeats the chain for the article's two practices and three neighbourhoods; any other size is in the Excel, R and Python code.

  Inputs and outputs: S_X, S_Y, P_1 to P_3 and d_1X to d_3Y are as in HE-FM-ACCS-002; beta is the friction coefficient; k the number of people the results are scaled to; V_X, V_Y are the demand potentials and A_1 to A_3 the scores in GPs per person; M and V are the population-weighted mean per k people and variance per k people, squared. In the code S, P and D are as in HE-CF-ACCS-001.

  Assumption: The assumptions of HE-FM-ACCS-002 hold (HE-AS-ACCS-004 to HE-AS-ACCS-006). Every travel time is above zero; the R code stops with a message and the Python code raises AssertionError when one is not, and Excel returns #DIV/0! for a zero or blank time when beta is above zero and #NUM! when beta is 0. The total population must be above zero; with every population at 0 the Python code raises ZeroDivisionError and the R code returns infinite scores and a NaN mean and variance.

  Worked example (Article's inputs at beta of 1): The scores are 4.440, 2.978 and 2.672 per 10,000, the mean 3.333 and the variance 0.483 (computed here). S_X = 8; S_Y = 4; P_1 = 10000; P_2 = 20000; P_3 = 6000; d_1X = 10; d_1Y = 40; d_2X = 20; d_2Y = 25; d_3X = 45; d_3Y = 15; beta = 1; k = 10000; V_X = 2133.333; V_Y = 1450; A_1 = 0.000443966; A_2 = 0.000297845; A_3 = 0.000267241; M = 3.3333; V = 0.4828

  Worked example (Article's inputs at beta of 2): The scores are 5.615, 2.293 and 2.997 per 10,000 and the variance 2.066 (computed here). S_X = 8; S_Y = 4; P_1 = 10000; P_2 = 20000; P_3 = 6000; d_1X = 10; d_1Y = 40; d_2X = 20; d_2Y = 25; d_3X = 45; d_3Y = 15; beta = 2; k = 10000; V_X = 152.963; V_Y = 64.917; A_1 = 0.000561513; A_2 = 0.000229339; A_3 = 0.000299683; M = 3.3333; V = 2.066

  Excel: With TravelTimes, PopColumn and SupplyRow laid out as for HE-CF-ACCS-001 and the coefficient in a cell named BetaCF, enter in a row of n cells named DemandRow =SUMPRODUCT(PopColumn,INDEX(TravelTimes,0,COLUMNS($A$1:A$1))^(-BetaCF)) and fill right, and in a column of m cells named AccessColumnG =SUMPRODUCT(SupplyRow/DemandRow,INDEX(TravelTimes,ROWS($A$1:A1),0)^(-BetaCF)) and fill down. A cell named MeanAccessG holds =SUMPRODUCT(PopColumn,AccessColumnG)/SUM(PopColumn) and a cell named VarAccessG holds =SUMPRODUCT(PopColumn,(AccessColumnG-MeanAccessG)^2)/SUM(PopColumn), in GPs per person.

  R: grav <- function(S, P, D, beta, k = 10000) { D <- matrix(D, nrow = length(P)); if (any(D <= 0)) stop("travel times must be above zero"); W <- D^(-beta); V <- colSums(W * P); A <- as.vector(W %*% (S / V)); M <- sum(P * A) / sum(P); list(V = V, A = k * A, mean = k * M, var = k^2 * sum(P * (A-M)^2) / sum(P)) } Called as grav(c(8, 4), c(10000, 20000, 6000), rbind(c(10, 40), c(20, 25), c(45, 15)), 1) it returns V of 2133.333 and 1450, A of 4.440, 2.978 and 2.672, mean 3.333 and variance 0.483.

  Python: def grav(S, P, D, beta, k=10000): assert min(min(r) for r in D) > 0, 'travel times must be above zero'; n, m = len(S), len(P); V = [sum(P[i] * D[i][j] ** -beta for i in range(m)) for j in range(n)]; A = [sum(S[j] * D[i][j] ** -beta / V[j] for j in range(n)) for i in range(m)]; T = sum(P); M = sum(p * a for p, a in zip(P, A)) / T; return V, [k * a for a in A], k * M, k * k * sum(p * (a-M) ** 2 for p, a in zip(P, A)) / T Standard library only; grav([8, 4], [10000, 20000, 6000], [[10, 40], [20, 25], [45, 15]], 2) returns A of 5.615, 2.293 and 2.997 and variance 2.066.

  Test (Supply conserved under gravity decay): The population-weighted sum of AccessColumnG equals the total of SupplyRow for any BetaCF. Expected result TRUE. With the article's inputs at BetaCF 1, a doubled AccessColumnG gives 24 against 12, a halved one 6, and scores without the division by DemandRow 22,866.7, all FALSE. Excel check =ABS(SUMPRODUCT(PopColumn,AccessColumnG)-SUM(SupplyRow))<1E-9

  Test (Article's neighbourhood B at beta of 1): The second cell of AccessColumnG times 10,000 is 2.9784 within 0.0001 for the article's inputs at BetaCF 1. Expected result TRUE. BetaCF entered as minus 1 gives 3.0880 and FALSE, and doubling (5.9569) or halving (1.4892) the cell also gives FALSE. Excel check =ABS(INDEX(AccessColumnG,2)*10000-2.9784)<0.0001

  Common error (Friction coefficient entered with its sign): Typing minus 1 for beta raises travel time to a positive power, so farther practices count more; for the article's inputs C, the neighbourhood farthest from the larger practice, then scores highest at 5.281 per 10,000 (computed here).

  Source: Luo W, Wang F. Measures of spatial accessibility to healthcare in a GIS environment: synthesis and a case study in Chicago region. Environment and Planning B: Planning and Design. 2003;30(6):865-884. Author manuscript read (HHS Author Manuscripts, PMC8238135); the published version, doi 10.1068/b29120, not seen. Section The Gravity-Based Method and a Synthesis, equation (6). Wang F. Measurement, optimization, and impact of health care accessibility: a methodological review. Annals of the Association of American Geographers. 2012;102(5):1104-1112. Author manuscript read (NIHPA Author Manuscripts, PMC3547595); the published version, doi 10.1080/00045608.2012.657146, not seen. Section Health Care Accessibility Measures Synthesized, equation (1).

  V_X = P_1 * d_1X^(-beta) + P_2 * d_2X^(-beta) + P_3 * d_3X^(-beta); V_Y = P_1 * d_1Y^(-beta) + P_2 * d_2Y^(-beta) + P_3 * d_3Y^(-beta); A_1 = S_X * d_1X^(-beta) / V_X + S_Y * d_1Y^(-beta) / V_Y; A_2 = S_X * d_2X^(-beta) / V_X + S_Y * d_2Y^(-beta) / V_Y; A_3 = S_X * d_3X^(-beta) / V_X + S_Y * d_3Y^(-beta) / V_Y; M = k * (P_1 * A_1 + P_2 * A_2 + P_3 * A_3) / (P_1 + P_2 + P_3); V = (P_1 * (k * A_1 - M)^2 + P_2 * (k * A_2 - M)^2 + P_3 * (k * A_3 - M)^2) / (P_1 + P_2 + P_3)

Try this function

Implementations

• Excel

  Gravity accessibility scores for two practices and three neighbourhoods from named cells

  With the inputs named as in the Excel implementation of HE-FM-ACCS-001 and the friction coefficient in a cell named FrictionBeta, the five formulas return the demand potentials in cells named DemandPotentialX and DemandPotentialY and the scores in GravityAccessA, GravityAccessB and GravityAccessC, in GPs per person. The brackets around the negative exponent keep Excel's order of operations right. For the article's inputs at FrictionBeta 1 they return 2133.333, 1450, 0.000443966, 0.000297845 and 0.000267241.

  =PopNbhdA*TimeAtoX^(-FrictionBeta)+PopNbhdB*TimeBtoX^(-FrictionBeta)+PopNbhdC*TimeCtoX^(-FrictionBeta); =PopNbhdA*TimeAtoY^(-FrictionBeta)+PopNbhdB*TimeBtoY^(-FrictionBeta)+PopNbhdC*TimeCtoY^(-FrictionBeta); =GPsPracticeX*TimeAtoX^(-FrictionBeta)/DemandPotentialX+GPsPracticeY*TimeAtoY^(-FrictionBeta)/DemandPotentialY; =GPsPracticeX*TimeBtoX^(-FrictionBeta)/DemandPotentialX+GPsPracticeY*TimeBtoY^(-FrictionBeta)/DemandPotentialY; =GPsPracticeX*TimeCtoX^(-FrictionBeta)/DemandPotentialX+GPsPracticeY*TimeCtoY^(-FrictionBeta)/DemandPotentialY

Assumptions

• Power decay with one friction coefficient for every trip in the gravity accessibility index

  Travel time discounts both supply and demand by the same power, d^(-beta), with one beta for every pair. Wang 2012 notes that beta requires additional data and work to define and might be region specific; Luo and Wang tested seven values of 1.0 to 2.2, citing 1.85 from an earlier commuting study in the same area, and found larger variances of the scores at higher values. At a beta of 1 two GPs 20 minutes away count as much as one GP 10 minutes away, a trade-off that Luo and Wang note may be considered questionable, particularly by people outside the field of geography.

• Positive travel times in one unit for the gravity accessibility index

  Every travel time must be above zero, because d^(-beta) is undefined at zero when beta is above zero (derived); a neighbourhood that contains a practice needs a positive time to it, such as the article's 10 minutes between A and X. All times are in one unit. Rescaling every time by the same factor leaves the scores unchanged, because the factor cancels between each numerator and the demand potential (derived, checked on a grid of 2,000 cases).

• No catchment edge, so every practice adds to every gravity accessibility score

  Every practice adds to every neighbourhood's score however far away it is, which removes the sharp edge of HE-FM-ACCS-001 but gives distant supply some weight. Luo and Wang found that the gravity form tends to give higher scores than the 2SFCA to areas with low accessibility and so could conceal local pockets of poor access. Like HE-FM-ACCS-001 it measures potential access for residents placed at one point per neighbourhood.

Worked examples

• Article's practices and neighbourhoods under gravity decay with beta of 1

  With the article's practices, populations and travel times and an illustrative beta of 1, the demand potentials are 2,133.3 for X and 1,450 for Y, and A, B and C score 4.440, 2.978 and 2.672 per 10,000 (computed here). The population-weighted mean is 3.333 per 10,000, as for the 2SFCA, but B now ranks below A, because B's practices are 20 and 25 minutes away while A's is 10 minutes away.

  S_X = 8; S_Y = 4; P_1 = 10000; P_2 = 20000; P_3 = 6000; d_1X = 10; d_1Y = 40; d_2X = 20; d_2Y = 25; d_3X = 45; d_3Y = 15; beta = 1; V_X = 2133.333; V_Y = 1450; A_1 = 0.000443966; A_2 = 0.000297845; A_3 = 0.000267241

• Article's practices and neighbourhoods under gravity decay with beta of 2

  With an illustrative beta of 2, A rises to 5.615 per 10,000, B falls to 2.293 and C rises to 2.997, so B and C change places, and the population-weighted variance of the scores rises to 2.066 (GPs per 10,000, squared) against 0.483 at a beta of 1 (computed here), in line with Luo and Wang's finding of larger variances at larger beta.

  S_X = 8; S_Y = 4; P_1 = 10000; P_2 = 20000; P_3 = 6000; d_1X = 10; d_1Y = 40; d_2X = 20; d_2Y = 25; d_3X = 45; d_3Y = 15; beta = 2; V_X = 152.963; V_Y = 64.917; A_1 = 0.000561513; A_2 = 0.000229339; A_3 = 0.000299683

• Gravity accessibility with no distance decay at beta of 0

  With beta of 0 every weight is 1, both demand potentials equal the 36,000 residents and every neighbourhood scores 12 / 36,000, or 3.333 per 10,000, the GPs per head of the whole area (derived). Travel time then plays no part, the same limit that HE-FM-ACCS-001 reaches with a threshold every travel time meets.

  S_X = 8; S_Y = 4; P_1 = 10000; P_2 = 20000; P_3 = 6000; d_1X = 10; d_1Y = 40; d_2X = 20; d_2Y = 25; d_3X = 45; d_3Y = 15; beta = 0; V_X = 36000; V_Y = 36000; A_1 = 0.000333333; A_2 = 0.000333333; A_3 = 0.000333333

Common errors

• Gravity sum computed without the population demand potential

  Leaving out V_X and V_Y gives the simple gravity form, which Guagliardo notes models supply only. At a beta of 1 it returns 0.9, 0.56 and 0.444 for A, B and C (computed here), GPs per minute rather than GPs per head, and the scores do not change when more people compete for a practice; Guagliardo's example is two providers at the same distance, one serving 1,000 people and one 5,000, which the simple form scores alike.

• Zero or token travel time for a neighbourhood that contains a practice

  Entering 0 minutes between A and X, where A contains practice X, makes 0 to the power minus beta undefined for any beta above zero, and Excel returns #DIV/0! in DemandPotentialX and in all three scores, which divide by it; a blank time cell counts as 0 and gives the same, and at a beta of 0 Excel returns #NUM! instead (checked in Excel). Entering 1 minute in place of the article's 10 raises A's score at a beta of 1 to 7.875 per 10,000 from 4.440 and lowers B's to 1.463 from 2.978 (computed here), so the time within a neighbourhood is an input that drives the result.

• Friction coefficient borrowed without checking how the gravity scores rank

  At a beta of 1 the article's neighbourhoods rank A, B, C (4.440, 2.978 and 2.672 per 10,000) and at a beta of 2 they rank A, C, B (5.615, 2.997 and 2.293), while the 2SFCA at 30 minutes ranks B first (computed here). The order that guides a siting or funding decision therefore depends on beta. Wang 2012 notes that beta might be region specific and that the best decay function cannot be settled without real utilisation data.

Sources

• Luo and Wang on the gravity-based accessibility measure

  Luo W, Wang F. Measures of spatial accessibility to healthcare in a GIS environment: synthesis and a case study in Chicago region. Environment and Planning B: Planning and Design. 2003;30(6):865-884. Author manuscript read (HHS Author Manuscripts, PMC8238135); the published version, doi 10.1068/b29120, not seen. Section The Gravity-Based Method and a Synthesis, equation (6), the gravity-based index as the sum over physician locations of S_j d_ij to the power minus beta divided by V_j, with V_j the sum over population locations of P_k d_kj to the power minus beta, after equation (5), Hansen's form without V_j; it states that the population-weighted average of the index equals the physician to population ratio of the whole area and shows the 2SFCA to be a special case by recoding travel times to 1 within d_0 and infinity beyond. Section The Case Study and Sensitivity Analysis tests seven values of beta of 1.0 to 2.2, cites 1.85 from a commuting study, reports larger variances at larger beta and higher gravity scores than 2SFCA scores in areas with low accessibility, and notes the trade-off between the number of physicians and travel time at a beta of 1.

  View source →

• Wang's review on the gravity-based accessibility model

  Wang F. Measurement, optimization, and impact of health care accessibility: a methodological review. Annals of the Association of American Geographers. 2012;102(5):1104-1112. Author manuscript read (NIHPA Author Manuscripts, PMC3547595); the published version, doi 10.1080/00045608.2012.657146, not seen. Section Health Care Accessibility Measures Synthesized, equation (1), with P_k, S_j, d, beta and the numbers of physician and population locations n and m defined, the description of the model as the ratio of supply to demand each discounted by a distance factor, the note that beta requires additional data and work to define and might be region specific, and the statement that the best decay function cannot be settled without analysing real utilisation. Section Optimization Models in Health Care Accessibility states the weighted mean property for equations 1 to 3.

  View source →

• Guagliardo on the simple and improved gravity models

  Guagliardo MF. Spatial accessibility of primary care: concepts, methods and challenges. International Journal of Health Geographics. 2004;3:3. Full text read (PMC394340), formula images included. Section Measuring spatial accessibility, the simplest gravity formula, A_i as the sum over j of S_j divided by d_ij to the power beta, the note that it models supply only with no adjustment for demand and the example of two providers serving 1,000 and 5,000 people, the population demand adjustment factor V_j of Joseph and Bantock as the sum over k of P_k divided by d_kj to the power beta, and the improved gravity model, A_i as the sum over j of S_j divided by d_ij to the power beta times V_j, all printed as unnumbered expressions, with the note that beta is usually unknown and must be estimated empirically.

  View source →

Canonical Identity