Nothing
Tutorial for TDLM
In TDLM: Systematic Comparison of Trip Distribution Laws and Models
knitr::opts_chunk$set(echo = TRUE, message = TRUE, warning = FALSE,
fig.width = 8, fig.height = 8)
# Packages --------------------------------------------------------------------
suppressPackageStartupMessages({
suppressWarnings({
library(TDLM)
library(sf)
})
})
options(tinytex.verbose = TRUE)
Introduction
This tutorial aims to describe the different features of the R package TDLM.
The main purpose of the TDLM package is to provide a rigorous framework for
fairly comparing trip distribution laws and models [@Lenormand2016]. This
general framework is based on a two-step approach for generating mobility flows
by separating the trip distribution law, such as gravity or intervening
opportunities, from the modeling approach used to generate flows from this law.
A short note on terminology
This framework is part of the four-step travel model. It corresponds to the
second step, called trip distribution, which aims to match trip origins with
trip destinations. The model used to generate trips or flows, and more
generally the degree of interaction between locations, is often referred to as
a spatial interaction model. Depending on the research area, a matrix or a
network formalism may be used to describe these spatial interactions.
Origin-Destination matrices (or trip tables) are common in geography or
transportation, while in statistical physics or complex systems studies, the
term mobility networks is usually preferred.
Origin–Destination matrix
The description of movements within a given area is represented by an
Origin-Destination (OD) matrix. The area of interest is divided into $n$
locations, and $T_{ij}$ represents the volume of flows between location $i$ and
location $j$. This volume typically corresponds to the number of trips or a
commuting flow (i.e., the number of individuals living in $i$ and working in
$j$). The OD matrix is square, contains only positive values, and may have a
zero diagonal (Figure 1).
Figure 1: Schematic representation of an Origin-Destination matrix.
Aggregated inputs information
Three categories of inputs are typically considered to simulate an OD matrix
(Figure 2). The masses and distances are the primary ingredients used to
generate a matrix of probabilities based on a given distribution law. Thus, the
probability $p_{ij}$ of observing a trip from location $i$ to location $j$
depends on the masses, the demand at the origin ($m_i$), and the offer at
the destination ($m_j$). Typically, population is used as a surrogate for
demand and offer. The probability of movement also depends on the costs,
which are based on the distance $d_{ij}$ between locations or the number of
opportunities $s_{ij}$ between locations, depending on the chosen law (more
details are provided in the next "Trip distribution laws" section). In general,
the effect of the cost can be adjusted with a parameter.
The margins are used to generate an OD matrix from the matrix of probabilities
while preserving the total number of trips ($N$), the number of outgoing trips
($O_i$), and/or the number of incoming trips ($D_j$) (more details are provided
in the "Constrained distribution models" section).
Figure 2: Schematic representation of the aggregated inputs information.
Trip distribution laws
The purpose of a trip distribution law is to estimate the probability $p_{ij}$
that, out of all the possible travels in the system, there is one between
location $i$ and location $j$. This probability is asymmetric in $i$ and $j$,
as are the flows themselves. It takes the form of a square matrix of probabilities.
This probability is normalized across all possible pairs of origins and
destinations, such that $\sum_{i,j=1}^n p_{ij} = 1$. Therefore, a matrix of
probabilities can be obtained by normalizing any OD matrix (Figure 3).
Figure 3: Schematic representation of the matrix of probabilities.
As mentioned in the previous section, most trip distribution laws depend
on the demand at the origin ($m_i$), the offer at the destination ($m_j$), and
a cost to move from $i$ to $j$. There are two major approaches for estimating
the probability matrix. The traditional gravity approach, in analogy with
Newton's law of gravitation, is based on the assumption that the amount of trips
between two locations is related to their populations and decays according to a
function of the distance $d_{ij}$ between locations. In contrast to the gravity
law, the laws of intervening opportunities hinge on the assumption that the
number of opportunities $s_{ij}$ between locations plays a more important role
than the distance [@Lenormand2016]. This fundamental difference between the two
schools of thought is illustrated in Figure 4.
Figure 4: Illustration of the fundamental difference between gravity and
intervening opportunity laws.
It is important to note that the effect of the cost between locations (distance
or number of opportunities) can usually be adjusted with a parameter, which can
be calibrated automatically or by comparing the simulated matrix with observed
data.
Constrained distribution models
The purpose of the trip distribution models is to generate an OD matrix
$\tilde{T}=(\tilde{T}{ij})$ by
drawing at random $N$ trips from the trip distribution law
$(p{ij})_{1 \leq i,j \leq n}$ respecting different level of constraints
according to the model. We considered four different types of models in this
package. As can be observed in Figure 5, the four models respect different level
of constraints from the total number of trips to the total number of out-going
and in-coming trips by locations (i.e. the margins).
Figure 5: Schematic representation of the constrained distribution models.
More specifically, the volume of flows $\tilde{T}_{ij}$ is generated from the
matrix of probability with multinomial random draws that will take different
forms according to the model used [@Lenormand2016]. Therefore, since the process
is stochastic, each simulated matrix is unique and composed of integers. Note
that it is also possible to generate an average matrix from the multinomial
trials.
Goodness-of-fit measures
Finally, the trip distribution laws can be calibrated, and both the trip
distribution laws and models can be evaluated by comparing a simulated matrix
$\tilde{T}$ with the observed matrix $T$. These comparisons rely on various
goodness-of-fit measures, which may or may not account for the distance
between locations. These measures are described above.
Common Part of Commuters (CPC)
$$\displaystyle CPC(T,\tilde{T}) = \frac{2\cdot\sum_{i,j=1}^n min(T_{ij},\tilde{T}_{ij})}{N + \tilde{N}}$$
[@Gargiulo2012;@Lenormand2012;@Lenormand2016]
Normalized Root Mean Square Error (NRMSE)
$$\displaystyle NRMSE(T,\tilde{T}) = \sqrt{\frac{\sum_{i,j=1}^n (T_{ij}-\tilde{T}_{ij})^2}{N}}$$
Kullback–Leibler divergence (KS)
$$\displaystyle KL(T,\tilde{T}) = \sum_{i,j=1}^n \frac{T_{ij}}{N}\log\left(\frac{T_{ij}}{N}\frac{\tilde{N}}{\tilde{T}_{ij}}\right)$$
[@Kullback1951]
Common Part of Links (CPL)
$$\displaystyle CPL(T,\tilde{T}) = \frac{2\cdot\sum_{i,j=1}^n 1_{T_{ij}>0} \cdot 1_{\tilde{T}{ij}>0}}{\sum{i,j=1}^n 1_{T_{ij}>0} + \sum_{i,j=1}^n 1_{\tilde{T}_{ij}>0}}$$
[@Lenormand2016]
Common Part of Commuters based on the disance (CPC_d)
Let us consider $N_k$ (and
$\tilde{N}_k$) the sum of observed (and simulated) flows at a distance comprised
in the bin $[\mbox{bin_size} \cdot k-\mbox{bin_size}, \mbox{bin_size} \cdot k[$.
$$\displaystyle CPC_d(T,\tilde{T}) = \frac{2\cdot\sum_{k=1}^{\infty} min(N_{k},\tilde{N}_{k})}{N+\tilde{N}}$$
[@Lenormand2016]
Kolmogorv-Smirnov statistic and p-value (KS).
These measures, described in @Massey1951, are based on the observed and
simulated flow distance distributions and are computed using the ks_test
function from the Ecume package.
Example of commuting in Kansas
Data
In this example, we will use commuting data from US Kansas in 2000 to illustrate
the main functions of the package. The dataset comprises three tables providing
information on commuting flows between the 105 US
Kansas counties in 2000. The observed OD matrix
od is a zero-diagonal square
matrix of integers. Each element of the matrix represents the number of
commuters between a pair of US Kansas counties.
data(od)
od[1:5, 1:5]
dim(od)
The aggregated data are composed of the
distance matrix,
data(distance)
distance[1:5, 1:5]
dim(distance)
and the masses and margins contained in the data.frame
mass.
data(mass)
mass[1:10,]
dim(mass)
mi <- as.numeric(mass[,1])
names(mi) <- rownames(mass)
mj <- mi
Oi <- as.numeric(mass[,2])
names(Oi) <- rownames(mass)
Dj <- as.numeric(mass[,3])
names(Dj) <- rownames(mass)
The data must always be based on the same number of locations sorted in the same
order. The function
check_format_names
can be used to verify the validity of all inputs before running the main
functions of the package.
check_format_names(vectors = list(mi = mi, mj = mj, Oi = Oi, Dj = Dj),
matrices = list(od = od, distance = distance),
check = "format_and_names")
Optional spatial information are also provided here.
county is a spatial
object containing the geometry of the 105 US Kansas counties in 2000.
library(sf)
data(county)
county[1:10,]
plot(county)
coords and
coords_xy are two
dataframes providing longitude/latitude and X/Y coordinates for each location,
respectively.
coords[1:10,]
coords_xy[1:10,]
Extract Additional Spatial Information
The functions
extract_distances,
extract_opportunities,
and
extract_spatial_information
can be used to extract matrices of distances and the number of intervening
opportunities.
The first function computes distances in kilometers between pairs of locations
based on geographical coordinates. It can calculate either great-circle
distances, using longitude/latitude coordinates and the Haversine formula
haversine_d <- extract_distances(coords = coords,
method = "Haversine")
haversine_d[1:5, 1:5]
distance[1:5, 1:5]
or Euclidean distances based on X/Y coordinates
xy_d <- extract_distances(coords = coords_xy,
method = "Euclidean")
oldmar <- par()$mar
par(mar = c(4.5, 6, 1, 1))
plot(haversine_d, xy_d, xlim=c(0,900), ylim=c(0,900),
type="p", pch=16, cex=2, lty=1, lwd=3,
col="steelblue3", axes=FALSE, xlab="", ylab="")
axis(1, cex.axis=1.2)
axis(2, cex.axis=1.2, las=1)
mtext("Haversine (km)", 1, line = 3.25, cex = 1.75)
mtext("Euclidean (km)", 2, line = 4, cex = 1.75)
box(lwd=1.5)
par(mar = oldmar)
The second function computes the number of opportunities between pairs of
locations. For a given pair of locations, the number of opportunities between
the origin and the destination is based on the number of opportunities within
a circle of radius equal to the distance between the origin and the destination,
centered at the origin. The number of opportunities at the origin and the
destination are not included. In our case, the number of inhabitants ($m_i$) is
used as a proxy for the number of opportunities.
sij <- extract_opportunities(opportunity = mi,
distance = distance,
check_names = TRUE)
sij[1:5, 1:5]
The last function takes as input a spatial object containing the geometry of
the locations that can be handled by the
sf package. It returns a matrix of
great-circle distances between locations (expressed in km). An optional id
can also be provided to be used as names for the outputs.
spi <- extract_spatial_information(county, id = "ID")
sp_d <- spi$distance
sp_d[1:5, 1:5]
distance[1:5, 1:5]
This function also allows extracting the surface area of each location
(in square kilometers), which can be useful to calibrate the trip distribution
laws' parameter value (see below).
mean(spi$surface)
Run functions
The main function of the package is
run_law_model. The
function has two sets of arguments, one for the law and another one for the
model. The inputs (described above) necessary to run this function depends on
the law (either the matrix of distances or number of opportunities can be used,
or neither of them for the uniform law) and on the model and its associated
constraints (number of trips, out-going trips and/or in-coming trips). The
example below will generate three simulated ODs with the normalized gravity law
with an exponential distance decay function [@Lenormand2016] and the Doubly
Constrained Model.
res <- run_law_model(law = "NGravExp",
mass_origin = mi,
mass_destination = mj,
distance = distance,
opportunity = NULL,
param = 0.01,
write_proba = TRUE,
model = "DCM",
nb_trips = NULL,
out_trips = Oi,
in_trips = Dj,
average = FALSE,
nbrep = 3)
The output is an object of class TDLM. In this case it is a list of matrices
composed of the three simulated matrices (replication_1, replication_2 and
replication_3), the matrix of probabilities (called proba) associated with
the law and returned only if write_proba = TRUE. The objects of class TDLM
contain a table info summarizing the simulation run.
print(res)
str(res)
This simulation run was based on one parameter value. It is possible to use a
vector instead of a scalar for the param argument.
res <- run_law_model(law = "NGravExp",
mass_origin = mi,
mass_destination = mj,
distance = distance,
opportunity = NULL,
param = c(0.01,0.02),
write_proba = TRUE,
model = "DCM",
nb_trips = NULL,
out_trips = Oi,
in_trips = Dj,
average = FALSE,
nbrep = 3)
In this case a list of list of matrices will be returned (one for each parameter
value).
print(res)
str(res)
It is also important to note that the radiation law and the uniform law are free
of parameter.
res <- run_law_model(law = "Rad",
mass_origin = mi,
mass_destination = mj,
distance = NULL,
opportunity = sij,
param = NULL,
write_proba = TRUE,
model = "DCM",
nb_trips = NULL,
out_trips = Oi,
in_trips = Dj,
average = FALSE,
nbrep = 3)
print(res)
The argument average can be used to generate an average matrix based on a
multinomial distribution (based on an infinite number of drawings). In this
case, the models' inputs can be either positive integer or real numbers and the
output (nbrep = 1 in this case) will be a matrix of positive real numbers.
res$replication_1[1:10,1:10]
res <- run_law_model(law = "Rad",
mass_origin = mi,
mass_destination = mj,
distance = NULL,
opportunity = sij,
param = NULL,
write_proba = TRUE,
model = "DCM",
nb_trips = NULL,
out_trips = Oi,
in_trips = Dj,
average = TRUE,
nbrep = 3)
print(res)
res$replication_1[1:10,1:10]
The functions run_law
and run_model have
been designed to run only one of the two components of the two-step approach.
They function the same as a
run_law_model, but
it is worth noting that only inter-location flows are considered for the
distribution laws, meaning that the matrix of probabilities (and associated
simulated OD matrices) generated by a given distribution law with
run_law_model or
run_law is a
zero-diagonal matrix. Nevertheless, it is possible to generate intra-location
flows with run_model
taking any kind of matrix of probabilities as input.
Parameters' calibration & models' evaluation
The package contains two function to help calibrating and evaluating the model.
The function gof computes
goodness-of-fit measures between observed and simulated OD matrices and the
function calib_param
that estimates the optimal parameter value for a given law and a given spatial
distribution of location based on the Figure 8 in [@Lenormand2016].
Let us illustrate the trip distribution laws and models' calibration with the
the normalized gravity law with an exponential distance decay function and the
Doubly Constrained Model. Based on the average surface area of the Kansas
counties (in square kilometers) it seems that the optimal parameter value of the
normalized gravity law with an exponential distance decay function (as described
in [@Lenormand2016]) for commuting in US Kansas counties is around 0.085.
print(calib_param(av_surf = mean(spi$surface), law = "NGravExp"))
This is just an estimation that help us to identify the potential range of
parameter value, so in order to rigorously calibrate and evaluate the
trip distribution law and model we need to compute the goodness-of-fit
measure for a wide range of parameter values.
res <- run_law_model(law = "NGravExp",
mass_origin = mi,
mass_destination = mj,
distance = distance,
opportunity = NULL,
param = seq(0.05,0.1,0.005),
write_proba = TRUE,
model = "DCM",
nb_trips = NULL,
out_trips = Oi,
in_trips = Dj,
average = FALSE,
nbrep = 3)
calib <- gof(sim = res, obs = od, measures = "all", distance = distance)
print(calib)
All the necessary information is stored in the object calib, most of the
goodness-of-fit measures agree on a parameter value of 0.075 in that case with
an associated average Common Part of Commuter equal to 85.6%.
cpc <- aggregate(calib$CPC, list(calib$Parameter_value), mean)[,2]
oldmar <- par()$mar
par(mar = c(4.5, 6, 1, 1))
plot(seq(0.05,0.1,0.005), cpc, type="b", pch=16, cex=2, lty=1, lwd=3,
col="steelblue3", axes=FALSE, xlab="", ylab="")
axis(1, cex.axis=1.2)
axis(2, cex.axis=1.2, las=1)
mtext("Parameter value", 1, line = 3.25, cex = 1.75)
mtext("Common Part of Commuters", 2, line = 4, cex = 1.75)
box(lwd=1.5)
par(mar = oldmar)
Reference
Try the TDLM package in your browser
Any scripts or data that you put into this service are public.
TDLM documentation built on June 8, 2025, 10:39 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo = TRUE, message = TRUE, warning = FALSE, fig.width = 8, fig.height = 8) # Packages -------------------------------------------------------------------- suppressPackageStartupMessages({ suppressWarnings({ library(TDLM) library(sf) }) }) options(tinytex.verbose = TRUE)
Introduction
This tutorial aims to describe the different features of the R package TDLM.
The main purpose of the TDLM package is to provide a rigorous framework for
fairly comparing trip distribution laws and models [@Lenormand2016]. This
general framework is based on a two-step approach for generating mobility flows
by separating the trip distribution law, such as gravity or intervening
opportunities, from the modeling approach used to generate flows from this law.
A short note on terminology
This framework is part of the four-step travel model. It corresponds to the second step, called trip distribution, which aims to match trip origins with trip destinations. The model used to generate trips or flows, and more generally the degree of interaction between locations, is often referred to as a spatial interaction model. Depending on the research area, a matrix or a network formalism may be used to describe these spatial interactions. Origin-Destination matrices (or trip tables) are common in geography or transportation, while in statistical physics or complex systems studies, the term mobility networks is usually preferred.
Origin–Destination matrix
The description of movements within a given area is represented by an Origin-Destination (OD) matrix. The area of interest is divided into $n$ locations, and $T_{ij}$ represents the volume of flows between location $i$ and location $j$. This volume typically corresponds to the number of trips or a commuting flow (i.e., the number of individuals living in $i$ and working in $j$). The OD matrix is square, contains only positive values, and may have a zero diagonal (Figure 1).
Aggregated inputs information
Three categories of inputs are typically considered to simulate an OD matrix (Figure 2). The masses and distances are the primary ingredients used to generate a matrix of probabilities based on a given distribution law. Thus, the probability $p_{ij}$ of observing a trip from location $i$ to location $j$ depends on the masses, the demand at the origin ($m_i$), and the offer at the destination ($m_j$). Typically, population is used as a surrogate for demand and offer. The probability of movement also depends on the costs, which are based on the distance $d_{ij}$ between locations or the number of opportunities $s_{ij}$ between locations, depending on the chosen law (more details are provided in the next "Trip distribution laws" section). In general, the effect of the cost can be adjusted with a parameter.
The margins are used to generate an OD matrix from the matrix of probabilities while preserving the total number of trips ($N$), the number of outgoing trips ($O_i$), and/or the number of incoming trips ($D_j$) (more details are provided in the "Constrained distribution models" section).
Trip distribution laws
The purpose of a trip distribution law is to estimate the probability $p_{ij}$ that, out of all the possible travels in the system, there is one between location $i$ and location $j$. This probability is asymmetric in $i$ and $j$, as are the flows themselves. It takes the form of a square matrix of probabilities. This probability is normalized across all possible pairs of origins and destinations, such that $\sum_{i,j=1}^n p_{ij} = 1$. Therefore, a matrix of probabilities can be obtained by normalizing any OD matrix (Figure 3).
As mentioned in the previous section, most trip distribution laws depend on the demand at the origin ($m_i$), the offer at the destination ($m_j$), and a cost to move from $i$ to $j$. There are two major approaches for estimating the probability matrix. The traditional gravity approach, in analogy with Newton's law of gravitation, is based on the assumption that the amount of trips between two locations is related to their populations and decays according to a function of the distance $d_{ij}$ between locations. In contrast to the gravity law, the laws of intervening opportunities hinge on the assumption that the number of opportunities $s_{ij}$ between locations plays a more important role than the distance [@Lenormand2016]. This fundamental difference between the two schools of thought is illustrated in Figure 4.
It is important to note that the effect of the cost between locations (distance or number of opportunities) can usually be adjusted with a parameter, which can be calibrated automatically or by comparing the simulated matrix with observed data.
Constrained distribution models
The purpose of the trip distribution models is to generate an OD matrix $\tilde{T}=(\tilde{T}{ij})$ by drawing at random $N$ trips from the trip distribution law $(p{ij})_{1 \leq i,j \leq n}$ respecting different level of constraints according to the model. We considered four different types of models in this package. As can be observed in Figure 5, the four models respect different level of constraints from the total number of trips to the total number of out-going and in-coming trips by locations (i.e. the margins).
More specifically, the volume of flows $\tilde{T}_{ij}$ is generated from the matrix of probability with multinomial random draws that will take different forms according to the model used [@Lenormand2016]. Therefore, since the process is stochastic, each simulated matrix is unique and composed of integers. Note that it is also possible to generate an average matrix from the multinomial trials.
Goodness-of-fit measures
Finally, the trip distribution laws can be calibrated, and both the trip distribution laws and models can be evaluated by comparing a simulated matrix $\tilde{T}$ with the observed matrix $T$. These comparisons rely on various goodness-of-fit measures, which may or may not account for the distance between locations. These measures are described above.
Common Part of Commuters (CPC)
$$\displaystyle CPC(T,\tilde{T}) = \frac{2\cdot\sum_{i,j=1}^n min(T_{ij},\tilde{T}_{ij})}{N + \tilde{N}}$$
[@Gargiulo2012;@Lenormand2012;@Lenormand2016]
Normalized Root Mean Square Error (NRMSE)
$$\displaystyle NRMSE(T,\tilde{T}) = \sqrt{\frac{\sum_{i,j=1}^n (T_{ij}-\tilde{T}_{ij})^2}{N}}$$
Kullback–Leibler divergence (KS)
$$\displaystyle KL(T,\tilde{T}) = \sum_{i,j=1}^n \frac{T_{ij}}{N}\log\left(\frac{T_{ij}}{N}\frac{\tilde{N}}{\tilde{T}_{ij}}\right)$$ [@Kullback1951]
Common Part of Links (CPL)
$$\displaystyle CPL(T,\tilde{T}) = \frac{2\cdot\sum_{i,j=1}^n 1_{T_{ij}>0} \cdot 1_{\tilde{T}{ij}>0}}{\sum{i,j=1}^n 1_{T_{ij}>0} + \sum_{i,j=1}^n 1_{\tilde{T}_{ij}>0}}$$ [@Lenormand2016]
Common Part of Commuters based on the disance (CPC_d)
Let us consider $N_k$ (and $\tilde{N}_k$) the sum of observed (and simulated) flows at a distance comprised in the bin $[\mbox{bin_size} \cdot k-\mbox{bin_size}, \mbox{bin_size} \cdot k[$.
$$\displaystyle CPC_d(T,\tilde{T}) = \frac{2\cdot\sum_{k=1}^{\infty} min(N_{k},\tilde{N}_{k})}{N+\tilde{N}}$$
[@Lenormand2016]
Kolmogorv-Smirnov statistic and p-value (KS).
These measures, described in @Massey1951, are based on the observed and simulated flow distance distributions and are computed using the ks_test function from the Ecume package.
Example of commuting in Kansas
Data
In this example, we will use commuting data from US Kansas in 2000 to illustrate
the main functions of the package. The dataset comprises three tables providing
information on commuting flows between the 105 US
Kansas counties in 2000. The observed OD matrix
od is a zero-diagonal square
matrix of integers. Each element of the matrix represents the number of
commuters between a pair of US Kansas counties.
data(od) od[1:5, 1:5] dim(od)
The aggregated data are composed of the distance matrix,
data(distance) distance[1:5, 1:5] dim(distance)
and the masses and margins contained in the data.frame mass.
data(mass) mass[1:10,] dim(mass) mi <- as.numeric(mass[,1]) names(mi) <- rownames(mass) mj <- mi Oi <- as.numeric(mass[,2]) names(Oi) <- rownames(mass) Dj <- as.numeric(mass[,3]) names(Dj) <- rownames(mass)
The data must always be based on the same number of locations sorted in the same
order. The function
check_format_names
can be used to verify the validity of all inputs before running the main
functions of the package.
check_format_names(vectors = list(mi = mi, mj = mj, Oi = Oi, Dj = Dj), matrices = list(od = od, distance = distance), check = "format_and_names")
Optional spatial information are also provided here. county is a spatial object containing the geometry of the 105 US Kansas counties in 2000.
library(sf) data(county) county[1:10,] plot(county)
coords and coords_xy are two dataframes providing longitude/latitude and X/Y coordinates for each location, respectively.
coords[1:10,] coords_xy[1:10,]
Extract Additional Spatial Information
The functions
extract_distances,
extract_opportunities,
and
extract_spatial_information
can be used to extract matrices of distances and the number of intervening
opportunities.
The first function computes distances in kilometers between pairs of locations based on geographical coordinates. It can calculate either great-circle distances, using longitude/latitude coordinates and the Haversine formula
haversine_d <- extract_distances(coords = coords, method = "Haversine") haversine_d[1:5, 1:5] distance[1:5, 1:5]
or Euclidean distances based on X/Y coordinates
xy_d <- extract_distances(coords = coords_xy, method = "Euclidean") oldmar <- par()$mar par(mar = c(4.5, 6, 1, 1)) plot(haversine_d, xy_d, xlim=c(0,900), ylim=c(0,900), type="p", pch=16, cex=2, lty=1, lwd=3, col="steelblue3", axes=FALSE, xlab="", ylab="") axis(1, cex.axis=1.2) axis(2, cex.axis=1.2, las=1) mtext("Haversine (km)", 1, line = 3.25, cex = 1.75) mtext("Euclidean (km)", 2, line = 4, cex = 1.75) box(lwd=1.5) par(mar = oldmar)
The second function computes the number of opportunities between pairs of locations. For a given pair of locations, the number of opportunities between the origin and the destination is based on the number of opportunities within a circle of radius equal to the distance between the origin and the destination, centered at the origin. The number of opportunities at the origin and the destination are not included. In our case, the number of inhabitants ($m_i$) is used as a proxy for the number of opportunities.
sij <- extract_opportunities(opportunity = mi, distance = distance, check_names = TRUE) sij[1:5, 1:5]
The last function takes as input a spatial object containing the geometry of
the locations that can be handled by the
sf package. It returns a matrix of
great-circle distances between locations (expressed in km). An optional id
can also be provided to be used as names for the outputs.
spi <- extract_spatial_information(county, id = "ID") sp_d <- spi$distance sp_d[1:5, 1:5] distance[1:5, 1:5]
This function also allows extracting the surface area of each location
(in square kilometers), which can be useful to calibrate the trip distribution
laws' parameter value (see below).
mean(spi$surface)
Run functions
The main function of the package is run_law_model. The function has two sets of arguments, one for the law and another one for the model. The inputs (described above) necessary to run this function depends on the law (either the matrix of distances or number of opportunities can be used, or neither of them for the uniform law) and on the model and its associated constraints (number of trips, out-going trips and/or in-coming trips). The example below will generate three simulated ODs with the normalized gravity law with an exponential distance decay function [@Lenormand2016] and the Doubly Constrained Model.
res <- run_law_model(law = "NGravExp", mass_origin = mi, mass_destination = mj, distance = distance, opportunity = NULL, param = 0.01, write_proba = TRUE, model = "DCM", nb_trips = NULL, out_trips = Oi, in_trips = Dj, average = FALSE, nbrep = 3)
The output is an object of class TDLM. In this case it is a list of matrices
composed of the three simulated matrices (replication_1, replication_2 and
replication_3), the matrix of probabilities (called proba) associated with
the law and returned only if write_proba = TRUE. The objects of class TDLM
contain a table info summarizing the simulation run.
print(res) str(res)
This simulation run was based on one parameter value. It is possible to use a
vector instead of a scalar for the param argument.
res <- run_law_model(law = "NGravExp", mass_origin = mi, mass_destination = mj, distance = distance, opportunity = NULL, param = c(0.01,0.02), write_proba = TRUE, model = "DCM", nb_trips = NULL, out_trips = Oi, in_trips = Dj, average = FALSE, nbrep = 3)
In this case a list of list of matrices will be returned (one for each parameter value).
print(res) str(res)
It is also important to note that the radiation law and the uniform law are free of parameter.
res <- run_law_model(law = "Rad", mass_origin = mi, mass_destination = mj, distance = NULL, opportunity = sij, param = NULL, write_proba = TRUE, model = "DCM", nb_trips = NULL, out_trips = Oi, in_trips = Dj, average = FALSE, nbrep = 3) print(res)
The argument average can be used to generate an average matrix based on a
multinomial distribution (based on an infinite number of drawings). In this
case, the models' inputs can be either positive integer or real numbers and the
output (nbrep = 1 in this case) will be a matrix of positive real numbers.
res$replication_1[1:10,1:10] res <- run_law_model(law = "Rad", mass_origin = mi, mass_destination = mj, distance = NULL, opportunity = sij, param = NULL, write_proba = TRUE, model = "DCM", nb_trips = NULL, out_trips = Oi, in_trips = Dj, average = TRUE, nbrep = 3) print(res) res$replication_1[1:10,1:10]
The functions run_law and run_model have been designed to run only one of the two components of the two-step approach. They function the same as a run_law_model, but it is worth noting that only inter-location flows are considered for the distribution laws, meaning that the matrix of probabilities (and associated simulated OD matrices) generated by a given distribution law with run_law_model or run_law is a zero-diagonal matrix. Nevertheless, it is possible to generate intra-location flows with run_model taking any kind of matrix of probabilities as input.
Parameters' calibration & models' evaluation
The package contains two function to help calibrating and evaluating the model. The function gof computes goodness-of-fit measures between observed and simulated OD matrices and the function calib_param that estimates the optimal parameter value for a given law and a given spatial distribution of location based on the Figure 8 in [@Lenormand2016].
Let us illustrate the trip distribution laws and models' calibration with the the normalized gravity law with an exponential distance decay function and the Doubly Constrained Model. Based on the average surface area of the Kansas counties (in square kilometers) it seems that the optimal parameter value of the normalized gravity law with an exponential distance decay function (as described in [@Lenormand2016]) for commuting in US Kansas counties is around 0.085.
print(calib_param(av_surf = mean(spi$surface), law = "NGravExp"))
This is just an estimation that help us to identify the potential range of parameter value, so in order to rigorously calibrate and evaluate the trip distribution law and model we need to compute the goodness-of-fit measure for a wide range of parameter values.
res <- run_law_model(law = "NGravExp", mass_origin = mi, mass_destination = mj, distance = distance, opportunity = NULL, param = seq(0.05,0.1,0.005), write_proba = TRUE, model = "DCM", nb_trips = NULL, out_trips = Oi, in_trips = Dj, average = FALSE, nbrep = 3) calib <- gof(sim = res, obs = od, measures = "all", distance = distance) print(calib)
All the necessary information is stored in the object calib, most of the goodness-of-fit measures agree on a parameter value of 0.075 in that case with an associated average Common Part of Commuter equal to 85.6%.
cpc <- aggregate(calib$CPC, list(calib$Parameter_value), mean)[,2] oldmar <- par()$mar par(mar = c(4.5, 6, 1, 1)) plot(seq(0.05,0.1,0.005), cpc, type="b", pch=16, cex=2, lty=1, lwd=3, col="steelblue3", axes=FALSE, xlab="", ylab="") axis(1, cex.axis=1.2) axis(2, cex.axis=1.2, las=1) mtext("Parameter value", 1, line = 3.25, cex = 1.75) mtext("Common Part of Commuters", 2, line = 4, cex = 1.75) box(lwd=1.5) par(mar = oldmar)
Reference
Try the TDLM package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.