Internship SIM Prep: Modelling Population Flows Using Spatial Interaction Models

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

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

Based on Adam Dennet’s Australian Population Studies

Modernized version of Adam’s accompanying : https://rpubs.com/adam_dennett/376877

Setting up some spatial data

As the name suggests, to run a spatial interaction model, you are going to need some spatial data and some data on interactions (flows). Let’s start with some spatial data:

• EPSG:4283 - Australia coordinates in degrees

library(sf)
# I/O for GeoJSON
library(geojsonio)

# Read geojson file as spatial object
aus <- geojson_read("https://www.dropbox.com/s/0fg80nzcxcsybii/GCCSA_2016_AUST_New.geojson?raw=1", what = "sp")




# Read data as simple features object and set coordinate reference system
aus_sf <- aus %>% 
  st_as_sf() %>% 
  st_set_crs(4283)

aus_sf %>% 
  head()

• sf::as() - Methods to coerce sf to sp –>Spatial* and Spatial*DataFrame objects

# Now you may have noticed that the code order is a bit weird, so let's fix that and reorder
aus_sf1 <- aus_sf %>% 
  arrange(GCCSA_CODE)

# Now let's create an 'sp' object from our new ordered SF object
aus <- as_Spatial(aus_sf1)

Check your boundaries have downloaded ok

• sf::st_make_valid: makes an invalid geometry valid – Required probably due to updates in sf packages. See: https://giters.com/RamiKrispin/coronavirus/issues/95

library(tmap)
tmap_mode("view")
#qtm(aus_sf)

aus_sf %>% 
  st_make_valid() %>% 
  tm_shape() +
  tm_borders(col = "black", lwd = 2, alpha = 0.4) +
  tm_fill(col = "gold", alpha = 0.1)

Calculating a distance matrix

In our spatial interaction model, space is one of the key predictor variables. In this example we will use a very simple Euclidean distance measure between the centroids of the Greater Capital City Statistical Areas as our measure of space.

Now, with some areas so huge, there are obvious potential issues with this (for example we could use the average distance to larger settlements in the noncity areas), however as this is just an example, we will proceed with a simple solution for now.

• EPSG:3112 – Projected Australia coordinates in metres

• sf::st_transform - transform coordinate in sf

library(sp)
# Use the spDists function to create a distance matrix
# First reproject into a projected (metres) coordinate system
aus_proj <- spTransform(aus, "+init=epsg:3112")

summary(aus_proj)

## Object of class SpatialPolygonsDataFrame
## Coordinates:
##        min      max
## x -2083066  2346598
## y -4973093 -1115948
## Is projected: TRUE 
## proj4string :
## [+proj=lcc +lat_0=0 +lon_0=134 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0
## +ellps=GRS80 +units=m +no_defs]
## Data attributes:
##   GCCSA_CODE         GCC_CODE16         GCCSA_NAME         STATE_CODE       
##  Length:15          Length:15          Length:15          Length:15         
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##   STATE_NAME          AREA_SQKM      
##  Length:15          Min.   :   1695  
##  Class :character   1st Qu.:   4838  
##  Mode  :character   Median :  15842  
##                     Mean   : 512525  
##                     3rd Qu.: 884729  
##                     Max.   :2520230

Now let’s calculate the distances

dist <- spDists(aus_proj) %>%
  as_tibble() %>% 
  mutate(n = row_number())

dist %>% 
  slice_head(n = 5)

# Probably do the same in sf?

#aus_projection_sf = st_transform(aus_sf, "+init=epsg:3112")
#sf_dist = st_distance(aus_projection_sf)

# Pivot longer the data and make distance km
dist <- dist %>% 
  pivot_longer(!n, names_to = "var2", values_to = "dist") %>%
  mutate(dist = dist/1000)

dist %>% 
  slice_head(n = 5)

These distances are notionally in km - although you may notice that they are not 100% accurate. This is not a big problem for now as this is just an example, but for real applications, more accurate distances may be used.

Flow Data

The data we are going to use to test our spatial interaction models with is migration data from the 2011 Australian Census. The Australian Census has usual address indicator on Census night (UAICP) and address one year ago and 5 years ago indicators. From these, one year and 5 year migration transitions can be recorded - here we will use the 5 year transitions.

As well as flow data, there are additional data on unemployment rates, weekly income and the percentage of people living in rented accommodation for each origin and destination. We will use these as destination attractiveness / mass term and origin emissiveness / mass term proxies in the models which follow.

These data can be read straight into R with the following command:

#read in your Australian Migration Data
mdata <- read_csv("https://www.dropbox.com/s/wi3zxlq5pff1yda/AusMig2011.csv?raw=1", show_col_types = FALSE)

mdata %>% 
  slice_head(n = 5)

Now to finish, we need to add in our distance data that we generated earlier and create a new column of total flows which excludes flows that occur within areas (we could keep the within-area (intra-area) flows in, but they can cause problems so for now we will just exclude them).

# First create a new total column which excludes intra-zone totals
# Sets them to a very small number for reasons we will see later
mdata <- mdata %>% 
  mutate(
    flow_no_intra = case_when(
      Orig_code == Dest_code ~ 0,
      TRUE ~ Flow),
    
    offset = case_when(
      Orig_code == Dest_code ~ 1e-10,
      TRUE ~ 1))

mdata %>% 
  slice_head(n = 5)

Now we ordered our spatial data earlier so that our zones are in their code order. We can now easily join these data together with our flow data as they are in the correct order.

# and while we are here, rather than setting the intra zonal distances to 0
# we should set them to something small (most intrazonal moves won't occur over 0 distance)

mdata <- mdata %>% 
  bind_cols(dist %>% select(dist)) %>% 
  mutate(dist = case_when(
    dist == 0 ~ 5,
    TRUE ~ dist))

mdata %>% 
  slice_head(n = 5)

And this is what those flows look like on a map - quick and dirty style… Although first we’ll remove the intra-zonal flows.

# remove intra-zonal flows
mdatasub <- mdata %>% 
  filter(Orig_code != Dest_code)

Now let’s create a flow-line object and weight the lines according to the flow volumes

• Work is needed to convert the OD data into ‘desire lines’. Desire lines are straight lines between the origin and destination and represent where people would go if they were not constrained by the route network (see Figure 3 from (Lovelace et al. 2017)).

• Origin, destination,flow columns

# use the od2line function from Robin Lovelace's excellent stplanr package - remove all but the origin, destination and flow columns
library(stplanr)
mdatasub_skinny <- mdatasub %>% 
  select(Orig_code, Dest_code, Flow)
# Desire lines
travel_network <- od2line(flow = mdatasub_skinny, zones = aus)

# Converts flows to WGS84
travel_networkwgs <- spTransform(travel_network, "+init=epsg:4326")

# Set the line widths to some sensible values according to flow
w <- mdatasub_skinny$Flow / max(mdatasub_skinny$Flow) * 10

Make a map of the linestrings representing movement from one place to another.

library(leaflet)
# Plot in leaflet
leaflet() %>% 
  addTiles() %>% 
  addPolylines(
    data = travel_networkwgs,
    weight = w)

Or you can view your flows as a matrix…

• rowwise: allows you to compute on a data frame a row-at-a-time

• c_across() is designed to work with rowwise() to make it easy to perform row-wise aggregations.

# Pivot wider and find the total flow from a specific area
mdatasubmat <- mdata %>% 
  select(Orig_code, Dest_code, flow_no_intra) %>% 
  pivot_wider(names_from = Dest_code, values_from = flow_no_intra) %>%
  rowwise(Orig_code) %>% 
  mutate(total_flow_frm_area = sum(c_across(where(is.numeric))) ) 

mdatasubmat

Okay, we’ve set everything up, now it’s…

Modellin’ Time!!

Population flows can be conceptualised as interactions between two entities – origins and destinations – which have different properties of emissivity and attractiveness (see Lee 1966 for the classic paper on this topic in relation to migration).

In explaining how to run and calibrate spatial interaction models in R, I will adopt the notation used by Taylor Oshan in his excellent primer for running spatial interation models in Python. The paper is well worth a read and can be found here: http://openjournals.wu.ac.at/region/paper_175/175.html

Below is the classic multiplicative gravity model:

1. \[T_{ij} = k \frac{V_i^\mu W_j^\alpha}{d_{ij}^\beta}\]

This gravity model can be written in the form more familiar from Wilson’s 1971 paper - http://journals.sagepub.com/doi/abs/10.1068/a030001

\(T_{ij} = k V_i^\mu W_j^\alpha d_{ij}^-\beta\)

This model just says that the flows between an origin and destination are proportional to the product of the mass of the origin and destination and inversely proportional to the distance between them.

As origin and destination masses increase, flows increase, but as distance increases, flows decrease, and vice versa.

• where TijTij is the transition or flow, TT, between origin ii (always the rows in a matrix) and destination jj (always the columns in a matrix). If you are not overly familiar with matrix notation, the ii and jj are just generic indexes to allow us to refer to any cell in the matrix more generally.

