Reilly’s law of retail gravitation

The standard Reilly’s law of retail gravitation

A classical gravity model can be expressed as

T_ij=k M_i M_j/d_ij²

where T_ij is the interaction (trade, commuting, etc.) between locations i and j; k a constant; M_i and M_j are masses (population, total income, etc.) of locations i and j respectively; and d_ij the distance between i and j.

According to the gravity model, a retail center’s attraction to its customers is M_i /d_i² or M_j /d_j².

When the attractions from the two centers are equal, there is M_i /d_i² = M_j /d_j²

Fig. 1 Market delineation between centers i and j

---

In Fig. 1, point X shows the location where attractions from the two centers, i and j, are equal. Thus, to the left of X, attraction from i is higher and customers shop at i. The opposite is true to the right of X. Clearly, X delineates market between centers i and j.

From M_i /d_i² = M_j /d_j², it can be shown that d _iX = d_ij/(1+ √M_j/M_i). Where d_iX is the distance from center i to the market boundary X. Thus the market boundary from center j to X is d_jX=d_ij-d_iX.

The purpose of Reilly’s law of retail gravilation is to find d_iX or d_jX while other conditions (d_ij, M_i, and M_j) are known.

Reilly’s law of retail gravitation in two dimensional space

The standard Reilly’s model is one dimensional. A meaningful extension is to apply the principle to two dimensional space. As a result, several factors need to be considered.

First, while the one dimensional linear market is also the travel path, a two dimensional space requires a road network.

Second, customers travel along the road network to a given set of centers to purchase goods or services. Suppose they prefer to purchase at the nearest center. There is a need to calculate the pairwise market divides where the attractions from different centers are equal.

Third, there should be visualization of the market areas on the two dimensional space based on the principle of Reilly’s model.

The work-flow is indicated in Fig. 2. While the writing states what needs to be done at each step, the red symbol suggests the software to use for that step.

Fig. 2 Work-flow for two dimensional Reilly’s model

---

Example 1: Reilly’s law of retail gravitation in two dimensional space

In this exmple, I used a point shapefile called places_2county (Fig. 3), which contains places and designated places defined by the U.S. Bureau of Census in the Madison and St. Clair two county region of Illinois; a line shapefile roads_2county (Fig. 3), contains roads within the bi-county area with about a 19 km (11 miles) buffer of roads; the third file is a polygon shapefile which is the outline of the Madison and St. Clair two county region. This last file is used for reilly( ) function.

Fig. 3 Centers and road network                     Fig. 4 A portion of the edge list, distance in meters

---

Fig. 5 A portion of the distance matrix, distance converted to km, place names added

---

Fig. 4 shows the edge list created in QGIS, which is converted into a distance matrix (Fig. 5) using R.

Input shapefiles and the distance matrix with necessary set-up parameters into the reilly( ) , the resultant market delineations are shown in Fig. 6 to Fig. 8.

Fig. 6 Raster map                                 Fig. 7 Vector map                               Fig. 8 Smoothed vector map

---

Fig. 6 shows a raster map created in R from package SpatialPosition ;Fig. 7 shows a vector map converted from the raster map in Fig 6; Fig. 8 shows the result of processing Fig. 7 using R packages raster, rmapshaper, and rgeos. It is clear that some centers are placed within the same market due to either short distance between them or one center is significantly larger or more attractive than the other.

Example 2: Hierarchies of central places and market areas

Application of Reilly’s model can be extended to incorporate concepts from central place theory. One such notion is the association of the size of places with their market functions with larger places offerring both basic and specialized goods and services and smaller places offerting only basic services.

Following this notion, we can divide the places into size categories and assume competition occurs at various functional levels. That is, competition for basic services occurs among centers of all sizes; competition for more specialized services occurs among larger centers, and competition for the most specialized services occurs among the largert centers.

For example, the centers in the study area are divided into three population size groups, small, medium, and large. Fig 9. shows centers of all sizes and they all compete for markets in basic goods and services. Fig.10 shows medium and large sized centers and they compete for markets in more specialized goods and services. Fig. 11 shows large sized centers and they compete for markets in most specialized goods and services. At the same time, the distance matrices are calculated separately for the three point files.

Fig. 9 All centers                                 Fig. 10 Medium and large centers               Fig. 11 Large centers

---

Running the work-flow separately with the three point files and distance matrices results in Fig. 12 to Fig. 14, each showing the market areas in basic, more specialized and most specialized goods and services.

Fig. 12 Markets for basic goods   Fig. 13 Markets for more specilized goods         Fig. 14 Markets for most specialized goods

---

Combining centers and/or markets for different functions and sizes gives rise to hierarchical cental places and market areas, as shown in Fig. 15 to Fig. 17.

Fig. 15 Central place hierarchy   Fig. 16 Market area hierarchy       Fig. 17 Central place and market area hierarchies

---

Fig. 15 shows the hierarchy of central places of various sizes and functions. Fig. 15 shows the hierarchy of market areas of various sizes and functions. Fig 17 is the combination of the hierarchies of central places and market areas from Fig. 15 and Fig. 16.

