diff --git a/DESCRIPTION b/DESCRIPTION
index 64d98ec2..72d28fa9 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -23,7 +23,7 @@ Depends:
     R (>= 3.4)
 Imports: 
     adehabitatHR (>= 0.4.21),
-    data.table (>= 1.10.5),
+    data.table (>= 1.15.0),
     igraph,
     sf,
     lwgeom,
@@ -43,3 +43,4 @@ Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.3.2
 SystemRequirements: GDAL (>= 2.0.1), GEOS (>= 3.4.0), PROJ (>= 4.8.0),
     sqlite3
+Remotes: r-quantities/units@889cf39
diff --git a/NAMESPACE b/NAMESPACE
index fad615cb..42873b26 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -13,6 +13,7 @@ export(direction_to_leader)
 export(distance_to_centroid)
 export(distance_to_leader)
 export(dyad_id)
+export(edge_delay)
 export(edge_dist)
 export(edge_nn)
 export(fusion_id)
diff --git a/R/edge_delay.R b/R/edge_delay.R
new file mode 100644
index 00000000..a8f70ba3
--- /dev/null
+++ b/R/edge_delay.R
@@ -0,0 +1,228 @@
+#' Directional correlation delay based edge lists
+#'
+#' \code{edge_delay} returns edge lists defined by the directional correlation
+#' delay between individuals. The function expects a distance based edge list
+#' generated by \code{edge_dist} or \code{edge_nn}, a \code{data.table} with
+#' relocation data, individual identifiers and a window argument. The window
+#' argument is used to specify the temporal window within which to measure the
+#' directional correlation delay. Relocation data should be in two columns
+#' representing the X and Y coordinates.
+#'
+#' The \code{edges} and \code{DT} must be \code{data.table}s. If your data is a
+#' \code{data.frame}, you can convert it by reference using
+#' \code{\link[data.table:setDT]{data.table::setDT}}.
+#'
+#' The \code{edges} argument expects a distance based edge list generated with
+#' \code{edge_nn} or \code{edge_dist}. The \code{DT} argument expects relocation
+#' data with a timegroup column generated with \code{group_times}.
+#'
+#' The rows in \code{edges} and \code{DT} are internally matched in
+#' \code{edge_delay} using the columns \code{timegroup} (from
+#' \code{group_times}) and \code{ID1} and \code{ID2} (in \code{edges}, from
+#' \code{dyad_id}) with \code{id} (in \code{DT}). This function expects a
+#' \code{fusionID} present, generated with the \code{fusion_id} function, and a
+#' \code{dyadID} present, generated with the \code{dyad_id} function. The
+#' \code{id}, and \code{direction} arguments expect the names of a column in
+#' \code{DT} which correspond to the id, and direction columns.
+#'
+#' @inheritParams centroid_fusion
+#' @inheritParams direction_group
+#' @param window temporal window in unit of timegroup column generated with
+#'   \code{group_times}, eg. \code{window = 4} corresponds to the 4 timegroups
+#'   before and after the focal observation
+#'
+#' @return \code{edge_delay} returns the input \code{edges} appended with
+#'   a 'dir_corr_delay' column indicating the temporal delay (in units of
+#'   timegroups) at which ID1's direction of movement is most similar to
+#'   ID2's direction of movement, within the temporal window defined. For
+#'   example, if focal individual 'A' moves in a 45 degree direction at time 2
+#'   and individual 'B' moves in a most similar direction within the window
+#'   at time 5, the directional correlation delay between A and B is 3. Positive
+#'   values of directional correlation delay indicate a directed leadership
+#'   edge from ID1 to ID2.
+#'
+#' @export
+#'
+#' @references
+#'
+#' The directional correlation delay is defined in Nagy et al. 2010
+#' (<https://doi.org/10.1038/nature08891>).
+#'
+#' See examples of measuring the directional correlation delay:
+#'  * <https://doi.org/10.1016/j.anbehav.2013.07.005>
+#'  * <https://doi.org/10.1073/pnas.1305552110>
+#'  * <https://doi.org/10.1111/jfb.15315>
+#'  * <https://doi.org/10.1371/journal.pcbi.1003446>
+#'
+#' @family Edge-list generation
+#' @family Direction functions
+#'
+#' @examples
+#' # Load data.table
+#' library(data.table)
+#' \dontshow{data.table::setDTthreads(1)}
+#'
+#' # Read example data
+#' DT <- fread(system.file("extdata", "DT.csv", package = "spatsoc"))
+#'
+#' # Select only individuals A, B, C for this example
+#' DT <- DT[ID %in% c('A', 'B', 'C')]
+#'
+#' # Cast the character column to POSIXct
+#' DT[, datetime := as.POSIXct(datetime, tz = 'UTC')]
+#'
+#' # Temporal grouping
+#' group_times(DT, datetime = 'datetime', threshold = '20 minutes')
+#'
+#' # Calculate direction
+#' direction_step(
+#'   DT = DT,
+#'   id = 'ID',
+#'   coords = c('X', 'Y'),
+#'   projection = 32736
+#' )
+#'
+#' # Distance based edge list generation
+#' edges <- edge_dist(
+#'   DT,
+#'   threshold = 100,
+#'   id = 'ID',
+#'   coords = c('X', 'Y'),
+#'   timegroup = 'timegroup',
+#'   returnDist = TRUE,
+#'   fillNA = FALSE
+#' )
+#'
+#' # Generate dyad id
+#' dyad_id(edges, id1 = 'ID1', id2 = 'ID2')
+#'
+#' # Generate fusion id
+#' fusion_id(edges, threshold = 100)
+#'
+#' # Directional correlation delay
+#' delay <- edge_delay(
+#'   edges = edges,
+#'   DT = DT,
+#'   window = 3,
+#'   id = 'ID'
+#' )
+#'
+#' delay[, mean(dir_corr_delay, na.rm = TRUE), by = .(ID1, ID2)][V1 > 0]
+edge_delay <- function(
+    edges,
+    DT,
+    window = NULL,
+    id = NULL,
+    direction = 'direction') {
+  # due to NSE notes in R CMD check
+  . <- timegroup <- fusionID <- min_timegroup <- max_timegroup <-
+    delay_timegroup <- ID1  <- ID2 <- dir_corr_delay <- NULL
+
+  if (is.null(DT)) {
+    stop('input DT required')
+  }
+
+  if (is.null(edges)) {
+    stop('input edges required')
+  }
+
+  if (is.null(id)) {
+    stop('id column name required')
+  }
+
+  check_cols_edges <- c('ID1', 'ID2', 'timegroup')
+  if (any(!(check_cols_edges %in% colnames(edges)))) {
+    stop(paste0(
+      as.character(paste(setdiff(
+        check_cols_edges,
+        colnames(edges)
+      ), collapse = ', ')),
+      ' field(s) provided are not present in input DT'
+    ))
+  }
+
+  check_cols_DT <- c(id, 'timegroup', direction)
+  if (any(!(check_cols_DT %in% colnames(DT)
+  ))) {
+    stop(paste0(
+      as.character(paste(setdiff(
+        check_cols_DT,
+        colnames(DT)
+      ), collapse = ', ')),
+      ' field(s) provided are not present in input DT'
+    ))
+  }
+
+  if (is.null(window)) {
+    stop('window is required')
+  }
+
+  if (!is.numeric(window)) {
+    stop('window should be a numeric, in the units of timegroup')
+  }
+
+  if (!is.integer(DT$timegroup) || !is.integer(edges$timegroup)) {
+    stop('timegroup should be an integer, did you use group_times?')
+  }
+
+  if (!'fusionID' %in% colnames(edges)) {
+    stop('fusionID field not present in edges, did you run fusion_id?')
+  }
+
+  if (!'dyadID' %in% colnames(edges)) {
+    stop('dyadID field not present in edges, did you run dyad_id?')
+  }
+
+  # "Forward": all edges ID1 -> ID2
+  forward <- edges[!is.na(fusionID),
+                   data.table::first(.SD),
+                   by = .(fusionID, timegroup)]
+
+  forward[, min_timegroup :=
+            data.table::fifelse(timegroup - window < min(timegroup),
+                                min(timegroup),
+                                timegroup - window),
+          by = fusionID,
+          env = list(window = window)]
+
+  forward[, max_timegroup :=
+            data.table::fifelse(timegroup + window > max(timegroup),
+                                max(timegroup),
+                                timegroup + window),
+          by = fusionID,
+          env = list(window = window)]
+
+  forward[, delay_timegroup := {
+    focal_direction <- DT[timegroup == .BY$timegroup &
+                            id == ID1, direction]
+    DT[between(timegroup, min_timegroup, max_timegroup) & id == ID2,
+       timegroup[which.min(diff_rad(focal_direction, direction))]]
+  },
+  by = c('timegroup',  'dyadID'),
+  env = list(id = id, direction = direction)]
+
+  forward[, dir_corr_delay := delay_timegroup - timegroup]
+
+  data.table::set(forward,
+                  j = c('min_timegroup', 'max_timegroup','delay_timegroup'),
+                  value = NULL)
+
+  # "Reverse": replicate forward but reverse direction ID1 <- ID2
+  reverse <- data.table::copy(forward)
+  setnames(reverse, c('ID1', 'ID2'), c('ID2', 'ID1'))
+  reverse[, dir_corr_delay := - dir_corr_delay]
+
+  out <- data.table::rbindlist(list(
+    forward,
+    reverse
+  ), use.names = TRUE)
+
+  data.table::setorder(out, timegroup)
+  data.table::setcolorder(
+    out,
+    c('timegroup', 'ID1', 'ID2', 'dyadID', 'fusionID')
+  )
+
+  return(out)
+}
+
diff --git a/R/internal.R b/R/internal.R
new file mode 100644
index 00000000..5d3fc73b
--- /dev/null
+++ b/R/internal.R
@@ -0,0 +1,63 @@
+#' Difference of two angles measured in radians
+#'
+#' Internal function
+#'
+#' @param x angle in radians
+#' @param y angle in radians
+#' @param signed boolean if signed difference should be returned, default FALSE
+#' @param return_units return difference with units = 'rad'
+#'
+#' @return Difference between x and y in radians. If signed is TRUE, the signed difference is returned. If signed is FALSE, the absolute difference is returned. Note: The difference is the smallest difference, eg.
+#' @references adapted from https://stackoverflow.com/a/7869457
+#'
+#' @examples
+#' # Load data.table
+#' library(data.table)
+#' \dontshow{data.table::setDTthreads(1)}
+#'
+#' # Read example data
+#' DT <- fread(system.file("extdata", "DT.csv", package = "spatsoc"))
+#'
+#' # Cast the character column to POSIXct
+#' DT[, datetime := as.POSIXct(datetime, tz = 'UTC')]
+#'
+#' # Set order using data.table::setorder
+#' setorder(DT, datetime)
+#'
+#' # Calculate direction
+#' direction_step(
+#'   DT = DT,
+#'   id = 'ID',
+#'   coords = c('X', 'Y'),
+#'   projection = 32736
+#' )
+#'
+#' # Differences
+#' spatsoc:::diff_rad(DT[1, direction], DT[2, direction])
+diff_rad <- function(x, y,  signed = FALSE, return_units = FALSE) {
+  if (!inherits(x, 'units') || units(x)$numerator != 'rad') {
+    stop('units(x) is not radians')
+  }
+  if (!inherits(y, 'units') || units(y)$numerator != 'rad') {
+    stop('units(y) is not radians')
+  }
+
+  d <- units::drop_units(y) - units::drop_units(x)
+  d <- ((d + pi) %% (2 * pi)) - pi
+
+  if (signed) {
+    out <- d
+  } else {
+    out <- abs(d)
+  }
+
+  if (return_units) {
+    return(units::as_units(out, 'rad'))
+  } else {
+    return(out)
+  }
+}
+
+
+# Requires version with this PR merged
+# remotes::install_github('https://github.com/r-quantities/units/pull/365')
diff --git a/codemeta.json b/codemeta.json
index d484277d..139051b0 100644
--- a/codemeta.json
+++ b/codemeta.json
@@ -138,7 +138,7 @@
       "@type": "SoftwareApplication",
       "identifier": "data.table",
       "name": "data.table",
