First merge of dosimetry metrics

.gitignore

@@ -50,3 +50,6 @@ rsconnect/
 .Rproj.user
 docs
 inst/doc
+
+# DS_Store file
+.DS_Store

DESCRIPTION

@@ -8,6 +8,10 @@ Authors@R: c(
     person("Manuel", "Spitschan", 
       email = "[email protected]", role = "aut",
       comment = c(ORCID = "0000-0002-8572-9268")),
+    person("Steffen", "Hartmeyer", 
+      email = "[email protected]", role = "aut",
+      comment = c(ORCID = "0000-0002-2813-2668")),
+    person("MeLiDos", role = "fnd"),
     person("EURAMET", role = "fnd", comment = "European Association of National Metrology Institutes. Website: www.euramet.org. Grant Number: 22NRM05 MeLiDos. Grant Statement: The project (22NRM05 MeLiDos) has received funding from the European Partnership on Metrology, co-financed from the European Union’s Horizon Europe Research and Innovation Programme and by the Participating States."),    
     person("European Union", role = "fnd", comment = "Co-funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or EURAMET. Neither the European Union nor the granting authority can be held responsible for them."),
     person("TSCN-Lab", comment = c(URL = "www.tscnlab.org"), role = "cph"))
@@ -43,6 +47,7 @@ Imports:
     rsconnect,
     scales,
     shiny,
+    slider,
     stats,
     stringr,
     tibble,

NAMESPACE

@@ -7,11 +7,13 @@ export(Brown_rec)
 export(Datetime_breaks)
 export(Datetime_limits)
 export(aggregate_Datetime)
+export(bright_dark_period)
 export(count_difftime)
 export(create_Timedata)
 export(cut_Datetime)
 export(data2reference)
 export(dominant_epoch)
+export(duration_above_threshold)
 export(filter_Date)
 export(filter_Datetime)
 export(filter_Datetime_multiple)
@@ -25,10 +27,13 @@ export(gg_overview)
 export(import)
 export(import_Dataset)
 export(import_Statechanges)
+export(interdaily_stability)
 export(interval2state)
+export(intradaily_variability)
 export(join_datasets)
 export(sc2interval)
 export(sleep_int2Brown)
 export(symlog_trans)
+export(timing_above_threshold)
 importFrom(magrittr,"%>%")
 importFrom(rlang,":=")

R/helper.R

@@ -86,3 +86,35 @@ pick.grouping.columns <- function(dataset) {
     dplyr::group_vars(dataset)
   )
 }
+
+# Compare with threshold
+compare_threshold <- function(Light.vector,
+                      threshold,
+                      comparison = c("above", "below"),
+                      na.replace = FALSE){
+
+  comparison = match.arg(comparison)
+
+  stopifnot(
+    "`Light.vector` must be numeric!" = is.numeric(Light.vector),
+    "`threshold` must be numeric!" = is.numeric(threshold),
+    "`threshold` must be either one or two values!" = length(threshold) %in% c(1, 2),
+    "`na.replace` must be logical!" = is.logical(na.replace)
+  )
+
+  if(length(threshold) == 1){
+    out <- switch(comparison,
+                  "above" = Light.vector >= threshold,
+                  "below" = Light.vector <= threshold)
+  }
+  else{
+    threshold <- sort(threshold)
+    out <- Light.vector >= threshold[1] & Light.vector <= threshold[2]
+  }
+
+  if(na.replace){
+    out <- tidyr::replace_na(out, FALSE)
+  }
+
+  return(out)
+}

R/metric_bright_dark_period.R

