First half of work on expanding functionality to process external function generated event-data

NEWS.md

@@ -1,3 +1,7 @@
+# CHANGES IN GGIR VERSION 3.1-?
+
+- Part 2: Code revisions in preparation for expansion of functionality to better facilitate external function produced event data. #653
+
 # CHANGES IN GGIR VERSION 3.1-6
 
 - Part 6:

R/MXLX.R

@@ -0,0 +1,71 @@
+MXLX = function(Y = NULL, X = 5, epochSize = 1, tseg = c(0, 24), resolutionMin = 10) {
+  # Y 24 hours of data
+  # tseg is vector of length 2 with the start and end time as hour in the day
+  # X is the window size in hours to look for
+
+  do.MXLX = mean(Y, na.rm = TRUE) > 0
+  Nwindows = (tseg[2] - X) - tseg[1]
+  if (length(do.MXLX) == 0 ||
+      is.na(do.MXLX) == TRUE ||
+      Nwindows < 1 ||
+      length(Y) < length(X) * (3600 / epochSize)) {
+    do.MXLX = FALSE
+  }
+
+  nStepsPerHour = 60 / resolutionMin # number of steps per hour
+  nEpochsPerStep = (resolutionMin * 60) / epochSize # number of epochs per step
+  if (do.MXLX == TRUE) { # only do the analysis if Y has values other than zero
+    Y = Y[((((tseg[1] - tseg[1]) * 3600)/epochSize) + 1):(((tseg[2] - tseg[1]) * 3600)/epochSize)]
+    Nwindows = Nwindows * nStepsPerHour
+    rollingMean = matrix(NA, Nwindows, 1)
+    for (hri in 1:Nwindows) { #e.g.9am-9pm
+      # start and end in terms of steps
+      indexStart = hri - 1 #e.g. 9am
+      indexEnd = indexStart + (X * nStepsPerHour) #e.g. 9am + 5 hrs
+      # start and end in terms of Y index
+      indexStart = (indexStart * nEpochsPerStep) + 1
+      indexEnd = indexEnd * nEpochsPerStep
+      einclude = indexStart:indexEnd
+      if (indexEnd <= length(Y)) {
+        rollingMean[hri,1] = mean(Y[einclude])
+      }
+    }
+    valid = which(is.na(rollingMean) == F)
+    LXvalue = min(rollingMean[valid], na.rm = T)
+    MXvalue = max(rollingMean[valid], na.rm = T)
+    LXhr = ((which(rollingMean == LXvalue & is.na(rollingMean) == F) - 1) / nStepsPerHour) + tseg[1]
+    MXhr = ((which(rollingMean == MXvalue & is.na(rollingMean) == F) - 1) / nStepsPerHour) + tseg[1]
+
+    #-------------------------------------
+    if (length(LXvalue) > 1) { LXvalue = sort(LXvalue)[ceiling(length(LXvalue)/2)] }
+    if (length(LXhr) > 1) { LXhr = sort(LXhr)[ceiling(length(LXhr)/2)] }
+    if (length(MXvalue) > 1) { MXvalue = sort(MXvalue)[ceiling(length(MXvalue)/2)] }
+    if (length(MXhr) > 1) { MXhr = sort(MXhr)[ceiling(length(MXhr)/2)] }
+
+    #Note it is + 1 because first epoch is one and not zero:
+    LXstart = ((LXhr - tseg[1]) * (3600 / epochSize)) + 1
+    LXend = LXstart + (X * (3600 / epochSize)) - 1
+    MXstart = ((MXhr - tseg[1]) * (3600 / epochSize)) + 1 
+    MXend = MXstart + (X * (3600 / epochSize)) - 1
+
+    MXLX = data.frame(LXhr = LXhr[1],
+                      LXvalue = LXvalue,
+                      LXindex0 = LXstart,
+                      LXindex1 = LXend,
+                      MXhr = MXhr[1],
+                      MXvalue = MXvalue,
+                      MXindex0 = MXstart,
+                      MXindex1 = MXend, stringsAsFactors = TRUE)
+  } else {
+    MXLX = data.frame(LXhr = NA, LXvalue = NA, 
+                      LXindex0 = NA,
+                      LXindex1 = NA,
+                      MXhr = NA, MXvalue = NA,
+                      MXindex0 = NA,
+                      MXindex1 = NA, stringsAsFactors = TRUE)
+  }
+  MXLXnames = c(paste0("L", X, "hr"), paste0("L", X), paste0("start_L", X),  paste0("end_L", X),
+                paste0("M", X, "hr"), paste0("M", X),  paste0("start_M", X),  paste0("end_M", X))
+  names(MXLX) = MXLXnames
+  return(MXLX)
+}