-      "version": ">= 1.10.5",
+      "version": ">= 1.15.0",
       "provider": {
         "@id": "https://cran.r-project.org",
         "@type": "Organization",
@@ -214,7 +214,7 @@
     },
     "SystemRequirements": "GDAL (>= 2.0.1), GEOS (>= 3.4.0), PROJ (>= 4.8.0),\n    sqlite3"
   },
-  "fileSize": "2023.039KB",
+  "fileSize": "2044.39KB",
   "citation": [
     {
       "@type": "ScholarlyArticle",
diff --git a/man/diff_rad.Rd b/man/diff_rad.Rd
new file mode 100644
index 00000000..a3b20235
--- /dev/null
+++ b/man/diff_rad.Rd
@@ -0,0 +1,51 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/internal.R
+\name{diff_rad}
+\alias{diff_rad}
+\title{Difference of two angles measured in radians}
+\usage{
+diff_rad(x, y, signed = FALSE, return_units = FALSE)
+}
+\arguments{
+\item{x}{angle in radians}
+
+\item{y}{angle in radians}
+
+\item{signed}{boolean if signed difference should be returned, default FALSE}
+
+\item{return_units}{return difference with units = 'rad'}
+}
+\value{
+Difference between x and y in radians. If signed is TRUE, the signed difference is returned. If signed is FALSE, the absolute difference is returned. Note: The difference is the smallest difference, eg.
+}
+\description{
+Internal function
+}
+\examples{
+# Load data.table
+library(data.table)
+\dontshow{data.table::setDTthreads(1)}
+
+# Read example data
+DT <- fread(system.file("extdata", "DT.csv", package = "spatsoc"))
+
+# Cast the character column to POSIXct
+DT[, datetime := as.POSIXct(datetime, tz = 'UTC')]
+
+# Set order using data.table::setorder
+setorder(DT, datetime)
+
+# Calculate direction
+direction_step(
+  DT = DT,
+  id = 'ID',
+  coords = c('X', 'Y'),
+  projection = 32736
+)
+
+# Differences
+spatsoc:::diff_rad(DT[1, direction], DT[2, direction])
+}
+\references{
+adapted from https://stackoverflow.com/a/7869457
+}
diff --git a/man/direction_group.Rd b/man/direction_group.Rd
index 560d441b..c70ff450 100644
--- a/man/direction_group.Rd
+++ b/man/direction_group.Rd
@@ -88,6 +88,7 @@ See examples of using mean group direction:
 Other Direction functions: 
 \code{\link{direction_polarization}()},
 \code{\link{direction_step}()},
-\code{\link{direction_to_leader}()}
+\code{\link{direction_to_leader}()},
+\code{\link{edge_delay}()}
 }
 \concept{Direction functions}
diff --git a/man/direction_polarization.Rd b/man/direction_polarization.Rd
index 4db6d66e..8c0b66ca 100644
--- a/man/direction_polarization.Rd
+++ b/man/direction_polarization.Rd
@@ -88,6 +88,7 @@ See examples of using polarization:
 Other Direction functions: 
 \code{\link{direction_group}()},
 \code{\link{direction_step}()},
-\code{\link{direction_to_leader}()}
+\code{\link{direction_to_leader}()},
+\code{\link{edge_delay}()}
 }
 \concept{Direction functions}
