wadpac · vincentvanhees · Dec 3, 2024 · Oct 29, 2024 · Oct 30, 2024 · Oct 30, 2024
diff --git a/NEWS.md b/NEWS.md
@@ -1,3 +1,7 @@
+# CHANGES IN GGIR VERSION 3.1-?
+
+- Part 2: Add parameters nonwearFiltermaxHours and nonwearFilterWindow to give user the option to filter short lasting nighttime nonwear #1218.
+
 # CHANGES IN GGIR VERSION 3.1-7
 
 - Part 3: Improved handling of DST, #1225

diff --git a/R/GGIR.R b/R/GGIR.R
@@ -360,7 +360,7 @@ GGIR = function(mode = 1:5, datadir = c(), outputdir = c(),
       if (length(f0) == 0) f0 = 1
       if (f1 == 0) f1 = N.files.ms2.out
       if (length(params_247[["qwindow"]]) > 2 |
-          is.character(params_247[["qwindow"]]) |
+          (is.character(params_247[["qwindow"]]) && length(grep(pattern = "onlyfilter|filteronly", x = params_247[["qwindow"]])) == 0) |
           (length(params_247[["qwindow"]]) == 2 & !all(c(0, 24) %in% params_247[["qwindow"]]))) {
         store.long = TRUE
       } else {

diff --git a/R/check_params.R b/R/check_params.R
@@ -104,7 +104,8 @@ check_params = function(params_sleep = c(), params_metrics = c(),
   if (length(params_cleaning) > 0) {
     numeric_params = c("includedaycrit", "ndayswindow", "data_masking_strategy", "maxdur", "hrs.del.start",
                        "hrs.del.end", "includedaycrit.part5", "minimum_MM_length.part5",
-                       "includenightcrit", "max_calendar_days", "includecrit.part6", "includenightcrit.part5")
+                       "includenightcrit", "max_calendar_days", "includecrit.part6", "includenightcrit.part5",
+                       "nonwearFiltermaxHours", "nonwearFilterWindow")
     boolean_params = c("excludefirstlast.part5", "do.imp", "excludefirstlast",
                        "excludefirst.part4", "excludelast.part4", "nonWearEdgeCorrection")
     character_params = c("data_cleaning_file", "TimeSegments2ZeroFile")
@@ -283,6 +284,27 @@ check_params = function(params_sleep = c(), params_metrics = c(),
                      " a fraction of the day between zero and one or the number ",
                      "of hours in a day."))
     }
+    if (!is.null(params_cleaning[["nonwearFiltermaxHours"]])) {
+      if (params_cleaning[["nonwearFiltermaxHours"]] < 0 ||
+          params_cleaning[["nonwearFiltermaxHours"]] > 12) {
+        stop("Parameters nonwearFiltermaxHours is expected to have a value > 0 and < 12")
+      }
+      if (!is.null(params_cleaning[["nonwearFilterWindow"]])) {
+        if (length(params_cleaning[["nonwearFilterWindow"]]) != 2) {
+          stop("Parameter nonwearFilterWindow does not have expected length of 2, please fix.", call. = FALSE)
+        }
+        if (params_cleaning[["nonwearFilterWindow"]][1] < params_cleaning[["nonwearFilterWindow"]][2] &&
+            params_cleaning[["nonwearFilterWindow"]][2] > 18 &&
+            params_cleaning[["nonwearFilterWindow"]][1] < 12) {
+          warning(paste0("The NonwearFilter applied to window starting at ", 
+                         params_cleaning[["nonwearFilterWindow"]][1], " and ending at ",
+                         params_cleaning[["nonwearFilterWindow"]][2], 
+                         " this is probably not the night, please check that order of",
+                         " values in nonwearFilterWindow is correct"), call. = FALSE)
+        }
+      }
+    }
+
     if (params_cleaning[["includenightcrit.part5"]] < 0) {
       stop("\nNegative value of includenightcrit.part5 is not allowed, please change.")
     } else if (params_cleaning[["includenightcrit.part5"]]  > 24) {

diff --git a/R/filterNonwearNight.R b/R/filterNonwearNight.R
@@ -0,0 +1,103 @@
+filterNonwearNight = function(r1, metalong, qwindowImp, desiredtz,
+                              params_cleaning, ws2) {
+  nonwearFiltermaxHours  = params_cleaning[["nonwearFiltermaxHours"]]
+  nonwearFilterWindow = params_cleaning[["nonwearFilterWindow"]]
+  nonwearEventsFiltered = nonwearHoursFiltered = 0
+  # Identify method to use
+  if (!is.null(nonwearFilterWindow)) {
+    filter_method = 1 # set window as provided by user
+  } else {
+    # If nonwearFilterWindow is not provided then
+    # look for qwindow is available
+    if (is.null(qwindowImp)) {
+      stop(paste0("Please specify parameter nonwearFilterWindow or ",
+                  "qwindow as diary with columns to define the window for ",
+                  "filtering short nonwear. See documentation for ",
+                  "parameter nonwearFiltermaxHours"), call. = FALSE)
+    }
+    if (inherits(qwindowImp, "data.frame")) {
+      filter_method = 2 # Use it as a dataframe
+    } else {
+      filter_method = 1 
+      nonwearFilterWindow = qwindowImp
+    }
+  }
+  # Get metalong timestamps
+  metalong$time_POSIX = iso8601chartime2POSIX(metalong$timestamp, tz = desiredtz)
+  metalong$hour = as.numeric(format(metalong$time_POSIX, "%H")) + as.numeric(format(metalong$time_POSIX, "%M")) / 60
+  # Mask short nonwear during the night
+  x = rle(as.numeric(r1))
+  if (length(x$values) != 1) {
+    r1B = as.data.frame(lapply(x, rep, times =  x$lengths))
+    r1B$hr = metalong$hour
+    r1B$filter_method = filter_method
+    r1B$lengths = (r1B$lengths * ws2) / 3600 # convert unit from epoch to hours
+    r1B$filterWindow = 0
+    # Define night window
+    if (filter_method == 1) { 
+      # set window same for all recordings
+      if (nonwearFilterWindow[1] > nonwearFilterWindow[2]) {
+        r1B$filterWindow[which(r1B$hr >= nonwearFilterWindow[1] |
+                                 r1B$hr < nonwearFilterWindow[2])] = 1
+      } else {
+        r1B$filterWindow[which(r1B$hr >= nonwearFilterWindow[1] &
+                                 r1B$hr < nonwearFilterWindow[2])] = 1
+      }
+    } else if (filter_method == 2) { 
+      # window defined by diary input from qwindow specific per participant
+      r1B$date = as.Date(metalong$time_POSIX)
+      for (qi in 1:nrow(qwindowImp)) {
+        qwindow_temp = qwindowImp$qwindow_values[[qi]]
+        # Only consider window when both start and end are reported
+        available = c(length(grep(pattern = "inbed|sleeponset|lightsout", x = qwindowImp$qwindow_names[[qi]])) > 0,
+                      length(grep(pattern = "outbed|wakeup|lightsoff", x = qwindowImp$qwindow_names[[qi]])) > 0)
+        # filter only SPT and time in bed reports
+        qwindow_temp = qwindow_temp[grep(pattern = "bed|wakeup|sleeponset|lights", x = qwindowImp$qwindow_names[[qi]])]
+        qwindow_temp = sort(qwindow_temp)
+        # Ignore if somehow only teh default window is available
+        isDefaultWindow = length(qwindow_temp) == 2 && qwindow_temp[1] == 0 && qwindow_temp[2] == 24
+        if (length(qwindow_temp) > 1 && isDefaultWindow == FALSE && all(available)) {
+          # convert to continuous scale to ease finding start and end
+          below18 = which(qwindow_temp < 18)
+          if (length(below18) > 0) {
+            qwindow_temp = qwindow_temp[below18] + 24
+          }
+          start = min(qwindow_temp)
+          end = max(qwindow_temp)
+          if (length(below18) > 0) {
+            start = ifelse(start >= 24, yes = start - 24, no = start)
+            end = ifelse(end >= 24, yes = end - 24, no = end)
+          }
+          if (start > end) {
+            r1B$filterWindow[which(r1B$date == qwindowImp$date[qi] & 
+                                     r1B$hr >= start |
+                                     r1B$hr < end)] = 1
+          } else {
+            r1B$filterWindow[which(r1B$date == qwindowImp$date[qi] & 
+                                     r1B$hr >= start &
+                                     r1B$hr < end)] = 1
+          }
+        } else {
+          # If diary has one date missing use default window
+          r1B$filter_method[which(r1B$hr >= 0 &
+                                    r1B$hr < 6)] = 4
+          # use midnight - 6am as fall back option
+          r1B$filterWindow[which(r1B$hr >= 0 &
+                                   r1B$hr < 6)] = 1
+        }
+      }
+    }
+    # Identify short nonwear during the night
+    short_nonwear_night = which(r1B$values == 1 &
+                                  r1B$lengths < nonwearFiltermaxHours &
+                                  r1B$filterWindow == 1)
+
+    if (length(short_nonwear_night) > 0) {
+      nonwearHoursFiltered = (length(which(diff(short_nonwear_night) == 1)) * ws2) / 3600
+      nonwearEventsFiltered = length(which(diff(short_nonwear_night) != 1)) + 1
+      r1[short_nonwear_night] = 0
+    }
+  }
+  invisible(list(r1 = r1, nonwearHoursFiltered = nonwearHoursFiltered,
+                 nonwearEventsFiltered = nonwearEventsFiltered))
+}
diff --git a/R/g.analyse.R b/R/g.analyse.R
@@ -46,14 +46,14 @@ g.analyse =  function(I, C, M, IMP, params_247 = c(), params_phyact = c(),
   metalong = M$metalong
   metashort = IMP$metashort
   rout = IMP$rout
+  nonwearHoursFiltered = IMP$nonwearHoursFiltered
+  nonwearEventsFiltered = IMP$nonwearEventsFiltered
   wdaycode = M$wday
   wdayname = M$wdayname
   if (length(params_phyact[["mvpadur"]]) > 0) params_phyact[["mvpadur"]] = sort(params_phyact[["mvpadur"]])
   LC2 = IMP$LC2
   LC = IMP$LC
   dcomplscore = IMP$dcomplscore
-  r1 = as.numeric(as.matrix(rout[,1]))
-  r2 = as.numeric(as.matrix(rout[,2]))
   r4 = as.numeric(as.matrix(rout[,4]))
   r5 = as.numeric(as.matrix(rout[,5]))
   ws3 = windowsizes[1]
@@ -282,7 +282,9 @@ g.analyse =  function(I, C, M, IMP, params_247 = c(), params_phyact = c(),
                                 meas_dur_dys =  LD/1440,
                                 dcomplscore = dcomplscore,
                                 meas_dur_def_proto_day = LMp / 1440,
-                                wear_dur_def_proto_day = LWp / 1440)
+                                wear_dur_def_proto_day = LWp / 1440,
+                                nonwearHoursFiltered = nonwearHoursFiltered,
+                                nonwearEventsFiltered = nonwearEventsFiltered)
   file_summary = data.frame(wdayname = wdayname,
                             deviceSerialNumber = deviceSerialNumber,
                             sensor.location = sensor.location,

diff --git a/R/g.analyse.perday.R b/R/g.analyse.perday.R
@@ -408,7 +408,7 @@ g.analyse.perday = function(ndays, firstmidnighti, time, nfeatures,
                 # Then also extract count metric
 
                 varnum_event = as.numeric(as.matrix(vari[,ExtFunColsi]))
-                if (NRV != length(averageday[, ExtFunColsi])) {
+                if (NRV < length(averageday[, ExtFunColsi])) {
                   if (di == 1) {
                     varnum_event = c(averageday[1:abs(deltaLength), ExtFunColsi], varnum_event)
                   } else {

diff --git a/R/g.analyse.perfile.R b/R/g.analyse.perfile.R
@@ -57,6 +57,12 @@ g.analyse.perfile = function(I, C, metrics_nav,
   s_names[vi:(vi + 4)] = c("clipping_score", "meas_dur_dys", "complete_24hcycle", 
                            "meas_dur_def_proto_day", "wear_dur_def_proto_day")
   vi = vi + 5
+  if (!is.null(params_cleaning[["nonwearFiltermaxHours"]])) {
+    filesummary[vi] = dataqual_summary$nonwearHoursFiltered
+    filesummary[vi + 1] = dataqual_summary$nonwearEventsFiltered
+    s_names[vi:(vi + 1)] = c("nonwear_hours_filtered", "nonwear_events_filtered")
+    vi = vi + 2
+  }
   # calibration error after auto-calibration
   if (length(C$cal.error.end) == 0)   C$cal.error.end = c(" ")
   filesummary[vi] = C$cal.error.end

diff --git a/R/g.conv.actlog.R b/R/g.conv.actlog.R
@@ -4,6 +4,7 @@ g.conv.actlog = function(qwindow, qwindow_dateformat="%d-%m-%Y", epochSize = 5)
   # local functions:
   time2numeric = function(x) {
     x = unlist(x)
+    x = x[grep(pattern = "impute", x = names(x), invert = TRUE)]
     c2t = function(x2) {
       tmp = as.numeric(unlist(strsplit(as.character(x2),":")))
       if (length(tmp) == 2) hourinday = tmp[1] + (tmp[2]/60)

diff --git a/R/g.impute.R b/R/g.impute.R
@@ -1,6 +1,6 @@
 g.impute = function(M, I, params_cleaning = c(), desiredtz = "",
                     dayborder = 0, TimeSegments2Zero = c(), acc.metric = "ENMO", 
-                    ID, ...) {
+                    ID, qwindowImp = c(), ...) {
 
   #get input variables
   input = list(...)
@@ -51,14 +51,18 @@ g.impute = function(M, I, params_cleaning = c(), desiredtz = "",
 
   #========================================
   # Extracting non-wear and clipping and make decision on which additional time needs to be considered non-wear
-  out = g.weardec(M, wearthreshold, ws2, nonWearEdgeCorrection = params_cleaning[["nonWearEdgeCorrection"]])
+  out = g.weardec(M, wearthreshold, ws2,
+                  params_cleaning = params_cleaning,
+                  desiredtz = desiredtz,
+                  qwindowImp = qwindowImp)
   r1 = out$r1 #non-wear
   r2 = out$r2 #clipping
   r3 = out$r3 #additional non-wear
   r4 = matrix(0,length(r3),1) #protocol based decisions on data removal
   LC = out$LC
   LC2 = out$LC2
-
+  nonwearHoursFiltered = out$nonwearHoursFiltered
+  nonwearEventsFiltered = out$nonwearEventsFiltered
   #========================================================
   # Check whether TimeSegments2Zero exist, because this means that the
   # user wants to ignore specific time windows. This feature is used
@@ -412,5 +416,6 @@ g.impute = function(M, I, params_cleaning = c(), desiredtz = "",
   invisible(list(metashort = metashort, rout = rout, r5long = r5long, dcomplscore = dcomplscore,
                  averageday = averageday, windowsizes = windowsizes, data_masking_strategy = params_cleaning[["data_masking_strategy"]],
                  LC = LC, LC2 = LC2, hrs.del.start = params_cleaning[["hrs.del.start"]], hrs.del.end = params_cleaning[["hrs.del.end"]],
-                 maxdur = params_cleaning[["maxdur"]]))
+                 maxdur = params_cleaning[["maxdur"]], nonwearHoursFiltered = nonwearHoursFiltered,
+                 nonwearEventsFiltered = nonwearEventsFiltered))
 }
diff --git a/R/g.part2.R b/R/g.part2.R
@@ -18,9 +18,16 @@ g.part2 = function(datadir = c(), metadatadir = c(), f0 = c(), f1 = c(),
   params_output = params$params_output
   params_general = params$params_general
   #-----------------------------
+  use_qwindow_as_diary = TRUE # If there is a diary specified via qwindow use it as qwindow
   if (is.numeric(params_247[["qwindow"]])) {
     params_247[["qwindow"]] = params_247[["qwindow"]][order(params_247[["qwindow"]])]
   } else if (is.character(params_247[["qwindow"]])) {
+    if (length(grep(pattern = "onlyfilter|filteronly", x = params_247[["qwindow"]])) > 0) {
+      # Do not use diary specified for qwindow if it has the word
+      # "onlyfilter" or "filteronly", but use it for filterning nighttime nonwear
+      # note that this filtering is only use if parameter nonwearFiltermaxHours is specified.
+      use_qwindow_as_diary = FALSE 
+    }
     params_247[["qwindow"]] = g.conv.actlog(params_247[["qwindow"]],
                                             params_247[["qwindow_dateformat"]],
                                             epochSize = params_general[["windowsizes"]][1])
@@ -83,7 +90,7 @@ g.part2 = function(datadir = c(), metadatadir = c(), f0 = c(), f1 = c(),
                         myfun=c(), params_cleaning = c(), params_247 = c(),
                         params_phyact = c(), params_output = c(), params_general = c(),
                         path, ms2.out, foldername, fullfilenames, folderstructure, referencefnames,
-                        daySUMMARY, pdffilenumb, pdfpagecount, csvfolder, cnt78, verbose) {
+                        daySUMMARY, pdffilenumb, pdfpagecount, csvfolder, cnt78, verbose, use_qwindow_as_diary) {
     Nappended = I_list = tail_expansion_log =  NULL
     if (length(ffdone) > 0) {
       if (length(which(ffdone == as.character(unlist(strsplit(fnames[i], "eta_"))[2]))) > 0) {
@@ -150,13 +157,27 @@ g.part2 = function(datadir = c(), metadatadir = c(), f0 = c(), f1 = c(),
             }
           }
         }
+        qwindowImp = params_247[["qwindow"]]
+        if (use_qwindow_as_diary == FALSE) {
+          # reset qwindow to default, because it is only used
+          # for filtering short nighttime nonwear 
+          params_247[["qwindow"]] = c(0, 24)
+        }
+        if (inherits(qwindowImp, "data.frame")) {
+          qwindowImp = qwindowImp[which(qwindowImp$ID == ID),]
+          if (nrow(qwindowImp) == 0) {
+            qwindowImp = c(0, 6) # If participant not present in diary
+          }
+        } else {
+          qwindowImp = NULL
+        }
         IMP = g.impute(M, I,
                        params_cleaning = params_cleaning,
                        dayborder = params_general[["dayborder"]],
                        desiredtz = params_general[["desiredtz"]],
                        TimeSegments2Zero = TimeSegments2Zero,
                        acc.metric = params_general[["acc.metric"]],
-                       ID = ID)
+                       ID = ID, qwindowImp = qwindowImp)
 
         if (params_cleaning[["do.imp"]] == FALSE) { #for those interested in sensisitivity analysis
           IMP$metashort = M$metashort
@@ -315,7 +336,8 @@ g.part2 = function(datadir = c(), metadatadir = c(), f0 = c(), f1 = c(),
                                                   params_phyact, params_output, params_general,
                                                   path, ms2.out, foldername, fullfilenames,
                                                   folderstructure, referencefnames,
-                                                  daySUMMARY, pdffilenumb, pdfpagecount, csvfolder, cnt78, verbose)
+                                                  daySUMMARY, pdffilenumb, pdfpagecount, 
+                                                  csvfolder, cnt78, verbose, use_qwindow_as_diary)
 
                                      })
                                      return(tryCatchResult)
@@ -335,7 +357,8 @@ g.part2 = function(datadir = c(), metadatadir = c(), f0 = c(), f1 = c(),
                  params_phyact, params_output, params_general,
                  path, ms2.out, foldername, fullfilenames,
                  folderstructure, referencefnames,
-                 daySUMMARY, pdffilenumb, pdfpagecount, csvfolder, cnt78, verbose)
+                 daySUMMARY, pdffilenumb, pdfpagecount,
+                 csvfolder, cnt78, verbose, use_qwindow_as_diary)
     }
   }
 }
diff --git a/R/g.part5.R b/R/g.part5.R
@@ -78,6 +78,21 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(),
   } else {
     sleeplog = logs_diaries = c()
   }
+  # Extract activity diary if applicable
+  if (is.character(params_247[["qwindow"]])) {
+    if (length(grep(pattern = "onlyfilter|filteronly", x = params_247[["qwindow"]])) == 0) {
+      epochSize_tmp = ifelse(params_general[["part5_agg2_60seconds"]], yes = 60, no = params_general[["windowsizes"]][1])
+      params_247[["qwindow"]] = g.conv.actlog(params_247[["qwindow"]],
+                                              params_247[["qwindow_dateformat"]],
+                                              epochSize = epochSize_tmp)
+      # This will be an object with numeric qwindow values for all individuals and days
+    } else {
+      # ignore the diary specified by qwindow because user only want to use
+      # it for filtering night time nonwear in part 2, but not as a way to
+      # do day segment analysis.
+      params_247[["qwindow"]] = c(0, 24)
+    }
+  }
   #------------------------------------------------
   # specify parameters
   ffdone = fnames.ms5 #ffdone is now a list of files that have already been processed by g.part5
@@ -453,12 +468,6 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(),
                     indjump = 1
                     qqq_backup = c()
                     add_one_day_to_next_date = FALSE
-                    if (is.character(params_247[["qwindow"]])) {
-                      params_247[["qwindow"]] = g.conv.actlog(params_247[["qwindow"]],
-                                                              params_247[["qwindow_dateformat"]],
-                                                              epochSize = ws3new)
-                      # This will be an object with numeric qwindow values for all individuals and days
-                    }
                     lastDay = ifelse(Nwindows > 0, yes = FALSE, no = TRUE) # skip while loop if there are no days to analyses
                     wi = 1
                     while (lastDay == FALSE) { #loop through windows