R/aggregateEvent.R

@@ -0,0 +1,167 @@
+aggregateEvent = function(metric_name, epochsize,
+                          daysummary,  ds_names, fi, di,
+                          vari, segmentInfo, myfun = NULL, params_247) {
+  anwi_nameindices = segmentInfo$anwi_nameindices
+  anwi_index = segmentInfo$anwi_index
+  anwi_t0 = segmentInfo$anwi_t0
+  anwi_t1 = segmentInfo$anwi_t1
+  if ("ilevels" %in% names(myfun) == FALSE) myfun$ilevels = c(0, 80)
+  if ("clevels" %in% names(myfun) == FALSE) myfun$clevels = c(0, 30)
+  if ("qlevels" %in% names(myfun) == FALSE) myfun$qlevels = c(0.25, 0.5, 0.75)
+  if (length(myfun$ilevels) == 0) myfun$ilevels = 0
+  if (length(myfun$clevels) == 0) myfun$clevels = 0
+  if (length(myfun$qlevels) == 0) myfun$qlevels = 0.5
+  acc.thresholds = myfun$ilevels
+  cadence.thresholds = myfun$clevels
+  qlevels = myfun$qlevels
+  fi_start = fi
+  cadence = vari[, metric_name] * (60/epochsize) # assumption now that cadence is also relevant if the event detected is not steps
+  # aggregate per window total
+  varnameevent = paste0("ExtFunEvent_tot_", metric_name, anwi_nameindices[anwi_index])
+  daysummary[di, fi] = sum(vari[, metric_name])
+  ds_names[fi] = varnameevent; fi = fi + 1
+
+  # cadence
+  varnameevent = paste0("ExtFunEvent_mn_cad", anwi_nameindices[anwi_index])
+  daysummary[di, fi] = mean(cadence, na.rm = TRUE)
+  ds_names[fi] = varnameevent; fi = fi + 1
+
+  # cadence percentiles
+  cadence_quantiles = quantile(cadence, probs = qlevels)
+  varnameevent = paste0("ExtFunEvent_cad_p", qlevels * 100, anwi_nameindices[anwi_index])
+  daysummary[di, fi:(fi + length(qlevels) - 1)] = as.numeric(cadence_quantiles)
+  ds_names[fi:(fi + length(qlevels) - 1)] = varnameevent; fi = fi + length(qlevels)
+
+  #========================================
+  # per acceleration level
+  cn_vari = colnames(vari)
+  acc.metrics = cn_vari[cn_vari %in% c("timestamp", "anglex", "angley", "anglez", metric_name) == FALSE]
+
+  for (ami in 1:length(acc.metrics)) {
+    vari[, acc.metrics[ami]] = as.numeric(vari[, acc.metrics[ami]])
+    for (ti in 1:length(acc.thresholds)) {
+      # step_count per acceleration level
+      if (ti < length(acc.thresholds)) {
+        acc_level_name = paste0(acc.thresholds[ti], "-", acc.thresholds[ti + 1], "mg",  "_",  acc.metrics[ami])
+        whereAccLevel = which(vari[, acc.metrics[ami]] >= (acc.thresholds[ti]/1000) &
+                                vari[, acc.metrics[ami]] < (acc.thresholds[ti + 1]/1000))
+      } else {
+        acc_level_name = paste0("atleast", acc.thresholds[ti], "mg",  "_",  acc.metrics[ami])
+        whereAccLevel = which(vari[,acc.metrics[ami]] >= (acc.thresholds[ti]/1000))
+      }
+      varnameevent = paste0("ExtFunEvent_tot_", metric_name, "_acc", acc_level_name, anwi_nameindices[anwi_index])
+      if (length(whereAccLevel) > 0)  {
+        daysummary[di, fi] = sum(vari[whereAccLevel, metric_name], na.rm = TRUE)
+      } else {
+        daysummary[di, fi] = 0
+      }
+      ds_names[fi] = varnameevent; fi = fi + 1
+
+      # cadence per acceleration level
+      varnameevent = paste0("ExtFunEvent_mn_cad_acc", 
+                            acc_level_name, anwi_nameindices[anwi_index])
+      if (length(whereAccLevel) > 0)  {
+        daysummary[di, fi] = mean(cadence[whereAccLevel], na.rm = TRUE)
+      } else {
+        daysummary[di, fi] = 0
+      }
+      ds_names[fi] = varnameevent; fi = fi + 1
+    }
+  }
+
+  #========================================
+  # per cadence level
+  for (ti in 1:length(cadence.thresholds)) {
+    # define cadence level
+    if (ti < length(cadence.thresholds)) {
+      cadence_level_name = paste0(cadence.thresholds[ti], "-", cadence.thresholds[ti + 1], "spm")
+      whereCadenceLevel = which(cadence >= cadence.thresholds[ti] &
+                                  cadence < cadence.thresholds[ti + 1])
+    } else {
+      cadence_level_name = paste0("atleast", cadence.thresholds[ti], "spm")
+      whereCadenceLevel = which(cadence >= cadence.thresholds[ti])
+    }
+    # cadence per cadence level
+    varnameevent = paste0("ExtFunEvent_mn_cad_cad", cadence_level_name, anwi_nameindices[anwi_index])
+    if (length(whereCadenceLevel) > 0)  {
+      daysummary[di, fi] = mean(cadence[whereCadenceLevel], na.rm = TRUE)
+    } else {
+      daysummary[di, fi] = 0
+    }
+    ds_names[fi] = varnameevent; fi = fi + 1
+
+    # step count per cadence level
+    varnameevent = paste0("ExtFunEvent_tot_", metric_name, "_cad", 
+                          cadence_level_name, anwi_nameindices[anwi_index])
+    if (length(whereCadenceLevel) > 0)  {
+      daysummary[di, fi] = sum(vari[whereCadenceLevel, metric_name], na.rm = TRUE)
+    } else {
+      daysummary[di, fi] = 0
+    }
+    ds_names[fi] = varnameevent; fi = fi + 1
+
+    # time per cadence level
+    varnameevent = paste0("ExtFunEvent_dur_cad", cadence_level_name, anwi_nameindices[anwi_index])
+    daysummary[di, fi] = length(whereCadenceLevel) / (60/epochsize)
+    ds_names[fi] = varnameevent; fi = fi + 1
+
+    for (ami in 1:length(acc.metrics)) {
+      # acceleration per cadence level
+      varnameevent = paste0("ExtFunEvent_mn_",  acc.metrics[ami],"_cad",
+                            cadence_level_name, anwi_nameindices[anwi_index])
+      if (length(whereCadenceLevel) > 0)  {
+        daysummary[di,fi] = mean(vari[whereCadenceLevel, acc.metrics[ami]], na.rm = TRUE) * 1000
+      } else {
+        daysummary[di, fi] = 0
+      }
+      ds_names[fi] = varnameevent; fi = fi + 1
+    }
+  }
+  if (length(params_247[["winhr"]]) > 0) { 
+    # Event based MXLX
+    tseg = c((anwi_t0[anwi_index] - 1) / (3600 / epochsize), anwi_t1[anwi_index] / (3600 / epochsize))
+    for (winhr_value in params_247[["winhr"]]) {
+      MXLXout = MXLX(Y = cadence, X = winhr_value, epochSize = epochsize,
+                     tseg = tseg,
+                     resolutionMin = params_247[["M5L5res"]])
+      # Describe for MX and LX:
+      LXMXwindow_name = anwi_nameindices[anwi_index]
+      for (WX in c("L", "M")) {
+        if (!all(is.na(MXLXout[, grep("M", names(MXLXout))]))) {
+          i0 = MXLXout[, grep(paste0("start_", WX, winhr_value), names(MXLXout))]
+          i1 = MXLXout[, grep(paste0("end_", WX, winhr_value), names(MXLXout))]
+          WXi = i0:i1
+          # # Code used to help check that MXLX makes sense with real study data
+          # x11()
+          # plot(cadence, type = "l")
+          # lines(WXi, cadence[WXi], type = "l", col = "red")
+
+          # Timing
+          daysummary[di, fi] = MXLXout[, which(names(MXLXout) == paste0(WX, winhr_value, "hr"))]
+          ds_names[fi] = paste0("ExtFunEvent_", WX, winhr_value,  "hr_cad", LXMXwindow_name); fi = fi + 1
+          # cadence in the form of mean, percentiles.
+          daysummary[di, fi] = mean(cadence[WXi])
+          ds_names[fi] = paste0("ExtFunEvent_", WX, winhr_value, "_cad_meancad", LXMXwindow_name); fi = fi + 1
+          if (length(params_247[["qM5L5"]]) > 0) {
+            for (qle in params_247[["qM5L5"]]) {
+              daysummary[di, fi] = quantile(x = cadence[WXi], probs = qle)
+              ds_names[fi] = paste0("ExtFunEvent_", WX, winhr_value, "_cad_q", qle*100, "cad", LXMXwindow_name); fi = fi + 1
+            }
+          }
+          for (ami in 1:length(acc.metrics)) {
+            # Acceleration in the form of mean, percentiles.
+            daysummary[di, fi] = mean(vari[WXi, acc.metrics[ami]]) * 1000
+            ds_names[fi] = paste0("ExtFunEvent_", WX, winhr_value, "_cad_mean_",  acc.metrics[ami], "_mg", LXMXwindow_name); fi = fi + 1
+            if (length(params_247[["qM5L5"]]) > 0) {
+              for (qle in params_247[["qM5L5"]]) {
+                daysummary[di, fi] = quantile(x = vari[WXi, acc.metrics[ami]], probs = qle)
+                ds_names[fi] = paste0("ExtFunEvent_", WX, winhr_value, "_cad_q", qle*100, "acc", LXMXwindow_name); fi = fi + 1
+              }
+            }
+          }
+        }
+      }
+    }
+  }
+  invisible(list(ds_names = ds_names, daysummary = daysummary, fi = fi))
+}

