Internship prep: Spatial Interaction Models

Based on Groningen Spatial Analysis Workshop - Session 6: Spatial Interaction Models By Chris Brunsdon and Lex Comber

What are spatial interaction models?

Spatial interaction models are mathematical models which , in general terms, deal with flows and movement between places, based on:

• their spatial separation

• their complementarity

• whatever other intervening opportunities or spatial structural elements serve to augment or diminish the expected flow

suppressWarnings(if(!require("pacman")) install.packages("pacman"))

## Loading required package: pacman

pacman::p_load('tidyverse', 'sf', 'tmap')

Download files, unzip them and read them in R.

• sf::st_read(): Read simple features or layers from file or database

library(here)
# Download file from the internet
#download.file("https://www.dropbox.com/s/8yz4r73v9zmugu7/groningen_polys.zip?dl=1","groningen_polys.zip")

# Extract files
unzip("groningen_polys.zip")

# Read simple features
gron_polys <- st_read(".", "groningen_polys")

# Download file
#download.file("https://www.dropbox.com/s/t6mjbycebknzjd4/shops_pts.zip?dl=1","shops_pts.zip")

# Extract files
unzip("shops_pts.zip")

# Read sf 
gron_shops <- st_read(".", "shops_pts")

#download.file("https://www.dropbox.com/s/frriyfou3ujj998/groningen_pts.zip?dl=1","groningen_pts.zip")
unzip("groningen_pts.zip")
gron_pts <- st_read(".","groningen_pts")

Plot administrative areas in Groningen

• The first element always is tm_shape(), which specified the input shape object. Next, map layers, additional map elements, and overall layout can be customized. — Great like ggplot!

# View data set
gron_polys %>% 
  slice_head(n = 5)

# Set tmap mode to interactive viewing
tmap_mode(mode = "view")

## tmap mode set to interactive viewing

# Make map
gron_polys %>% 
  tm_shape() +
  tm_borders(col = "black", lwd = 2, alpha = 0.7) +
  tm_fill(col = "midnightblue", alpha = 0.2)

Plot the stores and their respective floor areas

## Legend for symbol sizes not available in view mode.

A third data set relates to the administrative areas, but as a set of points - this is useful to compute ‘crow fly’ distances between the administrative areas and the shops. These data sets are called:

• gron_pts - label points for the administrative areas
• gron_polys - polygons for the administrative areas
• gron_shops - point locations for the stores

Let’s view the gron_pts data

gron_pts

The variables here are the population AANT_INW - and the potential spending for each administrative unit Money - this is a simple estimate of €20 per person. The polygons object contains the same information - it just uses polygon geometry instead of point.

Thinking about distances

The data we currently have tells you where the stores are, and where the population is located. It also tells you how much spending power is in each place, and how attractive each store is (measured by floor space). To predict interaction the final piece in the jigsaw is the distance between the stores and the population. In broad terms, the interaction between place \(i\) and store \(j\) (annotated here is \(I_{ij}\)) is modelled as

\[ I_{ij} = f(P_i,A_j,d_{ij}) \]

where \(P_i\) is the productiveness in place \(i\) (Money),

\(A_j\) is the attractiveness of store \(j\) (Floorspace)

and \(d_{ij}\) is the distance between the two.

So we need to compute the distances. One way of doing this is to extract the point coordinates from gron_pts and gron_shops and compute the straight line (crow fly) distances between each pair, using Pythagoras’ theorem. The following function does this.

• st_coordinates {sf}: retrieve coordinates in matrix form

• base::outer: The outer product of the arrays X and Y is the array A with dimension c(dim(X), dim(Y))

\[ d_{ji} = \sqrt (shop_jx - place_ix)^2 +(shop_jy - place_iy)^2 \]

# The first step would be to extract the distances from each df
# Retrieve coordinates in matrix form
shop_coord <- st_coordinates(gron_shops)
rownames(shop_coord) <- gron_shops$shopcode
area_coord <- st_coordinates(gron_pts)
rownames(area_coord) <- gron_pts$BU_CODE

# Compute euclidean distance between each shop coordinate and place in KM (x/1000)
dmat <- sqrt(outer(shop_coord[, 1], area_coord[, 1], FUN = "-")^2 + outer(shop_coord[, 2], area_coord[, 2], FUN = "-")^2)/1000


