Events recorded on a network often have a temporal dimension. In that context, one could estimate the density of events in both time and network spaces.

The spatio-temporal kernel is calculated as the product of the network kernel density and the time kernel density. For a sample point at location *l* and time *t*, the Temporal Network Kernel Density Estimate (TNKDE) is calculated as follows:

\[tnkde(l,t) = \frac{1}{bw_{net} * bw_{time}} * \sum^n_{i=1}(k_{net}(d(l,i_l),bw_{net})*k_{time}(d(t,i_t),bw_{time}))\]

with:

- \(k_{net}\) the network kernel density function
- \(k_{time}\) the time kernel density function
- \(bw_{net}\) the bandwidth used for the network kernel
- \(bw_{time}\) the bandwidth used for the time kernel
*n*the number of events*i*an event at location \(i_l\) and time \(i_t\)- \(d(a,b)\) the distance (in time or space) between
*a*and*b*

As specified in the previous formula, two bandwidths are necessary, one for space and one for time.

We give in this vignette a short example with the bike accidents data in 2016 in the central districts of Montreal.

We start here by exploring the density of events in time.

```
# first load data and packages
library(spNetwork)
library(tmap)
library(sf)
data(bike_accidents)
# converting the Date field to a numeric field (counting days)
$Time <- as.POSIXct(bike_accidents$Date, format = "%Y/%m/%d")
bike_accidents<- as.POSIXct("2016/01/01", format = "%Y/%m/%d")
start $Time <- difftime(bike_accidents$Time, start, units = "days")
bike_accidents$Time <- as.numeric(bike_accidents$Time)
bike_accidents
<- as.character(1:12)
months <- ifelse(nchar(months)==1, paste0("0", months), months)
months <- paste("2016/",months,"/01", sep = "")
months_starts_labs <- as.POSIXct(months_starts_labs, format = "%Y/%m/%d")
months_starts_num <- difftime(months_starts_num, start, units = "days")
months_starts_num <- as.numeric(months_starts_num)
months_starts_num <- gsub("2016/", "", months_starts_labs, fixed = TRUE)
months_starts_labs
ggplot(bike_accidents) +
geom_histogram(aes(x = Time), bins = 30, color = "white") +
scale_x_continuous(breaks = months_starts_num, labels = months_starts_labs)
```

It is not surprising to observe that most of the accidents occur during spring, summer, and fall. We will remove the lonely observations during the first three months of the year because they might influence the bandwidthsâ€™ size latter.

`<- subset(bike_accidents, bike_accidents$Time >= 90) bike_accidents `

We can now calculate the kernel density values in time for several bandwidths.

```
<- rep(1,nrow(bike_accidents))
w <- seq(0, max(bike_accidents$Time), 0.5)
samples
<- data.frame(
time_kernel_values bw_10 = tkde(bike_accidents$Time, w = w, samples = samples, bw = 10, kernel_name = "quartic"),
bw_20 = tkde(bike_accidents$Time, w = w, samples = samples, bw = 20, kernel_name = "quartic"),
bw_30 = tkde(bike_accidents$Time, w = w, samples = samples, bw = 30, kernel_name = "quartic"),
bw_40 = tkde(bike_accidents$Time, w = w, samples = samples, bw = 40, kernel_name = "quartic"),
bw_50 = tkde(bike_accidents$Time, w = w, samples = samples, bw = 50, kernel_name = "quartic"),
bw_60 = tkde(bike_accidents$Time, w = w, samples = samples, bw = 60, kernel_name = "quartic"),
time = samples
)
<- reshape2::melt(time_kernel_values,id.vars = "time")
df_time $variable <- as.factor(df_time$variable)
df_time
ggplot(data = df_time) +
geom_line(aes(x = time, y = value)) +
scale_x_continuous(breaks = months_starts_num, labels = months_starts_labs) +
facet_wrap(vars(variable), ncol=2, scales = "free") +
theme(axis.text = element_text(size = 5))
```