diff --git a/man/direction_step.Rd b/man/direction_step.Rd
index bfb5d1b4..20dcc3e5 100644
--- a/man/direction_step.Rd
+++ b/man/direction_step.Rd
@@ -112,6 +112,7 @@ example[, .(direction, units::set_units(direction, 'degree'))]
 Other Direction functions: 
 \code{\link{direction_group}()},
 \code{\link{direction_polarization}()},
-\code{\link{direction_to_leader}()}
+\code{\link{direction_to_leader}()},
+\code{\link{edge_delay}()}
 }
 \concept{Direction functions}
diff --git a/man/direction_to_leader.Rd b/man/direction_to_leader.Rd
index b200361f..65a02d58 100644
--- a/man/direction_to_leader.Rd
+++ b/man/direction_to_leader.Rd
@@ -105,6 +105,7 @@ See examples of using direction to leader and position within group:
 Other Direction functions: 
 \code{\link{direction_group}()},
 \code{\link{direction_polarization}()},
-\code{\link{direction_step}()}
+\code{\link{direction_step}()},
+\code{\link{edge_delay}()}
 }
 \concept{Direction functions}
diff --git a/man/edge_delay.Rd b/man/edge_delay.Rd
new file mode 100644
index 00000000..373a58cb
--- /dev/null
+++ b/man/edge_delay.Rd
@@ -0,0 +1,140 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/edge_delay.R
+\name{edge_delay}
+\alias{edge_delay}
+\title{Directional correlation delay based edge lists}
+\usage{
+edge_delay(edges, DT, window = NULL, id = NULL, direction = "direction")
+}
+\arguments{
+\item{edges}{edge list generated generated by \code{edge_dist} or
+\code{edge_nn}, with fusionID column generated by \code{fusion_id}}
+
+\item{DT}{input data.table with timegroup column generated with
+\code{group_times} matching the input data.table used to generate the edge
+list with \code{edge_nn} or \code{edge_dist}}
+
+\item{window}{temporal window in unit of timegroup column generated with
+\code{group_times}, eg. \code{window = 4} corresponds to the 4 timegroups
+before and after the focal observation}
+
+\item{id}{character string of ID column name}
+
+\item{direction}{character string of direction column name, default
+"direction"}
+}
+\value{
+\code{edge_delay} returns the input \code{edges} appended with
+a 'dir_corr_delay' column indicating the temporal delay (in units of
+timegroups) at which ID1's direction of movement is most similar to
+ID2's direction of movement, within the temporal window defined. For
+example, if focal individual 'A' moves in a 45 degree direction at time 2
+and individual 'B' moves in a most similar direction within the window
+at time 5, the directional correlation delay between A and B is 3. Positive
+values of directional correlation delay indicate a directed leadership
+edge from ID1 to ID2.
+}
+\description{
+\code{edge_delay} returns edge lists defined by the directional correlation
+delay between individuals. The function expects a distance based edge list
+generated by \code{edge_dist} or \code{edge_nn}, a \code{data.table} with
+relocation data, individual identifiers and a window argument. The window
+argument is used to specify the temporal window within which to measure the
+directional correlation delay. Relocation data should be in two columns
+representing the X and Y coordinates.
+}
+\details{
+The \code{edges} and \code{DT} must be \code{data.table}s. If your data is a
+\code{data.frame}, you can convert it by reference using
+\code{\link[data.table:setDT]{data.table::setDT}}.
+
+The \code{edges} argument expects a distance based edge list generated with
+\code{edge_nn} or \code{edge_dist}. The \code{DT} argument expects relocation
+data with a timegroup column generated with \code{group_times}.
+
+The rows in \code{edges} and \code{DT} are internally matched in
+\code{edge_delay} using the columns \code{timegroup} (from
+\code{group_times}) and \code{ID1} and \code{ID2} (in \code{edges}, from
+\code{dyad_id}) with \code{id} (in \code{DT}). This function expects a
+\code{fusionID} present, generated with the \code{fusion_id} function, and a
+\code{dyadID} present, generated with the \code{dyad_id} function. The
+\code{id}, and \code{direction} arguments expect the names of a column in
+\code{DT} which correspond to the id, and direction columns.
+}
+\examples{
+# Load data.table
+library(data.table)
+\dontshow{data.table::setDTthreads(1)}
+
+# Read example data
+DT <- fread(system.file("extdata", "DT.csv", package = "spatsoc"))
+
+# Select only individuals A, B, C for this example
+DT <- DT[ID \%in\% c('A', 'B', 'C')]
+
+# Cast the character column to POSIXct
+DT[, datetime := as.POSIXct(datetime, tz = 'UTC')]
+
+# Temporal grouping
+group_times(DT, datetime = 'datetime', threshold = '20 minutes')
+
+# Calculate direction
+direction_step(
+  DT = DT,
+  id = 'ID',
+  coords = c('X', 'Y'),
+  projection = 32736
+)
+
+# Distance based edge list generation
+edges <- edge_dist(
+  DT,
+  threshold = 100,
+  id = 'ID',
+  coords = c('X', 'Y'),
+  timegroup = 'timegroup',
+  returnDist = TRUE,
+  fillNA = FALSE
+)
+
+# Generate dyad id
+dyad_id(edges, id1 = 'ID1', id2 = 'ID2')
+
+# Generate fusion id
+fusion_id(edges, threshold = 100)
+
+# Directional correlation delay
+delay <- edge_delay(
+  edges = edges,
+  DT = DT,
+  window = 3,
+  id = 'ID'
+)
+
+delay[, mean(dir_corr_delay, na.rm = TRUE), by = .(ID1, ID2)][V1 > 0]
+}
+\references{
+The directional correlation delay is defined in Nagy et al. 2010
+(\url{https://doi.org/10.1038/nature08891}).
+
+See examples of measuring the directional correlation delay:
+\itemize{
+\item \url{https://doi.org/10.1016/j.anbehav.2013.07.005}
+\item \url{https://doi.org/10.1073/pnas.1305552110}
+\item \url{https://doi.org/10.1111/jfb.15315}
+\item \url{https://doi.org/10.1371/journal.pcbi.1003446}
+}
+}
+\seealso{
+Other Edge-list generation: 
+\code{\link{edge_dist}()},
+\code{\link{edge_nn}()}
+
+Other Direction functions: 
+\code{\link{direction_group}()},
+\code{\link{direction_polarization}()},
+\code{\link{direction_step}()},
+\code{\link{direction_to_leader}()}
+}
+\concept{Direction functions}
+\concept{Edge-list generation}
diff --git a/man/edge_dist.Rd b/man/edge_dist.Rd
index c0f78a64..34c48bb0 100644
--- a/man/edge_dist.Rd
+++ b/man/edge_dist.Rd
@@ -117,6 +117,7 @@ edges <- edge_dist(
 }
 \seealso{
 Other Edge-list generation: 
+\code{\link{edge_delay}()},
 \code{\link{edge_nn}()}
 }
 \concept{Edge-list generation}
diff --git a/man/edge_nn.Rd b/man/edge_nn.Rd
index d19c71a3..c17ce380 100644
--- a/man/edge_nn.Rd
+++ b/man/edge_nn.Rd
@@ -120,6 +120,7 @@ edge_nn(DT, id = 'ID', coords = c('X', 'Y'),
 }
 \seealso{
 Other Edge-list generation: 
+\code{\link{edge_delay}()},
 \code{\link{edge_dist}()}
 }
 \concept{Edge-list generation}
diff --git a/tests/testthat/test-diff-rad.R b/tests/testthat/test-diff-rad.R
new file mode 100644
index 00000000..633ac8d7
--- /dev/null
+++ b/tests/testthat/test-diff-rad.R
@@ -0,0 +1,34 @@
+# Test diff_rad
+context('test diff_rad')
+
+library(units)
+
+pt_01 <- as_units(0.1, 'rad')
+pt_02 <- as_units(0.2, 'rad')
+pi_rad <- as_units(pi, 'rad')
+
+# Adapted from: https://gist.github.com/bradphelan/7fe21ad8ebfcb43696b8
+expect_equal(diff_rad(pt_01, pt_02, TRUE, TRUE), pt_01)
+expect_equal(diff_rad(pt_01, pt_02 + 2 * pi_rad, TRUE, TRUE), pt_01)
+expect_equal(diff_rad(pt_01, pt_02 - 2 * pi_rad, TRUE, TRUE), pt_01)
+expect_equal(diff_rad(pt_01 + 2 * pi_rad, pt_02, TRUE, TRUE), pt_01)
+expect_equal(diff_rad(pt_01 - 2 * pi_rad, pt_02, TRUE, TRUE), pt_01)
+expect_equal(diff_rad(pt_02, pt_01, TRUE, TRUE), -pt_01)
+expect_equal(diff_rad(pt_02 + 2 * pi_rad, pt_01, TRUE, TRUE), -pt_01)
+expect_equal(diff_rad(pt_02 - 2 * pi_rad, pt_01, TRUE, TRUE), -pt_01)
+expect_equal(diff_rad(pt_02, pt_01 + 2 * pi_rad, TRUE, TRUE), -pt_01)
+expect_equal(diff_rad(pt_02, pt_01 - 2 * pi_rad, TRUE, TRUE), -pt_01)
+
+expect_equal(diff_rad(pt_01, pt_02, FALSE, TRUE), pt_01)
+expect_equal(diff_rad(pt_01, pt_02 + 2 * pi_rad, FALSE, TRUE), pt_01)
+expect_equal(diff_rad(pt_01, pt_02 - 2 * pi_rad, FALSE, TRUE), pt_01)
+expect_equal(diff_rad(pt_01 + 2 * pi_rad, pt_02, FALSE, TRUE), pt_01)
+expect_equal(diff_rad(pt_01 - 2 * pi_rad, pt_02, FALSE, TRUE), pt_01)
+expect_equal(diff_rad(pt_02, pt_01, FALSE, TRUE), pt_01)
+expect_equal(diff_rad(pt_02 + 2 * pi_rad, pt_01, FALSE, TRUE), pt_01)
+expect_equal(diff_rad(pt_02 - 2 * pi_rad, pt_01, FALSE, TRUE), pt_01)
+expect_equal(diff_rad(pt_02, pt_01 + 2 * pi_rad, FALSE, TRUE), pt_01)
+expect_equal(diff_rad(pt_02, pt_01 - 2 * pi_rad, FALSE, TRUE), pt_01)
+
+expect_error(diff_rad(pt_01, 0.2), 'rad')
+expect_error(diff_rad(0.1, pt_02), 'rad')
diff --git a/tests/testthat/test-edge-delay.R b/tests/testthat/test-edge-delay.R
new file mode 100644
index 00000000..0d68de7b
--- /dev/null
+++ b/tests/testthat/test-edge-delay.R
@@ -0,0 +1,182 @@
+# Test edge_delay
+context('test edge_delay')
+
+library(spatsoc)
+
+DT <- fread('../testdata/DT.csv')
+id <- 'ID'
+datetime <- 'datetime'
+timethreshold <- '20 minutes'
+threshold <- 50
+coords <- c('X', 'Y')
+timegroup <- 'timegroup'
+group <- 'group'
+projection <- 32736
+window <- 3
+
+
+DT[, datetime := as.POSIXct(datetime, tz = 'UTC')]
+group_times(DT, datetime = datetime, threshold = timethreshold)
+direction_step(DT, id, coords, projection)
+edges <- edge_dist(DT, threshold = threshold, id = id,
+                   coords = coords, timegroup = timegroup,
+                   returnDist = TRUE, fillNA = FALSE)
+dyad_id(edges, id1 = 'ID1', id2 = 'ID2')
+fusion_id(edges, threshold = threshold)
+
+clean_DT <- copy(DT)
+clean_edges <- copy(edges)
+
+# edge_delay(DT = DT, edges = edges, id = id, window = window)
+
+test_that('edges, DT are required', {
+  expect_error(edge_delay(edges, DT = NULL))
+  expect_error(edge_delay(edges = NULL, DT))
+})
+
+test_that('arguments required, otherwise error detected', {
+  expect_error(edge_delay(edges, DT, id = NULL),
+               'id column name required')
+  expect_error(edge_delay(edges, DT, id = id, window = NULL),
+               'window is required')
+})
+
+test_that('window is numeric, timegroup is integer', {
+  expect_error(edge_delay(edges, DT, id = id, window = 'potato'),
+               'numeric')
+  copy_edges <- copy(clean_edges)
+  copy_edges[, timegroup := as.character(timegroup)]
+  expect_error(edge_delay(copy_edges, DT, id = id, window = 2),
+               'integer')
+
+  copy_DT <- copy(clean_DT)
+  copy_DT[, timegroup := as.character(timegroup)]
+  expect_error(edge_delay(edges, copy_DT, id = id, window = 2),
+               'integer')
+})
+
+test_that('column names must exist in DT', {
+  expect_error(edge_delay(edges, DT, id = 'potato'),
+               'potato field')
+
+  copy_edges <- copy(clean_edges)
+  copy_edges[, timegroup := NULL]
+  expect_error(edge_delay(copy_edges, DT, id = id, window = window),
+               'timegroup field')
+
+  copy_edges <- copy(clean_edges)
+  copy_edges[, fusionID := NULL]
+  expect_error(edge_delay(copy_edges, DT, id = id, window = window),
+               'fusionID field')
+
+  copy_edges <- copy(clean_edges)
+  copy_edges[, dyadID := NULL]
+  expect_error(edge_delay(copy_edges, DT, id = id, window = window),
+               'dyadID field')
+
+  expect_error(edge_delay(edges, DT, id = id, window = window,
+                          direction = 'potato'),
+               'potato field')
+
+  copy_DT <- copy(clean_DT)
+  copy_DT[, timegroup := NULL]
+  expect_error(edge_delay(edges, copy_DT, id = id, window = window),
+               'timegroup field')
+})
+
+test_that('no rows are added to the result edges', {
+  expect_equal(nrow(edges),
+               nrow(edge_delay(edges, DT, id = id, window = window)))
+})
+
+test_that('column added to the result DT', {
+  copyEdges <- copy(edges)
+
+  expected_cols <- c(colnames(copyEdges), 'dir_corr_delay')
+
+  expect_setequal(expected_cols,
+                  colnames(edge_delay(edges, DT, id = id, window = window)))
+  expect_equal(ncol(edge_delay(edges, DT, id = id, window = window)),
+               length(expected_cols))
+
+})
+
+test_that('column added to the result DT is integer', {
+  expect_type(edge_delay(edges, DT, id = id, window = window)$dir_corr_delay, 'integer')
+})
+
+test_that('returns a data.table', {
+  expect_s3_class(edge_delay(edges, DT, id = id, window = window), 'data.table')
+})
+
+test_that('direction colname can different than default', {
+  copyDT <- copy(clean_DT)
+  setnames(copyDT, 'direction', 'potato')
+  expect_s3_class(edge_delay(edges, copyDT, id = id, window = window,
+                             direction = 'potato'),
+                  'data.table')
+})
+
+test_that('window column in edge, DT does not influence results', {
+  copyDT <- copy(clean_DT)
+  copyDT[, window := 42]
+
+  copyEdges <- copy(clean_edges)
+  copyEdges[, window := 42]
+
+  expect_equal(
+    edge_delay(edges, copyDT, id = id, window = window),
+    edge_delay(edges, DT, id = id, window = window)
+  )
+
+  expect_equal(
+    edge_delay(copyEdges, DT, id = id,
+               window = window)[, .SD, .SDcols = -'window'],
+    edge_delay(edges, DT, id = id, window = window)
+  )
+
+  expect_equal(
+    edge_delay(copyEdges, copyDT, id = id,
+               window = window)[, .SD, .SDcols = -'window'],
+    edge_delay(edges, DT, id = id, window = window)
+  )
+
+})
+
+
+N_id <- 5
+N_seq <- 10
+seq_xy <- c(seq(0, 5, length.out = N_seq / 2),
+            seq(5.1, 0, length.out = N_seq / 2))
+DT_expect <- data.table(
+  X = rep(seq_xy, each = N_id),
+  Y = rep(seq_xy, each = N_id),
+  ID = LETTERS[seq.int(N_id)]
+)
+DT_expect[, timegroup := seq.int(.GRP)[.GRP] + seq.int(.N), by = ID]
+setorder(DT_expect, timegroup)
+direction_step(DT_expect, id, coords, projection = 4326)
+DT_expect
+
+edges_expect <- edge_dist(DT_expect, threshold = 100, id, coords, timegroup,
+                          returnDist = TRUE)
+dyad_id(edges_expect, 'ID1', 'ID2')
+fusion_id(edges_expect)
+
+window <- 5
+delay <- edge_delay(edges_expect, DT_expect, window = window, id = id)
+
+test_that('expected results are returned', {
+  expect_lte(nrow(delay), nrow(edges_expect))
+  expect_lte(nrow(DT_expect), nrow(delay))
+
+  mean_delays <- delay[, mean(dir_corr_delay, na.rm = TRUE), by = ID1]
+
+  expect_equal(mean_delays[V1 == min(V1), ID1], LETTERS[N_id])
+  expect_equal(mean_delays[V1 == median(V1), ID1], LETTERS[ceiling(N_id / 2)])
+  expect_equal(mean_delays[V1 == max(V1), ID1], LETTERS[1])
+
+  mean_delays_wrt_A <- delay[ID1 == 'A', mean(dir_corr_delay, na.rm = TRUE),
+                             by = ID2]
+  expect_equal(mean_delays_wrt_A[V1 == max(V1), ID2], LETTERS[N_id])
+})