Sleep minor revisions

NEWS.md

@@ -1,8 +1,21 @@
+# CHANGES IN GGIR VERSION 3.1-?
+
+- Part 3: 
+
+  - Identification of daylight saving time days in the detection of the spt is improved as it had the risk to misclassify partial last days as daylight saving time.
+
+  - In days classified as daysleeper, the window over which the fraction_night_invalid is calculated now also shifts to 6pm-6pm, as it used to report the nonwear within 12pm-12pm.
+
+
 # CHANGES IN GGIR VERSION 3.1-8
 
-- Part 4: Parameter sib_must_fully_overlap_with_TimeInBed added to control whether sib should overlap fully with the start and/or end of time in bed to be considered sleep (default TRUE),
+- Part 4: 
+
+  - Parameter sib_must_fully_overlap_with_TimeInBed added to control whether sib should overlap fully with the start and/or end of time in bed to be considered sleep (default TRUE),
 this is consistent with functionality in the past. #1223
 
+  - Handle unexpected combinations of sleep diary column names. For example, if only inbed and wakeup are available then treat these as the time in bed period, while if only sleeponset and outbed are available then these are treated as SPT window.
+
 - Part 5: Expand functionality for exploring possibility of nap detection, this includes the addition of new parameter possible_nap_gap.
 
 - Visual report: Added new visualreport that is automatically generated when visualreport = TRUE and intended to eventually replace the problematic legacy report. Add parameter old_visualreport to turn off the old visualreport generation. #1173
@@ -14,8 +27,6 @@ this is consistent with functionality in the past. #1223
 - Part 5 and 6: part6_threshold_combi when not specified now defaults to first threshold as
 specified for light, moderate and vigorous intensity respectively.
 
-- Part 4: Handle unexpected combinations of sleep diary column names. For example, if only inbed and wakeup are available then treat these as the time in bed period, while if only sleeponset and outbed are available then these are treated as SPT window.
-
 # CHANGES IN GGIR VERSION 3.1-7
 
 - Part 3: Improved handling of DST, #1225

R/HASPT.R

@@ -185,4 +185,4 @@ HASPT = function(angle, sptblocksize = 30, spt_max_gap = 60, ws3 = 5,
   }
   invisible(list(SPTE_start = SPTE_start, SPTE_end = SPTE_end, tib.threshold = tib.threshold,
                  part3_guider = part3_guider))
-}
+}

R/g.sib.det.R

