Merge branch 'main' into dev

JZauner · web-flow · commit b05f38aabe17 · 2025-05-20T23:33:49.000+02:00
diff --git a/LightLogR.Rproj b/LightLogR.Rproj
@@ -1,5 +1,4 @@
 Version: 1.0
-ProjectId: 6b298e31-e359-4202-b3a9-ee65538e1c99
 
 RestoreWorkspace: Default
 SaveWorkspace: Default
diff --git a/R/metric_interdaily_stability.R b/R/metric_interdaily_stability.R
@@ -9,6 +9,9 @@
 #'
 #' @param Light.vector Numeric vector containing the light data.
 #' @param Datetime.vector Vector containing the time data. Must be POSIXct.
+#' @param use.samplevar Logical. Should the sample variance be used (divide by N-1)? 
+#'    By default (`FALSE`), the population variance (divide by N) is used, as described
+#'    in Van Someren et al. (1999).
 #' @param na.rm Logical. Should missing values be removed? Defaults to `FALSE`.
 #' @param as.df Logical. Should the output be returned as a data frame? If `TRUE`, a data
 #'    frame with a single column named `interdaily_stability` will be returned.
@@ -31,8 +34,7 @@
 #'   \doi{10.1177/14771535231170500}
 #'
 #' @examples
-#'
-#' set.seed(1)
+#'set.seed(1)
 #' N <- 24 * 7
 #' # Calculate metric for seven 24 h days with two measurements per hour
 #' dataset1 <-
@@ -45,17 +47,19 @@
 #'   dplyr::summarise(
 #'     "Interdaily stability" = interdaily_stability(MEDI, Datetime)
 #'   )
