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

@@ -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) {

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)

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

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

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

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

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{

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

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)

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