Here are a few interactive maps which show various layers within a herarchical system

Interactive map of hierarchy of central places

Interactive map of hierarchy of market areas

Interactive map of hierarchies of central places and market areas

Example 3: Calculation of market area boundaries

One limitation of R package SpatialPosition reilly( ) is that it does not output results of calculations, which means that you do not know the actual distance between a center and its market boundaries except seeing them on the map. You can do the calculation by youself using the formula d _iX = d_ij/(1+ √ M _j/M_i). The work-flow is shown in Fig. 18.

Fig. 18 Work-flow for market boundary calculation

---

Example 4: General function forms and their paratemers

So far, we have been using the standard gravity model and standard Reilly’s model in which the exponent of the masses is 1 and that of the distance is 2. The general gravity model with an inverse distance weight of distanc decay takes the form of

T_ij=k M_i ^aM_j^a/d_ij^b

where a is the exponent of the masses (attractions) and b the exponent of distance. Other terms are as before.

The general gravity model with a negative exponential weight of distanc decay takes the form of

T_ij=k M_i ^aM_j^a/exp(bd_ij)

As a result, the general form of Reilly’s model becomes

d_iX = d_ij/(1+ (M_j/M_i)^a/b) and d_iX = 0.5(d_ij-ln(M_j/M_i)^a/b) for the inverse distance and negative exponential distance function respectively.

The difference in functional forms leads to the distance decay parameter b with very different magnitudes.

Fig. 22 to Fig. 24 compare distance exponents in power and negative exponenetial functions. When exponent in power function (red curve) is -1, moving the slider between -0.1 to -0.5 to change the exponent of the negative exponential function (black curve). At around -0.25, a visual examination determines the negative exnonential function has comparable values to those of the power funtion.

Fig. 22 b_IDW =-1, b_Exp =-0.1                                 Fig. 23 b_IDW =-1, b_Exp=-0.25                   Fig. 24 b_IDW=-1, b_Exp=-0.5

---

In Fig. 25 to Fig. 27, The exponent in power function (red curve) is -2, the exponent in negative exponenetial function (black curve) varies from -0.1 to -0.75. At approximately -0.65, the resultant weight reaches comparable values to those the power funtion with exponent of -2.

Fig. 25 b_IDW=-2, b_Exp=-0.1                                 Fig. 26 b_IDW=-2, b_Exp=-0.65                   Fig. 24 b_IDW=-2, b_Exp=-0.75

---

Finally, In Fig. 28 to Fig. 30, the power function (red curve) exponent is -3, and the negative exponenetial function (black curve) exponent varies from -0.1 to -1.5. At about -0.9, the resultant weight reaches comparable values to those the power funtion with exponent of -3.

Fig. 28 b_IDW=-3, b_Exp=-0.1                                 Fig. 29 b_IDW=-3, b_Exp=-0.9                   Fig. 30 b_IDW=-3, b_Exp=-1.5

---

In above examples, when the exponent becomes more negative, the distance decay is more profound in both functional forms. The ratios of power function exponent to negative exponential function exponent are 4=1/0.25, 3.1=2/0.65, and 3.3=3/0.9. Clearly, there is no strict regularity here. However, it seems that the power function exponent is at least 3 to 4 times as large as the negative exponential function exponent.

The above empirical ratios based on visual comparison may not be very accurate. An optimization analysis is performed. The objective is to seek the optimal exponent value for the negative exponential function which minimzes the difference in the weights from the power function. The process experiments different values as the exponent of the negative exponential function to discount distance values from a distance dataset. The discounted distances are then compared with the discounted distances obtained from using a power function. The search process continues until an exponent value is found which generates the smallest difference in discounted distance values from the power function. For this purpose, the R package rgenoud is used which uses the genetic algorithm for optimization. The algorithm simulates the biological process where parents give births to children with genetic characteristics from the parents. The characteristics which help reach the optimization objective will be passed on generation after genereation. The selection process continues until a convergence takes place or a pre-set number of generations is reached. The distance data used in the optimization is the distance matrix obtained from the places_2county point file and roads_2county line file. The results are shown in Fig. 31 to Fig. 33.

Fig. 31 b_IDW=-1, b_Exp=-0.21                                 Fig. 32 b_IDW=-2, b_Exp=-0.62                   Fig. 33 b_IDW=-3, b_Exp=-1

---

As can be seen, the optimal exponent value for the negative exponential function which matches the power function with exponent =-1 is -0.21; the optimal exponent value for the negative exponential function which matches the power function with exponent =-2 is -0.62; and the optimal exponent value for the negative exponential function which matches the power function with exponent =-3 is -1. These exponent values largely conform to the visual examination results.

The above discussion only focues on the difference in magnitude between power and negative exponential distance decay exponents without attraction factors (masses). Knowing the differences helps in simulation studies. The actual distance functional forms and/or exponents are often calibrated from sample information.