R/geneabout_functions.R

In GENEAcore: Pre-Processing of 'GENEActiv' Data

#' Non-wear and Rest Coverage for Epochs/Events
#'
#' @param MPI Measurement Period Information.
#' @param aggregated_data Aggregated epochs or events.
#' @param start_time Time stamp to start the coverage assignment, defaults to first time stamp of aggregated data object passed.
#' @param end_time Time stamp to end the coverage assignment, defaults to end time stamp of aggregated data object passed.
#' @param nonwear_day Data frame of non-wear periods and their durations.
#' @keywords internal
#' @export
nonwear_rest_coverage <- function(MPI,
                                  aggregated_data,
                                  start_time = NULL,
                                  end_time = NULL,
                                  nonwear_day) {
  if (is.null(start_time) || is.null(end_time)) {
    start_time <- aggregated_data$TimeUTC[1]
    end_time <- aggregated_data$TimeUTC[nrow(aggregated_data)] + aggregated_data$Duration[nrow(aggregated_data)]
  }

  transitions <- data.frame(
    time_UTC = c(aggregated_data$TimeUTC, aggregated_data$TimeUTC[nrow(aggregated_data)] + aggregated_data$Duration[nrow(aggregated_data)]),
    index = c(1, cumsum(aggregated_data$Duration) + 1)
  )

  # create non-wear time by bout
  nonwear <- nonwear_day[nonwear_day$nonwear_start < end_time & nonwear_day$nonwear_end > start_time, , drop = FALSE]

  nonwear$nonwear_start <- pmax(nonwear$nonwear_start, start_time)
  nonwear$nonwear_end <- pmin(nonwear$nonwear_end, end_time)
  nonwear$MPI_nonwear_duration <- nonwear$nonwear_end - nonwear$nonwear_start

  nonwear <- nonwear[!apply(is.na(nonwear), 1, all), ]

  if (nrow(nonwear) != 0) {
    one_sequence <- data.frame(
      start_time = nonwear$nonwear_start,
      lengths = nonwear$MPI_nonwear_duration,
      values = 1
    )

    wear_start <- append(
      start_time,
      (
        nonwear$nonwear_start + nonwear$MPI_nonwear_duration
      )
    )
    wear_end <- append(
      nonwear$nonwear_start,
      end_time
    )

    zero_sequence <- data.frame(
      start_time = wear_start,
      lengths = wear_end - wear_start,
      values = 0
    )
    non_wear_sequence <- rbind(zero_sequence, one_sequence)
    non_wear_sequence <- non_wear_sequence[order(non_wear_sequence$start_time), ]
    non_wear_sequence <- subset(non_wear_sequence, select = -start_time)
    non_wear_sequence <- inverse.rle(non_wear_sequence)
  } else {
    non_wear_sequence <- rep(
      0,
      end_time - start_time
    )
  }

  restintervals <- MPI$non_movement$rest_intervals
  restintervals$end_time <- restintervals$start_time + restintervals$duration
  rest_intervals <- restintervals[restintervals$start_time < end_time & restintervals$end_time > start_time, , drop = FALSE]

  rest_intervals$start_time <- pmax(rest_intervals$start_time, start_time)
  rest_intervals$end_time <- pmin(rest_intervals$end_time, end_time)
  rest_intervals$duration <- rest_intervals$end_time - rest_intervals$start_time

  rest_intervals <- rest_intervals[!apply(is.na(rest_intervals), 1, all), ]

  if (nrow(rest_intervals) != 0) {
    one_sequence <- data.frame(
      start_time = rest_intervals$start_time,
      lengths = rest_intervals$duration,
      values = 1
    )
    move_start <- append(
      start_time,
      (
        rest_intervals$start_time + rest_intervals$duration
      )
    )
    move_end <- append(
      rest_intervals$start_time,
      end_time
    )
    zero_sequence <- data.frame(
      start_time = move_start,
      lengths = move_end - move_start,
      values = 0
    )
    rest_sequence <- rbind(zero_sequence, one_sequence)
    rest_sequence <- rest_sequence[order(rest_sequence$start_time), ]
    rest_sequence <- subset(rest_sequence, select = -start_time)
    rest_sequence <- inverse.rle(rest_sequence)
  } else {
    rest_sequence <- rep(
      0,
      end_time - start_time
    )
  }

  still_bouts <- MPI$non_movement$still_bouts
  still_bouts$end_time <- still_bouts$start_time + still_bouts$duration
  still_bouts <- still_bouts[still_bouts$start_time < end_time & still_bouts$end_time > start_time, , drop = FALSE]

  still_bouts$start_time <- pmax(still_bouts$start_time, start_time)
  still_bouts$end_time <- pmin(still_bouts$end_time, end_time)
  still_bouts$duration <- still_bouts$end_time - still_bouts$start_time

  still_bouts <- still_bouts[!apply(is.na(still_bouts), 1, all), ]

  if (nrow(still_bouts) != 0) {
    one_sequence <- data.frame(
      start_time = still_bouts$start_time,
      lengths = still_bouts$duration,
      values = 1
    )
    move_start <- append(
      start_time,
      (
        still_bouts$start_time + still_bouts$duration
      )
    )
    move_end <- append(
      still_bouts$start_time,
      end_time
    )
    zero_sequence <- data.frame(
      start_time = move_start,
      lengths = move_end - move_start,
      values = 0
    )
    still_sequence <- rbind(zero_sequence, one_sequence)
    still_sequence <- still_sequence[order(still_sequence$start_time), ]
    still_sequence <- subset(still_sequence, select = -start_time)
    still_sequence <- inverse.rle(still_sequence)
  } else {
    still_sequence <- rep(
      0,
      end_time - start_time
    )
  }

  # allocate non-wear & rest interval time to aggregated data
  bout_sequence <- data.frame(
    lengths = diff(transitions$index),
    values = seq(1, length(diff(
      transitions$index
    )))
  )

  if ((nrow(transitions) > 1) & (sum(bout_sequence$lengths) <= end_time - start_time)) {
    bout_sequence <- rbind(
      bout_sequence,
      c(
        (end_time - start_time - sum(bout_sequence$lengths)),
        nrow(transitions)
      )
    )

    bout_sequence <- inverse.rle(bout_sequence)
  } else {
    bout_sequence <- rep(
      1,
      end_time - start_time
    )
  }


  bout_sequence <- data.frame(bout_sequence, non_wear_sequence, rest_sequence, still_sequence)

  bout_still <- aggregate(
    data = bout_sequence,
    still_sequence ~ bout_sequence,
    FUN = sum
  )

  bout_non_wear <- aggregate(
    data = bout_sequence,
    non_wear_sequence ~ bout_sequence,
    FUN = sum
  )

  bout_rest <- aggregate(
    data = bout_sequence,
    rest_sequence ~ bout_sequence,
    FUN = sum
  )

  aggregated_data$NonWearSeconds <- bout_non_wear$non_wear_sequence
  aggregated_data$RestSeconds <- bout_rest$rest_sequence
  aggregated_data$StillSeconds <- bout_still$still_sequence

  # for partial epochs
  aggregated_data$NonWearSeconds <- pmin(aggregated_data$NonWearSeconds, aggregated_data$Duration)
  aggregated_data$RestSeconds <- pmin(aggregated_data$RestSeconds, aggregated_data$Duration)
  aggregated_data$StillSeconds <- pmin(aggregated_data$StillSeconds, aggregated_data$Duration)

  return(aggregated_data)
}