It seems that a bandwidth between 30 and 40 days capture the bimodal shape of the temporal dimension of bike accidents. The lower densities during July and August are likely caused by the lower traffic due to holidays.

It is also possible to use the classical functions from R to find a bandwidth with a data-driven approach.

```
<- bw.bcv(bike_accidents$Time, nb = 1000, lower = 1, upper = 80)
bw1 <- bw.ucv(bike_accidents$Time, nb = 1000, lower = 1, upper = 80)
bw2 <- bw.SJ(bike_accidents$Time, nb = 1000, lower = 1, upper = 80)
bw3
<- data.frame(
time_kernel_values bw_bcv = tkde(bike_accidents$Time, w = w, samples = samples, bw = bw1, kernel_name = "quartic"),
bw_ucv = tkde(bike_accidents$Time, w = w, samples = samples, bw = bw2, kernel_name = "quartic"),
bw_SJ = tkde(bike_accidents$Time, w = w, samples = samples, bw = bw3, kernel_name = "quartic"),
time = samples
)
<- reshape2::melt(time_kernel_values,id.vars = "time")
df_time $variable <- as.factor(df_time$variable)
df_time
ggplot(data = df_time) +
geom_line(aes(x = time, y = value)) +
scale_x_continuous(breaks = months_starts_num, labels = months_starts_labs) +
facet_wrap(vars(variable), ncol=2, scales = "free") +
theme(axis.text = element_text(size = 5))
```

In that case, the automatic bandwidth selection methods yield much noisier results with small bandwidths. Bike accidents are rare events and too small time bandwidths would create results with meaningless hot spots.

Before considering the spatio-temporal case, we can also investigate the spatial dimension.

```
# loading the road network
data(mtl_network)
tm_shape(mtl_network) +
tm_lines(col = "black") +
tm_shape(bike_accidents) +
tm_dots(col = "red", size = 0.1)
```

As suggested in the vignette NKDE, we will use an adaptive bandwidth of 450 metres for the discontinuous kernel.

```
# creating sample points
<- lixelize_lines(mtl_network, 50)
lixels <- lines_center(lixels)
sample_points
# calculating the densities
<- nkde(lines = mtl_network,
nkde_densities events = bike_accidents,
w = rep(1,nrow(bike_accidents)),
samples = sample_points,
kernel_name = "quartic",
bw = 450,
adaptive = TRUE, trim_bw = 900,
method = "discontinuous",
div = "bw",
max_depth = 10,
digits = 2, tol = 0.1, agg = 5,
grid_shape = c(1,1),
verbose = FALSE)
$density <- nkde_densities$k * 1000 sample_points
```

```
tm_shape(sample_points) +
tm_dots(col = "density", style = "kmeans", n = 8, palette = "viridis", size = 0.05) +
tm_layout(legend.outside = TRUE)
```

Several hot spots could already be identified, but we have no information about their location in time. The previous map represents the average density during the full period.

We can now estimate the kernel density in both space and time. Note that increasing the dimension leads to larger bandwidths.

```
<- bw_tnkde_cv_likelihood_calc(
cv_scores bw_net_range = c(100,1000),
bw_net_step = 100,
bw_time_range = c(10,60),
bw_time_step = 10,
lines = mtl_network,
events = bike_accidents,
time_field = "Time",
w = rep(1, length(bike_accidents)),
kernel_name = "quartic",
method = "discontinuous",
diggle_correction = FALSE,
study_area = NULL,
max_depth = 10,
digits = 2,
tol = 0.1,
agg = 15,
sparse=TRUE,
grid_shape=c(1,1),
sub_sample=1,
verbose = FALSE,
check = TRUE)
```

`::kable(cv_scores) knitr`

10 | 20 | 30 | 40 | 50 | 60 | |
---|---|---|---|---|---|---|

100 | -685.83959 | -653.23098 | -643.02769 | -622.70860 | -596.24179 | -559.70872 |

200 | -587.79223 | -489.93987 | -453.39790 | -408.76651 | -366.04846 | -339.61514 |

