R/g.sib.det.R

In GGIR: Raw Accelerometer Data Analysis

g.sib.det = function(M, IMP, I, twd = c(-12, 12),
                     acc.metric = "ENMO", desiredtz = "",
                     myfun=c(), sensor.location = "wrist", params_sleep = c(), zc.scale = 1, ...) {
  
  #get input variables
  input = list(...)
  if (length(input) > 0 || length(params_sleep) == 0) {
    # Extract and check parameters if user provides more arguments than just the parameter arguments,
    # or if params_[...] aren't specified (so need to be filled with defaults).
    # So, inside GGIR this will not be used, but it is used when g.sleep is used on its own
    # as if it was still the old g.sleep function
    params = extract_params(params_sleep = params_sleep,
                            input = input,
                            params2check = "sleep") # load default parameters
    params_sleep = params$params_sleep
  }
  #==============================================================
  perc = 10; spt_threshold = 15; sptblocksize = 30; spt_max_gap = 60 # default configurations (keep hardcoded for now
  # Abbreviaton SPTE = Sleep Period Time Estimate, although in case of HorAngle it is the Time in Bed estimate
  # but then we would have to come up with yet another term to represent the main sleep and/or time in bed window of the day
  # 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"))
    }
    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
      }
    }
    return(calc_SPTE_end)
  }
  #==================================================
  # get variables
  nD = nrow(IMP$metashort)
  mon = I$monn
  ws3 = M$windowsizes[1] #short epoch size
  ws2 = M$windowsizes[2] # medium length epoch size used for non-wear detection
  n_ws2_perday = (1440*60) / ws2
  n_ws3_perday = (1440*60) / ws3
  #--------------------
  # get indicator of non-wear periods
  rout = IMP$rout
  r5 = as.numeric(as.matrix(rout[,5]))
  r5long = matrix(0, length(r5), (ws2/ws3))
  r5long = replace(r5long, 1:length(r5long), r5)
  r5long = t(r5long)
  dim(r5long) = c((length(r5) * (ws2/ws3)),1)
  if (nD < length(r5long)) {
    invalid = r5long[1:nD]
  } else {
    invalid = c(r5long, rep(0, (nD - length(r5long))))
  }
  rm(r5,rout)
  # Deriving file characteristics from 15 min summary files
  ND = nD/n_ws3_perday #number of days
  if (ND > 0.2) {
    #========================================================================
    # timestamps
    time = format(IMP$metashort[,1])
    # # angle
    # if (length(which(colnames(IMP$metashort) == "anglez")) == 0 ) {
    #   cat("metric anglez was not extracted, please make sure that anglez  is extracted")
    # }
    fix_NA_invector = function(x){
      if (length(which(is.na(x) == TRUE)) > 0) {
        x[which(is.na(x) == T)] = 0
      }
      return(x)
    }
    anglez = as.numeric(as.matrix(IMP$metashort[,which(colnames(IMP$metashort) == "anglez")]))
    anglez = fix_NA_invector(anglez)
    
    
    anglex = angley = c()
    do.HASPT.hip = FALSE
    if (sensor.location == "hip" &
        "anglex" %in% colnames(IMP$metashort) &
        "angley" %in% colnames(IMP$metashort) &
        "anglez" %in% colnames(IMP$metashort)) {
      do.HASPT.hip = TRUE
      if (length(colnames(IMP$metashort) == "anglex") > 0) {
        anglex =  IMP$metashort[, which(colnames(IMP$metashort) == "anglex")]
        anglex = fix_NA_invector(anglex)
      }
      if (length(colnames(IMP$metashort) == "angley") > 0) {
        angley = IMP$metashort[,which(colnames(IMP$metashort) == "angley")]
        angley = fix_NA_invector(angley)
      }
    }
    if (acc.metric %in% colnames(IMP$metashort) == FALSE) {
      warning("Argument acc.metric is set to ",acc.metric," but not found in GGIR part 1 output data")
    }
    ACC = as.numeric(as.matrix(IMP$metashort[,which(colnames(IMP$metashort) == acc.metric)]))
    night = rep(0, length(ACC))
    if (params_sleep[["HASIB.algo"]] == "Sadeh1994" | 
        params_sleep[["HASIB.algo"]] == "Galland2012" |
        params_sleep[["HASIB.algo"]] == "ColeKripke1992") { # extract zeroCrossingCount
      zeroCrossingCount =  IMP$metashort[,which(colnames(IMP$metashort) == paste0("ZC", params_sleep[["Sadeh_axis"]]))]
      zeroCrossingCount = fix_NA_invector(zeroCrossingCount)
      zeroCrossingCount = zeroCrossingCount * zc.scale
      # always do zeroCrossingCount but optionally also add BrondCounts to output for comparison
      BrondCount_colname = paste0("BrondCount_", tolower(params_sleep[["Sadeh_axis"]]))
      if (BrondCount_colname %in% colnames(IMP$metashort)) {
        BrondCount =  IMP$metashort[, BrondCount_colname]
        BrondCount = fix_NA_invector(BrondCount)
      } else {
        BrondCount = c()
      }
      # optionally add NeishabouriCounts for comparison
      NeishabouriCount_colname = paste0("NeishabouriCount_", tolower(params_sleep[["Sadeh_axis"]]))
      if (NeishabouriCount_colname %in% colnames(IMP$metashort)) {
        NeishabouriCount =  IMP$metashort[, NeishabouriCount_colname]
        NeishabouriCount = fix_NA_invector(NeishabouriCount)
      } else {
        NeishabouriCount = c()
      }
    } else {
      zeroCrossingCount = c()
      BrondCount = c()
      NeishabouriCount = c()
    }
    #==================================================================
    # 'sleep' detection if sleep is not provided by external function.
    # Note that inside the code we call it sustained inactivity bouts
    # to emphasize that we know that this is not actually neurological sleep
    getSleepFromExternalFunction = FALSE
    if (length(myfun) != 0) {
      if ("wake_sleep" %in% myfun$colnames & myfun$outputtype == "character") {
        getSleepFromExternalFunction = TRUE
      }
    }
    if (getSleepFromExternalFunction == FALSE) {
      sleep = HASIB(HASIB.algo = params_sleep[["HASIB.algo"]], 
                    timethreshold = params_sleep[["timethreshold"]],
                    anglethreshold = params_sleep[["anglethreshold"]], 
                    time = time, anglez = anglez, ws3 = ws3,
                    zeroCrossingCount = zeroCrossingCount,
                    BrondCount = BrondCount,
                    NeishabouriCount = NeishabouriCount, activity = ACC)
    } else { # getSleepFromExternalFunction == TRUE
      # Code now uses the sleep estimates from the external function
      # So, the assumption is that the external function provides a 
      # character "Sleep" when it detects sleep
      sleep = matrix(0, nD, 1)
      sleep[which(M$metashort$wake_sleep == "Sleep")] = 1 
    }
    #-------------------------------------------------------------------
    # detect midnights
    detemout = g.detecmidnight(time, desiredtz, dayborder = 0) # ND, # for sleep dayborder is always 0, it is the summary of sleep that will be dayborder specific
    midnights = detemout$midnights
    midnightsi = detemout$midnightsi
    countmidn = length(midnightsi)
    tib.threshold = SPTE_end = SPTE_start = L5list = rep(NA,countmidn)
    if (countmidn != 0) {
      if (countmidn == 1) {
        tooshort = 1
        lastmidnight = midnights[length(midnights)]
        lastmidnighti = midnightsi[length(midnights)]
        firstmidnight = time[1]
        firstmidnighti = 1
      } else {
        cut = which(as.numeric(midnightsi) == 0)
        if (length(cut) > 0) {
          midnights = midnights[-cut]
          midnightsi = midnightsi[-cut]
        }
        lastmidnight = midnights[length(midnights)]
        lastmidnighti = midnightsi[length(midnights)]
        firstmidnight = midnights[1]
        firstmidnighti = midnightsi[1]
      }
      # if recording started before 4am, then also derive first awakening
      first4am = grep("04:00:00", time[1:pmin(nD, (n_ws3_perday + 1))])[1]
      # only do this if there is a 4am in the recording
      if (!is.na(first4am)) {
        if (first4am < firstmidnighti) { # this means recording started after midnight and before 4am
          midn_start = 0
        } else {
          midn_start = 1
        }
      } else {
        # since there's no 4am in the recording, then even if there is a midnight, skip this midnight. 
        # (basically treat this midnight like any other midnight without a preceding 4am)
        midn_start = countmidn
      }
      sptei = 0
      for (j in midn_start:(countmidn)) { #Looping over the midnight
        if (j == 0) {
          qqq1 = 1 # preceding noon (not available in recording)
          qqq2 = midnightsi[1] + (twd[1] * (3600 / ws3))# first noon in recording
        } else {
          qqq1 = midnightsi[j] + (twd[1] * (3600 / ws3)) #preceding noon
          qqq2 = midnightsi[j] + (twd[2] * (3600 / ws3)) #next noon
        }
        sptei = sptei + 1
        if (qqq2 > length(time))  qqq2 = length(time)
        if (qqq1 < 1)             qqq1 = 1
        night[qqq1:qqq2] = sptei
        detection.failed = FALSE
        #------------------------------------------------------------------
        # calculate L5 because this is used as back-up approach
        tmpACC = ACC[qqq1:qqq2]
        windowRL = round((3600/ws3)*5)
        if ((windowRL / 2) == round(windowRL / 2)) windowRL = windowRL + 1
        if (length(tmpACC) < windowRL) {  0 # added 4/4/2-17
          L5 = 0
        } else {
          ZRM = zoo::rollmean(x = c(tmpACC), k = windowRL, fill = "extend", align = "center") #
          L5 = which(ZRM == min(ZRM))[1]
          if (sd(ZRM) == 0) {
            L5 = 0
          } else {
            L5 = (L5  / (3600 / ws3)) + 12
          }
          if (length(L5) == 0) L5 = 0 #if there is no L5, because full they is zero
        }
        L5list[sptei] = L5
        # Estimate Sleep Period Time window, because this will be used by g.part4 if sleeplog is not available
        tmpANGLE = anglez[qqq1:qqq2]
        tmpTIME = time[qqq1:qqq2]
        daysleep_offset = 0
        if (do.HASPT.hip == TRUE & params_sleep[["HASPT.algo"]] != "NotWorn") {
          params_sleep[["HASPT.algo"]] = "HorAngle"
          if (length(params_sleep[["longitudinal_axis"]]) == 0) { #only estimate long axis if not provided by user
            count_updown = matrix(0,3,2)
            count_updown[1,] = sort(c(length(which(anglex[qqq1:qqq2] < 45)), length(which(anglex[qqq1:qqq2] > 45))))
            count_updown[2,] = sort(c(length(which(angley[qqq1:qqq2] < 45)), length(which(angley[qqq1:qqq2] > 45))))
            count_updown[3,] = sort(c(length(which(anglez[qqq1:qqq2] < 45)), length(which(anglez[qqq1:qqq2] > 45))))
            ratio_updown = count_updown[,1] / count_updown[,2]
            validval = which(abs(ratio_updown) != Inf & is.na(ratio_updown) == FALSE)
            if (length(validval) > 0) {
              ratio_updown[validval] = ratio_updown
              params_sleep[["longitudinal_axis"]] = which.max(ratio_updown)
            } else {
              params_sleep[["longitudinal_axis"]] = 2 # y-axis as fall back option if detection does not work
            }
          }
          if (params_sleep[["longitudinal_axis"]] == 1) {
            tmpANGLE = anglex[qqq1:qqq2]
          } else if (params_sleep[["longitudinal_axis"]] == 2) {
            tmpANGLE = angley[qqq1:qqq2]
          }
        }
        if (length(params_sleep[["def.noc.sleep"]]) == 1) {
          spt_estimate = HASPT(angle = tmpANGLE, ws3 = ws3,
                               constrain2range = params_sleep[["constrain2range"]],
                               perc = perc, spt_threshold = spt_threshold,
                               sptblocksize = sptblocksize, spt_max_gap = spt_max_gap,
                               HASPT.algo = params_sleep[["HASPT.algo"]],
                               invalid = invalid,
                               HASPT.ignore.invalid = params_sleep[["HASPT.ignore.invalid"]],
                               activity = tmpACC)
        } else {
          spt_estimate = list(SPTE_end = NULL, SPTE_start = NULL, tib.threshold = NULL)
        }
        if (length(spt_estimate$SPTE_end) != 0 & length(spt_estimate$SPTE_start) != 0) {
          if (spt_estimate$SPTE_end + qqq1 >= qqq2 - (1 * (3600 / ws3))) {
            # if estimated SPT ends within one hour of noon, re-run with larger window
            # to be able to detect sleep in daysleepers
            daysleep_offset = 6 # hours in which the window of data sent to SPTE is moved fwd from noon
            newqqq1 = qqq1 + (daysleep_offset * (3600 / ws3))
            newqqq2 = qqq2 + (daysleep_offset * (3600 / ws3))
            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)) ) {
              spt_estimate_tmp = HASPT(anglez[newqqq1:newqqq2], ws3 = ws3,
                                       constrain2range = params_sleep[["constrain2range"]],
                                       perc = perc, spt_threshold = spt_threshold, sptblocksize = sptblocksize,
                                       spt_max_gap = spt_max_gap,
                                       HASPT.algo = params_sleep[["HASPT.algo"]], invalid = invalid,
                                       HASPT.ignore.invalid = params_sleep[["HASPT.ignore.invalid"]],
                                       activity = tmpACC[newqqq1:newqqq2])
              if (length(spt_estimate_tmp$SPTE_start) > 0) {
                if (spt_estimate_tmp$SPTE_start + newqqq1 >= newqqq2) {
                  spt_estimate_tmp$SPTE_start = (newqqq2 - newqqq1) - 1
                  spt_estimate = spt_estimate_tmp
                } else {
                  daysleep_offset  = 0
                }
              } else {
                daysleep_offset  = 0
              }
            } else {
              daysleep_offset  = 0
            }
          } 
          if (qqq1 == 1) {  # 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))
            SPTE_end[sptei] = (spt_estimate$SPTE_end / (3600 / ws3)) + startTimeRecord + daysleep_offset
            SPTE_start[sptei] = (spt_estimate$SPTE_start / (3600 / ws3)) + startTimeRecord + daysleep_offset
          } else {
            SPTE_end[sptei] = (spt_estimate$SPTE_end / (3600 / ws3)) + 12 + daysleep_offset
            SPTE_start[sptei] = (spt_estimate$SPTE_start / (3600 / ws3)) + 12 + 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])
          tib.threshold[sptei] = spt_estimate$tib.threshold
        }
      }
      detection.failed = FALSE
    } else {
      cat("No midnights found")
      detection.failed = TRUE
    }
    metatmp = data.frame(time, invalid, night = night, sleep = sleep, stringsAsFactors = T)
  } else {
    metatmp = L5list = SPTE_end = SPTE_start = tib.threshold = c()
    detection.failed = TRUE
  }
  invisible(list(output = metatmp, detection.failed = detection.failed, L5list = L5list,
                 SPTE_end = SPTE_end, SPTE_start = SPTE_start,
                 tib.threshold = tib.threshold, longitudinal_axis = params_sleep[["longitudinal_axis"]]))
}