# Putting this in a function
get_dist <- function(gron_shops, gron_pts){
  shop_coord <- st_coordinates(gron_shops)
  rownames(shop_coord) <- gron_shops$shopcode
  area_coord <- st_coordinates(gron_pts)
  rownames(area_coord) <- gron_pts$BU_CODE
  dmat <- sqrt(outer(shop_coord[, 1], area_coord[, 1], FUN = "-")^2 + outer(shop_coord[, 2], area_coord[, 2], FUN = "-")^2)/1000
  
  return(dmat)
}

Now we have a distance matrix, but to merge it with the rest of the information, it would be better to have this in tidy format - essentially in the format of a data frame, with three columns - the store ID, the administrative area ID and the distance between the two. This is useful because it may then be used as a look-up table for a given store, or a given administrative area, to find distances. To do this, first convert the matrix object into a data frame, and add a variable column based on the row name - this gives the shop code. Here we use tidyverse-style pipelining to ‘mutate’ the data frame.

# Convert matrix to a tibble
dist_tab <- as_tibble(gron_shops) %>% 
  select(shopcode) %>% 
  bind_cols(as_tibble(dmat))

dist_tab %>% 
  slice_head(n = 5)

This is still not in the format mentioned earlier, as each origin has its own column. To re-shape the data into three columns the pivot_longer command is needed. This takes data spread over several columns (here all of the origins) and ‘gathers’ them into a single column. The column names get a new variable, and the distance values are stacked.

• Basically putting all the places columns under 1 column instead of them being spread out in 108 columns.

dist_tab %>% 
  pivot_longer(!shopcode, names_to = "BU_CODE", values_to = "Dist") -> dist_tab

# View dist mat
dist_tab %>% 
  slice_head(n = 5)

Next, we want to join the attractiveness and productiveness variables onto this list, using database join-type commands that are provided by the left_join function. First off, join the attractiveness information from gron_shops :

• By default use a left_join: includes all rows in x.

dist_tab %>% 
  left_join(gron_shops, by = "shopcode") %>% 
  left_join(gron_pts, by = "BU_CODE") %>% 
  select(-starts_with("geometry")) -> dist_tab

# View data with all model variables
dist_tab %>% 
  slice_head(n = 5)

So now dist_tab has a series of rows corresponding to each \(i,j\) pair, with the values of \(d_{ij}\) \(P_i\) and \(A_j\) among the variables -

• \(d_{ij}\) is Dist
• \(P_i\) is Money
• \(A_j\) is Floorspace

We can now actually carry out the modelling.

Fitting the model

Before we noted that the spatial interaction model takes the form \(I_{ij} = f(P_i,A_j,d_{ij})\) - here, we will use a specific, multiplicative format:

\[ I_{ij} = \textsf{const}_j \times A_i^\alpha g(d_{ij}) \] where \(g\) is a decaying desirability function of distance. Note that the productiveness does not appear explicitly here - in fact it features as a constraint - essentially for each administrative area, the total money going to all of the stores must add up to the total money for that area. This is achieved by multiplying the attractiveness and distance-deterrence terms by a constant for each area - the constant re-scales the amount spent so that the total for the area is equal to the variable Money for that area.

The form that \(g\) takes here will be that of exponential decay. Often this is written as \(e^{-\textsf{const} \times d_{ij}}\) but here we will use the form

\[ g(d_{ij}) = 2^{-d_{ij}/h} \] where \(h\) is the distance half-life. That is, for every \(h\) kilometers between the store and the home location, the desirability of making the trip halves. Although this can be written in more conventional exponential form as \(e^{-\log(2)d_{ij}/h}\), the parameter \(h\) is more easily interpreted here. Thus the full model is

\[ I_{ij} = \textsf{const}_j A_i^\alpha 2^{-d_{ij}/h} \]

Having arrived at the equation we now need to compute this and add a column to the dist_tab data frame to model the \(I_{ij}\) values. These will be the estimated expenditure for the population of area \(j\) at store \(i\). For a starting point, assume a fairly strong distance deterrence with a distance half-life of 200m. We can change this later.

The first part of the model (before applying the production constraint) is to compute \(A_i^\alpha 2^{-d_ij/h}\). \(\alpha\) here allows for extra attractiveness of the larger stores. In this case it is assumed to be quite close to 1 - say 1.05 - so that larger stores have a slightly disproportionate increase in attractiveness, but the effect is fairly weak.

dist_tab %>% 
  mutate(Flow_Factor=(Floorspace^1.05)*2^(-Dist/0.2)) -> sim_tab
head(sim_tab)