R/check_myfun.R

@@ -3,81 +3,87 @@ check_myfun = function(myfun, windowsizes) { # Function to check myfun object
   # check that myfun is a list:
   if (is.list(myfun) == F) {
     status = 1
-    stop("Error in check_myfun.R: Object myfun is not a list.") 
+    stop("Object myfun is not a list.", call. = FALSE)
   }
   # check that there are no foreign object:
-  foreignElements = which(names(myfun) %in% c("FUN", "parameters", "expected_sample_rate", "expected_unit", 
-                                              "colnames", "minlength", "outputres", "outputtype", "aggfunction",
-                                              "timestamp","reporttype") == FALSE)
+  foreignElements = which(names(myfun) %in% c("FUN", "parameters", "expected_sample_rate",
+                                              "expected_unit", "colnames",
+                                              "minlength", "outputres",
+                                              "outputtype", "aggfunction",
+                                              "timestamp","reporttype",
+                                              "clevels", "ilevels", "qlevels",
+                                              "ebout.dur", "ebout.th.cad", "ebout.th.acc",
+                                              "ebout.criter", "ebout.condition", "name") == FALSE)
   if (length(foreignElements) != 0) {
     status = 1
-    stop("Error in check_myfun.R: Object myfun has unexpected elements.")
+    stop("Object myfun has unexpected elements.", call. = FALSE)
   }
   # check that essential objects are included:
-  expectedElements = c("FUN", "parameters", "expected_sample_rate", "expected_unit", 
-                       "colnames", "minlength", "outputres") 
+  expectedElements = c("FUN", "parameters", "expected_sample_rate", "expected_unit",
+                       "colnames", "minlength", "outputres")
   missingElements = which(expectedElements %in% names(myfun) == FALSE)
   if (length(missingElements) != 0) {
     status = 1
-    stop(paste0("Error in check_myfun.R: Object myfun misses the following elements: ",
-                paste(expectedElements[missingElements],collapse=", "),"."), call. = FALSE)
+    stop(paste0("Object myfun misses the following elements: ",
+                paste(expectedElements[missingElements],collapse = ", "), "."), call. = FALSE)
   }
   # check that FUN is a function:
   if (is.function(myfun$FUN) ==  FALSE) {
     status = 1
-     stop("Error in check_myfun.R: Element FUN in myfun is not a function.", call. = FALSE)
+     stop("Element FUN in myfun is not a function.", call. = FALSE)
   }
   # check that expected_sample_rate is numeric
   if (is.numeric(myfun$expected_sample_rate) == F) {
     status = 1
-    stop("Error in check_myfun.R: Element expected_sample_rate in myfun is not numeric.", call. = FALSE)
+    stop("Element expected_sample_rate in myfun is not numeric.", call. = FALSE)
   }
   # check that unit is specified:
-  if (length(which(myfun$expected_unit %in%  c("mg","g", "ms2") == TRUE)) != 1) {
+  if (length(which(myfun$expected_unit %in%  c("mg", "g", "ms2") == TRUE)) != 1) {
     status = 1
-    stop("Error in check_myfun.R: Object myfun lacks a clear specification of the expected_unit.", call. = FALSE)
+    stop("Object myfun lacks a clear specification of the expected_unit.", call. = FALSE)
   }
   # check that colnames has at least one character value:
   if (length(myfun$colnames) == 0) {
     status = 1
-    stop("Error in check_myfun.R: Element colnames in myfun does not have a value.", call. = FALSE)
+    stop("Element colnames in myfun does not have a value.", call. = FALSE)
   }
   # Check that colnames is a character:
   if (is.character(myfun$colnames) == F) {
     status = 1
-    stop("Error in check_myfun.R: Element colnames in myfun does not hold a character value.", call. = FALSE)
+    stop("Element colnames in myfun does not hold a character value.", call. = FALSE)
   }
   # check that minlegnth has one value:
   if (length(myfun$minlength) != 1) {
     status = 1
-    stop("Error in check_myfun.R: Element minlength in myfun does not have one value.", call. = FALSE)
+    stop("Element minlength in myfun does not have one value.", call. = FALSE)
   }
   # Check that minlength is a number
   if (is.numeric(myfun$minlength) == F) {
     status = 1
-    stop("Error in check_myfun.R: Element minlength in myfun is not numeric.", call. = FALSE)
+    stop("Element minlength in myfun is not numeric.", call. = FALSE)
   }
   # check that outputres has one value:
   if (length(myfun$outputres) != 1) {
     status = 1
-    stop("Error in check_myfun.R: Element outputres in myfun does not have one value.", call. = FALSE)
+    stop("Element outputres in myfun does not have one value.", call. = FALSE)
   }
   # Check that outputres is a number:
   if (is.numeric(myfun$outputres) == F) {
     status = 1
-    stop("Error in check_myfun.R: Element outputres in myfun is not numeric.", call. = FALSE)
+    stop("Element outputres in myfun is not numeric.", call. = FALSE)
   }
   # check that it is a round number can add up to 900 seconds:
   if (myfun$outputres != round(myfun$outputres)) {
     status = 1
-    stop("Error in check_myfun.R: Element outputres in myfun should be a round number.", call. = FALSE)
+    stop("Element outputres in myfun should be a round number.", call. = FALSE)
   }
   # check that outputres is either a multitude of the epochs size or vica versa:
   if (myfun$outputres/windowsizes[1] != round(myfun$outputres/windowsizes[1]) &
       windowsizes[1]/myfun$outputres != round(windowsizes[1]/myfun$outputres)) {
     status = 1
-    stop(paste0("Error in check_myfun.R: Element outputres and the epochsize used",
-                " in GGIR (first element of windowsizes) are not a multitude of each other."), call. = FALSE)
+    stop(paste0("Element outputres and the epochsize used in GGIR",
+                " (first element of windowsizes) are not a multitude",
+                " of each other.", call. = FALSE))
   }
   if ("outputtype" %in% names(myfun)) { # If outputtype is available:
     if (is.character(myfun$outputtype) == F) {
@@ -92,6 +98,7 @@ check_myfun = function(myfun, windowsizes) { # Function to check myfun object
       stop("Error in check_myfun.R: Element aggfunction is not a function object.", call. = FALSE)
     }
   }
+
   # if ("timestamp" %in% names(myfun)) { # If timestamp is available:
   #   if (is.logical(myfun$timestamp) == F) {
   #     status = 1