Add day segment analyses to GGIR part5 for WW and OO defined days

NEWS.md

@@ -1,3 +1,7 @@
+# CHANGES IN GGIR VERSION 3.3-?
+
+- Part 5: Allow for day segment analysis for the WW and OO window definition. #1407
+
 # CHANGES IN GGIR VERSION 3.3-8
 
 - Updated documentation to point to the new web dashboard URL (#1509)

R/g.part5.R

@@ -503,16 +503,20 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(),
                           for (si in next_si:(next_si + length(segments) - 1)) {
                             fi = 1
                             current_segment_i = si - next_si + 1
-                            segStart = segments[[current_segment_i]][1]
-                            segEnd = segments[[current_segment_i]][2]
+                            Nindices = length(segments[[current_segment_i]])
+                            segStart = segments[[current_segment_i]][seq(1, Nindices, by = 2)]
+                            segEnd = segments[[current_segment_i]][seq(2, Nindices, by = 2)]
+                            Nsegments = pmin(length(segStart), length(segEnd))
                             extraRowsNeeded = max(c(si, di)) - nrow(dsummary)
                             if (extraRowsNeeded > 0) {
                               dsummary = rbind(dsummary, matrix(data = "", nrow = extraRowsNeeded, ncol = ncol(dsummary)))
                             }
-                            if (timewindowi == "MM" & si > 1) { # because first segment is always full window
+                            if (timewindowi %in% c("MM", "WW") & si > 1) { # because first segment is always full window
                               if (("segment" %in% colnames(ts)) == FALSE) ts$segment = NA
-                              if (!is.na(segStart) && !is.na(segEnd)) {
-                                ts$segment[segStart:segEnd] = si
+                              for (gi in 1:Nsegments) {
+                                if (!is.na(segStart[gi]) && !is.na(segEnd[gi])) {
+                                  ts$segment[segStart[gi]:segEnd[gi]] = si
+                                }
                               }
                             }
                             # Already store basic information about the file
@@ -788,7 +792,8 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(),
                            "g.part5.savetimeseries", "g.part5.wakesleepwindows",
                            "g.part5.onsetwaketiming", "g.part5_analyseSegment",
                            "g.part5_initialise_ts", "g.part5.analyseRest",
-                           "g.fragmentation", "g.intensitygradient")
+                           "g.fragmentation", "g.intensitygradient",
+                           "g.part4_extractid", "markerButtonForRest")
       errhand = 'stop'
     }
     i = 0 # declare i because foreach uses it, without declaring it

R/g.part5.definedays.R

@@ -64,61 +64,6 @@ g.part5.definedays = function(nightsi, wi, indjump, epochSize, qqq_backup = c(),
       qqq[2] = Nts
       lastDay = TRUE
     }
-    qqq_backup = qqq
-    # in MM, also define segments of the day based on qwindow
-    if (!is.na(qqq[1]) & !is.na(qqq[2])) {
-      segments_timing = NULL
-      if (qqq[2] > Nts) qqq[2] = Nts
-      fullQqq = qqq[1]:qqq[2]
-      firstepoch = format(ts$time[qqq[1]],  "%H:%M:%S")
-      lastepoch = format(ts$time[qqq[2]],  "%H:%M:%S")
-      qnames = NULL
-      if (is.data.frame(qwindow)) {
-        date_of_interest = substr(ts$time[qqq[1]], 1, 10)
-        qdate = which(qwindow$ID == ID & qwindow$date == date_of_interest)
-        if (length(qdate) == 1) { # if ID/date matched with activity log
-          qnames = unlist(qwindow$qwindow_names[qdate])
-          qwindow = unlist(qwindow$qwindow_values[qdate])
-          qwindow_order = order(qwindow)
-          qwindow = qwindow[qwindow_order]
-          qnames = qnames[qwindow_order]
-        } else { # if ID/date not correctly matched with activity log
-          qwindow = c(0, 24)
-        }
-      } else {
-        qwindow = sort(qwindow)
-        if (qwindow[1] != 0) qwindow = c(0, qwindow)
-        if (qwindow[length(qwindow)] != 24) qwindow = c(qwindow, 24)
-      }
-      # define segments timing in H:M:S format
-      breaks = qwindow2timestamp(qwindow, epochSize)
-      startOfSegments = breaks[-length(breaks)]
-      endOfSegments = subtractEpochFromTimeName(breaks[-1], epochSize)
-      if (length(startOfSegments) > 1) { # when qwindow segments are defined, add fullwindow at the beginning
-        startOfSegments = c(firstepoch, startOfSegments)
-        endOfSegments = c(lastepoch, endOfSegments)
-      } 
-      segments_timing = paste(startOfSegments, endOfSegments, sep = "-")
-      # define segment names based on qnames or segmentX
-      if (is.null(qnames)) {
-        segments_names = paste0("segment", 0:(length(segments_timing) - 1))
-        segments_names = gsub("segment0", "MM", segments_names)
-      } else {
-        segments_names = c("MM", paste(qnames[-length(qnames)], qnames[-1], sep = "-"))
-      }
-      # Get indices in ts for segments start and end limits
-      hms = format(ts$time[fullQqq], format = "%H:%M:%S")
-      segments = vector("list", length = length(segments_timing))
-      names(segments) = segments_timing
-      for (si in 1:length(segments_timing)) {
-        s0s1 = unlist(strsplit(segments_timing[si], split = "[-]"))
-        s0s1 = format(s0s1, format = "%H:%M:%S")
-        # tryCatch is needed in the case that the segment is not available in ts,
-        # then a no non-missing values warning would be triggered by the which function
-        segments[[si]] = tryCatch(range(fullQqq[which(hms >= s0s1[1] & hms <= s0s1[2])]), #segStart and segEnd
-                                  warning = function(w) rep(NA, 2))
-      }
-    }
   } else if (timewindowi == "WW" || timewindowi == "OO") {
     windowEdge = ifelse(timewindowi == "WW", yes = -1, no = 1)
     if (wi <= (Nwindows - 1)) { # all full windows
@@ -131,15 +76,87 @@ g.part5.definedays = function(nightsi, wi, indjump, epochSize, qqq_backup = c(),
       # not informative.
       qqq = c(NA, NA)
     }
-    # build up segments
-    if (!is.na(qqq[1]) & !is.na(qqq[2])) {
-      segments = list(qqq)
-      start = format(ts$time[qqq[1]], "%H:%M:%S")
-      end = format(ts$time[qqq[2]], "%H:%M:%S")
-      names(segments) = paste(start, end, sep = "-")
-      segments_names = timewindowi
+    if (wi == length(which(diff(ts$diur) == windowEdge)) - 1) {
+      lastDay = TRUE
+    }
+    if (qqq[2] >= Nts - 1) {
+      qqq[2] = Nts
+      lastDay = TRUE
+    }
+  }
+  qqq_backup = qqq
+  # in MM, also define segments of the day based on qwindow
+  if (!is.na(qqq[1]) & !is.na(qqq[2])) {
+    segments_timing = NULL
+    if (qqq[2] > Nts) qqq[2] = Nts
+    fullQqq = qqq[1]:qqq[2]
+    firstepoch = format(ts$time[qqq[1]],  "%H:%M:%S")
+    lastepoch = format(ts$time[qqq[2]],  "%H:%M:%S")
+    qnames = NULL
+    if (is.data.frame(qwindow)) {
+      date_of_interest = substr(ts$time[qqq[1]], 1, 10)
+      qdate = which(qwindow$ID == ID & qwindow$date == date_of_interest)
+      if (length(qdate) == 1) { # if ID/date matched with activity log
+        qnames = unlist(qwindow$qwindow_names[qdate])
+        qwindow = unlist(qwindow$qwindow_values[qdate])
+        qwindow_order = order(qwindow)
+        qwindow = qwindow[qwindow_order]
+        qnames = qnames[qwindow_order]
+      } else { # if ID/date not correctly matched with activity log
+        qwindow = c(0, 24)
+      }
+    } else {
+      qwindow = sort(qwindow)
+      if (qwindow[1] != 0) qwindow = c(0, qwindow)
+      if (qwindow[length(qwindow)] != 24) qwindow = c(qwindow, 24)
+    }
+    # define segments timing in H:M:S format
+    breaks = qwindow2timestamp(qwindow, epochSize)
+    startOfSegments = breaks[-length(breaks)]
+    endOfSegments = subtractEpochFromTimeName(breaks[-1], epochSize)
+    if (length(startOfSegments) > 1) { # when qwindow segments are defined, add fullwindow at the beginning
+      startOfSegments = c(firstepoch, startOfSegments)
+      endOfSegments = c(lastepoch, endOfSegments)
+    } 
+    segments_timing = paste(startOfSegments, endOfSegments, sep = "-")
+    # define segment names based on qnames or segmentX
+    if (is.null(qnames)) {
+      segments_names = paste0(paste0(timewindowi, "segment"), 0:(length(segments_timing) - 1))
+      segments_names = gsub(paste0(timewindowi, "segment0"), timewindowi, segments_names)
+    } else {
+      segments_names = c(timewindowi, paste(paste0(timewindowi, "segment"), 
+                                            qnames[-length(qnames)], qnames[-1], sep = "-"))
+    }
+    # Get indices in ts for segments start and end limits
+    hms = format(ts$time[fullQqq], format = "%H:%M:%S")
+    segments = vector("list", length = length(segments_timing))
+    names(segments) = segments_timing
+    for (si in 1:length(segments_timing)) {
+      s0s1 = unlist(strsplit(segments_timing[si], split = "[-]"))
+      s0s1 = format(s0s1, format = "%H:%M:%S")
+      # tryCatch is needed in the case that the segment is not available in ts,
+      # then a no non-missing values warning would be triggered by the which function
+      if (si == 1) {
+        test_seg_condition = NULL
+        segments[[si]] = range(fullQqq)
+      } else {
+        test_seg_condition = which(hms >= s0s1[1] & hms <= s0s1[2])
+        if (any(diff(test_seg_condition) > 1) & segments_names[si] %in% c("WW", "OO") == FALSE) {
+          # window is more than 24 hours and segment occurs twice
+          # we keep both index ranges
+          jump = which(diff(test_seg_condition) > 1)
+          test_seg_condition1 = test_seg_condition[1:jump]
+          test_seg_condition2 = test_seg_condition[(jump + 1):length(test_seg_condition)]
+          segments[[si]] = c(tryCatch(range(fullQqq[test_seg_condition1]), #segStart and segEnd
+                                      warning = function(w) rep(NA, 2)),
+                             tryCatch(range(fullQqq[test_seg_condition2]), #segStart and segEnd
+                                      warning = function(w) rep(NA, 2)))
+        } else {
+          segments[[si]] = tryCatch(range(fullQqq[test_seg_condition]), #segStart and segEnd
+                                    warning = function(w) rep(NA, 2))
+        }
+      }
     }
-    if (wi >= Nwindows) lastDay = TRUE
   }
   return(invisible(list(qqq = qqq, qqq_backup = qqq_backup, lastDay = lastDay,
                         segments = segments, segments_names = segments_names)))

R/g.part5.lux_persegment.R

@@ -24,14 +24,19 @@ g.part5.lux_persegment = function(ts, sse, LUX_day_segments, epochSize, desiredt
     sse = sse[1:Nepochperday]
   }
   # expanding first time segment
-  expected_N_seg1 = (rep(LUX_day_segments,2)[which(LUX_day_segments == first_hour_seg[1])+1] - LUX_day_segments[which(LUX_day_segments == first_hour_seg[1])]) * 60 * (60/epochSize)
-  actual_N_seg1 = length(which(first_hour_seg[1:round(Nepochperday*0.66)] == first_hour_seg[1]))
+  isFirstHour = which(LUX_day_segments == first_hour_seg[1])
+  expected_N_seg1 = (rep(LUX_day_segments, 2)[isFirstHour + 1] - 
+                       LUX_day_segments[isFirstHour]) * 60 * (60 / epochSize)
+  actual_N_seg1 = length(which(first_hour_seg[1:round(Nepochperday * 0.66)] == first_hour_seg[1]))
   missingN_seg1 = expected_N_seg1 - actual_N_seg1
-  if (missingN_seg1 > 0) {
-    extension = (sse[1]-missingN_seg1):(sse[1]-1)
-    ts$lightpeak[extension] = 0 # ensure that light during SPT is treated as zeros
-    sse = c(extension,sse)
-    first_hour_seg = c(rep(first_hour_seg[1], missingN_seg1), first_hour_seg)
+  if (length(sse) > 0) {
+    if (sse[1] - missingN_seg1 < 0) missingN_seg1 = missingN_seg1 + (sse[1] - missingN_seg1 - 1)
+    if (missingN_seg1 > 0) {
+      extension = (sse[1] - missingN_seg1):(sse[1] - 1)
+      ts$lightpeak[extension] = 0 # ensure that light during SPT is treated as zeros
+      sse = c(extension,sse)
+      first_hour_seg = c(rep(first_hour_seg[1], missingN_seg1), first_hour_seg)
+    }
   }
   # removing again the data after 24 hours to avoid overlap  Nepochperday = 24 * (3600 / ws3)
   if (length(first_hour_seg) > Nepochperday) {

R/g.part5.onsetwaketiming.R

@@ -16,6 +16,7 @@ g.part5.onsetwaketiming = function(qqq, ts, min, sec, hour, timewindowi) {
   if (timewindowi == "WW") {
     # For WW wake is by definition the start and end of the window
     wakei = qqq[2] + 1
+    if (wakei > length(hour)) wakei = length(hour) # in the case if MM use first wake-up time
   } else {
     # For OO and MM wake needs to be search in the window
     wakei = c(qqq[1]:qqq[2])[which(diff(ts$diur[qqq[1]:(qqq[2] - 1)]) == -1) + 1]

R/g.part5_analyseSegment.R

@@ -19,8 +19,14 @@ g.part5_analyseSegment = function(indexlog, timeList, levelList,
   qqq = indexlog$winStartEnd
   si = indexlog$segIndex1
   current_segment_i = indexlog$segIndex2
-  segStart = indexlog$segStartEnd[1]
-  segEnd = indexlog$segStartEnd[2]
+  Nsegments = length(indexlog$segStartEnd) / 2
+  if (Nsegments == 1) {
+    segStart = indexlog$segStartEnd[1]
+    segEnd = indexlog$segStartEnd[2]
+  } else {
+    segStart = indexlog$segStartEnd[1:2]
+    segEnd = indexlog$segStartEnd[3:4]
+  }
   fi = indexlog$columnIndex
 
   # timeList
@@ -114,9 +120,11 @@ g.part5_analyseSegment = function(indexlog, timeList, levelList,
                                  sumSleep$acc_available[dayofinterest])
     ds_names[fi:(fi + 5)] = c("night_number", "daysleeper", "cleaningcode",
                               "guider", "sleeplog_used", "acc_available");      fi = fi + 6
-    if (!is.na(segStart) & !is.na(segEnd)) {
-      # segment available in time series
-      ts$guider[segStart:segEnd] = sumSleep$guider[dayofinterest] # add guider also to timeseries
+    for (gi in 1:Nsegments) {
+      if (!is.na(segStart[gi]) & !is.na(segEnd[gi])) {
+        # segment available in time series
+        ts$guider[segStart[gi]:segEnd[gi]] = sumSleep$guider[dayofinterest] # add guider also to timeseries
+      }
     }
   } else {
     dsummary[si,fi:(fi + 5)] = rep(NA, 6)
@@ -139,12 +147,17 @@ g.part5_analyseSegment = function(indexlog, timeList, levelList,
   dsummary[si, fi:(fi + 2)] = c(TRLi, TRMi, TRVi)
   ds_names[fi:(fi + 2)] = c("TRLi", "TRMi", "TRVi")
   fi = fi + 3
-  wlih = ((qqq2 - qqq1) + 1)/((60/ws3new) * 60)
-  if (!is.na(qqq1)) {
-    if (qqq1 > length(LEVELS)) qqq1 = length(LEVELS)
-    sse = qqq1:qqq2
-  } else {
-    sse = NULL
+  for (gi in 1:Nsegments) {
+    if (!is.na(qqq1[gi])) {
+      if (qqq1[gi] > length(LEVELS)) qqq1[gi] = length(LEVELS)
+      if (gi == 1) {
+        sse = qqq1[gi]:qqq2[gi]
+      } else {
+        sse = c(sse, qqq1[gi]:qqq2[gi])
+      }
+    } else {
+      sse = NULL
+    }
   }
   doNext = FALSE
   if (length(sse) >= 1) { #next
@@ -288,11 +301,19 @@ g.part5_analyseSegment = function(indexlog, timeList, levelList,
     # }
     #===============================================
     # QUANTILES...
-    WLH = ((qqq2 - qqq1) + 1)/((60/ws3new) * 60)
-    if (WLH <= 1) WLH = 1.001
-    dsummary[si, fi] = quantile(ts$ACC[sse],probs = ((WLH - 1)/WLH), na.rm = TRUE)
+    if (Nsegments == 1) {
+      WLH = ((qqq2 - qqq1) + 1) / ((60 / ws3new) * 60)
+      if (WLH <= 1) WLH = 1.001
+    } else {
+      WLH = ((sum(qqq2 - qqq1)) + 1) / ((60 / ws3new) * 60)
+    }
+    if (segments_names[si] %in% c("WW", "OO", "MM")) {
+      dsummary[si, fi] = quantile(ts$ACC[sse], probs = ((WLH - 1)/WLH), na.rm = TRUE)
+    }
     ds_names[fi] = paste("quantile_mostactive60min_mg", sep = "");      fi = fi + 1
-    dsummary[si, fi] = quantile(ts$ACC[sse],probs = ((WLH - 0.5)/WLH), na.rm = TRUE)
+    if (segments_names[si] %in% c("WW", "OO", "MM")) {
+      dsummary[si, fi] = quantile(ts$ACC[sse],probs = ((WLH - 0.5)/WLH), na.rm = TRUE)
+    }
     ds_names[fi] = paste("quantile_mostactive30min_mg", sep = "");      fi = fi + 1
     #===============================================
     NANS = which(is.nan(dsummary[si,]) == TRUE) #average of no values will results in NaN
@@ -429,13 +450,18 @@ g.part5_analyseSegment = function(indexlog, timeList, levelList,
       }
       fi = fi + Nluxt
       # LUX per segment of the day
-      luxperseg = g.part5.lux_persegment(ts, sse,
-                                         LUX_day_segments = params_247[["LUX_day_segments"]],
-                                         epochSize = ws3new,
-                                         desiredtz = params_general[["desiredtz"]])
-      dsummary[si, fi:(fi + (length(luxperseg$values) - 1))] = luxperseg$values
-      ds_names[fi:(fi + (length(luxperseg$values) - 1))] = luxperseg$names
-      fi = fi + length(luxperseg$values)
+      if (timewindowi %in% c("WW", "OO", "MM")) {
+        # LUX per segment of the day
+        luxperseg = g.part5.lux_persegment(ts, sse,
+                                           LUX_day_segments = params_247[["LUX_day_segments"]],
+                                           epochSize = ws3new,
+                                           desiredtz = params_general[["desiredtz"]])
+        if (!is.null(unlist(luxperseg$values))) {
+          dsummary[si, fi:(fi + (length(luxperseg$values) - 1))] = luxperseg$values
+          ds_names[fi:(fi + (length(luxperseg$values) - 1))] = luxperseg$names
+        }
+        fi = fi + length(luxperseg$values)
+      }
     }
 
     #===============================================