This has the Flow_Factor variable proportional to the cash flow between \(i\) and \(j\) but as yet it hasn’t been scaled. To do this, we first need to rescale the flow factors for each administrative area so that they sum to one. Then, they are each multiplied by the Money variable for that area, so that the total modelled flow of cash equals the total money coming from that area. This is done by using the sum function. If the sim_tab data frame is first grouped by the administrative unit code, (using group_by(BU_CODE)) then the sum function for each row applies to the sum of values for its specific BU_CODE - rather than the sum for the entire column. Thus

sim_tab %>% 
  group_by(BU_CODE) %>% 
  mutate(Cash_Flow = (Flow_Factor/sum(Flow_Factor)) * Money) -> sim_tab

sim_tab %>% 
  slice_head()

Thus we now have a data frame whose rows are the details of each \(i\) and \(j\) combination, including the modelled cash flow. We can first check that the constraint has worked:

• unique(Money) could have worked too. Not sure why distinct() is not working

sim_tab %>% 
  group_by(BU_CODE) %>% 
  summarise(out_flow = sum(Cash_Flow), area_money = first(Money))

Here the summarise function creates a summary data frame, where statistics such as the sum or mean are computed for each group. Thus Out_Flow is the sum of all of the cash flows coming out of each administrative area. Note that the Money variable is the same for each member of the group in sim_tab so taking the first is just a way of summarising several repeated values with a single incidence of this value. Thus the output shows that the summed modelled cash flows from each adminstrative area is the same as the Money for that area. The constraint has worked. Phew!

Now we are happy with the model we may want to see what happens for different half-life distances. As with the data mungeing we can express the modelling process as a single pipeline. But to save time we can also put it into a function, with the half-life as an argument.

fit_sim <- function(h, dist_tab) {
  dist_tab %>% 
    mutate(Flow_Factor=(Floorspace^1.05)*2^(-Dist/h))  %>% 
    group_by(BU_CODE) %>% 
    mutate(Cash_Flow = Money*Flow_Factor/sum(Flow_Factor)) -> res
  return(res)
}

# Try out a halflife of 750m
fit_sim(h = 0.75, dist_tab = dist_tab) %>% 
  slice_head()

Mapping the results

We now have a data frame sim_tab that tabulates the results of the Spatial Interaction Model. The next step is to investigate the geographical patterns associated with this. Suppose we are mainly concerned with the AH/Albert Heijn group of stores. One starting point might be to discover the total cash flows predicted into each of the stores - this is a two-stage process - firstly filter out these stores (based on text patterns in the shop name) - and then group by the stores and sum over them using the summarise function. The pipeline is as below:

sim_tab %>% 
  filter(str_detect(shop, pattern = "AH|Heijn")) %>% 
  group_by(shopcode) %>%
  summarise(in_flow = sum(Cash_Flow)) %>% 
  arrange(desc(in_flow)) -> inflow_tab

inflow_tab

For now, it is enough to note that ‘AH|Heijn’ picks out shop names with either ‘AH’ or ‘Heijn’ in them. However, although this tabulates the influx of cash, it doesn’t show a map. To get this, firstly join this table to the gron_shops geographical object and then use tmap to create a map.

inflow_tab %>% 
  left_join(gron_shops, by = "shopcode") %>%
  # To make shop name appear on the hover, seems like the first column appears?
  relocate(shop) %>% 
  st_sf() %>% 
  tm_shape() +
  tm_dots(size = "in_flow", col = "darkorchid", alpha = 0.5)

## Legend for symbol sizes not available in view mode.

This depends on the distance half life. Once again we can encapsulate all of this in a function with \(h\) as the argument, making it easier to explore:

map_inflow <- function(h, dist_tab){
  fit_sim(h, dist_tab) %>% 
  filter(str_detect(shop, pattern = "AH|Heijn")) %>% 
  group_by(shopcode) %>%
  summarise(in_flow = sum(Cash_Flow)) %>% 
  arrange(desc(in_flow)) %>% 
  left_join(gron_shops, by = "shopcode") %>%
  relocate(shop) %>% 
  st_sf() %>% 
  tm_shape() +
  tm_dots(size = "in_flow", col = "darkorchid", alpha = 0.5)
    
}

map_inflow(h = 0.5, dist_tab = dist_tab)

## Legend for symbol sizes not available in view mode.

The above shows the effect if \(h\) is 500m - now check it for 1km:

map_inflow(h = 1, dist_tab = dist_tab)

