add option to require complete lastnight in part 5

NEWS.md

@@ -1,3 +1,7 @@
+# CHANGES IN GGIR VERSION 3.1-?
+
+- Part 5: Add parameters require_complete_lastnight_part5 to control whether last window is included if last night is incomplete. #1196
+
 # CHANGES IN GGIR VERSION 3.1-4
 
 - Part 3: Update threshold used for HorAngle to 60 degree, and auto-setting HASPT.ignore.invalid to NA when NotWorn guider is used. #1186

R/check_params.R

@@ -117,7 +117,7 @@ check_params = function(params_sleep = c(), params_metrics = c(),
     boolean_params = c("epochvalues2csv", "save_ms5rawlevels", "save_ms5raw_without_invalid",
                        "storefolderstructure", "dofirstpage", "visualreport", "week_weekend_aggregate.part5",
                        "do.part3.pdf", "outliers.only", "do.visual", "do.sibreport", "visualreport_without_invalid",
-                       "do.part2.pdf")
+                       "do.part2.pdf", "require_complete_lastnight_part5")
     character_params = c("save_ms5raw_format", "timewindow", "sep_reports", "sep_config",
                          "dec_reports", "dec_config")
     check_class("output", params = params_output, parnames = numeric_params, parclass = "numeric")

R/g.part5.R

@@ -600,7 +600,8 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(),
                                            includedaycrit.part5 = params_cleaning[["includedaycrit.part5"]],
                                            ID = ID,
                                            params_output = params_output,
-                                           params_247 = params_247)
+                                           params_247 = params_247,
+                                           timewindow = timewindowi)
                   }
                 }
               }
@@ -662,7 +663,9 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(),
                 if (length(GGIRversion) != 1) GGIRversion = sessionInfo()$otherPkgs$GGIR$Version
               }
               output$GGIRversion = GGIRversion
-              save(output, tail_expansion_log, GGIRversion,
+              # Capture final timestamp to ease filtering last window in g.report.part5
+              last_timestamp = time_POSIX[length(time_POSIX)] 
+              save(output, tail_expansion_log, GGIRversion, last_timestamp,
                    file = paste(metadatadir, ms5.out, "/", fnames.ms3[i], sep = ""))
             }
           }

R/g.part5.savetimeseries.R

