Enhance NotWorn algorithm to work with both count and raw data

.Rbuildignore

@@ -19,4 +19,5 @@ RELEASE_CYCLE.md
 codecov.yml
 index.md
 ^inst/doc/*.html
-^vignettes/HouseHoldCoanalysis.Rmd
+^vignettes/HouseHoldCoanalysis.Rmd
+^vignettes/Cookbook.Rmd

NEWS.md

@@ -4,6 +4,10 @@
 
 - Report part 5: fix bug that was introduced on 2024-Feb-19 in the calculation of wear percentage #1148
 
+- Part 3 and 4: Revise NotWorn algorithm to work with both count and raw data with varying degrees of nonwear. Further, when HASPT.algo is set to NotWorn then a second guider can optionally be specified in case the accelerometer unexpectedly worn by the participant. #1089
+
+- Visualreport: Improve handling of recordings where the accelerometer was not worn.
+
 - Part 1: Enable timegap imputation for ad-hoc csv data. #1154
 
 - Part 1: Reduce constraints on value for parameter chunksize #1155.
@@ -14,6 +18,8 @@
 
 - Vignette: Migrated many sections from main CRAN vignette to wadpac.github.io/GGIR/
 
+- Part 3: Fix handling of recordings with only 1 midnight that start before 4am, #1160
+
 # CHANGES IN GGIR VERSION 3.1-1
 
 - Part 2: Corrected calculation of LXhr and MXhr which had one hour offset when timing was after midnight, #1117

R/GGIR.R

@@ -335,7 +335,7 @@ GGIR = function(mode = 1:5, datadir = c(), outputdir = c(),
                           "GGIRversion",  "dupArgValues", "verbose", "is_GGIRread_installed", 
                           "is_read.gt3x_installed", "is_ActCR_installed", 
                           "is_actilifecounts_installed", "rawaccfiles", "is_readxl_installed", 
-                          "checkFormat", "getExt") == FALSE)]
+                          "checkFormat", "getExt", "rawaccfiles_formats") == FALSE)]
 
   config.parameters = mget(LS)
   config.matrix = as.data.frame(createConfigFile(config.parameters, GGIRversion))

R/HASIB.R

@@ -223,8 +223,17 @@ HASIB = function(HASIB.algo = "vanHees2015", timethreshold = c(), anglethreshold
     # around to avoid having to completely redesign GGIR just for those studies.
 
     sib_classification = as.data.frame(matrix(0, Nvalues, 1))
-    activityThreshold = as.numeric(quantile(x = activity, probs = 0.1) * 2)
-    activity2 = zoo::rollmax(x = activity, k = 3600 / epochsize, fill = 1)
+    activity2 = zoo::rollmax(x = activity, k = 300 / epochsize, fill = 1)
+    # ignore zeros because in ActiGraph with many zeros it skews the distribution
+    nonzero = which(activity2 != 0)
+    if (length(nonzero) > 0) {
+      activityThreshold = sd(activity2[nonzero], na.rm = TRUE) * 0.05
+      if (activityThreshold < min(activity)) {
+        activityThreshold = quantile(activity2[nonzero], probs = 0.1)
+      }
+    } else {
+      activityThreshold = 0
+    }
     zeroMovement = which(activity2 <= activityThreshold)
     if (length(zeroMovement) > 0) {
       sib_classification[zeroMovement, 1] = 1

R/HASPT.R

@@ -18,7 +18,7 @@ HASPT = function(angle, sptblocksize = 30, spt_max_gap = 60, ws3 = 5,
       # threshold = 10th percentile (constrained to 0.13-0.5 if required)
       medabsdi = function(angle) {
         #50th percentile, do not use mean because that will be outlier sensitive
-        angvar = stats::median(abs(diff(angle))) 
+        angvar = stats::median(abs(diff(angle)))
         return(angvar)
       }
       k1 = 5 * (60/ws3)
@@ -39,7 +39,7 @@ HASPT = function(angle, sptblocksize = 30, spt_max_gap = 60, ws3 = 5,
       # threshold = 45 degrees
       x = abs(angle)
       threshold = 45
-    } else if (HASPT.algo == "NotWorn") {  
+    } else if (HASPT.algo == "NotWorn") {
       # When protocol is to not wear sensor during the night,
       # and data is collected in count units we do not know angle
       # as needed for HorAngle and HDCZA.
@@ -52,15 +52,20 @@ HASPT = function(angle, sptblocksize = 30, spt_max_gap = 60, ws3 = 5,
       # smooth x to 5 minute rolling average to reduce sensitivity to sudden peaks
       ma <- function(x, n = 300 / ws3){stats::filter(x, rep(1 / n, n), sides = 2, circular = TRUE)}
       x = ma(x)
-      activityThreshold = sd(x, na.rm = TRUE) * 0.2
-      # For sensewear external data this will not work as it mostly has values of 1 and up.
-      if (activityThreshold < min(activity)) {
-        activityThreshold = quantile(x, probs = 0.1)  
+      nonzero = which(x != 0)
+      if (length(nonzero) > 0) {
+        activityThreshold = sd(x[nonzero], na.rm = TRUE) * 0.05
+        # For sensewear external data this will not work as it mostly has values of 1 and up.
+        if (activityThreshold < min(activity)) {
+          activityThreshold = quantile(x, probs = 0.1)
+        }
+      } else {
+        activityThreshold = 0
       }
       # this algorithm looked for x <= threshold, now a minimum quantity is added
       # to the threshold to allow for consistent definition of nomov below
       # i.e., x < threshold
-      threshold = activityThreshold + 0.001 
+      threshold = activityThreshold + 0.001
     }
     # Now define nomov periods with the selected strategy for invalid time
     nomov = rep(0,length(x)) # no movement
@@ -129,8 +134,8 @@ HASPT = function(angle, sptblocksize = 30, spt_max_gap = 60, ws3 = 5,
       if (SPTE_start == 0) SPTE_start = 1
       part3_guider = HASPT.algo
       if (is.na(HASPT.ignore.invalid)) {
-        # investigate if invalid time was included in the SPT definition, 
-        # and if so, keep track of that in the guider. This is needed in the 
+        # investigate if invalid time was included in the SPT definition,
+        # and if so, keep track of that in the guider. This is needed in the
         # case that sleeplog is used, to inform part 4 that it should
         # trust the sleeplog times for this specific night.
         spt_long = rep(0, length(invalid))
@@ -141,6 +146,28 @@ HASPT = function(angle, sptblocksize = 30, spt_max_gap = 60, ws3 = 5,
         }
       }
 
+      # # Code to help investigate classifications:
+      # plot(x, col = "black", type = "l")
+      # abline(v = SPTE_start, col = "green", lwd = 2)
+      # abline(v = SPTE_end, col = "red", lwd = 2)
+      # rect(xleft = s1, ybottom = rep(0, length(s1)),
+      #      xright = e1, ytop = rep(0.1, length(s1)),
+      #      col = rgb(0, 0, 255, max = 255, alpha = 50), border = NA)
+      # 
+      # rect(xleft = s5, ybottom = rep(0.1, length(s1)),
+      #      xright = e5, ytop = rep(1, length(s1)),
+      #      col = rgb(255, 0, 0, max = 255, alpha = 20), border = NA)
+      # lines(x, col = "black", type = "l")
+      # abline(h = threshold, col = "purple", lwd = 2)
+      # inva = which(invalid == 1)
+      # if (length(inva) > 0) {
+      #   lines(inva, rep(0.1, length(inva)),
+      #         type = "p", pch = 20, lwd = 4, col = "black")
+      # }
+      # lines(invalid* 0.05, type = "l", col = "red")
+      # # graphics.off()
+      # browser()
+
     } else {
       SPTE_end = c()
       SPTE_start = c()

R/check_params.R

@@ -176,7 +176,7 @@ check_params = function(params_sleep = c(), params_metrics = c(),
 
   if (length(params_sleep) > 0) {
     if (length(params_sleep[["def.noc.sleep"]]) != 2) {
-      if (params_sleep[["HASPT.algo"]] %in% c("HorAngle", "NotWorn") == FALSE) {
+      if (params_sleep[["HASPT.algo"]][1] %in% c("HorAngle", "NotWorn") == FALSE) {
         params_sleep[["HASPT.algo"]] = "HDCZA"
       }
     } else if (length(params_sleep[["def.noc.sleep"]]) == 2) {
@@ -200,13 +200,13 @@ check_params = function(params_sleep = c(), params_metrics = c(),
 
   if (length(params_general) > 0 & length(params_metrics) > 0 & length(params_sleep) > 0) {
     if (params_general[["sensor.location"]] == "hip" &&
-        params_sleep[["HASPT.algo"]] %in% c("notused", "NotWorn") == FALSE) {
+        params_sleep[["HASPT.algo"]][1] %in% c("notused", "NotWorn") == FALSE) {
       if (params_metrics[["do.anglex"]] == FALSE | params_metrics[["do.angley"]] == FALSE | params_metrics[["do.anglez"]] == FALSE) {
         warning(paste0("\nWhen working with hip data all three angle metrics are needed,",
                        "so GGIR now auto-sets arguments do.anglex, do.angley, and do.anglez to TRUE."), call. = FALSE)
         params_metrics[["do.anglex"]] = params_metrics[["do.angley"]] = params_metrics[["do.anglez"]] = TRUE
       }
-      if (params_sleep[["HASPT.algo"]] != "HorAngle") {
+      if (params_sleep[["HASPT.algo"]][1] != "HorAngle") {
         warning("\nChanging HASPT.algo value to HorAngle, because sensor.location is set as hip", call. = FALSE)
         params_sleep[["HASPT.algo"]] = "HorAngle"; params_sleep[["def.noc.sleep"]] = 1
       }
@@ -230,12 +230,12 @@ check_params = function(params_sleep = c(), params_metrics = c(),
       params_sleep[["def.noc.sleep"]] = 1
     }
 
-    if (params_sleep[["HASPT.algo"]] == "HorAngle" & params_sleep[["sleepwindowType"]] != "TimeInBed") {
+    if (params_sleep[["HASPT.algo"]][1] == "HorAngle" & params_sleep[["sleepwindowType"]] != "TimeInBed") {
       warning("\nHASPT.algo is set to HorAngle, therefore auto-updating sleepwindowType to TimeInBed", call. = FALSE)
       params_sleep[["sleepwindowType"]] = "TimeInBed"
     }
 
-    if (length(params_sleep[["loglocation"]]) == 0 & params_sleep[["HASPT.algo"]] != "HorAngle" & params_sleep[["sleepwindowType"]] != "SPT") {
+    if (length(params_sleep[["loglocation"]]) == 0 & params_sleep[["HASPT.algo"]][1] != "HorAngle" & params_sleep[["sleepwindowType"]] != "SPT") {
       warning("\nAuto-updating sleepwindowType to SPT because no sleeplog used and neither HASPT.algo HorAngle used.", call. = FALSE)
       params_sleep[["sleepwindowType"]] = "SPT"
     }
@@ -367,7 +367,7 @@ check_params = function(params_sleep = c(), params_metrics = c(),
   }
 
   if (length(params_general) > 0 & length(params_sleep) > 0) {
-    if (params_sleep[["HASPT.algo"]] == "HorAngle") {
+    if (params_sleep[["HASPT.algo"]][1] == "HorAngle") {
       # Not everywhere in the code we check that when HASPT.algo is HorAngle, sensor.location is hip.
       # However, the underlying assumption is that they are linked. Even if a study would
       # use a waist or chest worn sensor we refer to it as hip as the orientation and need

R/convertEpochData.R

@@ -522,6 +522,9 @@ convertEpochData = function(datadir = c(), metadatadir = c(),
         M$metashort = as.data.frame(cbind(time_shortEp_8601,
                                           D[1:length(time_shortEp_8601), ]))
         colnames(M$metashort) = c("timestamp", colnames(D))
+        for (ic in 2:ncol(M$metashort)) {
+          M$metashort[, ic] <- as.numeric(M$metashort[, ic])
+        }
       }
       if (length(which(is.na(M$metashort$ZCY) == TRUE)) > 0) {
         # impute missing data by zero