• VV is a vector (a 1 dimensional matrix - or, if you like, a single line of numbers) of origin attributes which relate to the emissiveness of all origins in the dataset, ii - in our sample dataset, we have a vector of origin populations (which I have called vi1_origpop) and a vector of origin average salaries (which I have called vi2_origsal) in 2001

• WW is a vector of desination of attributes relating to the attractivenss of all destinations in the dataset, jj - in our sample dataset, we have a vector of destination populations (which I have called wj1_destpop) and a vector of destination average salaries (which I have called wj2_destsal) in 2001

• dd is a matrix of costs relating to the flows between ii and jj - in our case the cost is distance and it is called ‘dist’ in our dataset.

• kk, μμ, αα and ββ are all model parameters to be estimated

kk is a constant of proportionality and leads to this particular model being more accurately described as a ‘total constrained’ model as all flows estimated by the model will sum to any observed flow data used to calibrate the parameters, where:

3. \[k = \frac{T}{\sum_i \sum_jV_i^\mu W_j^\alpha d_{ij}^-\beta }\]

and TT is the sum of our matrix of observed flows or:

4. \[T = \sum_i \sum_j T_{ij}\]

In plain language, this is just the sum of all observed flows divided by the sum of all of the other elements in the model.

Estimating Model Parameters

Now, it’s perfectly possible to produce some flow estimates by plugging some arbitrary or expected estimated values into our parameters. The parameters relate to the scaling effect / importance of the variables they are associated with.

Most simply, where the effects of origin and destination attributes on flows scale in a linear fashion (i.e. for a 1 unit increase in, say, population at origin, we might expect a 1 unit increase in flows of people from that origin, or for a halving in average salary at destination, we might expect a halving of commuters), μμ = 1 and αα = 1.

In Newton’s original gravity equation, ββ = -2 where the influence of distance on flows follows a power law - i.e. for a 1 unit increase in distance, we have a 1^-2 (1) unit decrease in flows, for a 2 unit increase in distance, we have 2^-2 (0.25 or 1/4) of the flows, for a 3 unit increase, 3^-2 (0.111) etc.

Let’s see if these parameters are a fair first guess

1. \[ d_{ij}^\beta \ with\ \beta=-2 \]

#First plot the commuter flows against distance and then fit a model line with a ^-2 parameter
theme_set(theme_light())
library(scales)

mdata %>% 
  ggplot(mapping = aes(x = dist, y = Flow)) +
  geom_point() +
  geom_function(fun = ~.x^-2, color = "red", lwd = 1) +
  scale_y_continuous(name = "Migration Flow", labels = comma)

2. \[ V_i^\mu \]

mdata %>% 
  ggplot(mapping = aes(x = vi1_origpop, y = Flow)) +
  geom_point() +
  geom_function(fun = ~.x^1, color = "red", lwd = 1) +
  scale_x_continuous(name = "Origin Population", labels = comma) +
  scale_y_continuous(name = "Migration Flow", labels = comma)

3. \[ W_j^\alpha \]

mdata %>% 
  ggplot(mapping = aes(x = wj3_destmedinc, y = Flow)) +
  geom_point() +
  geom_function(fun = ~.x^1, color = "red", lwd = 1) +
  scale_x_continuous(name = "Destination Median income", labels = comma) +
  scale_y_continuous(name = "Migration Flow", labels = comma)

OK, so it looks like we’re not far off (well, destination income doesn’t look too promising as a predictor, but we’ll see how we get on…), so let’s see what flow estimates with these starting parameters look like.

# Set up some variables to hold our parameter values in:
mu <- 1
alpha <- 1
beta <- -2
k <- 1
T2 <- sum(mdatasub %>% pull(Flow))

Now let’s create some flow estimates using Equation 2 (multiplicative gravity mode) above… Begin by applying the parameters to the variables:

# Emmisivity, Atrractiveness, Distance
vi1_mu <- (mdatasub %>% pull(vi1_origpop))^mu
wj3_alpha <- (mdatasub %>% pull(wj3_destmedinc))^alpha
dist_beta <- (mdatasub %>% pull(dist))^beta

T1 <- vi1_mu*wj3_alpha*dist_beta
k <- T2/sum(T1)
k

## [1] 3.112867

Then, just as in Equation 2 above, just multiply everything together to get your flow estimates:

• mdatasub : data without intra flows

# Run the model and store your flow estimates in a new column
mdatasub <- mdatasub %>% 
  mutate(unconstrainedEst1 = round(k*vi1_mu*wj3_alpha*dist_beta, 0))

# Check that the sum of these estimates makes sense
sum(mdatasub$unconstrainedEst1)

## [1] 1313520

Pivot wider the data and have a look at your handy work

mdatasubmat1 <- mdatasub %>% 
  select(Orig_code, Dest_code, unconstrainedEst1) %>% 
  pivot_wider(names_from = Dest_code, values_from = unconstrainedEst1, values_fill = 0) %>%
   relocate("1GSYD", .after = Orig_code) %>% 
  rowwise(Orig_code) %>% 
  mutate(total_flow_frm_area = sum(c_across(where(is.numeric))) ) 

mdatasubmat1

How does the flow compare with the original?

mdatasubmat

How good is my model?

So, looking at the two little matrices above you can see that in some cases the flow estimates aren’t too bad, but in others they are pretty rubbish. Whilst it’s OK to eyeball small flow matrices like this, when you have much larger matrices, we need another solution…

Testing the “goodness-of-fit”.

Yes, that’s what it’s called - I know, it doesn’t sound correct, but goodness-of-fit is the correct term for checking how well your model estimates match up with your observed flows.

So how do we do it?

Well… there are a number of ways but perhaps the two most common are to look at the coefficient of determination (r2r2) or the Square Root of Mean Squared Error (RMSE). You’ve probably come across r2r2 before if you have fitted a linear regression model, but you may not have come across RMSE. There are other methods and they all do more or less the same thing, which is essentially to compare the modelled estimates with the real data. r2r2 is popular as it is quite intuitive and can be compared across models. RMSE is less intuitive, but some argue is better for comparing changes to the same model. Here’s we’ll do both…

R-Squared

\(r^2\) is the square of the correlation coefficient, \(r\)

For our sample data, we can calculate this very easily using a little function

rsqrd <- function(truth, estimate){
  r = cor(truth, estimate)
  R2 = r^2
  return(R2)
}

rsqrd(truth = mdatasub$Flow, estimate = mdatasub$unconstrainedEst1)

## [1] 0.1953717

Using this function we get a value of 0.195 or around 20%. This tells us that our model accounts for about 20% of the variation of flows in the system. Not bad, but not brilliant either.

Root Mean Squared Error (RMSE)

We can use a similar simple function to calculate the RMSE for our data

rMSE <- function(truth, estimate){
  res = (truth - estimate)^2
 RMSE = sqrt(mean(res))
  return(round(RMSE, 3))
}

rMSE(truth = mdatasub$Flow, estimate = mdatasub$unconstrainedEst1)

## [1] 25858.11

This can be loosely interpreted as, on average, the flows are off by 25, 858.

The figure that is produced by the RMSE calculation is far less intuitive than the r2r2 value and this is mainly because it very much depends on things like the units the data are in and the volume of data. It can’t be used to compare different models run using different data sets. However, it is good for assessing whether changes to the model result in improvements. The closer to 0 the RMSE value, the better the model.

So how can we start to improve our fit…?

Improving our model: 1 - Calibrating parameters

Now, the model we have run above is probably the most simple spatial interaction model we could have run and the results aren’t terrible, but they’re not great either.

One way that we can improve the fit of the model is by calibrating the parameters on the flow data that we have.

The traditional way that this has been done computationally is by using the goodness-of-fit statistics. If you have the requisite programming skills, you can write a computer algorithm that iteratively adjusts each parameter, runs the model, checks the goodness-of-fit and then starts all over again until the goodness-of-fit statistic is maximised.

This is partly why spatial interaction modelling was the preserve of specialists for so long as acquiring the requisite skills to write such computer programmes can be challenging!

However, since the early days of spatial interaction modelling, a number of useful developments have occurred… For a more detailed explanation, read the accompanying paper, but I will skate over them again here.

The mathematically minded among you may have noticed that if you take the logarithms of both sides of Equation 2, you end up with the following equation:

5. \(\ln T_{ij} = k + \mu\ln V_i + \alpha\ln W_j - \beta \ln d_{ij}\)

Those of you who have played around with regression models in the past will realise that this is exactly that - a regression model.

And if you have played around with regression models you will be aware that there are various pieces of software available to run regressions (such as R) and calibrate the parameters for us, so we don’t have to be expert programmers to do this - yay!

Now, there are a couple of papers that are worth reading at this point. Perhaps the best is by Flowerdew and Aitkin (1982), titled “A METHOD OF FITTING THE GRAVITY MODEL BASED ON THE POISSON DISTRIBUTION” - the paper can be found here: http://onlinelibrary.wiley.com/doi/10.1111/j.1467-9787.1982.tb00744.x/abstract

