Merge pull request #70 from ropensci/feat/directional-correlation-delay

Feat/directional correlation delay

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

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)

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)
+}
+

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')

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",