MNO_uncertainty: Generation of network scenarios

This notebook provides supplementary material for my thesis project titled: “Tackling Uncertainty in the Spatial Density Estimation from Mobile Network Operator data”. It focuses on the first module of the MNO-simulator workflow, the generation of network scenarios. Its goal is to provide the complete code for reproduction and in-depth analysis of the four different network scenarios. Further information of this project can be found here.

As mentioned on page 4 of the thesis article, we need to create three data sources:

• "Operating area data: We assume the operating area to be discretized into a regular grid, e.g., 100m x 100m. Each square unit is termed as a tile, indexed in \(j = 1,2, \ldots, J\). \(u_j\) denotes the unknown non-negative number of mobile phones of the \(j\)th tile, stored in the column vector \(u := [u_{1}… u_{J}]^{T}\). This standardized unit helps us set a formalism framework that has been used in multiple studies beforehand [19, 22, 28].

• Aggregated network event data: Let the \(i\)th element \(c_{i}\) of the column vector \(c := [c_{1}… c_{I}]^{T}\) define the observed number of phones counted in cell \(i = 1,2, \ldots, I\). Furthermore, we formalize the total number of phones across all cells by \(C := \sum_{i = 1}^{I} c_{i} = 1_{I}^{T}c\), the assumed total number of individual phones in the operating area. As in [19] \(c\) describes a cell count vector, aggregating all estimable phones around a reference time \(t^*\). In this study we ignore the time dimension and focus only on the spatial density estimation of a certain point in time.

• Network topology data: This data source describes the configuration parameters of each relevant cell within the operating area. It can come in various levels of detail, ranging from only knowing the cell tower locations to having precise information on the azimuth direction, transmission power, beamwidth, etc. This data is used within MNO’s to improve the network via radio propagation modeling."

