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_outbreak.csv) 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_pts specifies the coordinate reference system (CRS) of the input coordinates (x and y), the default (crs_pts=NULL) sets the CRS to EPSG:4326, and input coordinates should be longtitude and lattitude. base_map(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, nbin = 50.
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. ngrid 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
calc.Rj() function calculate individual reproductive numbers given data (n-length vectors : t, x, y) and prespecified likelihood functions (lpdf_GI, lpdf_SP). adj.sp 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 calc.Rj(adj.sp = F) to calculate non-adjusted reproductive numbers, where only t and and lpdf_GI 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 image() 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. dt controls the time step size in the animation
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")
II. Polygon data with cases statistics
1. Preparation
Load demo dataset: COVID-19 outbreak
The demo data (covid_outbreak.csv) 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 (RandomPoints.csv). 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 three columns. date stands for onset data and TOWN stands for where case happened. SID is consist of TOWN and a number, to ensure each observations have an independent SID for sampling.
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. RndID is for independent RndID for sampling containing TOWN and a number, 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 (TPE_town) and case statistics (TOWN.df), 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 calc.Rj() for CRS of TWD97/TM2 to improve computing efficiency. Notice to set CRS_TWD97 to true.
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.4182 0.7415 0.9989 1.1766 17.7198 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.