## Legend for symbol sizes not available in view mode.

In a nutshell, as the distance half life gets larger, people are willing to travel further to stores, and the more attractive larger stores get a greater share of the money - people will travel to a more favourable store even when it is a long way away.

• In fact, money going to smaller stores decreases with increase in half life

Catchment Areas

The idea of a catchment area here is to identify, for each administrative unit, what is the most attractive store. This is done by considering the cash flow out of each unit to each of the stores and seeing which store gets the greatest share of the total outflow from the unit. In R (and in tidyverse terms) this is done using the rank function. Rank takes a variable and for each row, returns the rank of the value for that row, out of all of the rows. The ranks are in ascending order, so for \(n\) items the largest has the value \(n\).

To find the ranks of the outflow shares on an administrative unit by adminstrative unit basis, group_by is a helper function again. In this case it causes the rank function to return the relative ranks of the flow factors (\(A_i^\alpha 2^{-d_{ij}/h}\)) within each group. Finally, it would by useful to have ranks in descending order - you are interested in the maximum in each group, and so rather than trying to find ranks of \(n\) when you don’t know what \(n\) is, it is more useful to find ranks of 1. Changing the sign of the variable (by using desc()) given to the rank function achieves this - essentially it flips the orders of the variables. Check the results here:

sim_tab %>% 
  group_by(BU_CODE) %>% 
  mutate(rank = row_number(desc(Flow_Factor))) %>% 
  filter(rank == 1) %>% 
  arrange(desc(Cash_Flow)) -> rank_tab

rank_tab %>% 
  slice_head()

There is now a unique row for each administrative unit, and the shop also appearing in that row is the one where most money flows - and this states which catchment area the administrative unit is in. To map these, firstly join this table to the groningen_polys object:

rank_tab %>% 
  left_join(gron_polys) %>% 
  relocate(BU_CODE) %>%
  st_sf() %>% 
  tm_shape() +
  tm_polygons(col = "shopcode")

## Joining, by = c("BU_CODE", "AANT_INW", "Money")

We can further verify this by looking at the proportion which each shop serves:

rank_tab %>% 
  left_join(gron_polys) %>% 
  group_by(shopcode) %>% 
  summarise(n = n()) %>% 
  ungroup() %>% 
  arrange(desc(n)) 

## Joining, by = c("BU_CODE", "AANT_INW", "Money")

A further useful technique would be to create dissolved catchment polygons. This essentially involves removing the boundaries between administrative unit polygons assigned to the same store. This is done by applying a group_by operation to a simple features object and then a summarise - this merges all of the polygons in the same group. Here there is no summary variable, we just do this to make the merge:

rank_tab %>% 
  left_join(gron_polys) %>% 
  st_sf() %>% 
  group_by(shopcode) %>% 
  summarise() %>% 
  tm_shape() +
  tm_polygons(col = "shopcode")

## Joining, by = c("BU_CODE", "AANT_INW", "Money")

As before - all of this could be made into a small number of larger pipelines and encapsulated in a function of \(h\):

map_catchments <- function(h, dist_tab) {
  fit_sim(h, dist_tab) %>%
    filter(str_detect(shop, pattern = "AH|Heijn")) %>% 
    group_by(BU_CODE) %>% 
    mutate(rank = row_number(desc(Flow_Factor))) %>% 
    filter(rank == 1) -> rank_tab
  rank_tab %>% 
  left_join(gron_polys) %>% 
  st_sf() %>% 
  group_by(shopcode) %>% 
  summarise() %>% 
  tm_shape() +
  tm_polygons(col = "shopcode")
}

map_catchments(0.1, dist_tab)

## Joining, by = c("BU_CODE", "AANT_INW", "Money")

Now you can again experiment with different \(h\) values - for example 750m:

map_catchments(0.75, dist_tab)

## Joining, by = c("BU_CODE", "AANT_INW", "Money")

Assessing the impact of a proposed store.

Suppose now a new store was proposed with a floor space of 5000 square meters. Let’s see where it would range:

gron_shops %>% 
  slice_max(order_by = Floorspace, n = 5)

Therefore a fairly large store!

How could the effect of this be allowed for? The main way to see the impact is to add the new store to the gron_shops spatial layer, and to re-run the model with this extra store. This can then be compared to the original model.

To see this, firstly download a new spatial data object proposed:

#download.file("https://www.dropbox.com/s/qjzxl0cr7ehtrti/proposed.zip?dl=1","proposed.zip")
unzip("proposed.zip")
proposed_shop <- st_read(".","proposed")

