In dymium-org/dymiumCore: A Toolkit for Building a Dynamic Microsimulation Model for Integrated Urban Modelling

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>", 
  out.width = '100%',
  fig.align = "center", 
  fig.width =  7,
  fig.height =  5,
  dpi = 200,
  warning = FALSE
)

library(dymiumCore)
library(data.table)
library(sf)
library(ggplot2)
library(patchwork)
set.seed(728)

Introduction

Residential mobility or household relocation plays an important role in determining the urban traffic and the real estate market condition. In this example, we create a simple residential mobility model which comprises of two stages. The first stage is for households to decide whether to move from their current residence. If they decide to move then they will go onto the second stage otherwise they will remain where they are. The second stage is for the moving households to decide where they will relocate to.

knitr::include_graphics("residential-mobility-flowchart.png")

Required knowledge

This example requires the reader to know the grammar of data operations of the data.table package as it is quite different from the base R data.frame and also knows the R6 package to fully utilise the base classes of dymiumCore.

Assumptions

To make this model simple, the following assumptions were made.

1) All households have 5% chance to move within one year. 2) Moving households only consider the neighboring zones of their current zone. 3) Only interzone relocation is considered.

Here is the map of all the zones in this example.

# set seed for reproducible widget id
if (requireNamespace("htmltools", quietly = TRUE)) {
  htmlwidgets::setWidgetIdSeed(42)
}

library(leaflet)
library(sf)

centers <- toy_zones %>%
  sf::st_centroid()

leaflet() %>%
  addTiles() %>%
  setView(lat = -37.808, lng = 144.952, zoom = 12) %>%
  addPolygons(
    data = toy_zones,
    color = "#444444",
    weight = 1,
    smoothFactor = 0.5,
    opacity = 1.0,
    fillOpacity = 0.5,
    fillColor = ~ colorQuantile("BuPu", albers_sqm)(albers_sqm), 
    highlightOptions = highlightOptions(
      color = "white",
      weight = 2,
      bringToFront = TRUE
    )
  ) %>%
  addLabelOnlyMarkers(
    data = centers,
    label = ~ zid,
    labelOptions = labelOptions(
      noHide = TRUE,
      direction = 'top',
      # textOnly = TRUE,
      textsize = 5
    )
  )

Here is the list of queen contiguity of the zones.

st_queen <- function(a, b = a) st_relate(a, b, pattern = "F***T****", sparse = TRUE)
nb_list <- st_queen(toy_zones)

Since nb_list refers to the zones by their row numbers as listed in the toy_zones, we need to replace those row numbers with their respective zid.

# make the returned object from `st_queen` into a list by removing its class and attributes.
class(nb_list) <- NULL 
attributes(nb_list) <- NULL
# name all elements in the list by their associated `zid`.
names(nb_list) <- toy_zones$zid
# go inside each element of the list to replace row numbers with their `zid` in `toy_zones`.
nb_list <-
  lapply(nb_list, function(x) {
    as.integer(toy_zones$zid[x])
  })
nb_list

Create a microsimulation model

Prepare the initial state of the world

The provided toy datasets from the dymiumCore package are used in this example. We are going to create only Household agents.

new_toy_households <- merge(toy_households, toy_dwellings[, c("did", "zid")], by = "did")

# create a emptied World object to store agents and models
world <- World$new()
world$add(x = Household$new(new_toy_households, id_col = "hid"), name = "Household")

# calling the `world` object
world

You can see that there is only one entity inside the world object which Household agents. The Household object contains attributes and methods of its agents and can be accessed by calling the get() method of world

For example

world$get("Household")

A reference to any objects in world can be created by assigning the result from world$get() to a variable.

Hh <- world$get("Household")
Hh

To access the attributes of Individual agents you may use the get_data() method as following

Hh$get_data()

Prepare models

First we need to create a list that contains the assumed transition probabilities.

move_rate <- list(yes = 0.05, no = 0.95)
move_rate

As you can see, we have created a named list that contains two elements which are 'yes' and 'no' with numeric values that represent the associated probabilities to the actions. All households have 5% chance that they would move and 95% chance that they would remain where they are.

The next step is to store the move_rate model and the list that contains the indices of neighboring zones, nb_list, in side the world object as following.

world$add(x = move_rate, name = "move_rate")
world$add(x = nb_list, name = "nb_zones")

# it can be accessed using the get method as well. 
world$get("move_rate")

# the safer option is to use `get_model` method which only look for keys that contain model objects while the `get` method looks at all keys stored.
# world$get_model("move_rate")

