update NotWorn algorithm to work for both raw and count data and hand…

…le varying degrees of nonwear time

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()