One of the key points that Flowerdew and Aitkin make is that the model in Equation 5 (known as a log-normal model) has various problems associated with it which mean that the estimates produced might not be reliable. If you’d like to know more about these, read the paper (and also Wilson’s 1971 paper), but at this point it is worth just knowing that the way around many of these issues is to re-specify the model, not as a log-normal regression, but as a Poisson or negative binomial regression model.

Poisson regression

Again, I go into this in more detail in the accompanying paper, but the main theory (for non-experts like me anyway) behind the Poisson regression model is that the sorts of flows that spatial interaction models deal with (such as migration or commuting flows) relate to non-negative integer counts (you can’t have negative people moving between places and you can’t - normally, if they are alive - have fractions of people moving either).

As such, the continuous (normal) probabilty distributions which underpin standard regression models don’t hold. However, the discrtete probability distributions such as the Poisson distribution and the negative binomial distribution (of which the Poisson distribution is a special case - wikipedia it) do hold and so we can use these associations to model our flows.

At this point, it’s probably worth you looking at what a Poisson disribution looks like compared to a normal distribution, if you are not familiar.

Here’s a normal distribution:

# histogram with normal distribution of count 3000, a mean of 75 and a standard deviation of 5
tibble(ndist = rnorm(3000, mean = 10, sd = 5)) %>% 
  ggplot() +
  geom_histogram(aes(x = ndist), fill = "midnightblue", alpha = 0.7) +
  xlab("rnorm(3000, mean = 75, sd = 5)")

Now here’s a Poisson distribution with the same mean:

# histogram with a poisson distribution of count 3000, a mean of 75 
tibble(pdist = rpois(n = 3000, lambda = 75)) %>% 
  ggplot() +
  geom_histogram(aes(x = pdist), fill = "indianred", alpha = 0.7) +
  xlab("rpois(n = 3000, lambda = 75)")

Looks kind of similar doesn’t it! The thing with the Poisson distribution is, when the mean (λ - lambda) changes, so does the distribution. Normal distributions on the other hand retain their bell shape.

As the mean gets smaller (and this is often the case with flow data where small flows are very likely - have a look at the ‘Total’ column in your cdata dataframe, lots of small numbers aren’t there?) the distribution starts to look a lot more like a skewed or log-normal distrbution. They key thing is it’s not - it’s a Poisson distribution. Here’s a similar frequency distribution with a small mean:

# What about a lambda of 0.5?
tibble(pdist = rpois(n = 3000, lambda = 0.5)) %>% 
  ggplot() +
  geom_histogram(aes(x = pdist), fill = "indianred", alpha = 0.7, binwidth = 1) +
  xlab("rpois(n = 3000, lambda = 75)")

As far as we’re concerned, what this means is that if we are interested in all flows between all origins and destinations in our system, these flows will have a mean value of \(\lambda_{ij}\) and this will dictate the distribution. Here’s what the distrbution of our flows looks like:

# Distribution of flows
mdatasub %>% 
  ggplot(mapping = aes(x = Flow)) +
  geom_histogram(fill = "darkorange", alpha = 0.8)

Reveals a histogram which looks like a skewed normal or, more accurately, a Poisson distribution with a small mean.

So, what does all of this mean for our spatial interaction model?

Well the main thing it means is that Equation 5, for most sorts of spatial interaction models where we are modelling flows of people or whole things, is not correct.

By logging both sides of the equation in Equation 5, we are trying to get a situation where our \(T_{ij}\) flows can be modelled by using the values of our other variables such as distance, by using a straight line a bit like this:

mdatasub %>% 
  ggplot(mapping = aes(x = dist, y = log(Flow))) +
  geom_point() +
  geom_smooth(method = lm)

If you compare this graph with the graph above (the first scatter plot we drew in this practical exercise), it’s exactly the same data, but clearly by logging both the total and distance, we can get a bit closer to being able to fit a model estimate using a straight line.

What the Poisson distribution means is that the yy variable in our model is not logged as in the graph above, but it can still be modelled using something like the blue line - I hope that sort of makes sense. If not, don’t worry, just take it from me that this is good news.

The Poisson Regression Spatial Interaction Model

So, we can now re-specify Equation 5 as a Poisson Regression model. Instead of our independent variable being lnTij our dependent variable is now the mean of our Poisson distribution λij and the model becomes:

6. \(\lambda_{ij} = \exp(k + \mu\ln V_i + \alpha\ln W_j - \beta \ln d_{ij})\)

What this model says is \(\lambda_{ij}\) (our independent variable - the estimate of \(T_{ij}\)) is logarithmically linked to (or modelled by) a linear combination of the logged independent variables in the model.

Now we have Equation 6 at our disposal, we can use a Poisson regression model to produce estimates of \(k\), \(\mu\), \(\alpha\) and \(\beta\) - or put another way, we can use the regression model to calibrate our parameters.

So, let’s have a go at doing it!!

It’s very straight forward to run a Poisson regression model in R using the glm (Generalised Linear Models) function. In practical terms, running a GLM model is no different to running a standard regression model using lm. If you want to find out more about glm, try the R help system ?glm or google it to find the function details. If you delve far enough into the depths of what GLM does, you will find that the parameters are calibrated though an ‘iteratively re-weighted least squares’ algorithm. This algorithm does exaxtly the sort of job I described earlier, it fits lots of lines to the data, continually adjusting the parameters and then seeing if it can minimise the error between the observed and expected values useing some goodness-of-fit measure is maximised/minimised.

These sorts of algorithms have been around for years and are very well established so it makes sense to make use of them rather than trying to re-invent the wheel ourselves. So here we go…

• We are using glm to estimate the values of mu, alpha and beta

# Run the unconstrained model
uncon_sim <- glm(Flow ~ log(vi1_origpop) + log(wj3_destmedinc) + log(dist), family = poisson(link = "log"), na.action = na.exclude, data = mdatasub)

It’s a simple as that - runs in a matter of milliseconds. You should be able to see how the glm R code corresponds to Equation 6.

Total = \(T_{ij}\) = \(\lambda{ij}\)

~ means ‘is modelled by’

log(vi1_origpop) = \(ln{V_i}\)

log(wj2_destsal) = \(ln{W_j}\)

log(dist) = \(lnd_{ij}\)

family = poisson(link = "log") means that we are using a Poisson regression model (the link is always log with a Poisson model) where the left-hand side of the model equation is logarithmically linked to the variables on the right-hand side.

• Try and do this in Tidymodels

So what comes out of the other end?

Well, we can use the summary() function to have a look at the model parameters:

summary(uncon_sim)

## 
## Call:
## glm(formula = Flow ~ log(vi1_origpop) + log(wj3_destmedinc) + 
##     log(dist), family = poisson(link = "log"), data = mdatasub, 
##     na.action = na.exclude)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -177.78   -54.49   -24.50     9.21   470.11  
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          7.1953790  0.0248852  289.14   <2e-16 ***
## log(vi1_origpop)     0.5903363  0.0009232  639.42   <2e-16 ***
## log(wj3_destmedinc) -0.1671417  0.0033663  -49.65   <2e-16 ***
## log(dist)           -0.8119316  0.0010157 -799.41   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2750417  on 209  degrees of freedom
## Residual deviance: 1503573  on 206  degrees of freedom
## AIC: 1505580
## 
## Number of Fisher Scoring iterations: 5

We can see from the summary that the Poisson regression has calibrated all 4 parameters for us and these appear under the ‘estimate’ column:

$k$ (intercept) = 7.1953790

\(\mu\) = 0.5903363

\(\alpha\) = -0.1671417

and \(\beta\) = -0.8119316

We can also see from the other outputs that all variables are highly significant (***), with the z-scores revealing that distance has the most influence on the model (as we might have expected from the scatter plots we produced earlier which showed that distance had by far the strongest correlation with migration flows). — ?you take absolute values of z?

These parameters are not too far away from our initial guesses of μμ = 1, αα = 1 and ββ = -2, but how do the estimates compare?

One way to calculate the estimates is to plug all of the parameters back into Equation 6 like this:

mdatasub <- mdatasub %>% 
  mutate(unconstrainedEst2 = round(fitted(uncon_sim), 0)) 

sum(mdatasub$unconstrainedEst2)

## [1] 1313517

# Turn it into a little matrix and have a look at your handy work
mdatasubmat2 <- mdatasub %>% 
  select(Orig_code, Dest_code, unconstrainedEst2) %>% 
  pivot_wider(names_from = Dest_code, values_from = unconstrainedEst2, values_fill = 0) %>%
   relocate("1GSYD", .after = Orig_code) %>% 
  rowwise(Orig_code) %>% 
  mutate(total_flow_frm_area = sum(c_across(where(is.numeric))) ) 

mdatasubmat2

And the $1,000,000 question - has calibrating the parameters improved the model…?

# Evaluate R^2
rsqrd(truth = mdatasub$Flow, estimate = mdatasub$unconstrainedEst2)

## [1] 0.3245418

