Content
I. Point data with cases coordinates
II. Polygon data with cases statistics
- Preparation
- Descriptive plots
- Estimating reproductive numbers
- Visualizing town level reproductive numbers
I. Point data with cases coordinates
1. Preparation
Load demo dataset: dengue outbreak
The demo data (dengue) provided consists of the onset date and coordinates of 1,848 individual observations. The plot_epi() can draw the epidemic curve (incident cases over time) of this outbreak.
date long lat
1 2007-06-07 120.2165 23.02999
2 2007-06-09 120.2184 23.03016
3 2007-06-10 120.2187 23.03065
4 2007-06-10 120.2165 23.02999
5 2007-06-14 120.2081 22.91810
6 2007-06-15 120.2144 22.970042. Descriptive Plots
Several plotting tools are provided to visualize disease surveillance data.
Epidemic curve
Epidemic curve (time series of incident cases) is the most common way to characterize temporal pattern of an outbreak. Timestamp (here, the onset date) for each case is the only arguent needed for plot_epi().
Point pattern
Point pattern is the most direct and exact way to visualize spatial pattern of an outbreak. Coordinates of the cases are neccessary input for plot_points(). crs specifies the coordinate reference system (CRS) of the input coordinates (x and y), the default (crs = NULL) sets the CRS to EPSG:4326, and input coordinates should be longtitude and lattitude. basemap(optional) allows sf-multipolygon object to be inserted as base map.
Note that, here, we manually specify the boundary of the figure (or it will be decided by sporadic cases in the remote area) to get appropriate view of the main outbreak.
Hexagon-bining plot
However, point pattern can be overwhelmingly unreadable when the data points growing large. A solution “bining”, aggregate the point into specified bins. Here, we provide hexagon-bining function plot_hex() for the job. One can set number of bins on the x-axis (automatically determined on y-axis) to tune the resolution, for example, grid.n = 30.
Kernel density estimation
Also, one might need a smoothed overview of the clustering pattern of an outbreak. plot_kde() calculates the 2-d kernel density estimation, and plots the contour. grid.n controls the resolution in calculation.
3. Estimating Reproductive Numbers
Specify the parameters
Estimating spatial-adjusted reproductive numbers involves two important likelihood components that determine the transmission likelihood of each case-pair. Both of them should be chosen carefully based on prior knowledge about the disease. First is the probability distribution of generation interval, defined as the time interval between the infection (or onset) time of infector and infectee. In the case of dengue, we apply gamma distribution with mean = 20 days and variance = 9 day^2. We need to specify this distribution as a single-argument, log probability densiy function for the subsequent calculation for spatial-adjusted reproductive numbers. Another key component is the spatial weighting function. Here, we use an exponential distribution with mean tranmission distance = 125 m.
Visualize these two probability density functions.
Calculate spatial-adjusted reproductive numbers
Rj() function calculate individual reproductive numbers given data (n-length vectors : t, x, y) and prespecified likelihood functions (GI.pdf, SW.pdf). adjSP controls whether to calculate spatial-adjusted reproductive number (default = TRUE). This function returns a list containing the resulting individual reproductive numbers (n-length vector) and pair-wise transmission probability, an n*n upper triangular matrix (a sparse matrix of class “dtCMatrix”).
res_adj = Rj(t = dengue$date, x = dengue$long, y = dengue$lat,GI.pdf = lpdf_GI(),SW.pdf = lpdf_SW(), adjSP = T)
summary(res_adj$Rj) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.1631 0.5237 0.9995 1.2499 30.2761 Calculate non-adjusted reproductive numbers
We can use the same function Rj(adjSP = F) to calculate non-adjusted reproductive numbers, where only t and and GI.pdf are required
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.4659 0.8925 0.9995 1.2606 20.2298 Visualize pair-wise transmission likelihood
The pair-wise transmission probabilities can be visualized by plot_trans() method. Comparing two method, we can find that spatial-adjusted method produced irregular pattern for accounting the distance of each case-pair. On the other hand, non-adjusted produced very smooth pattern that is solely governed by the distribution of generation interval.
4. Visualizing Individual Reproductive Numbers
Time-varying reproductive numbers
plot_Rt() converts individual reproductive numbers into time-varying reproductive number which is frequently used to characterize the evolution in transmision potential of a whole outbreak. The shaded bars (only present in spatial-adjusted method) represent the variance in the transmissibility among individuals.
Spatial pattern
plot_points() and plot_hex() allow estimated individual reproductive numbers to be mapped to the color attribute of previous plots. Just specify the estimates in the Rj argument. Because the original scle of individual reproductive numbers is usually large, plot_points() will convert it to percentiles which is then mapped to the color and point size in the figure.
dengue$Rj = res_adj$Rj
dengue =dengue[dengue$date > as.Date("2007-09-15") & dengue$date < as.Date("2007-10-15"),]
plot_points(x = dengue$long, y= dengue$lat, Rj = dengue$Rj, basemap = Taiwan)plot_hex(), on the other hand, calculates the average individual reproductive numbers whithin each heaxgon, and maps them to color.
Spatio-temporal pattern
We use animation to visualize the spatio-temporal pattern of an outbreak’s transmision potential. animate_points() generation animation of evolving point pattern while animate_hex() is used for Hexagon-bining pattern.
animate_points(t = dengue$date, x = dengue$long, y= dengue$lat, Rj = dengue$Rj, interval = "week",basemap = Taiwan)
animate_grid(x=dengue$long, y=dengue$lat, t=dengue$date, Rj = dengue$Rj, grid.n=20, hex = T, basemap = Taiwan, interval="week")
animate_duration(x = dengue$long, y = dengue$lat,t = dengue$date,dt = 7,Rj = dengue$Rj, basemap = Taiwan)
II. Polygon data with cases statistics
1. Preparation
Load demo dataset: COVID-19 outbreak
The demo data (covid19) provided consists of the onset date and administrative district (TOWN) of 11,348 individual observations. This data contains COVID-19 local cases in Taipei Metropolitan Area between May to June from Taiwan CDC. Due to the spatial resolution of this data is town (polygon data), the method in the previous part cannot be used directly.
Therefore, our first goal is to convert polygon data to point coordinates. We import the administrative district data, and create random points according to populations in advance (you can use more detailed polygon census data to get more accurate population distribution), to treat as the real population distribution (RndPts). Further, through random selection, we can generate a set of case coordinates. If you are worried about the distortion of the sampling, you can use Monte Carlo method repeating sampling multiple times to get stable results.
Random sampling
covid contains two columns. date stands for onset data and region stands for where case happened.
region date
1 台北市萬華區 2021-05-01
2 台北市萬華區 2021-05-01
3 台北市萬華區 2021-05-01
4 台北市萬華區 2021-05-01
5 台北市信義區 2021-05-01
6 台北市中正區 2021-05-01RndPts contains four columns. region is for spatial region, x and y stand for coordinates of each random points. Here, due to processing data in Taiwan, these coordinates is based on CRS of TWD97/TM2 zone 121 (EPSG:3826).
region x y
1 基隆市信義區 328726.9 2780933
2 基隆市信義區 328941.2 2780944
3 基隆市信義區 328840.1 2780856
4 基隆市信義區 328840.9 2780815
5 基隆市信義區 328820.6 2780807
6 基隆市信義區 328797.6 2780919Calculate the number of cases and random points in each town. Through random sampling to get case point coordinates.
After getting the point coordinates, we can use previous method to calculate reproductive numbers.
2. Descriptive Plots
Choropleth map
Furthermore, joining the original town data (Taipei) and case statistics (covid), we can plot choropleth map of case distribution with the functions in tmap package. For more information about tmap.
3. Estimating Reproductive Numbers
Specify the parameters
Calculate spatial-adjusted reproductive numbers
In the same way, through specifying the parameters of generation interval and spatial weighting function, to calculate spatial-adjusted reproductive numbers. Here, we modified Rj() for CRS of TWD97/TM2 to improve computing efficiency. Notice to set unit_coord = "meter" to correspond CRS.
covid.R_adj=Rj(t = covid$date, x = covid$x, y = covid$y, GI.pdf = lpdf_GI(5.2,1.5), unit_coord = "meter")
summary(covid.R_adj$Rj) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.4218 0.7446 0.9989 1.1767 17.1406 Calculate non-adjusted reproductive numbers
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.6786 0.8106 0.9989 0.9967 6.3785 And then, we can sum up individual reproductive numbers to town level through ReduceMatrix() function.
4. Visualizing Town Level Reproductive Numbers
Time-varying reproductive numbers
plot_Rt() converts individual reproductive numbers into time-varying reproductive number which is frequently used to characterize the evolution in transmision potential of a whole outbreak. The shaded bars (only present in spatial-adjusted method) represent the variance in the transmissibility among individuals.
After summing up reproductive numbers to town level, the results mean the transitivity between towns. We can draw origin-destination (OD) map by following steps to visualize the transmission.