#' Determine Daily Recommended Primary Rest Interval
#'
#' @param start_expansion_percent Percentage of the still period duration used to extend its start time for overlap detection.
#' @param end_expansion_percent Percentage of the still period duration used to extend its end time for overlap detection.
#' @param cut_time_24hr Time in 24h to split days up by.
#' @param MPI Measurement Period Information.
#' @keywords internal
#' @export
find_rest_intervals <- function(start_expansion_percent,
                                end_expansion_percent,
                                cut_time_24hr,
                                MPI) {
  still_bouts <- MPI[["non_movement"]][["still_bouts"]]
  non_wear <- MPI[["non_movement"]][["non_wear"]]
  non_wear$end_time <- non_wear$start_time + non_wear$duration
  still_bouts$end_time <- still_bouts$start_time + still_bouts$duration
  still_bouts$start_adj <- still_bouts$start_time - floor(still_bouts$duration * start_expansion_percent / 100)
  still_bouts$end_adj <- still_bouts$end_time + floor(still_bouts$duration * end_expansion_percent / 100)

  result <- data.frame()
  cut_times <- get_cut_times(cut_time_24hr, MPI)

  if (nrow(still_bouts) == 0) {
    MPI$non_movement[["rest_intervals"]] <- result
    return(MPI)
  } else if (nrow(still_bouts) == 1) {
    result <- data.frame(
      start_time = still_bouts$start_time,
      duration = still_bouts$duration,
      day_number = 1
    )
    MPI$non_movement[["rest_intervals"]] <- result
    return(MPI)
  } else {
    is_in_non_wear <- function(start_time, end_time, non_wear) {
      any(
        (start_time >= non_wear$start_time & start_time <= non_wear$end_time) |
          (end_time >= non_wear$start_time & end_time <= non_wear$end_time) |
          (start_time <= non_wear$start_time & end_time >= non_wear$end_time)
      )
    }

    still_bouts$group <- 1
    for (row in 2:nrow(still_bouts)) {
      current_start <- still_bouts$start_adj[row]
      current_end <- still_bouts$end_adj[row]
      prev_end <- still_bouts$end_adj[row - 1]

      # Check if current interval touches a non_wear period
      if (is_in_non_wear(still_bouts$start_time[row], still_bouts$end_time[row], non_wear)) {
        # Start a new group due to overlap with non_wear
        still_bouts$group[row] <- still_bouts$group[row - 1] + 1
      } else if (current_start <= prev_end) {
        # Overlaps with previous bout and no non_wear, place in same group
        still_bouts$group[row] <- still_bouts$group[row - 1]
        still_bouts$end_adj[row] <- max(current_end, prev_end)
      } else {
        # No overlap, place in new group
        still_bouts$group[row] <- still_bouts$group[row - 1] + 1
      }
    }

    group_starts <- aggregate(start_time ~ group, data = still_bouts, FUN = min)
    group_ends <- aggregate(end_time ~ group, data = still_bouts, FUN = max)
    rest_bouts <- data.frame(
      start = group_starts$start_time,
      end = group_ends$end_time,
      duration = as.numeric(group_ends$end_time - group_starts$start_time, units = "secs")
    )

    for (day_number in 1:(length(cut_times) - 1)) {
      day_start <- cut_times[day_number]
      day_end <- cut_times[day_number + 1]

      day_rest_bouts <- rest_bouts[rest_bouts$start < day_end & rest_bouts$end > day_start, , drop = FALSE]
      true_day_bouts <- day_rest_bouts[day_rest_bouts$start >= day_start, , drop = FALSE]

      if (nrow(true_day_bouts) > 0) {
        primary <- true_day_bouts[which.max(true_day_bouts$duration), , drop = FALSE]
        result <- rbind(
          result,
          data.frame(
            start_time = primary$start,
            duration = primary$duration,
            day_number = day_number
          )
        )
      }
    }

    MPI$non_movement[["rest_intervals"]] <- result
    return(MPI)
  }
}