TransitionClassification allows a probabilistic model such as the one we have just created, move_rate, to be used for simulation (Monte Carlo).

The steps of the residential mobility event are inside the following function.

event_ltmobility_residentialRelocation <- function(world) {
  # Create a reference to the household object stored in world
  Hh <- world$get("Household")  

  # simulate the first stage - move or not move? ----
  TransMove <- TransitionClassification$new(x = Hh, model = world$get("move_rate")$get())

  # get_result returns a data.table with two columns `id` and `response`.
  moving_household_ids <- TransMove$get_result()[response == "yes", id]

  # simulate the second stage - where to? ----
  if (length(moving_household_ids) > 0) {
    # create a choiceset for each moving household
    # get only `hid` and `zid` of the moving households by creating a copy of 
    # their attribute data.
    mover_choicemodel <- Hh$get_data(ids = moving_household_ids)[, c("hid", "zid")]

    # create emptied columns to store choicesets and their probabilities
    mover_choicemodel[, `:=`(choices = list(), probs = list())] 

    # get the list of neighboring zones. Note that the second `$get()` is required 
    # to get the list instead of not a model object. 
    nb_zones <- world$get("nb_zones")$get() 

    for (i in 1:nrow(mover_choicemodel)) {
      self_zid <- mover_choicemodel[i, zid] # get `zid` in row `i`
      choices <- unlist(nb_zones[names(nb_zones) %in% self_zid]) # get the neighboring zones
      probs <- rep(x = 1 / length(choices), times = length(choices)) # assign equal probabilities to all choices
      set(mover_choicemodel, i = i, j = "choices", value = choices)
      set(mover_choicemodel, i = i, j = "probs", value = probs) 
    }

    # remove `zid` since the TransitionClassification assumes that all columns that 
    # are not named as `probs` and `choices` are supposed to be used as matching columns  
    mover_choicemodel[, zid := NULL]

    # simulate location selection
    TransWhereTo <- TransitionClassification$new(
      x = Hh,
      model = mover_choicemodel,
      targeted_agents = mover_choicemodel[, hid],
      model_by_id = TRUE
    )$update_agents(attr = "zid")
  }

  return(world)
}

Create a microsimulation pipeline

To simulate the residential mobility event for multiple iterations we can simply use a flow control statement such as a 'for' loop. The chunk of codes below runs the residential mobility event that we have just created for 10 iterations. At the end of each iteration the total number of households in each zone will be stored inside count_ls which will be indexed by its iteration number.

start_iter <- 1
end_iter <- 10
count_ls <- list() 
for (i in start_iter:end_iter) {
  world$start_iter(i, unit = "year") %>%
    event_ltmobility_residentialRelocation(.)

  # record the count of households by zone in each iteration
  count_ls[[i]] <-
    world$get("Household")$get_data() %>% # returns the attribute data of the households
    .[, .(count = .N, iteration = i), by = zid] # count the total number of households by zone
}

Visualise the simulation result

To visualise the spatial dynamic of households we use the gganimate and ggplot packages.

count <- rbindlist(count_ls)
count <- merge(count, toy_zones, by = "zid") %>%
  st_as_sf()

As a spatial plot..

p_before <-
  ggplot(data = subset(count, iteration == 1)) +
    geom_sf(aes(fill = count)) +
    geom_sf_label(aes(label = sa2_name11), label.size = 0.7) +
    scale_fill_viridis_c(option = "B", limits = c(0,30)) +
    labs(title = paste0("Iteration: ", 1))

p_after <- 
   ggplot(data = subset(count, iteration == 10)) +
    geom_sf(aes(fill = count)) +
    geom_sf_label(aes(label = sa2_name11), label.size = 0.7) +
    scale_fill_viridis_c(option = "B", limits = c(0,30)) +
    labs(title = paste0("Iteration: ", 10))

p_before / p_after +
  plot_layout(guides = 'collect') &
  plot_annotation(title = 'Spatial distribution of households')

As a time series plot..

ggplot(count, aes(iteration, count, group = sa2_name11)) +
  geom_line() +
  geom_point(size = 1.5) +
  coord_cartesian(clip = 'off') +
  theme_bw() +
  scale_x_continuous(breaks = scales::pretty_breaks()) +
  labs(title = "The total number of households by SA2 zone") +
  theme(plot.margin = margin(5.5, 40, 5.5, 5.5)) +
  facet_wrap(~ sa2_name11)

dymium-org/dymiumCore documentation built on July 18, 2021, 5:10 p.m.