300 | -491.89552 | -367.76376 | -302.90136 | -240.13889 | -199.50585 | -159.13773 |

400 | -398.35016 | -233.82688 | -183.05089 | -134.51681 | -102.11291 | -79.89925 |

500 | -270.20408 | -134.31060 | -83.69580 | -59.46859 | -47.27205 | -31.17958 |

600 | -190.88013 | -99.59822 | -63.26827 | -41.08603 | -33.00040 | -28.99820 |

700 | -132.14032 | -71.17085 | -45.00551 | -32.92649 | -28.89854 | -24.96029 |

800 | -105.43875 | -52.89760 | -34.79517 | -28.83948 | -28.87864 | -24.96563 |

900 | -73.14577 | -42.66894 | -30.69225 | -26.80958 | -26.90015 | -25.01240 |

1000 | -66.75177 | -38.62339 | -28.70405 | -26.83371 | -26.94606 | -25.07264 |

According to the â€śleave one out cross validationâ€ť method, the optimal set of bandwidths is 700 metres and 60 days. As expected, larger bandwidths are required because the density of the events are spread both in space and time.

```
# choosing sample in times (every 10 days)
<- seq(0, max(bike_accidents$Time), 10)
sample_time
# calculating densities
<- tnkde(lines = mtl_network,
tnkde_densities events = bike_accidents,
time_field = "Time",
w = rep(1, nrow(bike_accidents)),
samples_loc = sample_points,
samples_time = sample_time,
kernel_name = "quartic",
bw_net = 700, bw_time = 60,
adaptive = TRUE,
trim_bw_net = 900,
trim_bw_time = 80,
method = "discontinuous",
div = "bw", max_depth = 10,
digits = 2, tol = 0.01,
agg = 15, grid_shape = c(1,1),
verbose = FALSE)
# creating a color palette for all the densities
library(classInt)
library(viridis)
<- c(tnkde_densities$k)
all_densities <- classIntervals(all_densities, n = 10, style = "kmeans")
color_breaks
# generating a map at each sample time
<- lapply(1:ncol(tnkde_densities$k), function(i){
all_maps <- sample_time[[i]]
time <- as.Date(start) + time
date
$density <- tnkde_densities$k[,i]
sample_points<- tm_shape(sample_points) +
map1 tm_dots(col = "density", size = 0.01,
breaks = color_breaks$brks, palette = viridis(10)) +
tm_layout(legend.show=FALSE, main.title = as.character(date), main.title.size = 0.5)
return(map1)
})
# creating a gif with all the maps
tmap_animation(all_maps, filename = "images/animated_map.gif",
width = 1000, height = 1000, dpi = 300, delay = 50)
```

`::include_graphics("images/animated_map.gif") knitr`

The locations of the hot spots are changing during the year. This information was hidden when only a spatial NKDE was used.

One can decide to use adaptive bandwidth for the TNKDE. They are calculated using the Abramsonâ€™s smoothing regimen. Two approaches are available: separated and simultaneous. In the first case, the network and temporal density at each event location are calculated separately and used to calculate the local bandwidth. Isolated events on the network will obtain larger network bandwidth and isolated events in time will obtain larger time bandwidth. The second method calculates the temporal-network density at each event location and uses it to adjust both bandwidths simultaneously. This method will always lead to stronger variation in the local bandwidths.

```
<- tnkde(lines = mtl_network,
tnkde_densities events = bike_accidents,
time_field = "Time",
w = rep(1, nrow(bike_accidents)),
samples_loc = sample_points,
samples_time = sample_time,
kernel_name = "quartic",
bw_net = 700, bw_time = 60,
adaptive = TRUE,
adaptive_separate = FALSE,
trim_bw_net = 900,
trim_bw_time = 80,
method = "discontinuous",
div = "bw", max_depth = 10,
digits = 2, tol = 0.01,
agg = 15, grid_shape = c(1,1),
verbose = FALSE)
```