## Reading layer `proposed' from data source 
##   `C:\Users\ADMIN\Desktop\AI Kenya\Compressed\Internship_prep' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 235494.2 ymin: 580729.9 xmax: 235494.2 ymax: 580729.9
## Projected CRS: Double_Stereographic

Then, have a look at its location - the proposed store is in red.

gron_shops %>% 
  tm_shape() +
  tm_dots(size = "Floorspace", col = "darkorchid", alpha = 0.5) +
  tm_shape(proposed_shop) +
  tm_dots(size = "Floorspace", col = "indianred", alpha = 0.7)

## Legend for symbol sizes not available in view mode.

To re-run the model, again using a function is useful. Below is a function that encapsulates all of the steps in producing a model, from the provision of a shops object and an administrative area object as well as \(h\). The function looks intimidating, but really it is just a combination of the operations done above - the data munging and model fitting.

new_model <- function(admin,shops,h) {
  # Get the distance
  dmat <- get_dist(shops,admin)
  
  # Pivot longer the data
  as_tibble(shops) %>% 
  select(shopcode) %>% 
  bind_cols(as_tibble(dmat)) %>% 
  pivot_longer(!shopcode, names_to = "BU_CODE", values_to = "Dist") %>% 
  left_join(shops, by = "shopcode") %>% 
  left_join(admin, by = "BU_CODE") %>% 
  select(-starts_with("geometry")) -> dist_tab
  
  # Find flow factor
  dist_tab %>% 
  mutate(Flow_Factor=(Floorspace^1.05)*2^(-Dist/h)) %>% 
  group_by(BU_CODE) %>% 
  mutate(Cash_Flow = Money*Flow_Factor/sum(Flow_Factor)) -> res
  
  
  return(res) 
}

test <- new_model(gron_pts,gron_shops,0.3)
head(test)

Here the function takes the inputs and returns a data frame with the fitted model values, in the format of sim_tab used earlier.

To re-fit thew model with a new store, firstly create a new scenario by combining the new proposed store with the existing ones in Groningen:

all_shops <- proposed_shop %>% 
  bind_rows(gron_shops)

all_shops %>% 
  slice_head(n = 5)

Then update the simulation table:

new_sim_tab <- new_model(gron_pts, all_shops, 0.2)

new_sim_tab %>% 
  slice_head()

We can now compute the inflows to all of the stores.

new_sim_tab %>% 
  filter(str_detect(shop, "AH|Heijn|New")) %>% 
  group_by(shopcode) %>% 
  summarise(in_flow = sum(Cash_Flow)) %>% 
  arrange(desc(in_flow)) -> new_inflow_tab


new_inflow_tab %>% 
  slice_head(n = 5)

and then map the inflows:

new_inflow_tab %>% 
  left_join(all_shops, by = "shopcode") %>% 
  st_sf() %>% 
  tm_shape() +
  tm_dots(size = "in_flow", col = "darkgreen", alpha = 0.5)

## Legend for symbol sizes not available in view mode.

If we join the original inflow_tab to the new_inflow_tab using the shop code as a merge key we obtain a table with new and old cash inflow estimates, respectively called inflow.x and inflow.y.

The proposed store, and so doesn’t have an original inflow. Here it makes sense to set this as zero:

new_inflow_tab %>% 
  left_join(inflow_tab, by = "shopcode") %>%
  # Replace a 0 anywhere there may be an NA
  mutate(across(everything(), ~replace_na(.x, 0))) %>% 
  # Compute difference
  mutate(Diff = (in_flow.x - in_flow.y)) -> diff_inflow_tab


# Which stores have been impacted most by the new store
diff_inflow_tab %>% 
  slice_max(order_by = abs(Diff), n = 5)

Finally to map this as proportional symbols, since some values are positive and others are negative, create an absolute value of change column - and then map this

diff_inflow_tab %>% 
  mutate(abs_diff = abs(Diff)) %>% 
  left_join(all_shops, by = "shopcode") %>% 
  relocate(shop) %>% 
  st_sf() %>% 
  tm_shape() +
  tm_dots(size = "abs_diff", col = "Diff")

## Variable(s) "Diff" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

## Legend for symbol sizes not available in view mode.

This helps to identify the stores that will experience the most impact if a new store is opened as proposed.

Resources

• Making maps with R

• Geocomputation with R

• Elegant and informative maps with tmap