@@ -2,7 +2,8 @@ g.part5.savetimeseries = function(ts, LEVELS, desiredtz, rawlevels_fname,
                                   DaCleanFile = NULL,
                                   includedaycrit.part5 = 2/3, ID = NULL,
                                   params_output,
-                                  params_247 = NULL) {
+                                  params_247 = NULL,
+                                  timewindow = NULL) {
 
   ms5rawlevels = data.frame(date_time = ts$time, class_id = LEVELS,
                             # class_name = rep("",Nts),
@@ -23,6 +24,17 @@ g.part5.savetimeseries = function(ts, LEVELS, desiredtz, rawlevels_fname,
   names(mdat)[which(names(mdat) == "nonwear")] = "invalidepoch"
   names(mdat)[which(names(mdat) == "diur")] = "SleepPeriodTime"
   mdat = mdat[,-which(names(mdat) == "date_time")]
+  if ("require_complete_lastnight_part5" %in% names(params_output) &&
+      params_output[["require_complete_lastnight_part5"]] == TRUE) {
+    last_timestamp = as.numeric(format(mdat$timestamp[length(mdat$timestamp)], "%H"))
+    if ((timewindow == "MM" || timewindow == "OO") && last_timestamp < 9) {
+      mdat$window[which(mdat$window == max(mdat$window))] = 0 
+    }
+    if (timewindow == "WW" && last_timestamp < 15) {
+      mdat$window[which(mdat$window == max(mdat$window))] = 0 
+    }
+  }
+
   # Add invalid day indicator
   mdat$invalid_wakinghours = mdat$invalid_sleepperiod =  mdat$invalid_fullwindow = 100
   wakeup = which(diff(c(mdat$SleepPeriodTime,0)) == -1) + 1 # first epoch of each day

R/g.report.part5.R

@@ -47,21 +47,35 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c
 
       x$wear_min_day = (1 - (x$nonwear_perc_day / 100)) * x$dur_day_min #valid minute during waking hours
       x$wear_perc_day = 100 - x$nonwear_perc_day #wear percentage during waking hours
+      x$lasttimestamp = as.numeric(x$lasttimestamp)
 
       minimumValidMinutesMM = 0 # default
       if (length(params_cleaning[["includedaycrit"]]) == 2) {
         minimumValidMinutesMM = params_cleaning[["includedaycrit"]][2] * 60
       }
+      if (params_output[["require_complete_lastnight_part5"]] == FALSE) {
+        x$lastnight = FALSE
+      } else {
+        x$lastnight = x$window_number == max(x$window_number)
+      }
       if (window == "WW" | window == "OO") {
         indices = which(x$wear_perc_day >= includeday_wearPercentage &
                           x$wear_min_day >= includeday_absolute &
                           x$dur_spt_min > 0 & x$dur_day_min > 0 &
+                          ((x$lastnight == TRUE & x$lasttimestamp >= 15 & window == "WW") |
+                             (x$lastnight == TRUE & x$lasttimestamp >= 9 & window == "OO") |
+                             x$lastnight == FALSE) &
                           include_window == TRUE &
                           x$wear_min_day_spt >= minimumValidMinutesMM)
       } else if (window == "MM") {
+        # Note: For MM it would also be good to put requirements on the 
+        # timing of the last timestamp to include/exlcude the last window
+        # but not done yet
         indices = which(x$wear_perc_day >= includeday_wearPercentage &
                           x$wear_min_day >= includeday_absolute &
                           x$dur_spt_min > 0 & x$dur_day_min > 0 &
+                          ((x$lastnight == TRUE & x$lasttimestamp > 9) |
+                             x$lastnight == FALSE) &
                           x$dur_day_spt_min >= (params_cleaning[["minimum_MM_length.part5"]] * 60) &
                           include_window == TRUE &
                           x$wear_min_day_spt >= minimumValidMinutesMM)
@@ -118,12 +132,17 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c
                  " take a few minutes\n"))
     }
     myfun = function(x, expectedCols = c()) {
-      tail_expansion_log = NULL
+      tail_expansion_log = last_timestamp = output = NULL
       load(file = x)
       cut = which(output[, 1] == "")
       if (length(cut) > 0 & length(cut) < nrow(output)) {
         output = output[-cut, which(colnames(output) != "")]
       }
+      if (exists("last_timestamp") == TRUE) {
+        output$lasttimestamp = as.numeric(format(last_timestamp, "%H"))
+      } else {
+        output$lasttimestamp = Inf # use dummy value
+      }
       out = as.matrix(output)
       if (length(expectedCols) > 0) {
         tmp = as.data.frame(matrix(0, 0, length(expectedCols)))
@@ -269,8 +288,11 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c
                   sep = params_output[["sep_reports"]],
                   dec = params_output[["dec_reports"]])
                 # store all summaries in csv files with cleaning criteria
-                validdaysi = getValidDayIndices(x = OF3, window = uwi[j],
+                validdaysi = getValidDayIndices(x = OF3_clean, window = uwi[j],
                                                 params_cleaning = params_cleaning)
+                if ("lasttimestamp" %in% colnames(OF3_clean)) {
+                  OF3_clean = OF3_clean[, -which(colnames(OF3_clean) == "lasttimestamp")]
+                }
                 if (length(validdaysi) > 0) {
                   data.table::fwrite(
                     OF3_clean[validdaysi, ],

R/load_params.R

@@ -115,7 +115,8 @@ load_params = function(topic = c("sleep", "metrics", "rawdata",
                          do.sibreport = FALSE, do.part2.pdf = TRUE,
                          sep_reports = ",", sep_config = ",", 
                          dec_reports = ".", dec_config = ".", 
-                         visualreport_without_invalid = TRUE)
+                         visualreport_without_invalid = TRUE,
+                         require_complete_lastnight_part5 = FALSE)
 
   }
   if ("general" %in% topic) {

man/GGIR.Rd

@@ -1539,6 +1539,13 @@ GGIR(mode = 1:5,
           If TRUE, then reports generated with \code{visualreport = TRUE} only show
           the windows with sufficiently valid data according to \code{includedaycrit}
           when viewingwindow = 1 or \code{includenightcrit} when viewingwindow = 2}
+      \item{require_complete_lastnight_part5}{    
+          Boolean (default = FALSE).
+          When set to TRUE: The last WW window is excluded if recording ends 
+          between midnight and 3pm; The last OO and MM window are excluded if 
+          recording ends between midnight and 9am. This to avoid risk that recording
+          end biases the sleep estimates for the last night.
+      }
     }
   }
 }

man/g.part5.savetimeseries.Rd

@@ -12,7 +12,8 @@
                               DaCleanFile = NULL,
                               includedaycrit.part5 = 2/3,
                               ID = NULL, params_output,
-                              params_247 = NULL)
+                              params_247 = NULL,
+                              timewindow = NULL)
 }
 \arguments{
   \item{ts}{
@@ -48,6 +49,9 @@
   \item{params_247}{
     See \link{GGIR}
   }
+  \item{timewindow}{
+    See \link{GGIR}
+  }
 }
 \value{
   Function does not provide output, it only prepare data for saving

tests/testthat/test_load_check_params.R

@@ -17,7 +17,7 @@ test_that("load_params can load parameters", {
   expect_equal(length(params$params_247), 22)
   expect_equal(length(params$params_cleaning), 24)
   expect_equal(length(params$params_phyact), 14)
-  expect_equal(length(params$params_output), 21)
+  expect_equal(length(params$params_output), 22)
   expect_equal(length(params$params_general), 17)
 
   params_sleep = params$params_sleep

tests/testthat/test_recordingEndSleepHour.R

@@ -3,7 +3,7 @@ context("recordingEndSleepHour")
 test_that("recordingEndSleepHour works as expected", {
   skip_on_cran()
   #=======================
-  create_test_acc_csv(Nmin = 2.5 * 1440) # ends at 20:45
+  create_test_acc_csv(Nmin = 3.5 * 1440) # ends at 20:45
   fn = "123A_testaccfile.csv"
   tz = "Europe/Amsterdam"
   # this should NOT trigger data expansion
@@ -13,7 +13,7 @@ test_that("recordingEndSleepHour works as expected", {
        visualreport = FALSE, do.report = c(), verbose = FALSE)
   rn = dir("output_test/meta/basic/", full.names = TRUE)
   load(rn[1])
-  expect_true(nrow(M$metashort) == 43020)
+  expect_true(nrow(M$metashort) == 60300)
   # expect_true(M$metashort$timestamp[nrow(M$metashort)] == "2016-06-25T20:44:55+0200")
 
   # errors and warnings work properly
@@ -46,25 +46,25 @@ test_that("recordingEndSleepHour works as expected", {
        minimum_MM_length.part5 = 6, verbose = FALSE)
   rn = dir("output_test/meta/basic/", full.names = TRUE)
   load(rn[1])
-  expect_true(nrow(M$metashort) > 43020) # metashort is expanded
+  expect_true(nrow(M$metashort) > 60300) # metashort is expanded
 
   # expanded time is not reports
   p2 = read.csv("output_test/results/part2_daysummary.csv")
   expect_equal(p2$N.hours[nrow(p2)], 20.75)  # N hours in last day does not include the expanded time
 
   p4 = read.csv("output_test/results/part4_nightsummary_sleep_cleaned.csv")
   p4full = read.csv("output_test/results/QC/part4_nightsummary_sleep_full.csv")
-  expect_equal(nrow(p4), 2) # Night 3 is not in the part 4 reports
-  expect_equal(nrow(p4full), 2) # Night 3 is not in the part 4 reports
-  expect_equal(sum(p4$sleeponset), 41.848)
-  expect_equal(sum(p4$wakeup), 62.336)
-  expect_equal(sum(p4$guider_onset), 41.851)
-  expect_equal(sum(p4$guider_wakeup), 62.34)
-  expect_equal(sum(p4$number_sib_sleepperiod), 73)
+  expect_equal(nrow(p4), 3) # Night 3 is not in the part 4 reports
+  expect_equal(nrow(p4full), 3) # Night 3 is not in the part 4 reports
+  expect_equal(sum(p4$sleeponset), 62.77)
+  expect_equal(sum(p4$wakeup), 93.505)
+  expect_equal(sum(p4$guider_onset), 62.778)
+  expect_equal(sum(p4$guider_wakeup), 93.511)
+  expect_equal(sum(p4$number_sib_sleepperiod), 114)
   expect_true(all(is.na(p4$longitudinal_axis)))
 
   p5 = read.csv("output_test/results/part5_daysummary_MM_L40M100V400_T5A5.csv")
-  expect_equal(nrow(p5), 3) # expanded day appears in MM report
+  expect_equal(nrow(p5), 4) # expanded day appears in MM report
   expect_true(p5$dur_day_spt_min[nrow(p5)] < 23*60) # but expanded time is not accounted for in estimates
 
   #================================================================
@@ -83,26 +83,26 @@ test_that("recordingEndSleepHour works as expected", {
        minimum_MM_length.part5 = 6, verbose = FALSE)
   rn = dir("output_test/meta/basic/", full.names = TRUE)
   load(rn[1])
-  expect_true(nrow(M$metashort) > 43020) # metashort is expanded
+  expect_true(nrow(M$metashort) > 60300) # metashort is expanded
 
   # expanded time is not reports
   p2 = read.csv("output_test/results/part2_daysummary.csv")
   expect_equal(p2$N.hours[nrow(p2)], 20.75)  # N hours in last day does not include the expanded time
 
   p4 = read.csv("output_test/results/part4_nightsummary_sleep_cleaned.csv")
   p4full = read.csv("output_test/results/QC/part4_nightsummary_sleep_full.csv")
-  expect_equal(nrow(p4), 2) # Night 3 is not in the part 4 reports
-  expect_equal(nrow(p4full), 2) # Night 3 is not in the part 4 reports
-  expect_equal(sum(p4$sleeponset), 41.581)
-  expect_equal(sum(p4$wakeup), 58.924)
-  expect_equal(sum(p4$guider_inbedStart), 41.17)
-  expect_equal(sum(p4$guider_inbedEnd), 64.427)
-  expect_equal(sum(p4$number_sib_sleepperiod), 18)
-  expect_equal(sum(p4$sleepefficiency), 0.358)
-  expect_equal(sum(p4$longitudinal_axis), 6)
+  expect_equal(nrow(p4), 3) # Night 3 is not in the part 4 reports
+  expect_equal(nrow(p4full), 3) # Night 3 is not in the part 4 reports
+  expect_equal(sum(p4$sleeponset), 62.1)
+  expect_equal(sum(p4$wakeup), 81.096)
+  expect_equal(sum(p4$guider_inbedStart), 53.171)
+  expect_equal(sum(p4$guider_inbedEnd), 86.853)
+  expect_equal(sum(p4$number_sib_sleepperiod), 24)
+  expect_equal(sum(p4$sleepefficiency), 0.501)
+  expect_equal(sum(p4$longitudinal_axis), 9)
 
   p5 = read.csv("output_test/results/part5_daysummary_MM_L40M100V400_T5A5.csv")
-  expect_equal(nrow(p5), 2) # expanded day appears in MM report
+  expect_equal(nrow(p5), 3) # expanded day appears in MM report
   expect_true(p5$dur_day_spt_min[nrow(p5)] < 23*60) # but expanded time is not accounted for in estimates
 
   # test that time series output has all expected columns including the multiple angle columns
@@ -116,6 +116,19 @@ test_that("recordingEndSleepHour works as expected", {
                     "class_id", "invalid_fullwindow", "invalid_sleepperiod",
                     "invalid_wakinghours", "timestamp") %in% names(mdat)))
 
+  #========================================================
+  # Test require_complete_lastnight_part5 = TRUE
+  create_test_acc_csv(Nmin = 3 * 1440) # ends at 8:45
+  # delete config.csv such that GGIR uses its own default parameter values
+  if (file.exists("output_test/config.csv"))  unlink("output_test/config.csv", recursive = TRUE)
+  GGIR(datadir = fn, outputdir = getwd(), mode = 1:5,
+       studyname = "test", overwrite = TRUE, desiredtz = tz, 
+       verbose = FALSE, require_complete_lastnight_part5 = TRUE,
+       do.report = 5)
+  p5 = read.csv("output_test/results/part5_daysummary_MM_L40M100V400_T5A5.csv")
+  expect_equal(nrow(p5), 2) # last window is ignored from because the night ends at 8:45
+
+
   if (file.exists("output_test"))  unlink("output_test", recursive = TRUE)
   if (file.exists(fn)) file.remove(fn)
 })

vignettes/GGIRParameters.Rmd

@@ -215,6 +215,7 @@ find a description and default value for all the arguments.
 | save_ms5raw_format          | 5                 | params_output            |
 | save_ms5raw_without_invalid | 5                 | params_output            |
 | do.sibreport                | 5                 | params_output            |
+| require_complete_lastnight_part5 | 5                | params_output            |
 | visualreport_without_invalid| visualreport      | params_output            |
 | dofirstpage                 | visualreport      | params_output            |
 | visualreport                | visualreport      | params_output            |

vignettes/chapter12_TimeUseAnalysis.Rmd

@@ -45,7 +45,9 @@ GGIR part 5 facilitates the following time window definitions, which can be sele
 | “WW”       | waking up to waking up     | Each day is defined as the time from the participant wakes up a given day to the time they wake up the next day.                             |
 | “OO”       | sleep onset to sleep onset | Each day is defined as the time from the sleep onset a given day to the sleep onset of the next day.                                         |
 
-: For "WW" and "OO", the onset and waking times are guided by the estimates from part 4, but if they are missing, part 5 will attempt to retrieve the estimate from the guider method. Note that the parameter `timewindow` can consist of one of the options beforementioned or any combination of them, for example, the default value is `timewindow = c("MM", "WW")`.
+For "WW" and "OO", the onset and waking times are guided by the estimates from part 4, but if they are missing, part 5 will attempt to retrieve the estimate from the guider method. Note that the parameter `timewindow` can consist of one of the options beforementioned or any combination of them, for example, the default value is `timewindow = c("MM", "WW")`.
+
+When recordings end in the night or early morning the sleep estimates for the night are likely affected. For example, if a recording ends at 10am we cannot be sure that the participant did not sleep until after 10am, and if a recording ends at 2am we cannot be sure that the sleep onset time was reliably estimated. To handle this and ignore the final window in the data, set parameter `require_complete_lastnight_part5 = TRUE` (not default).
 
 ### Defining segments within the MM window