Prerequisites
You need to have installed
R and a suitable editor such as RStudio
or VSCode with R
plugin. See the package’s README for
instructions on installing the {simodels} package.
library(simodels)
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
library(ggplot2)
library(sf)
#> Linking to GEOS 3.10.2, GDAL 3.4.3, PROJ 8.2.0; sf_use_s2() is TRUE
Preparing a SIM
Key to the workings of the {simodels} package is the conversion of
geographic objects representing origins and destinations into an OD
dataset. In this case, we already have an OD dataset, so this step is
less relevant. However, we will take this step in any case because many
SIMs start without such a comprehensive OD dataset as we have in this
case.
Prepare the OD dataset as follows:
od_sim = si_to_od(origins = zones, destinations = zones)
#> Maximum distance is > 100km. The 'cheap' measure is inaccurate over such
#> large distances, you'd likely be better using a different 'measure'.
names(od_sim)
#> [1] "O" "D" "distance_euclidean"
#> [4] "origin_GCCSA_NAME" "origin_AREA_SQKM" "destination_GCCSA_NAME"
#> [7] "destination_AREA_SQKM" "geometry"
Note that the output has duplicate columns: si_to_od()
joins data from the origin and destination objects into the resulting OD
object.
An unconstrained
SIM
A simplistic SIM - in this case an inverse power distance decay
function (negative exponential is another commonly used decay function)
- can be created just based on the distance between points:
si_power = function(d, beta) (d / 1000)^beta
od_calculated = si_calculate(
od_sim,
fun = si_power,
d = distance_euclidean,
beta = -0.8
)
plot(od_calculated["interaction"], logz = TRUE)
#> Warning in classInt::classIntervals(v0, min(nbreaks, n.unq), breaks, warnSmallN
#> = FALSE): var has infinite values, omitted in finding classes
This approach, ignoring all variables at the level of trip origins
and destinations, results in flow estimates with no units. Before
learning how to run constrained SIMs, let’s scale the result by the
total flow and see how far we are from reality, just focussing on the
interzonal OD pairs:
od_interzonal = od %>%
filter(Orig_code != Dest_code)
od_calculated_interzonal = od_calculated %>%
filter(O != D)
scale_factor = sum(od_interzonal$Flow) /
sum(od_calculated_interzonal$interaction)
od_calculated_interzonal = od_calculated_interzonal %>%
mutate(interaction_scaled = interaction * scale_factor)
od_joined = inner_join(
od_calculated_interzonal,
od %>% rename(O = Orig_code, D = Dest_code)
)
#> Joining, by = c("O", "D")
od_joined %>%
ggplot() +
geom_point(aes(Flow, interaction_scaled))
cor(od_joined$Flow, od_joined$interaction_scaled)^2
#> [1] 0.1195504
The results show that a simple unconstrained model, without any
parameter fitting, can explain less than 20% of the variability in
flows. We can do better!
od_joined %>%
mutate(decay = distance_euclidean^-0.8) %>%
mutate(decay = decay * (sum(Flow) / sum(decay))) %>%
ggplot() +
geom_point(aes(distance_euclidean, Flow)) +
geom_line(aes(distance_euclidean, decay), colour = "red")
A production
constrained SIM
The first logical way to improve model fit is to run a production
constrained model. To do that, we’ll first calculate the total number of
people leaving each zone and then use the
constraint_production argument:
od_originating = od_joined %>%
group_by(O) %>%
mutate(originating_per_zone = sum(Flow)) %>%
ungroup()
od_constrained_p = si_calculate(
od_originating,
fun = si_power,
d = distance_euclidean,
beta = -0.8,
constraint_production = originating_per_zone
)
od_constrained_p %>%
ggplot() +
geom_point(aes(Flow, interaction))
cor(od_constrained_p$Flow, od_constrained_p$interaction)^2
#> [1] 0.3127995
Progress! We have more than doubled the predictive ability of our
model by using a ‘production constrained’ SIM, as defined mathematically
in the simodels-first-principles
vignette.
Training a SIM
An advantage of the flow data used in this example is that we already
know the interaction. (This raises the question of why a SIM is needed,
answer: to test our models and demonstrate the techniques.)
We can do this using the nls() function as follows:
library(minpack.lm)
f = Flow ~ a * (distance_euclidean)^b
m = nlsLM(
formula = f,
data = od_originating,
)
#> Warning in nlsLM(formula = f, data = od_originating, ): No starting values specified for some parameters.
#> Initializing 'a', 'b' to '1.'.
#> Consider specifying 'start' or using a selfStart model
#> Warning in nls.lm(par = start, fn = FCT, jac = jac, control = control, lower = lower, : lmdif: info = -1. Number of iterations has reached `maxiter' == 50.
m
#> Nonlinear regression model
#> model: Flow ~ a * (distance_euclidean)^b
#> data: od_originating
#> a b
#> 2.936e+07 -5.779e-01
#> residual sum-of-squares: 1.482e+10
#>
#> Number of iterations till stop: 50
#> Achieved convergence tolerance: 1.49e-08
#> Reason stopped: Number of iterations has reached `maxiter' == 50.
# Nonlinear regression model
# model: Flow ~ a * (distance_euclidean)^b
# data: od_originating
# a b
# 2.182e+07 -5.801e-01
od_joined %>%
mutate(decay = distance_euclidean^-5.801e-01) %>%
mutate(decay = decay * 2.182e+07) %>%
ggplot() +
geom_point(aes(distance_euclidean, Flow)) +
geom_line(aes(distance_euclidean, decay), colour = "red")
od_pred = si_predict(od_originating, model = m)
cor(od_pred$Flow, od_pred$interaction)^2
#> [1] 0.1234891
od_pred_const = si_predict(od_originating, model = m,
constraint_production = originating_per_zone)
cor(od_pred_const$Flow, od_pred_const$interaction)^2
#> [1] 0.2701773
library(tmap)
ttm()
#> tmap mode set to interactive viewing
tm_shape(od_pred_const) +
tm_lines("interaction_scaled", palette = "viridis")