@@ -0,0 +1,145 @@
+#' Brightest or darkest continuous period
+#'
+#' This function finds the brightest or darkest continuous period of a given
+#' timespan and calculates its `mean` light level, as well as the timing of the period's
+#' `onset`, `midpoint`, and `offset`. It is defined as the period with the maximum
+#' or minimum mean light level. Note that the data need to be regularly spaced 
+#' (i.e., no gaps) for correct results.
+#'
+#' @param Light.vector Numeric vector containing the light data.
+#' @param Time.vector Vector containing the time data. Can be HMS or POSIXct.
+#' @param period String indicating the type of period to look for. Can be either
+#'  `"brightest"`(the default) or `"darkest"`.
+#' @param timespan The timespan across which to calculate. Can be either a
+#'  `lubridate::duration()` or a `lubridate::duration()` string, e.g.,
+#'  `"1 day"` or `"10 sec"`.
+#' @param epoch The epoch at which the data was sampled. Can be either a
+#'  `lubridate::duration()` or a string. If it is a string, it needs to be
+#'  either `"dominant.epoch"` (the default) for a guess based on the data, or a valid
+#'  `lubridate::duration()` string, e.g., `"1 day"` or `"10 sec"`.
+#' @param loop Logical. Should the data be looped? If `TRUE`, a full copy of the data 
+#'    will be concatenated at the end of the data. Makes only sense for 24 h data.
+#'    Defaults to `FALSE`.
+#' @param na.rm Logical. Should missing values be removed for the calculation?
+#'    Defaults to `FALSE`.
+#' @param as.df Logical. Should the output be returned as a data frame? Defaults
+#'    to `TRUE`.
+#'
+#' @return A named list with the `mean`, `onset`, `midpoint`, and `offset` of the
+#'    calculated brightest or darkest period, or if `as.df == TRUE` a data frame 
+#'    with columns named `{period}_{timespan}_{metric}`. The output type corresponds
+#'    to the type of `Time.vector`, e.g., if `Time.vector` is HMS, the timing metrics 
+#'    will be also HMS, and vice versa for POSIXct and numeric. 
+#'
+#' @details Assumes regular 24h light data. Otherwise, results may not be
+#'    meaningful. Looping the data is recommended for finding the darkest period.
+#'
+#' @references 
+#'   Hartmeyer, S.L., Andersen, M. (2023). Towards a framework for light-dosimetry studies:
+#'   Quantification metrics. \emph{Lighting Research & Technology}. 
+#'   \url{https://doi.org/10.1177/14771535231170500}
+#'
+#' @export
+#' 
+#' @family metrics
+#'
+#' @examples
+# Dataset with light > 250lx between 06:00 and 18:00
+#' dataset1 <-
+#'   tibble::tibble(
+#'     Id = rep("A", 24),
+#'     Datetime = lubridate::as_datetime(0) + lubridate::hours(0:23),
+#'     MEDI = c(rep(1, 6), rep(250, 13), rep(1, 5))
+#'   )
+#' 
+#' dataset1 %>%
+#'   dplyr::reframe(bright_dark_period(MEDI, Datetime, "brightest", "10 hours",
+#'     as.df = TRUE))
+#' dataset1 %>%
+#'   dplyr::reframe(bright_dark_period(MEDI, Datetime, "darkest", "5 hours",
+#'     loop = TRUE, as.df = TRUE))
+
+bright_dark_period <- function(Light.vector,
+                               Time.vector,
+                               period = c("brightest", "darkest"),
+                               timespan = "10 hours",
+                               epoch = "dominant.epoch",
+                               loop = FALSE,
+                               na.rm = FALSE,
+                               as.df = FALSE) {
+  # Match arguments
+  period <- match.arg(period)
+
+  # Perform argument checks
+  stopifnot(
+    "`Light.vector` must be numeric!" = is.numeric(Light.vector),
+    "`Time.vector` must be POSIXct or HMS" = lubridate::is.POSIXct(Time.vector) | 
+      hms::is_hms(Time.vector),
+    "`epoch` must either be a duration or a string" =
+      lubridate::is.duration(epoch) | is.character(epoch),
+    "`timespan` must either be a duration or a string" =
+      lubridate::is.duration(timespan) | is.character(timespan),
+    "`na.rm` must be logical!" = is.logical(na.rm),
+    "`as.df` must be logical!" = is.logical(as.df)
+  )
+
+  # Check whether time series is regularly spaced
+  if (length(unique(diff(Time.vector))) > 1) {
+    warning("`Time.vector` is not regularly spaced. Calculated results may be incorrect!")
+  }
+
+  # Get the epochs based on the data
+  if (epoch == "dominant.epoch") {
+    epoch <- count_difftime(tibble::tibble(Datetime = Time.vector))$difftime[1]
+  }
+  # If the user specified an epoch, use that instead
+  epoch <- lubridate::as.duration(epoch)
+
+  # Convert timespan to seconds
+  timespan <- lubridate::as.duration(timespan)
+
+  # Check if timespan longer than Time.vector
+  time.total <- dplyr::last(Time.vector) - dplyr::first(Time.vector)
+  stopifnot("Timespan must be shorter than length of `Time.vector` interval!" = 
+              timespan < time.total)
+
+  # Loop data
+  if (loop) {
+    Light.vector <- c(Light.vector, Light.vector)
+    span <- dplyr::last(Time.vector) - Time.vector[1]
+    Time.vector <- c(Time.vector, Time.vector + span + epoch)
+  }
+
+  # Calculate window size
+  window <- floor(as.numeric(timespan) / as.numeric(epoch))
+  if (window %% 2 != 0) {
+    window <- window + 1
+  }
+
+  # Calculate rolling means
+  means <- slider::slide_vec(Light.vector, .f = mean, na.rm = na.rm,
+                         .before = window/2-1, .after = window/2, 
+                         .complete = TRUE)
+
+  # Find maximum/minimum mean value
+  center <- switch(period,
+    "brightest" = which(means == max(means, na.rm = TRUE))[1],
+    "darkest" = which(means == min(means, na.rm = TRUE))[1]
+  )
+
+  # Prepare output
+  out <- list(
+    "mean" = means[center],
+    "midpoint" = Time.vector[center],
+    "onset" = Time.vector[center - (window / 2 - 1)],
+    "offset" = Time.vector[center + (window / 2)]
+  )
+
+  # Return as data frame or numeric matrix
+  if (as.df) {
+    ts <- paste0(as.numeric(timespan, unit = "hours"), "h")
+    out <- tibble::as_tibble(out) %>%
+      dplyr::rename_with(~paste(period, ts, .x, sep = "_"))
+  }
+  return(out)
+}