-#'
 interdaily_stability <- function(Light.vector,
                                  Datetime.vector,
+                                 use.samplevar = FALSE,
                                  na.rm = FALSE,
-                                 as.df = FALSE) {
+                                 as.df = FALSE
+                                 ) {
   # Initial checks
   stopifnot(
     "`Light.vector` must be numeric!" = is.numeric(Light.vector),
     "`Datetime.vector` must be POSIXct!" = lubridate::is.POSIXct(Datetime.vector),
     "`Light.vector` and `Datetime.vector` must be same length!" = 
       length(Light.vector) == length(Datetime.vector),
+    "`use.samplevar` must be logical!" = is.logical(use.samplevar),
     "`na.rm` must be logical!" = is.logical(na.rm),
     "`as.df` must be logical!" = is.logical(as.df)
   )
@@ -114,19 +118,34 @@ interdaily_stability <- function(Light.vector,
   if (na.rm) {
     df <- df %>% tidyr::drop_na(Light)
   }
+  
+  # Subtract 1 if `use.samplevar == TRUE`
+  c = as.numeric(use.samplevar)
 
   # Hourly averages for each day
-  total_hourly <- df %>%
+  hourly_data <- df %>%
     dplyr::group_by(Datetime = lubridate::floor_date(Datetime, unit = "1 hour")) %>%
     dplyr::summarise(Light = mean(Light))
+  
+  # N hourly data
+  n <- length(hourly_data$Light)
+  
+  # Overall variance
+  var_total <- sum((hourly_data$Light-mean(hourly_data$Light))^2) / (n-c)
 
   # Hourly average across all days
-  avg_hourly <- total_hourly %>%
+  avg_hourly <- hourly_data %>%
     dplyr::group_by(Hour = lubridate::hour(Datetime)) %>%
     dplyr::summarise(Light = mean(Light))
+  
+  # N per day
+  p <- length(avg_hourly$Light)
+  
+  # Variance across average day
+  var_avg_day <- sum((avg_hourly$Light-mean(hourly_data$Light))^2) / (p-c)
 
   # Variance across average day / variance across all days
-  is <- stats::var(avg_hourly$Light) / stats::var(total_hourly$Light)
+  is <- var_avg_day / var_total
 
   # Return data frame or numeric vector
   if (as.df) {
diff --git a/R/metric_intradaily_variability.R b/R/metric_intradaily_variability.R
@@ -7,6 +7,9 @@
 #'
 #' @param Light.vector Numeric vector containing the light data.
 #' @param Datetime.vector Vector containing the time data. Must be POSIXct.
+#' @param use.samplevar Logical. Should the sample variance be used (divide by N-1)? 
+#'    By default (`FALSE`), the population variance (divide by N) is used, as described
+#'    in Van Someren et al. (1999).
 #' @param na.rm Logical. Should missing values be removed? Defaults to `FALSE`.
 #' @param as.df Logical. Should the output be returned as a data frame? If `TRUE`, a data
 #'    frame with a single column named `intradaily_variability` will be returned.
@@ -46,6 +49,7 @@
 #'
 intradaily_variability <- function(Light.vector,
                                    Datetime.vector,
+                                   use.samplevar = FALSE,
                                    na.rm = FALSE,
                                    as.df = FALSE) {
   # Initial checks
@@ -54,6 +58,7 @@ intradaily_variability <- function(Light.vector,
     "`Datetime.vector` must be POSIXct!" = lubridate::is.POSIXct(Datetime.vector),
     "`Light.vector` and `Datetime.vector` must be same length!" = 
       length(Light.vector) == length(Datetime.vector),
+    "`use.samplevar` must be logical!" = is.logical(use.samplevar),
     "`na.rm` must be logical!" = is.logical(na.rm),
     "`as.df` must be logical!" = is.logical(as.df)
   )
@@ -112,18 +117,26 @@ intradaily_variability <- function(Light.vector,
   if (na.rm) {
     df <- df %>% tidyr::drop_na(Light)
   }
+  
+  # Subtract 1 if `use.samplevar == TRUE`
+  c = as.numeric(use.samplevar)
 
   # Hourly averages for each day
-  total_hourly <- df %>%
+  hourly_data <- df %>%
     dplyr::group_by(lubridate::floor_date(Datetime, unit = "1 hour")) %>%
     dplyr::summarise(Light = mean(Light))
+  
+  # N hourly data
+  n <- length(hourly_data$Light)
 
-  # Variance of consecutive hourly differences
-  var_hourly_diff <-
-    sum(diff(total_hourly$Light)^2) / (length(total_hourly$Light) - 1)
+  # Overall variance
+  var_total <- sum((hourly_data$Light-mean(hourly_data$Light))^2) / (n-c)
+  
+  # Sum of squares of consecutive hourly differences
+  var_hourly_diff <- sum(diff(hourly_data$Light)^2) / (n-1)
 
-  # Variance of consecutive differences / variance across all days
-  iv <- var_hourly_diff / stats::var(total_hourly$Light)
+  # Sum of squares of consecutive differences / variance across all days
+  iv <- var_hourly_diff / var_total
 
   # Return data frame or numeric vector
   if (as.df) {
diff --git a/man/import_Dataset.Rd b/man/import_Dataset.Rd
diff --git a/man/interdaily_stability.Rd b/man/interdaily_stability.Rd
diff --git a/man/intradaily_variability.Rd b/man/intradaily_variability.Rd
diff --git a/tests/testthat/helper.R b/tests/testthat/helper.R
@@ -0,0 +1 @@
+trunc_decimal = function(x, accuracy){base::trunc(x * 10^accuracy) / 10^accuracy}
diff --git a/tests/testthat/test-metric_duration_above_threshold.R b/tests/testthat/test-metric_duration_above_threshold.R
@@ -12,7 +12,7 @@ test_that("Calculation works", {
     duration_above_threshold(MEDI, datetime, "below", 100), lubridate::dminutes(7)
   )
   expect_equal(
-    duration_above_threshold(MEDI, datetime, threshold = c(0,99)), lubridate::dminutes(4)
+    duration_above_threshold(MEDI, datetime, threshold = c(0,100)), lubridate::dminutes(5)
   )
 })
 
diff --git a/tests/testthat/test-metric_interdaily_stability.R b/tests/testthat/test-metric_interdaily_stability.R
@@ -0,0 +1,43 @@
+test_that("Population variance works", {
+  MEDI = sample.data.environment$MEDI
+  Datetime = sample.data.environment$Datetime
+  expect_equal(trunc_decimal(interdaily_stability(MEDI, Datetime),2), 0.77)
+})
+
+test_that("Sample variance works", {
+  MEDI = sample.data.environment$MEDI
+  Datetime = sample.data.environment$Datetime
+  expect_equal(trunc_decimal(interdaily_stability(MEDI, Datetime, use.samplevar = TRUE),2), 0.79)
+})
+
+test_that("IS for Gaussian noise near 2", {
+  N = 24*1000
+  MEDI = rnorm(N)
+  Datetime = lubridate::as_datetime(lubridate::dhours(0:(N-1)), tz = "UTC")
+  expect_equal(round(interdaily_stability(MEDI, Datetime),1), 0)
+})
+
+test_that("IS equals 1 for perfect IS", {
+  MEDI = rep(c(rep(1, 6), rep(250, 4), rep(500, 5), rep(750, 4), rep(1, 5)),7)
+  Datetime = lubridate::as_datetime(lubridate::dhours(0:(24*7-1)), tz = "UTC")
+  expect_equal(interdaily_stability(MEDI, Datetime), 1)
+})
+
+test_that("Handling missing values works", {
+  MEDI = rep(c(as.numeric(NA), rep(1, 5), rep(250, 4), rep(500, 5), rep(750, 4), rep(1, 5)),7)
+  Datetime = lubridate::as_datetime(lubridate::dhours(0:(24*7-1)), tz = "UTC")
+  expect_equal(suppressWarnings({interdaily_stability(MEDI, Datetime)}), as.numeric(NA))
+  expect_equal(suppressWarnings({interdaily_stability(MEDI, Datetime, na.rm=TRUE)}), 1)
+})
+
+test_that("Return data frame works", {
+  MEDI = sample.data.environment$MEDI
+  Datetime = sample.data.environment$Datetime
+  expect_equal(
+    round(interdaily_stability(MEDI, Datetime, as.df = TRUE), 1), 
+    tibble::tibble("interdaily_stability" = 0.8)
+  )
+})
+
+
+
diff --git a/tests/testthat/test-metric_intradaily_variability.R b/tests/testthat/test-metric_intradaily_variability.R
@@ -0,0 +1,40 @@
+test_that("Population variance works", {
+  MEDI = sample.data.environment$MEDI
+  Datetime = sample.data.environment$Datetime
+  expect_equal(trunc_decimal(intradaily_variability(MEDI, Datetime),3), 0.166)
+})
+
+test_that("Sample variance works", {
+  MEDI = sample.data.environment$MEDI
+  Datetime = sample.data.environment$Datetime
+  expect_equal(trunc_decimal(intradaily_variability(MEDI, Datetime, use.samplevar = TRUE),3), 0.165)
+})
+
+test_that("IV for Gaussian noise near 2", {
+  N = 24*1000
+  MEDI = rnorm(N)
+  Datetime = lubridate::as_datetime(lubridate::dhours(0:(N-1)), tz = "UTC")
+  expect_equal(round(intradaily_variability(MEDI, Datetime)), 2)
+})
+
+test_that("IV for Sine wave near 0", {
+  MEDI = sin(seq(0,7*pi,length.out = 24*7))
+  Datetime = lubridate::as_datetime(lubridate::dhours(0:(24*7-1)), tz = "UTC")
+  expect_equal(round(intradaily_variability(MEDI, Datetime),1), 0)
+})
+
+test_that("Handling missing values works", {
+  MEDI = rep(c(as.numeric(NA), rep(1, 5), rep(250, 4), rep(500, 5), rep(750, 4), rep(1, 5)),7)
+  Datetime = lubridate::as_datetime(lubridate::dhours(0:(24*7-1)), tz = "UTC")
+  expect_equal(suppressWarnings({intradaily_variability(MEDI, Datetime)}),as.numeric(NA))
+  expect_equal(suppressWarnings({round(intradaily_variability(MEDI, Datetime, na.rm=TRUE),1)}), 0.4)
+})
+
+test_that("Return data frame works", {
+  MEDI = sample.data.environment$MEDI
+  Datetime = sample.data.environment$Datetime
+  expect_equal(
+    round(intradaily_variability(MEDI, Datetime, as.df = TRUE), 2), 
+    tibble::tibble("intradaily_variability" = 0.17)
+  )
+})
diff --git a/tests/testthat/test-metric_period_above_threshold.R b/tests/testthat/test-metric_period_above_threshold.R
@@ -9,7 +9,7 @@ test_that("Calculation works", {
     period_above_threshold(MEDI, datetime, "below", 10), lubridate::dminutes(5)
   )
   expect_equal(
-    period_above_threshold(MEDI, datetime, threshold = c(12,13)), lubridate::dminutes(1)
+    period_above_threshold(MEDI, datetime, threshold = c(11,13)), lubridate::dminutes(2)
   )
   expect_equal(
     period_above_threshold(MEDI, datetime, threshold = 10, loop = TRUE), lubridate::dminutes(7)
diff --git a/tests/testthat/test-metric_timing_above_threshold.R b/tests/testthat/test-metric_timing_above_threshold.R
@@ -1,4 +1,4 @@
-test_that("Calculation works", {
+test_that("Above threshold works", {
   MEDI = c(rep(1, 6), rep(250, 4), rep(500, 5), rep(750, 4), rep(1, 5))
   Datetime = lubridate::as_datetime(lubridate::dhours(0:23), tz = "UTC")
   expect_equal(
@@ -9,6 +9,11 @@ test_that("Calculation works", {
       "last" = lubridate::as_datetime(lubridate::dhours(18), tz = "UTC")
     )
   )
+})
+
+test_that("Below threshold works", {
+  MEDI = c(rep(1, 6), rep(250, 4), rep(500, 5), rep(750, 4), rep(1, 5))
+  Datetime = lubridate::as_datetime(lubridate::dhours(0:23), tz = "UTC")
   expect_equal(
     timing_above_threshold(MEDI, Datetime, comparison = "below", threshold = 250), 
     list(
@@ -17,6 +22,11 @@ test_that("Calculation works", {
       "last" = lubridate::as_datetime(lubridate::dhours(23), tz = "UTC")
     )
   )
+})
+
+test_that("Between thresholds works", {
+  MEDI = c(rep(1, 6), rep(250, 4), rep(400, 5), rep(750, 4), rep(1, 5))
+  Datetime = lubridate::as_datetime(lubridate::dhours(0:23), tz = "UTC")
   expect_equal(
     timing_above_threshold(MEDI, Datetime, threshold = c(250, 500)), 
     list(

Original file line number	Diff line number	Diff line change
`@@ -1,5 +1,4 @@`
`1`	`1`	`Version: 1.0`
`2`		`-ProjectId: 6b298e31-e359-4202-b3a9-ee65538e1c99`
`3`	`2`
`4`	`3`	`RestoreWorkspace: Default`
`5`	`4`	`SaveWorkspace: Default`
Original file line number	Diff line number	Diff line change
`@@ -0,0 +1 @@`
	`1`	`+trunc_decimal = function(x, accuracy){base::trunc(x * 10^accuracy) / 10^accuracy}`
Original file line number	Diff line number	Diff line change
`@@ -12,7 +12,7 @@ test_that("Calculation works", {`
`12`	`12`	`duration_above_threshold(MEDI, datetime, "below", 100), lubridate::dminutes(7)`
`13`	`13`	`)`
`14`	`14`	`expect_equal(`
`15`		`- duration_above_threshold(MEDI, datetime, threshold = c(0,99)), lubridate::dminutes(4)`
	`15`	`+ duration_above_threshold(MEDI, datetime, threshold = c(0,100)), lubridate::dminutes(5)`
`16`	`16`	`)`
`17`	`17`	`})`
`18`	`18`
Original file line number	Diff line number	Diff line change
`@@ -9,7 +9,7 @@ test_that("Calculation works", {`
`9`	`9`	`period_above_threshold(MEDI, datetime, "below", 10), lubridate::dminutes(5)`
`10`	`10`	`)`
`11`	`11`	`expect_equal(`
`12`		`- period_above_threshold(MEDI, datetime, threshold = c(12,13)), lubridate::dminutes(1)`
	`12`	`+ period_above_threshold(MEDI, datetime, threshold = c(11,13)), lubridate::dminutes(2)`
`13`	`13`	`)`
`14`	`14`	`expect_equal(`
`15`	`15`	`period_above_threshold(MEDI, datetime, threshold = 10, loop = TRUE), lubridate::dminutes(7)`