From 2763b9f947a1ae086f90dfd40228339e02f893e3 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Mon, 2 Sep 2024 15:33:24 +0200 Subject: [PATCH] improvements to nap detection code #967 - added parameter passible_nap_gap - existing parameters now actually used - midnight better handled in sibreport --- R/check_params.R | 15 +++++++- R/g.part5.analyseRest.R | 47 ++++++++++++++++++++++--- R/g.part5_analyseSegment.R | 2 +- R/g.sibreport.R | 6 ++++ R/load_params.R | 3 +- man/GGIR.Rd | 5 +++ man/g.part5.analyseRest.Rd | 11 ++---- tests/testthat/test_chainof5parts.R | 6 +++- tests/testthat/test_load_check_params.R | 2 +- tests/testthat/test_part5_analyseRest.R | 23 ++++++------ 10 files changed, 90 insertions(+), 30 deletions(-) diff --git a/R/check_params.R b/R/check_params.R index f5d8c17be..2f31c9430 100644 --- a/R/check_params.R +++ b/R/check_params.R @@ -32,7 +32,8 @@ check_params = function(params_sleep = c(), params_metrics = c(), numeric_params = c("anglethreshold", "timethreshold", "longitudinal_axis", "possible_nap_window", "possible_nap_dur", "colid", "coln1", "def.noc.sleep", "nnights", - "sleepefficiency.metric", "possible_nap_edge_acc", "HDCZA_threshold") + "sleepefficiency.metric", "possible_nap_edge_acc", "HDCZA_threshold", + "possible_nap_gap") boolean_params = c("ignorenonwear", "HASPT.ignore.invalid", "relyonguider", "sleeplogidnum") character_params = c("HASPT.algo", "HASIB.algo", "Sadeh_axis", "nap_model", @@ -187,6 +188,18 @@ check_params = function(params_sleep = c(), params_metrics = c(), } else if (length(params_sleep[["def.noc.sleep"]]) == 2) { params_sleep[["HASPT.algo"]] = "notused" } + if (length(params_sleep[["possible_nap_gap"]]) != 1) { + stop(paste0("Parameter possible_nap_gap has length ", length(params_sleep[["possible_nap_gap"]]), + " while length 1 is expected"), call. = FALSE) + } + if (length(params_sleep[["possible_nap_window"]]) != 2) { + stop(paste0("Parameter possible_nap_window has length ", length(params_sleep[["possible_nap_window"]]), + " while length 2 is expected"), call. = FALSE) + } + if (length(params_sleep[["possible_nap_dur"]]) != 2) { + stop(paste0("Parameter possible_nap_dur has length ", length(params_sleep[["possible_nap_dur"]]), + " while length 2 is expected"), call. = FALSE) + } } if (length(params_metrics) > 0 & length(params_sleep) > 0) { diff --git a/R/g.part5.analyseRest.R b/R/g.part5.analyseRest.R index 33328d2ea..e3b088009 100644 --- a/R/g.part5.analyseRest.R +++ b/R/g.part5.analyseRest.R @@ -1,7 +1,7 @@ g.part5.analyseRest = function(sibreport = NULL, dsummary = NULL, ds_names = NULL, fi = NULL, di = NULL, - ts = NULL, tz = NULL, possible_nap_dur = 0, - possible_nap_edge_acc = Inf) { + ts = NULL, tz = NULL, + params_sleep = NULL) { # define function to summarise overlap between selfreported behaviour and sibs summarise_overlap = function(srep_tmp, X, Y, xi, yi, name = "", sumobject = NULL) { # X: column name in srep_temp to reflect overlap SIB with Selfreport @@ -55,16 +55,53 @@ g.part5.analyseRest = function(sibreport = NULL, dsummary = NULL, sibreport$end = as.POSIXct(sibreport$end, tz = tz) sibreport$start = as.POSIXct(sibreport$start, tz = tz) + #--------------------------------------- + # merge sibs when gap is shorter than possible_nap_gap + if (params_sleep[["possible_nap_gap"]] > 0) { + sibreport$gap2next = NA + Nrow = nrow(sibreport) + sibreport$gap2next[1:(Nrow - 1)] = as.numeric(sibreport$start[2:Nrow]) - as.numeric(sibreport$end[1:(Nrow - 1)]) + sibreport$gap2next[which(sibreport$type != "sib" | sibreport$gap2next < 0)] = NA + iter = 1 + while (iter < nrow(sibreport)) { + if (!is.na(sibreport$gap2next[iter]) && + sibreport$gap2next[iter] < params_sleep[["possible_nap_gap"]]) { + sibreport$end[iter] = sibreport$end[iter + 1] + sibreport$mean_acc_1min_after[iter] = sibreport$mean_acc_1min_after[iter + 1] + sibreport = sibreport[-(iter + 1),] + sibreport$gap2next[iter] = as.numeric(sibreport$start[iter + 1]) - as.numeric(sibreport$end[iter]) + # no need to increment iter, because by merging the sib blocks + # the current iter now refers to the next gap + } else { + iter = iter + 1 + } + if (iter > nrow(sibreport) - 1) { + break() + } + } + sibreport$duration = as.numeric(difftime(sibreport$end, sibreport$start, units = "mins")) + } + # Only consider sib episodes with minimum duration if (length(grep(pattern = "mean_acc_1min", x = colnames(sibreport))) > 0) { sibreport$acc_edge = pmax(sibreport$mean_acc_1min_before, sibreport$mean_acc_1min_after) } else { sibreport$acc_edge = 0 } + sibreport$startHour = as.numeric(format(sibreport$start, "%H")) + sibreport$endHour = as.numeric(format(sibreport$end, "%H")) + + overlapMidnight = which(sibreport$endHour < sibreport$startHour) + if (length(overlapMidnight) > 0) { + sibreport$endHour[overlapMidnight] = sibreport$endHour[overlapMidnight] + 24 + } + longboutsi = which((sibreport$type == "sib" & - sibreport$duration >= possible_nap_dur[1] & - sibreport$duration < possible_nap_dur[2] & - sibreport$acc_edge <= possible_nap_edge_acc) | + sibreport$duration >= params_sleep[["possible_nap_dur"]][1] & + sibreport$duration < params_sleep[["possible_nap_dur"]][2] & + sibreport$acc_edge <= params_sleep[["possible_nap_edge_acc"]] & + sibreport$startHour >= params_sleep[["possible_nap_window"]][1] & + sibreport$endHour < params_sleep[["possible_nap_window"]][2]) | (sibreport$type != "sib" & sibreport$duration >= 1)) # for qc purposes: dsummary[di,fi] = length(longboutsi) diff --git a/R/g.part5_analyseSegment.R b/R/g.part5_analyseSegment.R index 3760b1c1c..058f4ad88 100644 --- a/R/g.part5_analyseSegment.R +++ b/R/g.part5_analyseSegment.R @@ -505,7 +505,7 @@ g.part5_analyseSegment = function(indexlog, timeList, levelList, ds_names = ds_names, fi = fi, di = si, ts = ts[sse[ts$diur[sse] == 0], ], tz = params_general[["desiredtz"]], - possible_nap_dur = params_sleep[["possible_nap_dur"]]) + params_sleep = params_sleep) fi = restAnalyses$fi si = restAnalyses$di ts[sse[ts$diur[sse] == 0], ] = restAnalyses$ts diff --git a/R/g.sibreport.R b/R/g.sibreport.R index f6b15e9d4..4bca7a6cf 100644 --- a/R/g.sibreport.R +++ b/R/g.sibreport.R @@ -101,6 +101,12 @@ g.sibreport = function(ts, ID, epochlength, logs_diaries=c(), desiredtz="") { tt2 = as.POSIXlt(timestamps[(bi * 2)], tz = desiredtz) logreport_tmp$start[bi] = format(tt1) logreport_tmp$end[bi] = format(tt2) + if (length(unlist(strsplit(logreport_tmp$start[bi], " "))) == 1) { + logreport_tmp$start[bi] = paste0(logreport_tmp$start[bi], " 00:00:00") + } + if (length(unlist(strsplit(logreport_tmp$end[bi], " "))) == 1) { + logreport_tmp$end[bi] = paste0(logreport_tmp$end[bi], " 00:00:00") + } logreport_tmp$duration[bi] = abs(as.numeric(difftime(time1 = tt1, time2 = tt2, units = "mins"))) } } diff --git a/R/load_params.R b/R/load_params.R index 3366f01d5..2ee0b01b4 100644 --- a/R/load_params.R +++ b/R/load_params.R @@ -24,8 +24,9 @@ load_params = function(topic = c("sleep", "metrics", "rawdata", sleeplogsep = NULL, sleepwindowType = "SPT", possible_nap_window = c(9, 18), possible_nap_dur = c(15, 240), - nap_model = c(), sleepefficiency.metric = 1, + possible_nap_gap = 0, possible_nap_edge_acc = Inf, + nap_model = c(), sleepefficiency.metric = 1, HDCZA_threshold = c()) } if ("metrics" %in% topic) { diff --git a/man/GGIR.Rd b/man/GGIR.Rd index c49399018..bca8593dd 100755 --- a/man/GGIR.Rd +++ b/man/GGIR.Rd @@ -1253,6 +1253,11 @@ GGIR(mode = 1:5, Maximum acceleration before or after the SIB for the nap to be considered. By default this will allow all possible naps. } + \item{possible_nap_gap}{ + Numeric (default = 0). + Time gap expressed in seconds that is allowed between the sustained + inactivity bouts that form the naps. + } } } diff --git a/man/g.part5.analyseRest.Rd b/man/g.part5.analyseRest.Rd index 4347c0ab4..a05f7bb77 100644 --- a/man/g.part5.analyseRest.Rd +++ b/man/g.part5.analyseRest.Rd @@ -11,8 +11,7 @@ g.part5.analyseRest(sibreport = NULL, dsummary = NULL, ds_names = NULL, fi = NULL, di = NULL, ts = NULL, tz = NULL, - possible_nap_dur = 0, - possible_nap_edge_acc = Inf) + params_sleep = NULL) } \arguments{ \item{sibreport}{ @@ -37,12 +36,8 @@ \item{tz}{ Timezone database name } - \item{possible_nap_dur}{ - Minimum sib duration to be considered. For self-reported naps/nonwear - there is a minimum duration of 1 epoch. - } - \item{possible_nap_edge_acc}{ - Maximum acceleration before or after the SIB for it to be considered. + \item{params_sleep}{ + See \link{GGIR} } } \value{ diff --git a/tests/testthat/test_chainof5parts.R b/tests/testthat/test_chainof5parts.R index 1a33b1069..25856c758 100644 --- a/tests/testthat/test_chainof5parts.R +++ b/tests/testthat/test_chainof5parts.R @@ -204,7 +204,11 @@ test_that("chainof5parts", { overwrite = TRUE, excludefirstlast = FALSE, do.parallel = do.parallel, frag.metrics = "all", save_ms5rawlevels = TRUE, part5_agg2_60seconds = TRUE, do.sibreport = TRUE, nap_model = "hip3yr", - iglevels = 1, timewindow = c("MM", "WW", "OO")) + iglevels = 1, timewindow = c("MM", "WW", "OO"), + possible_nap_window = c(0, 24), + possible_nap_dur = c(0, 240), + possible_nap_edge_acc = Inf, + possible_nap_gap = 0) sibreport_dirname = "output_test/meta/ms5.outraw/sib.reports" expect_true(dir.exists(sibreport_dirname)) expect_true(file.exists(paste0(sibreport_dirname, "/sib_report_123A_testaccfile_T5A5.csv"))) diff --git a/tests/testthat/test_load_check_params.R b/tests/testthat/test_load_check_params.R index d7f41d5b3..d56f011cc 100644 --- a/tests/testthat/test_load_check_params.R +++ b/tests/testthat/test_load_check_params.R @@ -11,7 +11,7 @@ test_that("load_params can load parameters", { # Test length of objects expect_equal(length(params), 8) - expect_equal(length(params$params_sleep), 22) + expect_equal(length(params$params_sleep), 23) expect_equal(length(params$params_metrics), 41) expect_equal(length(params$params_rawdata), 39) expect_equal(length(params$params_247), 23) diff --git a/tests/testthat/test_part5_analyseRest.R b/tests/testthat/test_part5_analyseRest.R index 8ea2be3e3..8ec08f988 100644 --- a/tests/testthat/test_part5_analyseRest.R +++ b/tests/testthat/test_part5_analyseRest.R @@ -16,12 +16,12 @@ test_that("Overlap 1 nap and 1 sib", { end = c("2022-06-02 14:20:00", "2022-06-02 14:20:00")) sibreport$duration = as.numeric(difftime(time1 = sibreport$end, time2 = sibreport$start, units = "mins", tz = tz)) + params_sleep = load_params()$params_sleep + params_sleep[["possible_nap_dur"]] = c(0, 240) restAnalyses = g.part5.analyseRest(sibreport = sibreport, dsummary = dsummary, ds_names = ds_names, fi = fi, di = di, - ts = ts, - possible_nap_dur = c(0, 240), - possible_nap_edge_acc = Inf, - tz = tz) + ts = ts, tz = tz, + params_sleep = params_sleep) fi = restAnalyses$fi di = restAnalyses$di dsummary = restAnalyses$dsummary @@ -47,12 +47,11 @@ test_that("Overlap 1 nonwear and 1 sib", { start = c("2022-06-02 14:00:00", "2022-06-02 14:05:00"), end = c("2022-06-02 14:20:00", "2022-06-02 14:20:00")) sibreport$duration = as.numeric(difftime(time1 = sibreport$end, time2 = sibreport$start, units = "mins", tz = tz)) + params_sleep = load_params()$params_sleep + params_sleep[["possible_nap_dur"]] = c(0, 240) restAnalyses = g.part5.analyseRest(sibreport = sibreport, dsummary = dsummary, ds_names = ds_names, fi = fi, di = di, - ts = ts, - possible_nap_dur = c(0, 240), - possible_nap_edge_acc = Inf, - tz = tz) + ts = ts, tz = tz,params_sleep = params_sleep) fi = restAnalyses$fi di = restAnalyses$di dsummary = restAnalyses$dsummary @@ -79,12 +78,12 @@ test_that("No overlap 1 nonwear, 1 nap, and 1 sib", { start = c("2022-06-02 12:00:00", "2022-06-02 13:00:00", "2022-06-02 15:00:00"), end = c("2022-06-02 12:20:00", "2022-06-02 13:20:00", "2022-06-02 15:20:00")) sibreport$duration = as.numeric(difftime(time1 = sibreport$end, time2 = sibreport$start, units = "mins", tz = tz)) + params_sleep = load_params()$params_sleep + params_sleep[["possible_nap_dur"]] = c(0, 240) + restAnalyses = g.part5.analyseRest(sibreport = sibreport, dsummary = dsummary, ds_names = ds_names, fi = fi, di = di, - ts = ts, - possible_nap_dur = c(0, 240), - possible_nap_edge_acc = Inf, - tz = tz) + ts = ts, tz = tz, params_sleep = params_sleep) fi = restAnalyses$fi di = restAnalyses$di dsummary = restAnalyses$dsummary