#' Bouts Decision Tree
#'
#' @param bouts Aggregated events after non-wear and rest coverage have been assigned.
#' @param SDduration_threshold Threshold for sleep classification. Uses the average x,y,z standard deviation over duration metric.
#' @param AGSA_threshold Threshold for active classification.
#' @param running_threshold Threshold for run classification.
#' @returns Data frame of bouts with classification.
#' @keywords internal
#' @export
bouts_decision_tree <- function(bouts,
                                SDduration_threshold = c(1.38e-4, 5.7e-5),
                                AGSA_threshold = 0.0625,
                                running_threshold = 0.407) {
  bouts$sd_over_dur <- (bouts$xSD + bouts$ySD + bouts$zSD) / (3 * bouts$Duration)

  # decision tree thresholds
  bouts$non_wear <- ifelse(bouts$NonWearSeconds > bouts$Duration / 2, TRUE, FALSE)
  bouts$rest <- ifelse(bouts$RestSeconds > bouts$Duration / 2, TRUE, FALSE)
  bouts$active <- ifelse(bouts$AGSAMean > AGSA_threshold, TRUE, FALSE)
  bouts$ambulatory <- ifelse((bouts$Duration > 9) & (bouts$StepMean >= 45), TRUE, FALSE)
  bouts$vigourous <- ifelse(bouts$AGSAMean > running_threshold, TRUE, FALSE)
  bouts$fastwalk <- ifelse(bouts$StepMean >= 80, TRUE, FALSE)
  sleep_threshold <- ifelse(bouts$rest, SDduration_threshold[1], SDduration_threshold[2])
  bouts$sleep <- !is.na(bouts$sd_over_dur) & bouts$sd_over_dur < sleep_threshold

  # bout classification
  bouts$classification <- ifelse(bouts$non_wear, "non_wear",
    ifelse(bouts$active,
      ifelse(bouts$ambulatory,
        ifelse(bouts$vigourous, "run",
          ifelse(bouts$fastwalk, "fast_walk",
            "slow_walk"
          )
        ),
        "active"
      ),
      ifelse(bouts$sleep, "sleep", "sedentary")
    )
  )

  # reclassify non-wear bouts
  groups <- cumsum(c(TRUE, bouts$classification[-1] != bouts$classification[-nrow(bouts)]))
  reduced_bouts <- do.call(
    rbind,
    lapply(split(bouts, groups), function(x) {
      data.frame(
        TimeUTC = x$TimeUTC[1],
        Duration = sum(x$Duration),
        Classification = x$classification[1]
      )
    })
  )

  if (length(reduced_bouts$Classification == "non_wear") > 0) {
    nw_idx <- which(reduced_bouts$Classification == "non_wear")
    for (idx in nw_idx) {
      class_prior <- reduced_bouts$Classification[idx - 1]
      class_post <- reduced_bouts$Classification[idx + 1]
      if (isTRUE(class_prior == "sleep") && isTRUE(class_post == "sleep")) {
        nw_start_time <- reduced_bouts$TimeUTC[idx]
        nw_end_time <- reduced_bouts$TimeUTC[idx] + reduced_bouts$Duration[idx]
        bouts$classification[bouts$TimeUTC >= nw_start_time & bouts$TimeUTC < nw_end_time] <- "sleep"
      }
    }
  }

  bouts <- bouts[, !(names(bouts) %in% c("sd_over_dur"))]
  new_names <- c(
    non_wear = "NonWear",
    rest = "Rest",
    active = "Active",
    ambulatory = "Ambulatory",
    vigourous = "Vigorous",
    fastwalk = "FastWalk",
    sleep = "Sleep",
    classification = "Classification"
  )

  names(bouts)[names(bouts) %in% names(new_names)] <-
    new_names[names(bouts)[names(bouts) %in% names(new_names)]]

  return(bouts)
}