@@ -21,37 +21,27 @@ g.sib.det = function(M, IMP, I, twd = c(-12, 12),
   # So, out of convenience I keep the object name SPTE. 
   dstime_handling_check = function(tmpTIME = tmpTIME, spt_estimate = spt_estimate, tz = c(),
                                    calc_SPTE_end = c(), calc_SPTE_start = c()) {
-    time_SPTE_start = iso8601chartime2POSIX(tmpTIME[spt_estimate$SPTE_start], tz = tz)
-    time_SPTE_end = iso8601chartime2POSIX(tmpTIME[spt_estimate$SPTE_end], tz = tz)
-    time_SPTE_end_hr = as.numeric(format(time_SPTE_end, format = "%H"))
-    time_SPTE_start_hr = as.numeric(format(time_SPTE_start, format = "%H"))
-    t0 = as.numeric(format(iso8601chartime2POSIX(tmpTIME[1], tz = tz), format = "%H"))
-    t1 = as.numeric(format(iso8601chartime2POSIX(tmpTIME[length(tmpTIME)], tz = tz), format = "%H"))
-    if (is.na(t1) == TRUE | is.na(t0) == TRUE) {
-      # if timestamp is stored in POSIX format (older versions of GGIR prior 2017)
-      # for example when milestone data generated by old packages is now processed
-      # then this part of the code is to facilitate consistency
-      time_SPTE_start = as.POSIXlt(tmpTIME[spt_estimate$SPTE_start], origin = "1970-01-01", tz = tz)
-      time_SPTE_end = as.POSIXlt(tmpTIME[spt_estimate$SPTE_end], origin = "1970-01-01", tz = tz)
-      time_SPTE_end_hr = as.numeric(format(time_SPTE_end, format = "%H"))
-      time_SPTE_start_hr = as.numeric(format(time_SPTE_start, format = "%H"))
-      t0 = as.numeric(format(as.POSIXlt(tmpTIME[1],origin = "1970-01-01", tz = tz), format = "%H"))
-      t1 = as.numeric(format(as.POSIXlt(tmpTIME[length(tmpTIME)], origin = "1970-01-01", tz = tz), format = "%H"))
+    # Function to readjust estimated SPTE_start and SPTE_end in DST days
+    t = tmpTIME[c(1, spt_estimate$SPTE_start, spt_estimate$SPTE_end)]
+    timezone = format(GGIR::iso8601chartime2POSIX(t, tz = tz), "%z")
+    if (length(unique(timezone)) == 1) {
+      return(c(calc_SPTE_start, calc_SPTE_end))
+    } else {
+      sign = ifelse(substr(timezone, 1, 1) == "-", -1, 1)
+      hours = as.numeric(substr(timezone, 2, 3))
+      minutes = as.numeric(substr(timezone, 4, 5))
+      offset = sign * (hours + minutes / 60)
+      offset = offset[2:3] - offset[1]
+      return(c(calc_SPTE_start + offset[1], calc_SPTE_end + offset[2]))
     }
-    Nhoursdata = floor(length(tmpTIME) / (3600/ws3))
-    delta_t1_t0 = (t1 + 24) - t0
-    if (length(time_SPTE_end_hr) > 0 & length(time_SPTE_start_hr) > 0) {
-      if (Nhoursdata > (delta_t1_t0 + 0.1) &  # time has moved backward (autumn) so SPTE_end also needs to move backward
-          (time_SPTE_end_hr > 1 | time_SPTE_end_hr < 12) &
-          (time_SPTE_start_hr < 1 | time_SPTE_start_hr > 12)) { #extra DST hour not recognized
-        calc_SPTE_end = calc_SPTE_end - 1
-      } else if (Nhoursdata + 0.1 < delta_t1_t0 &  # time has moved forward (spring) so SPTE_end also needs to move forward
-                 (time_SPTE_end_hr > 1 | time_SPTE_end_hr < 12) &
-                 (time_SPTE_start_hr < 1 | time_SPTE_start_hr > 12)) { #missing DST hour not recognized
-        calc_SPTE_end = calc_SPTE_end + 1
-      }
+  }
+  decide_guider = function(HASPT.algo, nonwear_percentage) {
+    # Function to decide the guider to use in the case that HASPT.algo is of length 2
+    if (HASPT.algo[1] == "NotWorn" && nonwear_percentage < 25 && length(HASPT.algo) == 2) {
+      return(2)
+    } else {
+      return(1)
     }
-    return(calc_SPTE_end)
   }
   #==================================================
   # get variables
@@ -249,13 +239,7 @@ g.sib.det = function(M, IMP, I, twd = c(-12, 12),
         detection.failed = FALSE
         # Calculate nonwear percentage for this window
         nonwear_percentage = (length(which(invalid[qqq1:qqq2] == 1)) /  (qqq2 - qqq1 + 1)) * 100
-        guider_to_use = 1
-        if (params_sleep[["HASPT.algo"]][1] == "NotWorn" &&
-            nonwear_percentage < 25 &&
-            length(params_sleep[["HASPT.algo"]]) == 2) {
-          # Nonwear percentage was low, so use alternative guider specified as second element
-          guider_to_use = 2
-        }
+        guider_to_use = decide_guider(params_sleep[["HASPT.algo"]], nonwear_percentage)
         #------------------------------------------------------------------
         # calculate L5 because this is used as back-up approach
         tmpACC = ACC[qqq1:qqq2]
@@ -310,8 +294,12 @@ g.sib.det = function(M, IMP, I, twd = c(-12, 12),
             if (newqqq2 > length(anglez)) newqqq2 = length(anglez)
             # only try to extract SPT again if it is possible to extract a window of more than 23 hour
             if (newqqq2 < length(anglez) & (newqqq2 - newqqq1) > (23*(3600/ws3)) ) {
+              # Recalculate nonwear percentage for new window (6pm to 6pm)
+              nonwear_percentage = (length(which(invalid[newqqq1:newqqq2] == 1)) /  (newqqq2 - newqqq1 + 1)) * 100
+              guider_to_use = decide_guider(params_sleep[["HASPT.algo"]], nonwear_percentage)
+              # new TIME (6pm tp 6pm)
               tmpTIME = time[newqqq1:newqqq2]
-              if (params_sleep[["HASPT.algo"]][1] != "NotWorn") {
+              if (params_sleep[["HASPT.algo"]][guider_to_use] != "NotWorn") {
                 tmpANGLE = anglez[newqqq1:newqqq2]
                 if (do.HASPT.hip == TRUE) {
                   if (params_sleep[["longitudinal_axis"]] == 1) {
@@ -323,7 +311,7 @@ g.sib.det = function(M, IMP, I, twd = c(-12, 12),
               }
               spt_estimate_tmp = HASPT(angle = tmpANGLE, ws3 = ws3,
                                        sptblocksize = sptblocksize, spt_max_gap = spt_max_gap,
-                                       HASPT.algo = params_sleep[["HASPT.algo"]][1],
+                                       HASPT.algo = params_sleep[["HASPT.algo"]][guider_to_use],
                                        invalid = invalid[newqqq1:newqqq2],
                                        HDCZA_threshold = params_sleep[["HDCZA_threshold"]],
                                        HASPT.ignore.invalid = params_sleep[["HASPT.ignore.invalid"]],
@@ -342,7 +330,8 @@ g.sib.det = function(M, IMP, I, twd = c(-12, 12),
               daysleep_offset  = 0
             }
           }
-          if (qqq1 == 1 && partialFirstDay == TRUE) {  # only use startTimeRecord if the start of the block send into SPTE was after noon
+          if (qqq1 == 1 && partialFirstDay == TRUE) {  
+            # only use startTimeRecord if the start of the block send into SPTE was after noon
             startTimeRecord = unlist(iso8601chartime2POSIX(IMP$metashort$timestamp[1], tz = desiredtz))
             startTimeRecord = sum(as.numeric(startTimeRecord[c("hour", "min", "sec")]) / c(1, 60, 3600))
             daysleep_offset = daysleep_offset + startTimeRecord
@@ -351,9 +340,12 @@ g.sib.det = function(M, IMP, I, twd = c(-12, 12),
           }
           SPTE_end[sptei] = (spt_estimate$SPTE_end / (3600 / ws3)) + daysleep_offset
           SPTE_start[sptei] = (spt_estimate$SPTE_start / (3600 / ws3)) + daysleep_offset
-          SPTE_end[sptei] = dstime_handling_check(tmpTIME = tmpTIME, spt_estimate = spt_estimate,
-                                                  tz = desiredtz, calc_SPTE_end = SPTE_end[sptei],
-                                                  calc_SPTE_start = SPTE_start[sptei])
+          SPTE_dst = dstime_handling_check(tmpTIME = tmpTIME, spt_estimate = spt_estimate,
+                                           tz = desiredtz, 
+                                           calc_SPTE_end = SPTE_end[sptei],
+                                           calc_SPTE_start = SPTE_start[sptei])
+          SPTE_start[sptei] = SPTE_dst[1]
+          SPTE_end[sptei] = SPTE_dst[2]
           tib.threshold[sptei] = spt_estimate$tib.threshold
           part3_guider[sptei] = spt_estimate$part3_guider
         }
@@ -372,4 +364,4 @@ g.sib.det = function(M, IMP, I, twd = c(-12, 12),
                  SPTE_end = SPTE_end, SPTE_start = SPTE_start,
                  tib.threshold = tib.threshold, longitudinal_axis = params_sleep[["longitudinal_axis"]],
                  part3_guider = part3_guider))
-}
+}