# Evaluate RMSE
rMSE(truth = mdatasub$Flow, estimate = mdatasub$unconstrainedEst2)

## [1] 10789.17

Yes indeedy do!!

The \(r^2\) has improved from 0.20 to 0.32 and the RMSE has reduced from 25858.11 to 10789.17 so by calibrating our parameters using the Poisson Regression Model, we have markedly improved our model fit.

But we can do even better. We have just been playing with the unconstrained model, by adding constraints into the model we can both improve our fit further AND start to do cool things like esimate transport trip distributions from know information about people leaving an area, or in different contexts estimate the amount of money a shop is going to make from the available money that people in the surrounding area have to spend, or guess the number of migrants travelling between specific countries where we only know how many people in total leave one country and arrive in another.

Section 2 - Constrained Models

If we return to Alan Wilson’s 1971 paper, he introduces a full family of spatial interaction models of which the unconstrained model is just the start. And indeed since then, there have been all number of incremental advances and alternatives (such as Stewart Fotheringham’s Competing Destinations models, Pooler’s production/attraction/cost relaxed models, Stillwell’s origin/destination parameter specific models and mine and Alan’s own multi-level model (to name just a few).

In this section we will explore the rest of Wilson’s family - the Production (origin) Constrained Model; the Attraction (destination) constrained model; and the Doubly Constrained Model.

We will see how we can, again, use a Poisson regression model in R to calibrate these models and how, once calibrated, we can use the models in different contexts, such as Land Use Transportation Interaction (LUTI) modelling, retail modelling and migration modelling.

1. Production and Attraction Constrained Models

Wilson’s real contribution to the field was in noticing that the unconstrained gravity model was sub-optimal as it did not make use of all of the available data in the system we are studying.

If we recall the estimates from our unconstrained model, none of the estimates summed to the observed in and out-flow totals:

Our estimates did sum to the grand total of flows, but this is because we were really fitting a ‘total constrained’ model which used kk - our constant of proportionality - to ensure everything sort of added up (to within 1 commuter).

Where we have a full flow matrix to calibrate parameters, we can incorporate the row (origin) totals, column (destination) totals or both origina and destination totals to constrain our flow estimates to these known values.

As I outline in the accompanying paper, there are various reasons for wanting to do this, for example:

1. If We are interested in flows of money into businesses or customers into shops, might have information on the amount of disposable income and shopping habits of the people living in different areas from loyalty card data. This is known information about our origins and so we could constrain our spatial interaction model to this known information - we can make the assumption that this level of disposable income remains the same. We can then use other information about the attractiveness of places these people might like to shop in (store size, variety / specialism of goods etc.), to estimate how much money a new store opening in the area might make, or if a new out-of-town shopping centre opens, how much it might affect the business of shops in the town centre. This is what is known in the literature as the ‘retail model’ and is perhaps the most common example of a Production (orign) Constrained Spatial Interaction Model

2. We might be interested in understanding the impact of a large new employer in an area on the flows of traffic in the vicinity or on the demand for new worker accommodation nearby. A good example of where this might be the case is with large new infrastructure developments like new airports. For example, before the go-ahead for the new third runway at Heathrow was given, one option being considered was a new runway in the Thames Estuary. If a new airport was built here, what would be the potential impact on transport flows in the area and where might workers commute from? This sort of scenario could be tested with an Attraction (destination) Constrained Spatial Interaction Model where the number of new jobs in a destination is known (as well as jobs in the surrounding area) and the model could be used to estimate where the workers will be drawn from (and their likely travel-to-work patterns). This model is exactly the sort of model Land Use Transport Interaction (LUTI) model that was constructed by the Mechanicity Team in CASA - details here if you are interested…

3. We might be interested in understanding the changing patterns of commuting or migration over time. Data from the Census allows us to know an accurate snap-shot of migrating and commuting patterns every 10 years. In these full data matrices, we know both the numbers of commuters/migrants leaving origins and arriving at destinations as well as the interactions between them. If we constrain our model estimates to this known information at origin and destination, we can examine various things, including:

   1. the ways that the patterns of commuting/migration differ from the model predictions - where we might get more migrant/commuter flows than we would expect

   2. how the model parameters vary over time - for example how does distance / cost of travel affect flows over time? Are people prepared to travel further or less far than before?

2. Production-constrained Model

1. \(T_{ij} = A_i O_i W_j^\alpha d_{ij}^-\beta\)

where

\(O_{i} = \sum_j T_{ij}\)

and

\(A_i = \frac{1}{\sum_j W_j^\alpha d_{ij}^-\beta}\)

In the production-constrained model, OiOi does not have a parameter as it is a known constraint. AiAi is known as a balancing factor and is a vector of values which relate to each origin, ii, which do the equivalent job to kk in the unconstrained/total constrained model but ensure that flow estimates from each origin sum to the know totals, OiOi rather than just the overall total.

Now at this point, we could calculate all of the OiOis and the AiAis by hand for our sample system and then set about guessing/estimating the parameter values for the rest of the model, but as you might have already suspected from last time, we can use R and glm to make it really easy and do all of that for us - woo hoo!

We set about re-specifying the Production-Constrained model as a Poisson regression model in exactly the same way as we did before. We need to take logs of the right-hand side of equation and assume that these are logarithmially linked to the Poisson distributed mean (λijλij) of the TijTij variable. As such, Equation (1) becomes:

4. \(\lambda_{ij} = exp(\mu_{i} + \alpha \ln W_j - \beta \ln d_{ij})\)

In Equation (4) μiμi is the equivalent of the vector of balancing factors AiAi, but in regression / log-linear modelling terminology can also be described as either dummy variables or fixed effects. In practical terms, what this means is that in our regression model, μiμi is modelled as a categorical predictor and therefore in the Poisson regression model, we don’t use the numeric values of OiOi, we use a categorical identifier for the origin. In terms of the origin/destination migration matrix shown in Table 3, rather than the flow of 204,828 migrants leaving Sydney (row 1) being used as a predictor, simply the code ‘1GSYD’ is used as a dummy variable.

Before giving it a whirl, it’s important to note in the code example below the use of ‘-1’ after the distance variable (thanks to Hadrien Salat in CASA for bringing this to my attention). The -1 serves the purpose of removing the intercept that by default, GLM will insert into the model. As was mentioned earlier, the vector of origin parameters will replace the intercept in this model. It also serves the purpose

# Run the production constrained SIM (the "-1" indicates no intercept in the regression model)
prodSim <- glm(Flow ~ Orig_code + log(wj3_destmedinc) + log(dist) - 1, na.action = na.exclude, family = poisson(link = "log"), data = mdatasub)

# Summary of model
summary(prodSim)

## 
## Call:
## glm(formula = Flow ~ Orig_code + log(wj3_destmedinc) + log(dist) - 
##     1, family = poisson(link = "log"), data = mdatasub, na.action = na.exclude)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -225.71   -54.10   -15.94    20.45   374.27  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## Orig_code1GSYD      19.541851   0.023767  822.22   <2e-16 ***
## Orig_code1RNSW      19.425497   0.023913  812.35   <2e-16 ***
## Orig_code2GMEL      18.875763   0.023243  812.12   <2e-16 ***
## Orig_code2RVIC      18.335242   0.022996  797.31   <2e-16 ***
## Orig_code3GBRI      19.856564   0.024063  825.20   <2e-16 ***
## Orig_code3RQLD      20.094898   0.024300  826.94   <2e-16 ***
## Orig_code4GADE      18.747938   0.023966  782.28   <2e-16 ***
## Orig_code4RSAU      18.324029   0.024407  750.75   <2e-16 ***
## Orig_code5GPER      20.010551   0.024631  812.43   <2e-16 ***
## Orig_code5RWAU      19.392751   0.024611  787.96   <2e-16 ***
## Orig_code6GHOB      16.802016   0.024282  691.97   <2e-16 ***
## Orig_code6RTAS      17.013981   0.023587  721.33   <2e-16 ***
## Orig_code7GDAR      18.607483   0.025012  743.93   <2e-16 ***
## Orig_code7RNTE      17.798856   0.025704  692.45   <2e-16 ***
## Orig_code8ACTE      17.796693   0.023895  744.79   <2e-16 ***
## log(wj3_destmedinc) -0.272640   0.003383  -80.59   <2e-16 ***
## log(dist)           -1.227679   0.001400 -876.71   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 23087017  on 210  degrees of freedom
## Residual deviance:  1207394  on 193  degrees of freedom
## AIC: 1209427
## 
## Number of Fisher Scoring iterations: 6

So, what do we have?

Well, there are the elements of the model output that should be familiar from the unconstrained model:

the α parameter related to the destination attractiveness: -0.272640

the β distance decay parameter: -1.227679

We can see from the standard outputs from the model that all of the explanatory variables are statistically significant (***) and the z-scores indicate that the destination salary is having the most influence on the model, with distance following closely behind. And then we have a series of parameters which are the vector of μi values associated with our origin constraints.

2.1 Model Estimates

Now at this point you will be wanting to know what affect the constraints have had on the estimates produced by the model, so let’s plug the parameters back into Equation 4 and take a look…

Create some Oi and Dj columns and store the total in and out flow matrix margins in them.

• \(O_i\) represents the total flows from a given origin

• \(D_j\) represents the total flows going to a particular destination

# Create some Oi and Dj columns in the dataframe and store row and column totals in them:
mdatasub <- mdatasub %>% 
  group_by(Orig_code) %>% 
  summarise(O_i = sum(Flow)) %>% 
  right_join(mdatasub) 

# Destination flows
mdatasub <-  mdatasub %>% 
  group_by(Dest_code) %>% 
  summarise(D_j = sum(Flow)) %>% 
  right_join(mdatasub) %>% 
  arrange(Orig_code) %>% 
  relocate(D_j, .after = O_i)

pull out the coefficient values for μi and store them back in the dataframe along with Oi and Dj

# Extract coefficients
alpha <- prodSim$coefficients[16]
beta <- prodSim$coefficients[17]

mdatasub <- prodSim %>% 
  pluck("coefficients") %>% 
  tibble(mu_i = .) %>% 
  slice(1:(n() - 2)) %>% 
  bind_cols(distinct(mdatasub, Orig_code)) %>% 
  right_join(mdatasub) %>% 
  relocate(mu_i, .after = D_j)

slice_head(mdatasub)

Now lets’s add the estimated flows by the Production constrained model

mdatasub <- mdatasub %>% 
  mutate(prodsimFitted = round(fitted(prodSim), 0))

2.2 Assessing the model output

So what do the outputs from our Production Constrained Model look like? How has the goodness-of-fit improved and how can we start to use this a bit like a retail model and assess the likely impacts of changing destination attractiveness etc.?

2.2.1 The flow matrix

Now we can create pivot table to turn paired list into matrix (and compute the margins as well)

mdatasubmat3 <- mdatasub %>% 
  select(Orig_code, Dest_code, prodsimFitted) %>% 
  pivot_wider(names_from = Dest_code, values_from = prodsimFitted, values_fill = 0) %>% 
  relocate("1GSYD", .after = Orig_code) %>% 
  rowwise(Orig_code) %>% 
  mutate(total_flow_frm_area = sum(c_across(where(is.numeric))) ) 

mdatasubmat3

And compared with the original observed data?

mdatasubmat

Here it is very easy to see the Origin Constraints working. The sum across all destinations for each origin in the estimated matrix is exactly the same as the same sum across the observed matrix - \(\sum_j T_{ij} = \sum_j \lambda_{ij} = O_i\), but clearly, the same is not true when you sum across all origins for each destination - \(\sum_i T_{ij} \neq \sum_i \lambda_{ij} \neq D_j\)

2.2.2 How do the fits compare with the unconstrained model from last time?

# Evaluate R^2
rsqrd(truth = mdatasub$Flow, estimate = mdatasub$prodsimFitted)

## [1] 0.4345011

# Evaluate RMSE 
rMSE(truth = mdatasub$Flow, estimate = mdatasub$prodsimFitted)

## [1] 9872.693

Clearly by constraining our model estimates to known origin totals, the fit of the model has improved quite considerably - from around 0.32 in the unconstrained model to around 0.43 in this model. The RMSE has also dropped quite noticably.

2.2.3 A ‘what if…’ scenario

Now that we have calibrated our parameters and produced some estimates, we can start to play around with some what-if scenarios.

For example - What if the government invested lots of money in Tasmainia and average weekly salaries increased from 540.45 to 800.50 dollars a week? A far fetched scenario, but one that could make a good experiment.

First create create a new variable with these altered salaries:

mdatasub <- mdatasub %>% 
  mutate(wj3_destmedincScenario = case_when(
    wj3_destmedinc == 540.45 ~ 800.50,
    TRUE ~ wj3_destmedinc)
  )

Now let’s plug these new values into the model and see how this changes the flows in the system… \(\lambda_{ij} = exp(\mu_{i} + \alpha \ln W_j - \beta \ln d_{ij})\)

mdatasub <- mdatasub %>% 
  mutate(prodsimest2 = exp(mu_i + (alpha*log(wj3_destmedincScenario) + (beta*log(dist)))) %>% round(., 0))

#now we can create pivot table to turn paired list into matrix (and compute the margins as well)
mdatasubmat4 <- mdatasub %>% 
  select(Orig_code, Dest_code, prodsimest2) %>% 
  pivot_wider(names_from = Dest_code, values_from = prodsimest2, values_fill = 0) %>% 
  relocate("1GSYD", .after = Orig_code) %>% 
  rowwise(Orig_code) %>% 
  mutate(total_flow_frm_area = sum(c_across(where(is.numeric))) )

mdatasubmat4

You will notice that by increasing the average salary in the rest of Tazmania, we’ve reduced the flows into this area (yes, I know, counterintuitively, but this is just an example), but have not reduced the flows into other zones - the original constraints are still working on the other zones. One way to get around this, now that we have calibrated our parameters, is to return to the multiplicative model in Equation 1 and run this model after calculating our own \(A_i\) balancing factors.

\(T_{ij} = A_i O_i W_j^\alpha d_{ij}^-\beta\)

# Calculate new wj^alpha and dij=j^beta
wj2_alpha <- mdatasub$wj3_destmedinc^alpha
dist_beta <- mdatasub$dist^beta

# Calculate the first stage  of the Ai values
mdatasub <- mdatasub %>% 
  mutate(Ai1 = wj2_alpha*dist_beta) %>% 
  group_by(Orig_code) %>% 
  # Find constraint of a given location
  summarize(A_i = sum(Ai1)) %>% 
  mutate(A_i = 1/A_i) %>%
  right_join(mdatasub)

So that is it for calculating your \(A_i\) values. Now that you have these, it’s very simple to plug everything back into Equation 1 and generate some estimates.

# To check everything works, recreate the original estimates
mdatasub <- mdatasub %>% 
  mutate(prodsimest3 = round(A_i*O_i*wj3_destmedinc^alpha*dist^beta, 0))

You should see that your new estimates are exactly the same as your first estimates. If they are not then something has gone wrong. Now we have this though, we can keep messing around with some new estimates and keep the constraints.

Remember, though, that you will need to recalculate \(A_i\) each time you want to create a new set of estimates. Let’s try with our new values for the destination salary in the rest of Tazmania:

wj3_alpha <- mdatasub$wj3_destmedincScenario^alpha
dist_beta <- mdatasub$dist^beta

# Calculate first stages of Ai values
mdatasub <- mdatasub %>% 
  mutate(Ai1 = wj3_alpha*dist_beta) %>% 
  group_by(Orig_code) %>% 
  # Find constraint of a given location
  summarize(A_i = sum(Ai1)) %>% 
  mutate(A_i = 1/A_i) %>%
  right_join(mdatasub %>% select(!starts_with("A")))

Now we have some new \(A_i\)s, let’s generate some new scenario flow estimates.

mdatasubmat5 <- mdatasub %>% 
  mutate(prodsimest4_scenario = round(A_i*O_i*wj3_alpha*dist_beta)) %>%
  select(Orig_code, Dest_code, prodsimest4_scenario) %>% 
  pivot_wider(names_from = "Dest_code", values_from = "prodsimest4_scenario", values_fill = 0) %>% 
  relocate("1GSYD", .after = Orig_code) %>% 
  rowwise(Orig_code) %>% 
  mutate(total_flow_frm_area = sum(c_across(where(is.numeric))))

There are a number of things to note here. Firstly, flows into Tazmania have reduced, while flows into other regions have increased. Secondly, yes, I know this was a bad example, but try with some of the other variables for yourself.

Thirdly, Our origin constraints are now holding again.

3. Attraction-Constrained Model

The attraction constrained Model is virtually the same as the Production constrained model:

5. \(T_{ij} = D_j B_j V_i^\mu d_{ij}^-\beta\)

where

6. \(D_{j} = \sum_i T_{ij}\)

and

7. \(B_j = \frac{1}{\sum_i V_i^\mu d_{ij}^-\beta}\)

I won’t dwell on the attraction constrained model, except to say that it can be run in R as you would expect:

8. \(\lambda_{ij} = exp(\mu \ln V_i + \alpha_{i} - \beta \ln d_{ij})\)

or in R:

attrSim <- glm(Flow ~ Dest_code + log(vi1_origpop) + log(dist) - 1, family = poisson(link = "log"), data = mdatasub, na.action = na.exclude)

summary(attrSim)

## 
## Call:
## glm(formula = Flow ~ Dest_code + log(vi1_origpop) + log(dist) - 
##     1, family = poisson(link = "log"), data = mdatasub, na.action = na.exclude)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -138.69   -33.38   -10.47    11.72   293.39  
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## Dest_code1GSYD    8.8262922  0.0176638   499.7   <2e-16 ***
## Dest_code1RNSW    9.1809447  0.0178316   514.9   <2e-16 ***
## Dest_code2GMEL    8.6716196  0.0170155   509.6   <2e-16 ***
## Dest_code2RVIC    8.0861927  0.0173840   465.1   <2e-16 ***
## Dest_code3GBRI    9.5462594  0.0183631   519.9   <2e-16 ***
## Dest_code3RQLD   10.1295722  0.0184672   548.5   <2e-16 ***
## Dest_code4GADE    8.3051406  0.0184018   451.3   <2e-16 ***
## Dest_code4RSAU    8.1438651  0.0188772   431.4   <2e-16 ***
## Dest_code5GPER    9.9664486  0.0190008   524.5   <2e-16 ***
## Dest_code5RWAU    9.3061908  0.0190006   489.8   <2e-16 ***
## Dest_code6GHOB    6.9737562  0.0186288   374.4   <2e-16 ***
## Dest_code6RTAS    7.1546249  0.0183673   389.5   <2e-16 ***
## Dest_code7GDAR    8.3972440  0.0199735   420.4   <2e-16 ***
## Dest_code7RNTE    7.4521232  0.0206128   361.5   <2e-16 ***
## Dest_code8ACTE    7.3585270  0.0181823   404.7   <2e-16 ***
## log(vi1_origpop)  0.5828662  0.0009556   610.0   <2e-16 ***
## log(dist)        -1.1820013  0.0015267  -774.2   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 23087017  on 210  degrees of freedom
## Residual deviance:   665984  on 193  degrees of freedom
## AIC: 668017
## 
## Number of Fisher Scoring iterations: 5

we can examine how the constraints hold for destinations this time:

# First round of the estimates
mdatasub %>% 
  mutate(attrsimFitted = round(fitted(attrSim))) -> mdatasub

# Pivot the table to turn paired list into matrix (and compute the margins as well)
mdatasubmat6 <- mdatasub %>% 
  select(Orig_code, Dest_code, attrsimFitted) %>% 
  pivot_wider(names_from = "Dest_code", values_from = "attrsimFitted", values_fill = 0) %>% 
  relocate("1GSYD", .after = Orig_code) %>% 
  rowwise(Orig_code) %>% 
  mutate(total_flow_frm_area = sum(c_across(where(is.numeric))))

mdatasubmat6

compared to …

mdatasubmat

and we can test the goodess of fit in exactly the same way as before:

# Evaluate R^2
rsqrd(truth = mdatasub$Flow, estimate = mdatasub$attrsimFitted)

## [1] 0.6550357

# Evaluate RMSE 
rMSE(truth = mdatasub$Flow, estimate = mdatasub$attrsimFitted)

## [1] 7714.627

OK, that’s where I’ll leave singly constrained models for now. There are, of course, plenty of things you could try out. For example:

• You could try mapping the coefficients or the residual values from the model to see if there is any patterning in either the over or under prediction of flows.

• You could try running your own version of a LUTI model by first calibrating the model parameters and plugging these into a multiplicative version of the model, adjusting the destination constraints to see which origins are likely to generate more trips.

4. Doubly Constrained Model

Now, the model in the family you have all been waiting for - the big boss, the daddy, the doubly constrained model!

Let’s begin with the formula:

9. \(T_{ij} = A_i O_i B_j D_j d_{ij}^-\beta\)

where

10. \(O_{i} = \sum_j T_{ij}\)

11. \(D_{j} = \sum_i T_{ij}\)

and

12. \(A_i = \frac{1}{\sum_j B_j D_j d_{ij}^-\beta}\)

13. \(B_j = \frac{1}{\sum_i A_i O_i d_{ij}^-\beta}\)

Now, the astute will have noticed that the calculation of Ai relies on knowing Bj and the calculation of Bj relies on knowing Ai. A conundrum!! If I don’t know Ai how can I calcuate Bj and then in turn Ai and then Bj ad infinitum???!!

Well, I wrestled with that for a while until I came across this paper by Martyn Senior where he sketches out a very useful algorithm for iteratively arriving at values for Ai and Bj by setting each to equal 1 initially and then continuing to calculate each in turn until the difference between each value is small enough not to matter.

We will return to this later, but for now, we will once again use the awesome power of R to deal with all of this difficulty for us!

We can run the doubly constrained model in exactly the same way as we ran the singly constrained models:

15. \(\lambda_{ij} = exp(\mu_i + \alpha_{i} - \beta \ln d_{ij})\)

The code below has changed a litte from the singly constrained models I have removed the ‘-1’ which means that an intercept will appear in the model again. This is not because I want an intercept as it makes the origin and destination coefficients harder to interpret - reference categories zones will appear and the coefficients will need to be compared with the intercept - rather the ‘-1’ cheat for removing the intercept only works with one factor level - here we have two (origins and destinations). For full details and an explanation for alternative ways for dealing with this, please visit here - https://stats.stackexchange.com/questions/215779/removing-intercept-from-glm-for-multiple-factorial-predictors-only-works-for-fir - for ease, here we will just continue with the intercept.

The reference level means that the origin and destination coefficients need to be interpreted in relation to a reference category. In this example, the first zone in the system is used (Sydney), with the direction and size of the coefficients referring to whether another origin or destination zone has a greater or lesser positive or negative effect on migration flows in the system when compared to Sydney. The estimates produced by the Doubly Constrained Model are the most accurate in the Wilson family of models (in this example, an R2 value of 0.87). However, there is a loss of flexibility when compared to the singly constrained models, as only alternatives to origin/destination interaction explanatory variables such as historic flows or something other than distance can be experimented with. The double constraints mean that origin- and destination-specific explanatory variables cannot be used.

# Run a doubly constrained model
doubSim <- glm(Flow ~ Orig_code + Dest_code + log(dist), family = poisson(link = log), na.action = na.exclude, data = mdatasub)

# Summary of model
summary(doubSim)

## 
## Call:
## glm(formula = Flow ~ Orig_code + Dest_code + log(dist), family = poisson(link = log), 
##     data = mdatasub, na.action = na.exclude)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -93.018  -26.703    0.021   19.046  184.179  
## 
## Coefficients:
##                 Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)    20.208178   0.011308 1786.999   <2e-16 ***
## Orig_code1RNSW -0.122417   0.003463  -35.353   <2e-16 ***
## Orig_code2GMEL -0.455872   0.003741 -121.852   <2e-16 ***
## Orig_code2RVIC -1.434386   0.004511 -317.969   <2e-16 ***
## Orig_code3GBRI  0.241303   0.003597   67.091   <2e-16 ***
## Orig_code3RQLD  0.772753   0.003599  214.700   <2e-16 ***
## Orig_code4GADE -0.674261   0.004527 -148.936   <2e-16 ***
## Orig_code4RSAU -1.248974   0.005889 -212.091   <2e-16 ***
## Orig_code5GPER  0.742687   0.004668  159.118   <2e-16 ***
## Orig_code5RWAU -0.317806   0.005131  -61.943   <2e-16 ***
## Orig_code6GHOB -2.270736   0.008576 -264.767   <2e-16 ***
## Orig_code6RTAS -1.988784   0.007477 -265.981   <2e-16 ***
## Orig_code7GDAR -0.797620   0.007089 -112.513   <2e-16 ***
## Orig_code7RNTE -1.893522   0.008806 -215.022   <2e-16 ***
## Orig_code8ACTE -1.921309   0.005511 -348.631   <2e-16 ***
## Dest_code1RNSW  0.389478   0.003899   99.894   <2e-16 ***
## Dest_code2GMEL -0.007616   0.004244   -1.794   0.0727 .  
## Dest_code2RVIC -0.781258   0.004654 -167.854   <2e-16 ***
## Dest_code3GBRI  0.795909   0.004037  197.178   <2e-16 ***
## Dest_code3RQLD  1.516186   0.003918  386.955   <2e-16 ***
## Dest_code4GADE -0.331189   0.005232  -63.304   <2e-16 ***
## Dest_code4RSAU -0.627202   0.006032 -103.980   <2e-16 ***
## Dest_code5GPER  1.390114   0.005022  276.811   <2e-16 ***
## Dest_code5RWAU  0.367314   0.005362   68.509   <2e-16 ***
## Dest_code6GHOB -1.685934   0.008478 -198.859   <2e-16 ***
## Dest_code6RTAS -1.454819   0.007612 -191.112   <2e-16 ***
## Dest_code7GDAR -0.308516   0.007716  -39.986   <2e-16 ***
## Dest_code7RNTE -1.462020   0.009743 -150.060   <2e-16 ***
## Dest_code8ACTE -1.506283   0.005709 -263.866   <2e-16 ***
## log(dist)      -1.589102   0.001685 -942.842   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2750417  on 209  degrees of freedom
## Residual deviance:  335759  on 180  degrees of freedom
## AIC: 337818
## 
## Number of Fisher Scoring iterations: 6

And the various flows and goodness-of-fit statistics?

mdatasub <- mdatasub %>% 
  mutate(doubsimFitted = round(fitted(doubSim)))

#now we can create pivot table to turn paired list into matrix (and compute the margins as well)
mdatasubmat7 <- mdatasub %>% 
  select(Orig_code, Dest_code, doubsimFitted) %>% 
  pivot_wider(names_from = Dest_code, values_from = doubsimFitted, values_fill = 0) %>% 
  relocate("1GSYD", .after = "Orig_code") %>% 
  rowwise(Orig_code) %>% 
  mutate(Flow_from_zone = sum(c_across(where(is.numeric))))

mdatasubmat7

compared to ..

mdatasubmat

and we can test the goodness-of-fit in exactly the same way as before:

# Evaluate R^2
rsqrd(truth = mdatasub$Flow, estimate = mdatasub$doubsimFitted)

## [1] 0.8662571

# Evaluate RMSE 
rMSE(truth = mdatasub$Flow, estimate = mdatasub$doubsimFitted)

## [1] 4877.799

So the goodness of fit has shot up and we can clearly see the origin and destination constraints working, and for most sets of flows, the model is now producing some good estimates. However, there are still some errors in the flows.

Is there anything more we can do? Yes, of course there is.

4.1 Tweaking our Models

4.1.1 Distance Decay

Now, all of the way through these practicals, we have assumed that the distance decay parameter follows a negative power law. Well, it doesn’t need to.

In Wilson’s original paper, he generalised the distance decay parameter to:

\(f(d_{ij})\)

Where $f$ represents some function of distance describing the rate at which the flow interactions change as distance increases. Lots of people have written about this, including Tayor (1971) and more recently Robin Lovelace in a transport context, here.

For the inverse power law that we have been using is one possible function of distance, the other common one that is used is the negative exponential function:

\(exp(-\beta d_{ij})\)

We can get a feel for how different distance decay parameters work by plotting some sample data (try different parameters):

xdist <-  seq(1, 20)
invpower2 <-  xdist^-2
negexp0.3 <- exp(-0.3*xdist) 

df <- tibble(xdist, invpower2, negexp0.3) %>% 
  pivot_longer(!xdist, names_to = "decay", values_to = "values")

df %>% 
  ggplot(mapping = aes(x = xdist, y = values, color = decay)) +
  geom_line(size = 1)

In this particular example with these parameters and 𝛽𝛽 values, the inverse power function has a far more rapid distance decay effect than the negative exponential function. In real life, what this means is that if the observed interactions drop off very rapidly with distance, then they might be more likely to follow an inverse power law. This might be the case when looking at trips to the local convenience store by walking, for example. On the other hand, if the effect of distance is less severe – for example, migration across the country for a new job – then the negative exponential function with a small value of \(𝛽\) function might be more appropriate. There is no hard and fast rule as to which function to pick. It will just come down to which fits the data better.

As Tayor Oshan points out in his excellent Primer what this means in our Poisson regression model is that we simply substitute − \(βln⁡dij\) for$ −βdi$ in our model:

# Run a doubly constrained SIM
doubSim1 <- glm(Flow ~ Orig_code + Dest_code + dist, family = poisson(link = "log"), na.action = na.exclude, data = mdatasub)

summary(doubSim1)

## 
## Call:
## glm(formula = Flow ~ Orig_code + Dest_code + dist, family = poisson(link = "log"), 
##     data = mdatasub, na.action = na.exclude)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -127.953   -31.964    -4.223    22.000   224.899  
## 
## Coefficients:
##                  Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)     1.103e+01  4.029e-03 2737.291   <2e-16 ***
## Orig_code1RNSW -2.329e-01  3.365e-03  -69.212   <2e-16 ***
## Orig_code2GMEL -1.690e-01  3.437e-03  -49.177   <2e-16 ***
## Orig_code2RVIC -8.485e-01  4.000e-03 -212.105   <2e-16 ***
## Orig_code3GBRI  3.146e-02  3.528e-03    8.916   <2e-16 ***
## Orig_code3RQLD  5.454e-01  3.553e-03  153.532   <2e-16 ***
## Orig_code4GADE -9.537e-01  4.453e-03 -214.143   <2e-16 ***
## Orig_code4RSAU -1.525e+00  5.851e-03 -260.707   <2e-16 ***
## Orig_code5GPER  1.018e+00  5.249e-03  193.958   <2e-16 ***
## Orig_code5RWAU -7.128e-01  5.632e-03 -126.548   <2e-16 ***
## Orig_code6GHOB -2.065e+00  8.093e-03 -255.208   <2e-16 ***
## Orig_code6RTAS -1.876e+00  7.043e-03 -266.303   <2e-16 ***
## Orig_code7GDAR -6.602e-01  7.197e-03  -91.730   <2e-16 ***
## Orig_code7RNTE -2.069e+00  8.824e-03 -234.480   <2e-16 ***
## Orig_code8ACTE -1.771e+00  5.423e-03 -326.524   <2e-16 ***
## Dest_code1RNSW  3.090e-01  3.793e-03   81.473   <2e-16 ***
## Dest_code2GMEL  2.154e-01  4.000e-03   53.842   <2e-16 ***
## Dest_code2RVIC -2.193e-01  4.196e-03  -52.272   <2e-16 ***
## Dest_code3GBRI  5.674e-01  3.962e-03  143.209   <2e-16 ***
## Dest_code3RQLD  1.270e+00  3.853e-03  329.606   <2e-16 ***
## Dest_code4GADE -6.360e-01  5.150e-03 -123.501   <2e-16 ***
## Dest_code4RSAU -9.072e-01  5.982e-03 -151.663   <2e-16 ***
## Dest_code5GPER  1.692e+00  5.567e-03  304.000   <2e-16 ***
## Dest_code5RWAU  9.058e-03  5.795e-03    1.563    0.118    
## Dest_code6GHOB -1.543e+00  8.007e-03 -192.767   <2e-16 ***
## Dest_code6RTAS -1.408e+00  7.211e-03 -195.206   <2e-16 ***
## Dest_code7GDAR -1.465e-01  7.831e-03  -18.711   <2e-16 ***
## Dest_code7RNTE -1.631e+00  9.753e-03 -167.270   <2e-16 ***
## Dest_code8ACTE -1.281e+00  5.597e-03 -228.950   <2e-16 ***
## dist           -1.392e-03  1.858e-06 -749.415   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2750417  on 209  degrees of freedom
## Residual deviance:  521907  on 180  degrees of freedom
## AIC: 523966
## 
## Number of Fisher Scoring iterations: 6

mdatasub <- mdatasub %>% 
  mutate(doubsimFitted1 = round(fitted(doubSim1)))

# Evaluate R^2
rsqrd(truth = mdatasub$Flow, estimate = mdatasub$doubsimFitted1)

## [1] 0.7591672

# Evaluate RMSE 
rMSE(truth = mdatasub$Flow, estimate = mdatasub$doubsimFitted1)

## [1] 6577.102

So, it would appear that in this case using a negative exponential function in our model results in a worse outcome than the initial inverse power law - this may not always be the case, so it is worth experimenting.

4.1.2 Distance Decay

Yes, the nice thing about doing all of this in a regression modelling framework is we can just keep adding predictor variables into the mix and seeing whether they have an effect.

You can’t add origin or destination specific predictors into a doubly constrained model like this, however, switching back to the singly constrained models, as many different origin or destination predictor variables can be added as seems reasonable (subject to the usual restrictions on high correlation)

In addition to variables relating to median income, there are variables on unemployment rate and the percentage of households living in rented accommodation. Experiment with these variables for origins and destinations to see whether the singly constrained models can be improved in any way.

kitchensinkSIM <- glm(Flow ~ Dest_code + vi1_origpop + vi2_origunemp + vi3_origmedinc + vi4_origpctrent -1, na.action = na.exclude, family = poisson(link = "log"), data = mdatasub)
#let's have a look at it's summary...
summary(kitchensinkSIM)

## 
## Call:
## glm(formula = Flow ~ Dest_code + vi1_origpop + vi2_origunemp + 
##     vi3_origmedinc + vi4_origpctrent - 1, family = poisson(link = "log"), 
##     data = mdatasub, na.action = na.exclude)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -157.58   -50.78   -25.36    -1.61   376.17  
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## Dest_code1GSYD   8.508e+00  9.273e-03  917.53   <2e-16 ***
## Dest_code1RNSW   8.851e+00  9.302e-03  951.47   <2e-16 ***
## Dest_code2GMEL   8.661e+00  9.576e-03  904.37   <2e-16 ***
## Dest_code2RVIC   8.315e+00  9.551e-03  870.53   <2e-16 ***
## Dest_code3GBRI   8.651e+00  9.294e-03  930.72   <2e-16 ***
## Dest_code3RQLD   8.931e+00  9.276e-03  962.84   <2e-16 ***
## Dest_code4GADE   7.603e+00  9.972e-03  762.41   <2e-16 ***
## Dest_code4RSAU   7.185e+00  1.040e-02  691.08   <2e-16 ***
## Dest_code5GPER   8.116e+00  9.629e-03  842.84   <2e-16 ***
## Dest_code5RWAU   7.706e+00  9.868e-03  780.83   <2e-16 ***
## Dest_code6GHOB   6.444e+00  1.169e-02  551.42   <2e-16 ***
## Dest_code6RTAS   6.698e+00  1.112e-02  602.39   <2e-16 ***
## Dest_code7GDAR   6.593e+00  1.133e-02  581.73   <2e-16 ***
## Dest_code7RNTE   6.010e+00  1.280e-02  469.58   <2e-16 ***
## Dest_code8ACTE   7.340e+00  1.022e-02  718.15   <2e-16 ***
## vi1_origpop      4.240e-07  5.802e-10  730.74   <2e-16 ***
## vi2_origunemp    7.089e-02  1.375e-03   51.55   <2e-16 ***
## vi3_origmedinc   2.932e-04  7.138e-06   41.07   <2e-16 ***
## vi4_origpctrent -2.239e-02  2.265e-04  -98.84   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 23087017  on 210  degrees of freedom
## Residual deviance:  1441432  on 191  degrees of freedom
## AIC: 1443469
## 
## Number of Fisher Scoring iterations: 6

# First round of the estimates
mdatasub %>% 
  mutate(attrsimFitted2 = round(fitted(kitchensinkSIM))) -> mdatasub

and we can test the goodess of fit in exactly the same way as before:

# Evaluate R^2
rsqrd(truth = mdatasub$Flow, estimate = mdatasub$attrsimFitted2)

## [1] 0.3126191

# Evaluate RMSE 
rMSE(truth = mdatasub$Flow, estimate = mdatasub$attrsimFitted2)

## [1] 10883.35

5. Conclusions, further notes and ideas for additional activities

Hopefully you have now seen how it is extremely straight-forward to run and calibrate Wilson’s full family of Spatial Interaction Models in R using GLM and Poisson Regression.

5.1 Some Further Notes

Now might be the time to mention that despite everything I’ve shown you, there has been some discussion in the literature as to whether the Poisson Model is actually a misspecification, especially for modelling migration flows. If you have the stomach for it, this paper by Congdon goes into a lot of detail.

The issue is a thing called ‘overdispersion’ which, translated, essentially relates to the model not being able to capture all of the things that could be explaining the flows in the independent variables that are supplied to the model. The details are tedious and only really intelligible to those with a statistics background. If you want a starter, try here, but in practical terms, we can get around this problem by fitting a very similar sort of regression model called the negative binomial regression model.

If you wish, you can read up and experiment with this model - you can fit it in exactly the same way as the glm model but using a function called glm.nb which is part of the mass package. The negative binomial model has an extra parameter in the model for overdispersion. You you do try this, you will almost certainly discover that your results barely change - but hell, you might keep a pedantic reviewer at bay if you submit this to a journal (not that I’m speaking from experience or anything).

And some more comments

Another thing to note is that the example we used here had quite neat data. You will almost certainly run into problems if you have sparse data or predictors with 0s in them. If this happens, then you might need to either drop some rows in your data (if populated with 0s) or substitute 0s for very small numbers, much less than 1, but greater than 0 (this is because you can’t take the log of 0)

And another thing to note is that the models in this Australian example assumed that the flow data and predictors were all in and around the same order or magnitude. This is not necessarily the case, particularly with a 5 year migration transition and some large metropolitan areas. Where data that (such as population masses at origins and destinations) that are an order of magnitude different (i.e. populations about ten times larger in different locations) then the model estimates might be biased. Fortunately, there are packages available to help us with these problems as well. The robustbase package features a function called glmrob which will deal with this issues (again, your results probably won’t change much, but worth knowing).

Further Reading https://www.researchgate.net/publication/255576515_DEVELOPING_THE_SINGLY_CONSTRAINED_GRAVITY_MODEL_FOR_APPLICATION_IN_DEVELOPING_COUNTRIES

• Abel G J (2010) Estimation of international migration flow tables in Europe: international migration flow tables. Journal of the Royal Statistical Society. Series A, Statistics in Society 173(4): 797–825.

• Congdon P (1993) Approaches to modelling overdispersion in the analysis of migration. Environment and Planning A: Economy and Space 25(1): 1481–1510.

• Congdon P (1988) Modelling migration flows between areas: an analysis for London using the census and OPCS Longitudinal Study. Regional Studies 23(2): 87–103.

• Crymble A, Dennett A and Hitchcock T (2017) Modelling regional imbalances in English plebeian migration to late eighteenth-century London. The Economic History Review 71(3): 747-771.

• Dennett A and Wilson A (2013) A multi-level spatial interaction modelling framework for estimating inter-regional migration in Europe. Environment and Planning A: Economy and Space 45(6): 1491–1507.

• Erhardt G D and Dennett A (2017) Understanding the role and relevance of the census in a changing transportation data landscape. Paper presented at the Transportation Research Board Conference on Applying Census Data for Transportation, Kansas City, Missouri.

• Flowerdew R (2010) Modelling migration with Poisson regression. In: Stillwell J, Duke-Williams O and Dennett A (eds) Technologies for Migration and Commuting Analysis: Spatial Interaction Data Applications. Hershey PA: IGI Global.

• Flowerdew R (1982) Fitting the lognormal gravity model to heteroscedastic data. Geographical Analysis 14(3): 263–267.

• Flowerdew R and Aitkin M (1982) A method of fitting the gravity model based on the Poisson distribution. Journal of Regional Science 22(2): 191–202.

• Fotheringham A S (1983) A new set of spatial-interaction models: the theory of competing destinations.

• Environment and Planning A: Economy and Space 15(1): 15–36. Fotheringham A S, Nakaya T, Yano K, Openshaw S and Ishikawa Y (2001) Hierarchical destination choice and spatial interaction modelling: a simulation experiment. Environment and Planning A: Economy and Space 33(5): 901–920.

• Kim K and Cohen J E (2010) Determinants of international migration flows to and from industrialized countries: a panel data approach beyond gravity. International Migration Review 44(4): 899–932.

• Lee E S (1966) A theory of migration. Demography 3(1): 47–57.

• Lomax N and Norman P (2016) Estimating population attribute values in a table: “get me started in” iterative proportional fitting. The Professional Geographer 68(3): 451–461.

• Lovelace R (2015) Estimating distance decay for the national propensity to cycle tool. https://www.slideshare.net/ITSLeeds/estimating-distance-decay-for-the-national-propensityto-cycle-tool.

• Oshan T M (2016) A primer for working with the Spatial Interaction modeling (SpInt) module in the python spatial analysis library (PySAL). REGION 3: R11–R23.

• Pooler J (1994) An extended family of spatial interaction models. Progress in Human Geography 18(1): 17–39.

• Pooler J (1987) Modeling interprovincial migration using entropy-maximizing methods. The Canadian Geographer 31(1): 57–64.

• Raymer J (2007) The estimation of international migration flows: a general technique focused on the origin–destination association structure. Environment and Planning A: Economy and Space 39(4): 985–995.

• Raymer J and Abel G (2008) Methods to improve estimates of migration flows – the MIMOSA model for estimating international migration flows in the European Union. UNECE/Eurostat work session on migration statistics, Working Paper 8, Geneva.

• Raymer J, Abel G and Smith P W (2007) Combining census and registration data to estimate detailed elderly migration flows in England and Wales. Journal of the Royal Statistical Society. Series A, Statistics in Society 170(4): 891–908.

• Raymer J and Giulietti C (2010) Analysing structures of interregional migration in England. In: Stillwell J, Duke-Williams O and Dennett A (eds) Technologies for Migration and Commuting Analysis: Spatial Interaction Data Applications. Hershey PA: IGI Global.

• Rees P (1977) The measurement of migration from census data and other sources. Environment and Planning A: Economy and Space 9(3): 257–280.

• Rogers A and Raymer J (1998) The spatial focus of US interstate migration flows. Population, Space and Place 4(1): 63–80.

• Senior M L (1979) From gravity modelling to entropy maximizing: a pedagogic guide. Progress in Human Geography 3(2): 175–210.

• Shen J (2017) Modelling interregional migration in China in 2005–2010: the roles of regional attributes and spatial interaction effects in modelling error. Population, Space and Place 23(3): e2014.

• Shen J (2015) Explaining interregional migration changes in China, 1985–2000, using a decomposition approach. Regional Studies 49(7): 1176–1192.

• Stillwell J (1978) Interzonal migration: some historical tests of spatial-interaction models.

• Environment and Planning A: Economy and Space 10(1): 1187–1200.

• Taylor P J (1983) Distance decay in spatial interactions. Norwich: Geo Books.

• Willekens F (1999) Modeling approaches to the indirect estimation of migration flows: from entropy to EM. Mathematical Population Studies 7(3): 239–278.

• Wilson A (1971) A family of spatial interaction models, and associated developments. Environment and Planning A: Economy and Space 3(1): 1–32.

• Zipf G K (1946) The P1 P2 / D hypothesis: on the intercity movement of persons. American Sociological Review 11(6): 677–686