R/metric_duration_above_threshold.R

@@ -0,0 +1,113 @@
+#' #' Time above/below threshold or within threshold range
+#'
+#' This function calculates the time spent above/below a specified threshold
+#' light level or within a specified range of light levels.
+#'
+#' @param Light.vector Numeric vector containing the light data.
+#' @param Time.vector Vector containing the time data. Can be numeric, HMS or POSIXct.
+#' @param comparison String specifying whether the time above or below threshold
+#'    should be calculated. Can be either `"above"` (the default) or `"below"`. If
+#'    two values are provided for `threshold`, this argument will be ignored.
+#' @param threshold Single numeric value or two numeric values specifying the
+#'    threshold light level(s) to compare with. If a vector with two values is provided,
+#'    the time within the two thresholds will be calculated.
+#' @param epoch The epoch at which the data was sampled. Can be either a
+#'  `lubridate::duration()` or a string. If it is a string, it needs to be
+#'  either `"dominant.epoch"` (the default) for a guess based on the data, or a valid
+#'  `lubridate::duration()` string, e.g., `"1 day"` or `"10 sec"`.
+#' @param na.rm Logical. Should missing values (NA) be removed for the calculation?
+#'    Defaults to `FALSE`.
+#' @param as.df Logical. Should a data frame with be returned? If `TRUE`, a data
+#'    frame with a single column named `TAT_{threshold}` will be returned.
+#'    Defaults to `FALSE`.
+#'
+#' @return A duration object (see \code{\link[lubridate]{duration}}) as single value,
+#'    or single column data frame.
+#'    
+#' @references 
+#'   Hartmeyer, S.L., Andersen, M. (2023). Towards a framework for light-dosimetry studies:
+#'   Quantification metrics. \emph{Lighting Research & Technology}. 
+#'   \url{https://doi.org/10.1177/14771535231170500}
+#'
+#' @export
+#' 
+#' @family metrics
+#'
+#' @examples
+#' N <- 50
+#' # Dataset with epoch = 1min
+#' dataset1 <-
+#'   tibble::tibble(
+#'     Id = rep("A", N),
+#'     Datetime = lubridate::as_datetime(0) + lubridate::minutes(1:N),
+#'     MEDI = sample(c(sample(1:249, N / 2), sample(250:1000, N / 2))),
+#'   )
+#' # Dataset with epoch = 30s
+#' dataset2 <-
+#'   tibble::tibble(
+#'     Id = rep("B", N),
+#'     Datetime = lubridate::as_datetime(0) + lubridate::seconds(seq(30, N * 30, 30)),
+#'     MEDI = sample(c(sample(1:249, N / 2), sample(250:1000, N / 2))),
+#'   )
+#' dataset.combined <- rbind(dataset1, dataset2)
+#' 
+#' dataset1 %>%
+#'   dplyr::reframe("TAT >250lx" = duration_above_threshold(MEDI, Datetime, threshold = 250))
+#' 
+#' dataset1 %>%
+#'   dplyr::reframe(duration_above_threshold(MEDI, Datetime, threshold = 250, as.df = TRUE))
+#' 
+#' # Group by Id to account for different epochs
+#' dataset.combined %>%
+#'   dplyr::group_by(Id) %>%
+#'   dplyr::reframe("TAT >250lx" = duration_above_threshold(MEDI, Datetime, threshold = 250))
+#' 
+duration_above_threshold <- function(Light.vector,
+                                     Time.vector,
+                                     comparison = c("above", "below"),
+                                     threshold,
+                                     epoch = "dominant.epoch",
+                                     na.rm = FALSE,
+                                     as.df = FALSE) {
+  # Match input arguments
+  comparison <- match.arg(comparison)
+
+  # Perform argument checks
+  stopifnot(
+    "`Light.vector` must be numeric!" = is.numeric(Light.vector),
+    "`Time.vector` must be numeric, HMS, or POSIXct" =
+      is.numeric(Time.vector) | hms::is_hms(Time.vector) | lubridate::is.POSIXct(Time.vector),
+    "`threshold` must be numeric!" = is.numeric(threshold),
+    "`threshold` must be either one or two values!" = length(threshold) %in% c(1, 2),
+    "`epoch` must either be a duration or a string" =
+      lubridate::is.duration(epoch) | is.character(epoch),
+    "`na.rm` must be logical!" = is.logical(na.rm),
+    "`as.df` must be logical!" = is.logical(as.df)
+  )
+
+  # Get the epochs based on the data
+  if (epoch == "dominant.epoch") {
+    epoch <- count_difftime(tibble::tibble(Datetime = Time.vector))$difftime[1]
+    epoch <- ifelse(hms::is_hms(Time.vector) | lubridate::is.POSIXct(Time.vector),
+      lubridate::as.duration(epoch), epoch
+    )
+  }
+  # If the user specified an epoch, use that instead
+  else {
+    epoch <- lubridate::as.duration(epoch)
+  }
+
+  # Calculate TAT
+  tat <- sum(compare_threshold(Light.vector, threshold, comparison, na.rm)) * as.numeric(epoch)
+
+  # As duration object
+  tat <- lubridate::as.duration(tat)
+
+  # Return data frame or numeric value
+  if (as.df) {
+    threshold <- stringr::str_flatten(sort(threshold), collapse = "-")
+    return(tibble::tibble("duration_{comparison}_{threshold}" := tat))
+  } else {
+    return(tat)
+  }
+}