We plan to create one operating area, which resembles our ground truth population (GTP. Then, we create four network scenarios, which are based respectively on four network topology datasets. This in turn also creates four aggregated network event vectors (i.e., \(c\)). The remaining notebook is structured in the following way:

• General Setup

• Description of the operating area

• Description of the general network topology parameters

• Description of each network scenario

1 General Setup: Loading Packages and Custom Functions

# Data manipulation
library(tidyverse)
library(data.table) 

# Spatial operations
library(sf)
library(raster)
library(stars)

# Matrix operations
library(Matrix)

# MNO data handling and propagation model setup
# Credits to Prof. Martijn Tennekes https://github.com/mtennekes/mobloc
library(mobloc)

# Comparison of 2d histograms (Kantorovitch Wasserstein distance a.k.a. Earth Movers distance)
# Credits to Prof. Stefano Gualandi https://cran.r-project.org/web/packages/SpatialKWD/SpatialKWD.pdf
library(SpatialKWD)

# Output organisation and plotting support
library(ggthemes)
library(viridis)
library(ggrepel)
library(ggpointdensity)
library(scattermore)
library(grid)
library(gridExtra)
library(knitr)
library(DT)

# seed for reproducibility
set.seed(42)


# Loading Custom functions
source("pipeline functions.R")

2 Generation of the operating area

Our operating area is based on a semi-synthetic data generated process. For this, census data from Germany on a 100m*100m regular grid has been used, which can be downloaded here. Each element in this grid is expressed as a tile. For computation purposes only a small area of Germany was used for the toyworld, namely the area of Munich and its near surroundings. This focus area includes 160,000 tiles. The code for clipping this specific area can be found here, which is also part of this research repository. For a mobile phone population the regular census population values are used. To mimic the mobile phone population of one mobile network operator (MNO) the population is reduced to about a third.

2.1 Mobile phone population

In the following chunk we load in the focus area data and create the area object, which contains all information concerning the operating area, e.g., geographical information, mobile phone population, prior information, etc. Each object has its purpose in the remaining workflow.

# data read in
munich.raw <- readRDS("Data/munich.rds")

# define raster object from focus area and set the coordinate reference system (crs)
munich.raster <- rasterFromXYZ(munich.raw, crs = st_crs(3035)$proj4string)

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum Unknown based on GRS80 ellipsoid in CRS definition

# define empty list object where all GTP objects will be stored
munich <- NULL

# define sf version of raster object
munich$area.sf.raw <- munich.raster %>%
  st_as_stars() %>%
  st_as_sf() %>%
  mutate(tile.id.chr = as.character(tile.id)) %>% 
  mutate(tile.id.fac = factor(tile.id.chr)) %>% 
  mutate(tile.id.num = as.numeric(tile.id.fac)) %>% 
  mutate(type = "NA") %>%  # only necessary if different tile types can be defined (urban, rural, etc...)
  mutate(prior.1 = 1, # uniform prior
         prior.2 = case_when(pop <= 1 ~ 0.1,
                             pop > 1 & pop <= 50 ~ 1,
                             pop > 50 ~ 100)) %>% # 3 category prior
  dplyr::select(contains("tile.id"), pop, type, elevation, contains("prior"))

# define sf version with only geometry centroids of tiles
munich$area.sf <- munich$area.sf.raw %>%
  mutate(centroid.geometry = st_centroid(.$geometry)) %>% 
  mutate(X.centroid = unlist(map(.$centroid.geometry, 1)),
         Y.centroid = unlist(map(.$centroid.geometry, 2))) 
  

# regular dataframe version
munich$area.df <- munich$area.sf %>%
  st_drop_geometry()

# prior value reduced dataframe version
munich$area.prior <- munich$area.df %>% 
  dplyr::select(tile.id.chr, contains("prior."))

# variable reduced dataframe version
munich$area.reduced.df <- munich$area.df %>% 
  dplyr::select(tile.id.chr, contains("centroid"))

# unionized version of focus area
munich$area.union <- munich$area.sf %>%
  st_union()

# bounding box coordinates of focus area
munich$area.bbox <- munich$area.union %>%
  st_bbox(crs = sf::st_crs(3035))

# specify raster object and tile id number
munich$area.raster <- munich.raster %>%
  raster(., layer = "tile.id")

# specify raster object and elevation value of each tile (here considered as constant)
munich$area.elevation <- munich.raster %>%
  raster(., layer = "elevation")

# number of tiles
munich$area.params[["tile.num"]] <- length(munich$area.df$tile.id)

# size of tiles
munich$area.params[["base.tile.size"]] <- as.numeric(sqrt(st_area(munich$area.sf[1,])))

# storing everything in area object
area <- munich

saveRDS(area, "Data/area.rds")

In the following, we can see the spatial density of the GTP (see Figure 4 in the paper).

# adjustable break points for map categories
breaks <- c(0, 2, 5, 10, 20, 50, 100, 200, 350, Inf)
# plot map and print
(GTP.spat.dens.plot <- area$area.sf %>% 
  mutate(pop.cat = cut(pop, breaks = breaks, dig.lab = 7, right = F)) %>% 
  map_density(data = ., var = "pop.cat", label = "GTP"))

Figure 2.1: Spatial density of the ground truth population

# ggsave("Final Plots/GTP.png", GTP.spat.dens.plot, device = "png", bg = "transparent")

In the following, we can see the empirical cumulative complementary distribution (ECCDF) of the GTP (see Figure 4 in the paper). The capped empirical cumulative distribution (ECDF) is added as an insert.

# ECCDF and ECDF of GTP
(GTP.spat.dens.plot <- density_plots(area$area.df))

Figure 2.2: ECCDF and ECDF of the ground truth population, (Insert ECDF capped at 30)

# ggsave("Final Plots/ECCDF.png", GTP.spat.dens.plot, device = "png")

# warning can be ignored

2.2 Radio network

This subsection refers to the development of a radio network within our operating area. Developing the radio network is heavily dependent on the mobloc package, which is promoted through the European Statistical System. However, we have adjusted and extended some of the functions to our needs, e.g. concerning the generation of multiple cell layers.

In general, the mobloc package allows to define many parameters, however, if they are not defined, default parameters, set by the package, are used. This makes it very easy to implement as much / as little information one has on a certain network and always getting it to work. The following adjustable parameters and default values are provided through the package for creating a radio network and modelling the signal strength in a focus area:

# possible mobloc parameters
mobloc_param()

## $W
## [1] 10
## 
## $W_small
## [1] 5
## 
## $ple
## [1] 3.75
## 
## $ple_small
## [1] 6
## 
## $ple_0
## [1] 3.5
## 
## $ple_1
## [1] 4
## 
## $midpoint
## [1] -92.5
## 
## $steepness
## [1] 0.2
## 
## $range
## [1] 10000
## 
## $range_small
## [1] 1000
## 
## $height
## [1] 30
## 
## $height_small
## [1] 8
## 
## $tilt
## [1] 5
## 
## $beam_v
## [1] 9
## 
## $beam_h
## [1] 65
## 
## $azim_dB_back
## [1] -30
## 
## $elev_dB_back
## [1] -30
## 
## $sig_d_th
## [1] 0.005
## 
## $max_overlapping_cells
## [1] 100
## 
## $TA_step
## [1] 78.12
## 
## $TA_max
## [1] 1282
## 
## $TA_buffer
## [1] 1
## 
## attr(,"class")
## [1] "mobloc_param"

For our study we construct four networks. The first two networks are one-layer networks and differ in cell density, while the other two networks are multi-layer networks. We initialize each network with the number of layers and each layer refers to a certain kind of cells (Macro, Meso and Micro). The development of each layer starts with a hexagonal grid in which the points define tower locations. The hexagons have different sizes (i.e. tower distance) dependent on the layer and may contain some randomness within the layer (jitter). Each layer spans over the complete focus area.

On each tower three directional cells are placed that are directed in a 120° angle to each other. All layers contain a rotation parameter (azimuth) to prevent cells of different layers broadcasting into the exact same direction, in reference to the operating area. In this study, we do not implement omni-directional cells, therefore, all mobloc parameters with the suffix “_small” are not used.

The cells are specified with layer specific parameters (e.g. height, power, path loss exponent, etc.). All cells of one network scenario are specified in a so-called cellplan. This cellplan is a data frame on the cell level, which contains all information on the parameters. When the cellplan is completed the cell-specific broadcasting profile (i.e. cell profile) is estimated and projected onto the operating area. The function compute_sig_strength() computes the distance, signal strength and signal dominance between any tile and any cell. Furthermore, a minimum parameter is implemented that defines the minimum signal dominance value an cell-tile relationship needs to have in order to be considered “covered”. For the generation of a network scenario, we assure that all tiles are sufficiently covered, i.e. each mobile phones within each tile have at least one cell that covers them with a signal dominance value greater than the minimum threshold, in this case 0.05.

In the following, each cellplan will be presented respectively. First, we will report some important theoretical network parameters. Then, different coverage maps will be visualized. Finally, a histogram of

2.2.1 Cellplan.1: One layer, sparse network

# specify parameters of each cell
ME.cell.param.mobloc <- mobloc_param(W = 50, range = 8000, ple = 3.7, height = 10,
                                     midpoint = -85, steepness = 0.3, sig_d_th = 0.05)

layer.list <- list(ME.cell.param.mobloc)

# create dataframe for theoretical signal strength distribution
param.df <- map_dfr(layer.list, rbind.data.frame) %>% 
  mutate(cell.kind = c("ME"),
         label = c("Meso")) %>% 
  dplyr::select(cell.kind, label, everything(), dominance.th = sig_d_th)

# reduced data frame of theoretical signal strength distribution
param.df.reduced <- param.df %>% 
  dplyr::select(cell.kind, dominance.th)

# theoretical signal strength parameter plots
sig.pram.plots <- sig_param_plots(param.df = param.df, range.max = 20000, base_size = 11)

set.seed(100)

layer.params.ext <- list(
  ME = list(tower.dist = 4100,
            rotation.deg = 35,
            jitter = 3,
            subscript = "ME",
            seed = 7,
            mobloc.params = ME.cell.param.mobloc)
  
)


cellplan.1.layer.1 <- complete_cellplan_gen(area = area,
                                            layer.params.ext = layer.params.ext,
                                            param.df = param.df)



# how many tiles are not sufficiently covered (Check)
paste0("Number of tiles which do not reach the signal dominance threshold of: " , sum(cellplan.1.layer.1$signal.dom$missing))

## [1] "Number of tiles which do not reach the signal dominance threshold of: 0"

Here we can see some important theoretical network parameters.

(cellplan.1.theo.param <- ggpubr::as_ggplot(sig.pram.plots$final))

Here we can see some important spatial network parameters.

celplan.1.spat.param <- cell.spat.diag(area = area, cellplan = cellplan.1.layer.1, label.name = "Cellplan.1")

celplan.1.spat.param$tile.coverage.hist

celplan.1.spat.param$coverage.map.complete

2.2.2 Cellplan.2: One layer, dense network

# specify parameters of each cell
ME.cell.param.mobloc <- mobloc_param(W = 50, range = 8000, ple = 3.7, height = 10,
                                     midpoint = -85, steepness = 0.3, sig_d_th = 0.05)

layer.list <- list(ME.cell.param.mobloc)

# create dataframe for theoretical signal strength distribution
param.df <- map_dfr(layer.list, rbind.data.frame) %>% 
  mutate(cell.kind = c("ME"),
         label = c("Meso")) %>% 
  dplyr::select(cell.kind, label, everything(), dominance.th = sig_d_th)

# reduced data frame of theoretical signal strength distribution
param.df.reduced <- param.df %>% 
  dplyr::select(cell.kind, dominance.th)

# theoretical signal strength parameter plots
cellplan.2.sig.pram.plots <- sig_param_plots(param.df = param.df, range.max = 20000, base_size = 11)

set.seed(100)

layer.params.ext <- list(
  ME = list(tower.dist = 3000,
            rotation.deg = 35,
            jitter = 1,
            subscript = "ME",
            seed = 7,
            mobloc.params = ME.cell.param.mobloc)
  
)


cellplan.1.layer.2 <- complete_cellplan_gen(area = area,
                                            layer.params.ext = layer.params.ext,
                                            param.df = param.df)



# how many tiles are not sufficiently covered
paste0("Number of tiles which do not reach the signal dominance threshold of: " , sum(cellplan.1.layer.2$signal.dom$missing))

## [1] "Number of tiles which do not reach the signal dominance threshold of: 0"

Here we can see some important theoretical network parameters.

(cellplan.2.theo.param <- ggpubr::as_ggplot(cellplan.2.sig.pram.plots$final))

Here we can see some important spatial network parameters.

celplan.2.spat.param <- cell.spat.diag(area = area, cellplan = cellplan.1.layer.2, label.name = "Cellplan.2")

celplan.2.spat.param$tile.coverage.hist

celplan.2.spat.param$coverage.map.complete

2.2.3 Cellplan.3: Three layer network

# specify parameters of each cell
MA.cell.param.mobloc <- mobloc_param(W = 5, # Power in Watts
                                     range = 10000, # maximum coverage range
                                     ple = 3.4, # Path loss exponent
                                     height = 10, # height of the antenna
                                     midpoint = -85, # midpoint parameter of the logistic function for signal dominance
                                     steepness = 0.15, # steepness parameter of the logistic function for signal dominance
                                     sig_d_th = 0.05) # dominance minimum threshold 

ME.cell.param.mobloc <- mobloc_param(W = 50, range = 3500, ple = 3.8, height = 10,
                                     midpoint = -85, steepness = 0.3, sig_d_th = 0.05)

MI.cell.param.mobloc <- mobloc_param(W = 1, range = 3500, ple = 4, height = 6,
                                     midpoint = -85, steepness = 0.4, sig_d_th = 0.05)


layer.list <- list(MA.cell.param.mobloc, ME.cell.param.mobloc, MI.cell.param.mobloc)

# create dataframe for theoretical signal strength distribution
param.df <- map_dfr(layer.list, rbind.data.frame) %>% 
  mutate(cell.kind = c("MA", "ME", "MI"),
         label = c("Macro", "Meso", "Micro")) %>% 
  dplyr::select(cell.kind, label, everything(), dominance.th = sig_d_th)

# reduced data frame of theoretical signal strength distribution
param.df.reduced <- param.df %>% 
  dplyr::select(cell.kind, dominance.th)

# theoretical signal strength parameter plots
cellplan.3.sig.pram.plots <- sig_param_plots(param.df = param.df, range.max = 20000, base_size = 11)

## Warning: Using alpha for a discrete variable is not advised.

## Warning: Using alpha for a discrete variable is not advised.

set.seed(100)

layer.params.ext <- list(
  MA = list(tower.dist = 8500,
            rotation.deg = 0,
            jitter = 1000,
            subscript = "MA",
            seed = 3,
            mobloc.params = MA.cell.param.mobloc),
  ME = list(tower.dist = 3500,
            rotation.deg = 35,
            jitter = 700,
            subscript = "ME",
            seed = 7,
            mobloc.params = ME.cell.param.mobloc),
  MI = list(tower.dist = 10000,
            rotation.deg = 60,
            jitter = 2000,
            subscript = "MI",
            seed = 10,
            mobloc.params = MI.cell.param.mobloc)
  
)


cellplan.3.layer <- complete_cellplan_gen(area = area,
                                          layer.params.ext = layer.params.ext,
                                          param.df = param.df)

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

## Warning: Unknown or uninitialised column: `height`.

## Warning: Unknown or uninitialised column: `W`.

## Warning: Unknown or uninitialised column: `range`.

## Warning: Unknown or uninitialised column: `ple`.

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

## Warning: Unknown or uninitialised column: `height`.

## Warning: Unknown or uninitialised column: `W`.

## Warning: Unknown or uninitialised column: `range`.

## Warning: Unknown or uninitialised column: `ple`.

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

## Warning: Unknown or uninitialised column: `height`.

## Warning: Unknown or uninitialised column: `W`.

## Warning: Unknown or uninitialised column: `range`.

## Warning: Unknown or uninitialised column: `ple`.

# how many tiles are not sufficiently covered
paste0("Number of tiles which do not reach the signal dominance threshold of: " , sum(cellplan.3.layer$signal.dom$missing))

## [1] "Number of tiles which do not reach the signal dominance threshold of: 0"

Here we can see some important theoretical network parameters.

(cellplan.3.theo.param <- ggpubr::as_ggplot(cellplan.3.sig.pram.plots$final))

Here we can see some important spatial network parameters.

celplan.3.spat.param <- cell.spat.diag(area = area, cellplan = cellplan.3.layer, label.name = "Cellplan.3")

celplan.3.spat.param$tile.coverage.hist

celplan.3.spat.param$coverage.map.complete

2.2.4 Cellplan.4: Two layer network

# specify parameters of each cell
ME.cell.param.mobloc <- mobloc_param(W = 50, range = 8000, ple = 3.7, height = 10,
                                     midpoint = -85, steepness = 0.3, sig_d_th = 0.05)

MI.cell.param.mobloc <- mobloc_param(W = 1, range = 3500, ple = 4, height = 6,
                                     midpoint = -85, steepness = 0.4, sig_d_th = 0.05)

layer.list <- list(ME.cell.param.mobloc, MI.cell.param.mobloc)

# create dataframe for theoretical signal strength distribution
param.df <- map_dfr(layer.list, rbind.data.frame) %>% 
  mutate(cell.kind = c("ME", "MI"),
         label = c("Meso", "Micro")) %>% 
  dplyr::select(cell.kind, label, everything(), dominance.th = sig_d_th)

# reduced data frame of theoretical signal strength distribution
param.df.reduced <- param.df %>% 
  dplyr::select(cell.kind, dominance.th)

# theoretical signal strength parameter plots
cellplan.4.sig.pram.plots <- sig_param_plots(param.df = param.df, range.max = 20000, base_size = 11)

## Warning: Using alpha for a discrete variable is not advised.

## Warning: Using alpha for a discrete variable is not advised.

set.seed(100)

layer.params.ext <- list(
    ME = list(tower.dist = 4000,
            rotation.deg = 35,
            jitter = 700,
            subscript = "ME",
            seed = 7,
            mobloc.params = ME.cell.param.mobloc),
    MI = list(tower.dist = 10000,
            rotation.deg = 60,
            jitter = 2000,
            subscript = "MI",
            seed = 10,
            mobloc.params = MI.cell.param.mobloc)
  
)


cellplan.2.layer <- complete_cellplan_gen(area = area,
                                          layer.params.ext = layer.params.ext,
                                          param.df = param.df)

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

## Warning: Unknown or uninitialised column: `height`.

## Warning: Unknown or uninitialised column: `W`.

## Warning: Unknown or uninitialised column: `range`.

## Warning: Unknown or uninitialised column: `ple`.

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

## Warning: Unknown or uninitialised column: `height`.

## Warning: Unknown or uninitialised column: `W`.

## Warning: Unknown or uninitialised column: `range`.

## Warning: Unknown or uninitialised column: `ple`.

# how many tiles are not sufficiently covered
paste0("Number of tiles which do not reach the signal dominance threshold of: " , sum(cellplan.2.layer$signal.dom$missing))

## [1] "Number of tiles which do not reach the signal dominance threshold of: 0"

Here we can see some important theoretical network parameters.

(cellplan.4.theo.param <- ggpubr::as_ggplot(cellplan.4.sig.pram.plots$final))

Here we can see some important spatial network parameters.

celplan.4.spat.param <- cell.spat.diag(area = area, cellplan = cellplan.2.layer, label.name = "Cellplan.4")

celplan.4.spat.param$tile.coverage.hist

celplan.4.spat.param$coverage.map.complete

After this presentation of each cellplan we will perform now the phone-to-cell association, which is responsible for creating each vector \(c\).

2.3 Device-to-cell association (Generative model)

Here we bring the mobile phone population with each network together and stochastically assign the mobile phone population in each tile to the respective cells. The result is the so-called c-vector (i.e., \(c\)), which describes the number of mobile phones assigned to each cell. The basis parameter for this assignment is the signal dominance which is normalized in the form of a conditional probability. These conditional probabilities (emission probabilities) describe the elements of \(P\).

First, we combine and save all cellplan objects into the list object cellplans.list, which helps with automatizing each phone-to-cell association task.

cellplans.list <- list(
  cellplan.1 = cellplan.1.layer.1,
  cellplan.2 = cellplan.1.layer.2,
  cellplan.3 = cellplan.3.layer,
  cellplan.4 = cellplan.2.layer
)

# saveRDS(cellplans.list, "Data/cellplans.list.rds")

Then, we perform the phone-to-cell association. Finally, we save the object gen.model.objects, which will contain three objects for further use in the workflow:

• C.vec.df.list: This object contains the c-vector for each network scenario.

• P.long.complete.df: This object contains the for each network scenario the tile-cell relationships (edge list), which are characterized by the variables emission probabilities (generative model), signal strength and signal dominance values. Note, this object contains also the “zero” elements, meaning tile-cell relationships that have an emission probability of 0. Also, it contains certain cellplan parameters, making it a very large dataframe for each network scenario.

• P.long.df: This object is a reduced version of P.long.complete.df, without zero elements and no cellplan parameters.

# specify the option of differing parameters for sig_d_th and max_overlapping_cells depending on the cell type in custom create_strength_llh function

# Workaround: securing that sig_d_th and max_overlapping_cells are the same for each layer
signal.strength.llh.param <- list(sig_d_th = 0.05,
                                  max_overlapping_cells = 100)

# define connection llh and classify tiles depending on coverage status
connection.llh.list <- map(cellplans.list, ~create_strength_llh_custom(signal.strength.comb.dt = .x$signal.strength.comb.dt,
                                                                       signal.strength.llh.param = signal.strength.llh.param,
                                                                       smart.rounding.digits = 3,
                                                                       area.df = area$area.df))

## Number of tiles which are unsufficiently covered: 0Number of tiles which are unsufficiently covered: 0Number of tiles which are unsufficiently covered: 0Number of tiles which are unsufficiently covered: 0

# create c-vector
C.vec.df.list <- map(connection.llh.list, ~create_c_vector(signal.strength.llh.combined = .x$signal.strength.llh.combined))

# develop long format of P matrix which also contains certain cellplan parameters and "zero elements"
# develop different versions of id variables (tile and antennas) for easier joining
P.long.complete.df <- map2(connection.llh.list, C.vec.df.list, ~full_join(.x$signal.strength.llh.combined, .y, by = "cell")) %>%
  map(~dplyr::select(., tile.id.chr,, tile.id.fac, tile.id.num, pop, cell, type, dist, pij, phones.sum)) %>%
  map(~mutate(., cell.chr = as.character(cell))) %>%
  map2(cellplans.list, ~mutate(.x, cell.fac = factor(cell.chr, levels = fct_unique(.y$cellplan.combined.df$cell.fac)))) %>%
  map(~mutate(., cell.num = as.numeric(cell.fac)))

# Long format of P matrix with minimal variables and and unique rows
P.long.df <- P.long.complete.df %>%
  map(~dplyr::select(., tile.id.chr, tile.id.fac, tile.id.num, cell.chr, cell.fac, cell.num, pij)) %>%
  map(~distinct(.))


gen.model.objects <- list(C.vec.df.list = C.vec.df.list,
                          P.long.complete.df = P.long.complete.df,
                          P.long.df = P.long.df)

# saveRDS(gen.model.objects, "Data/gen.model.objects.rds")

sessionInfo()

## R version 4.0.5 (2021-03-31)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19043)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252   
## [3] LC_MONETARY=German_Germany.1252 LC_NUMERIC=C                   
## [5] LC_TIME=German_Germany.1252    
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] DT_0.14              knitr_1.34           gridExtra_2.3       
##  [4] scattermore_0.7      ggpointdensity_0.1.0 ggrepel_0.8.2       
##  [7] viridis_0.5.1        viridisLite_0.3.0    ggthemes_4.2.0      
## [10] SpatialKWD_0.4.0     mobloc_0.5-1         Matrix_1.3-2        
## [13] stars_0.5-2          abind_1.4-5          raster_3.3-7        
## [16] sp_1.4-2             sf_1.0-3             data.table_1.13.2   
## [19] forcats_0.5.1        stringr_1.4.0        dplyr_1.0.5         
## [22] purrr_0.3.4          readr_1.4.0          tidyr_1.1.3         
## [25] tibble_3.1.2         ggplot2_3.3.3        tidyverse_1.3.1     
## 
## loaded via a namespace (and not attached):
##  [1] fs_1.5.0           lubridate_1.7.10   doParallel_1.0.16  httr_1.4.2        
##  [5] tools_4.0.5        backports_1.1.7    utf8_1.2.1         rgdal_1.5-12      
##  [9] R6_2.5.0           KernSmooth_2.23-18 DBI_1.1.1          colorspace_1.4-1  
## [13] withr_2.4.1        tidyselect_1.1.0   curl_4.3           compiler_4.0.5    
## [17] cli_2.5.0          rvest_1.0.0        xml2_1.3.2         labeling_0.4.2    
## [21] bookdown_0.21      scales_1.1.1       classInt_0.4-3     digest_0.6.25     
## [25] foreign_0.8-81     rmarkdown_2.5      rio_0.5.16         pkgconfig_2.0.3   
## [29] htmltools_0.5.1.1  highr_0.8          dbplyr_2.1.1       htmlwidgets_1.5.3 
## [33] rlang_0.4.10       readxl_1.3.1       rstudioapi_0.13    generics_0.1.0    
## [37] farver_2.0.3       jsonlite_1.7.2     zip_2.2.0          car_3.0-8         
## [41] magrittr_2.0.1     Rcpp_1.0.7         munsell_0.5.0      fansi_0.4.1       
## [45] lifecycle_1.0.0    stringi_1.4.6      yaml_2.2.1         carData_3.0-4     
## [49] parallel_4.0.5     crayon_1.4.1       lattice_0.20-41    cowplot_1.0.0     
## [53] haven_2.3.1        hms_1.1.0          pillar_1.6.1       ggpubr_0.4.0      
## [57] ggsignif_0.6.0     codetools_0.2-18   reprex_2.0.0       glue_1.4.1        
## [61] evaluate_0.14      modelr_0.1.8       vctrs_0.3.8        foreach_1.5.1     
## [65] cellranger_1.1.0   gtable_0.3.0       assertthat_0.2.1   openxlsx_4.1.5    
## [69] xfun_0.26          lwgeom_0.2-5       broom_0.7.6        e1071_1.7-3       
## [73] rstatix_0.6.0      class_7.3-18       iterators_1.0.13   units_0.6-7       
## [77] ellipsis_0.3.2