| Title: | Estimate the spread-rate of epidemics |
|---|---|
| Description: | Estimate the local velocity of propagation of an epidemic event, given dates and locations of observed cases. The method estimates the surface of first date of invasion by interpolation of the earliest observations in a neighbourhood and derives the local spread rate as the inverse slope of the surface. Quantify the estimation uncertainty by a Monte Carlo approach. Visualise and simulate the spatio-temporal progress of epidemics. |
| Authors: | Facundo Muñoz [aut, cre] |
| Maintainer: | Facundo Muñoz <[email protected]> |
| License: | GPL-3 | file LICENSE |
| Version: | 0.1.13 |
| Built: | 2024-10-05 04:26:40 UTC |
| Source: | https://forgemia.inra.fr/umr-astre/spreadrate |
Compute an approximate angle that correspond to a given distance in both longitude and latitude.
dist2arc(x, dm)dist2arc(x, dm)
x |
A |
dm |
A target distance in m. |
The distance corresponding to an angle in longitude and latitude depends on the location in the earth and is typically different for both directions.
This function computes the angle that yields approximately the target distance when added to the longitude or the latitude by minimising the squared errors.
The average error is returned as an attribute.
A numeric value of angle in decimal degrees, such that a displacement in longitude and latitude by that magnitude corresponds to approximately the given distance.
require(sf) ## What angle yields 1 km at the equator? (and with which precision?) x0 <- st_sfc(st_point(c(0, 0)), crs = 4326) dist2arc(x0, 1e3) ## Note that this is independent of the longitude x1 <- st_sfc(st_point(c(90, 0)), crs = 4326) dist2arc(x1, 1e3) ## But is strongly dependent of the latitude x2 <- st_sfc(st_point(c(0, 45)), crs = 4326) x3 <- st_sfc(st_point(c(0, 75)), crs = 4326) dist2arc(x2, 1e3) dist2arc(x3, 1e3)require(sf) ## What angle yields 1 km at the equator? (and with which precision?) x0 <- st_sfc(st_point(c(0, 0)), crs = 4326) dist2arc(x0, 1e3) ## Note that this is independent of the longitude x1 <- st_sfc(st_point(c(90, 0)), crs = 4326) dist2arc(x1, 1e3) ## But is strongly dependent of the latitude x2 <- st_sfc(st_point(c(0, 45)), crs = 4326) x3 <- st_sfc(st_point(c(0, 75)), crs = 4326) dist2arc(x2, 1e3) dist2arc(x3, 1e3)
Given a set of observations, define a estimation template.
estimation_mask(x, buffer_size, res = buffer_size/8)estimation_mask(x, buffer_size, res = buffer_size/8)
x |
An object of class |
buffer_size |
Numeric. Expansion of the extent around the
convex-hull of |
res |
Numeric vector of length 1 or 2. Resolution of the raster template. Size(s) of the pixel, in metres. |
## Some random points in the unit square pp <- st_multipoint(matrix(runif(30), ncol = 2)) (em <- estimation_mask(pp, buffer_size = .3, res = .1)) plot(em) points(st_coordinates(pp))## Some random points in the unit square pp <- st_multipoint(matrix(runif(30), ncol = 2)) (em <- estimation_mask(pp, buffer_size = .3, res = .1)) plot(em) points(st_coordinates(pp))
Select a subset of the observation points at least tol m far
apart as representative of the emphneighbourhood. Assign the
earliest times observed in the neighbourhood to them.
filter_earliest_neigh(x)filter_earliest_neigh(x)
x |
sr_obs object. |
Note that the neighbouring distance must be expressed in m, even when the coordinates are geographical.
This function reproduces the operation for all the Monte Carlo replicates if any.
For the original dataset, the neighbouring tolerance parameter
is the mean value of the interval specified in its uq.
Another sr_obs object with a subset of the points and the earliest times observed in the neighbourhood of each.
d <- data.frame(lon = runif(30), lat = runif(30), date = 1:30) ## Use between 10 and 20 % of data diameter as neighbouring ## tolerance. The main result will use exactly 15% while the ## \code{mc} replicates will use random Gaussian values from ## that interval (at 99.9%) sruq <- sr_uq(10, 0, 1, neigh_tol = c(-10, -20)) sro <- sr_obs(d, "date", uq = sruq) srf <- filter_earliest_neigh(sro)d <- data.frame(lon = runif(30), lat = runif(30), date = 1:30) ## Use between 10 and 20 % of data diameter as neighbouring ## tolerance. The main result will use exactly 15% while the ## \code{mc} replicates will use random Gaussian values from ## that interval (at 99.9%) sruq <- sr_uq(10, 0, 1, neigh_tol = c(-10, -20)) sro <- sr_obs(d, "date", uq = sruq) srf <- filter_earliest_neigh(sro)
Fit a Thin-plate splines model to a set of observations
fit_surface(x)fit_surface(x)
x |
A |
An object of class Krig and Tps which includes
the fitted values and the model residuals. See
Tps. It also allows interpolation using
interpolate.
d <- data.frame(lon = runif(30), lat = runif(30)) d <- transform(d, date = round(10 * sqrt(lon**2 + lat**2))) sro <- sr_obs(d, "date") fm <- fit_surface(sro) summary(fm) plot(fm) # Some model diagnostics plots fields::surface(fm) # Quick image/contour plot of predicted surface.d <- data.frame(lon = runif(30), lat = runif(30)) d <- transform(d, date = round(10 * sqrt(lon**2 + lat**2))) sro <- sr_obs(d, "date") fm <- fit_surface(sro) summary(fm) plot(fm) # Some model diagnostics plots fields::surface(fm) # Quick image/contour plot of predicted surface.
Given a data.frame with two variables with longitudes and latitudes, return the variable names.
get_lonlatvars(x)get_lonlatvars(x)
x |
data.frame |
This function explores the names of the variables for the patterns 'lon' or 'long' and 'lat' ignoring case. If exactly one of each patterns is detected, it verifies that their magnitudes are within -365 and 365 for the longitude and within -90.1 and 90.1 for the latitude.
A character vector with the variable names for longitude and latitude, or an error.
d <- data.frame(lon = 1:10, lat = 1:10) get_lonlatvars(d)d <- data.frame(lon = 1:10, lat = 1:10) get_lonlatvars(d)
This function is used to derive local spread rate from a first-date-of-invasion surface
invslope(x, rmtop, bnd)invslope(x, rmtop, bnd)
x |
|
rmtop |
Numeric, between 0 and 100. Discard the |
bnd |
Numeric vector of length 2. Expected range of values. Everything beyond this range is filtered. |
To control for spurious numerical divergences, the user can remove
the largest rmtop
values beyond a reasonable interval bnd.
An object of the same class Raster* as the input.
## A flat surface with constant slope of 1 r_values <- replicate(100, seq(0.01, 1, length = 100)) r <- raster(r_values, crs = sp::CRS("+init=epsg:3857")) plot(r) invslope(r)## A flat surface with constant slope of 1 r_values <- replicate(100, seq(0.01, 1, length = 100)) r <- raster(r_values, crs = sp::CRS("+init=epsg:3857")) plot(r) invslope(r)
Shuffle observations in space and time.
mc_sample(x)mc_sample(x)
x |
An object of class |
Add an attribute mc with a list of datasets with some noise
added to the original coordinates and dates. The number of Monte
Carlo samples and the variablity of the noise is given by the
sr_uq attribute in the object.
d <- data.frame(lon = 1, lat = 1, date = 1) sr_obs(d, "date", uq = sr_uq(nsim = 3, space = 1, time = 1))d <- data.frame(lon = 1, lat = 1, date = 1) sr_obs(d, "date", uq = sr_uq(nsim = 3, space = 1, time = 1))
Get the value or interval of the Neighbourhood Tolerance Parameter
of a sr_obs object.
neigh_tol(x)neigh_tol(x)
x |
A |
This function retrieves the values set during creation with
sr_obs in absolute terms, performing the convertion from
values relative to the diameter of the dataset if necessary. See
sr_uq.
A numeric vector of size 1 or 2. Either a constant value for this parameter or a interval with distance values in m.
d <- data.frame(lon = 1:3, lat = 1:3, date = 1:3) ## Here, the neighbourhood tolerance parameter is set between ## 2.5 % of the dataset diameter and 1 km. ## Note that the resulting distances are in m, even if the ## original coordinates are geographic. obs <- sr_obs(d, "date", uq = sr_uq(neigh_tol = c(-2.5, 1e4))) neigh_tol(obs)d <- data.frame(lon = 1:3, lat = 1:3, date = 1:3) ## Here, the neighbourhood tolerance parameter is set between ## 2.5 % of the dataset diameter and 1 km. ## Note that the resulting distances are in m, even if the ## original coordinates are geographic. obs <- sr_obs(d, "date", uq = sr_uq(neigh_tol = c(-2.5, 1e4))) neigh_tol(obs)
Select a subset of points at a minimum tolerance distance from each other.
representative_points(x, dTolerance)representative_points(x, dTolerance)
x |
Object of class |
dTolerance |
Numeric. Distance tolerance in m. |
Object of class sf with a subset of the original
locations in x.
## 50 points in the unit square require(sf) x <- st_sfc(st_multipoint(matrix(runif(100), ncol = 2)), crs = 3857) representative_points(x, .5)## 50 points in the unit square require(sf) x <- st_sfc(st_multipoint(matrix(runif(100), ncol = 2)), crs = 3857) representative_points(x, .5)
Estimate the local spread rate of an epidemiological invasion.
sr(x, r = estimation_mask(x, buffer_size = max(st_distance(x))/10))sr(x, r = estimation_mask(x, buffer_size = max(st_distance(x))/10))
x |
A |
r |
A |
This function will compute spread-rate estimates for all of the
Monte Carlo samples in the dataset. It uses internally the function
future_map which will take advantage of
multiple processors or cluster access if you set a proper
plan beforehand. See examples.
require(sf) ## Randomly sample 100 point in the unit square x_mat <- matrix(runif(100 * 2), ncol = 2) ## Make it spatial point in a projected CRS x_sfc <- st_cast(st_sfc(st_multipoint(x_mat), crs = 3857), "POINT") ## Assign dates proportional to the squared distance to the origin ## So that spread-rate decreases linearly obs_dates <- sqrt(rowSums(x_mat**2)) sro <- sr_obs(st_sf(x_sfc, date = obs_dates), "date") ## Alternative setting with uncertainty quantification and an ## absolute specification of the neighbouring tolerance parameter. # sro2 <- sr_obs(st_sf(x_sfc, date = obs_dates), "date", # uq = sr_uq(6, .2, .2, .3)) ## Estimate local spread-rate ## in the units of the coordinates divided by the units of time sre <- sr(sro) plot(sre, col = hcl.colors(12)) points(st_coordinates(sro), pch = 19)require(sf) ## Randomly sample 100 point in the unit square x_mat <- matrix(runif(100 * 2), ncol = 2) ## Make it spatial point in a projected CRS x_sfc <- st_cast(st_sfc(st_multipoint(x_mat), crs = 3857), "POINT") ## Assign dates proportional to the squared distance to the origin ## So that spread-rate decreases linearly obs_dates <- sqrt(rowSums(x_mat**2)) sro <- sr_obs(st_sf(x_sfc, date = obs_dates), "date") ## Alternative setting with uncertainty quantification and an ## absolute specification of the neighbouring tolerance parameter. # sro2 <- sr_obs(st_sf(x_sfc, date = obs_dates), "date", # uq = sr_uq(6, .2, .2, .3)) ## Estimate local spread-rate ## in the units of the coordinates divided by the units of time sre <- sr(sro) plot(sre, col = hcl.colors(12)) points(st_coordinates(sro), pch = 19)
Estimate a surface of first-date of invasion, given a sparse set of observation locations and times.
sr_fdoi(x, estimation_mask)sr_fdoi(x, estimation_mask)
x |
A |
estimation_mask |
A |
This function first takes the earliest observed case in a neighbourhood,
as defined by the parameter neigh_tol in the definition of
x. Next, it interpolates the earliest observations with a
smooth surface, using a Thin-Plate splines model.
A RasterLayer with estimated dates of invasion at each pixel.
d <- data.frame(lon = runif(30), lat = runif(30)) d <- transform(d, date = round(10 * sqrt(lon**2 + lat**2))) sro <- sr_obs(d, "date") r <- raster(st_sf(sro), vals = 1) fdoi <- sr_fdoi(sro, r) plot(fdoi, col = hcl.colors(12)) points(d)d <- data.frame(lon = runif(30), lat = runif(30)) d <- transform(d, date = round(10 * sqrt(lon**2 + lat**2))) sro <- sr_obs(d, "date") r <- raster(st_sf(sro), vals = 1) fdoi <- sr_fdoi(sro, r) plot(fdoi, col = hcl.colors(12)) points(d)
Provide a dataset with the necessary meta-data for spread-rate estimation.
sr_obs(x, timevar, uq)sr_obs(x, timevar, uq)
x |
|
timevar |
Character. Variable name with observation times or dates. |
uq |
Object created with |
If x is a data.frame with geographical coordinates
(longitude and latitude) these are automatically recognised and
validated. See get_lonlatvars,
If your dataset is projected (no lon/lat variables) then first
build the sf object yourself using the corresponding
Coordinate Reference System. See example.
If uq is missing and x is a data.frame or
sf, a default value is generated with sr_uq() which
does not perform Uncertainty Quantification. See sr_uq.
If uq is missing and x is a sr_obs object,
This function will produce the Monte Carlo samples from the
dataset and the uq object, if necessary. It uses internally
the function future_map which will take advantage
of multiple processors or cluster access if you set a proper
plan beforehand. See examples.
Object of class sr_obs, which is a sf of type
POINT with complementary meta-data.
d_geo <- data.frame(lon = 1, lat = 1, date = 1) sr_obs(d_geo, "date") ## Projected coordinates d_prj <- data.frame(x = 1, y = 1, date = 1) d.sf <- st_as_sf(d_prj, coords = c("x", "y"), crs = 27561) sr_obs(d.sf, "date") ## Not run: ## Perform Monte Carlo samples in parallel library(furrr) plan(multiprocess) uq <- sr_uq(nsim = 3, space = 1, time = 1) sr_obs(d_geo, "date", uq = uq) ## End(Not run)d_geo <- data.frame(lon = 1, lat = 1, date = 1) sr_obs(d_geo, "date") ## Projected coordinates d_prj <- data.frame(x = 1, y = 1, date = 1) d.sf <- st_as_sf(d_prj, coords = c("x", "y"), crs = 27561) sr_obs(d.sf, "date") ## Not run: ## Perform Monte Carlo samples in parallel library(furrr) plan(multiprocess) uq <- sr_uq(nsim = 3, space = 1, time = 1) sr_obs(d_geo, "date", uq = uq) ## End(Not run)
Define the parameters for quantifying the uncertainty in the estimation of the spread-rate.
sr_uq(nsim = 1L, space = 0, time = 0, neigh_tol = -4.5)sr_uq(nsim = 1L, space = 0, time = 0, neigh_tol = -4.5)
nsim |
Integer > 0. Number of Monte Carlo replicates. |
space |
Numeric, non-negative. Uncertainty of spatial coordinates. |
time |
Numeric, non-negative. Uncertainty of temporal values. |
neigh_tol |
Number or interval (numeric vector of length 2). Neighouring-tolerance parameter in the units of the dataset coordinates if positive, or relative to the dataset diameter if negative. |
The actual observations will be The number of Monte Carlo replicates
The space and time arguments are expressed in the
data spatial and temporal units.
For each one of the nsim replicates, the observations will
be randomly shifted in space within a circle of radious
space, while the dates will be shifted within the interval
+- time with a discrete approximation to a Gaussian
curve.
The neighbouring-tolerance parameter neigh_tol is used to
filter the earliest observed cases in a neighbourhood. If
two points are closer than this value, they are considered roughly
in the same neighbourhood. The date of invasion of a location is
the earliest observed case in the neighbourhood. It can be defined
as a numeric value in the units of the coordinates, or as a range
(i.e. a numeric vector of length 2). Negative values will be
interpreted as a percentage of the diameter of the dataset.
An object of class sr_uq which is a list with the
given or default values, after some sanity checks.
## Use the default values: work with current observations ## without quantifying uncertainty uq_def <- sr_uq() ## Set the neighbouring-tolerance parameter only uq1 <- sr_uq(neigh_tol = 800) ## 10 MC replicates, shift locations within a circle of radious ## 1 km, shift dates by zero/one day up or down, consider points ## closer than 5% of the diameter of the dataset as in the same ## neighbourhood uq <- sr_uq(nsim = 10, space = 1e3, time = 1, neigh_tol = -5)## Use the default values: work with current observations ## without quantifying uncertainty uq_def <- sr_uq() ## Set the neighbouring-tolerance parameter only uq1 <- sr_uq(neigh_tol = 800) ## 10 MC replicates, shift locations within a circle of radious ## 1 km, shift dates by zero/one day up or down, consider points ## closer than 5% of the diameter of the dataset as in the same ## neighbourhood uq <- sr_uq(nsim = 10, space = 1e3, time = 1, neigh_tol = -5)
Transform the underlying sf object and all the MC replicas,
if any.
## S3 method for class 'sr_obs' st_transform(x, crs, ...)## S3 method for class 'sr_obs' st_transform(x, crs, ...)
x |
A |
crs |
Coordinate reference system: inteer with the EPSG code or character with proj4string |
... |
Ignored |
A sr_obs object in the target projection.
## Some observed data and MC replicates in longitude/latitude x <- sr_obs(data.frame(lat = 1, lon = 1, day = 1), "day", uq = sr_uq(2, 1, 1)) st_crs(x) lapply(attr(x, "mc"), st_crs) xt <- st_transform(x, 3857) st_crs(xt) lapply(attr(xt, "mc"), st_crs)## Some observed data and MC replicates in longitude/latitude x <- sr_obs(data.frame(lat = 1, lon = 1, day = 1), "day", uq = sr_uq(2, 1, 1)) st_crs(x) lapply(attr(x, "mc"), st_crs) xt <- st_transform(x, 3857) st_crs(xt) lapply(attr(xt, "mc"), st_crs)
Modify a vector of dates by a random quantity in [-e, e]. The distribution decreases exponentially, with rate 1/2
time_jitter(x, e)time_jitter(x, e)
x |
Numeric or Date vector. |
e |
Numeric. Semi-width of the jittering interval. |
Vector of the same size and type as x with jittered
values.
spreadrate:::time_jitter(1:10, 1) # Numeric input spreadrate:::time_jitter(Sys.Date() + 1:10, 1) # Date input ## Distribution of jittering x <- seq.int(1e3) xj <- spreadrate:::time_jitter(x, 5) barplot(table(x - xj))spreadrate:::time_jitter(1:10, 1) # Numeric input spreadrate:::time_jitter(Sys.Date() + 1:10, 1) # Date input ## Distribution of jittering x <- seq.int(1e3) xj <- spreadrate:::time_jitter(x, 5) barplot(table(x - xj))