Merge pull request #1219 from wadpac/work839_955_1028_1129

Revise and expand circadian rhythm and fragmentation analyses based on recent publications

DESCRIPTION

@@ -21,10 +21,16 @@ Authors@R: c(person("Vincent T","van Hees",role=c("aut","cre"),
                     comment = c(ORCID = "0000-0002-6033-960X")),
              person(given = "Taren",family = "Sanders", role = "ctb"),
              person("Chenxuan","Zhao",role="ctb"),
+             person("Ian","Meneghel Danilevicz", role="ctb",
+                    email="ian.meneghel-danilevicz @ inserm.fr",
+                    comment = c(ORCID = "0000-0003-4541-0524")),
+             person("Victor","Barreto Mesquita", role="ctb",
+                    email="[email protected]"),
              person("Gaia","Segantin",role="ctb"),
              person("Medical Research Council UK",  role = c("cph", "fnd")),
              person("Accelting",  role = c("cph", "fnd")),
-             person("French National Research Agency",  role = c("cph", "fnd")))
+             person("French National Research Agency",  role = c("cph", "fnd"))
+             )
 Maintainer: Vincent T van Hees <[email protected]>
 Description: A tool to process and analyse data collected with wearable raw acceleration sensors as described in Migueles and colleagues (JMPB 2019), and van Hees and colleagues (JApplPhysiol 2014; PLoSONE 2015). The package has been developed and tested for binary data from 'GENEActiv' <https://activinsights.com/>, binary (.gt3x) and .csv-export data from  'Actigraph' <https://theactigraph.com> devices, and binary (.cwa) and .csv-export data from 'Axivity' <https://axivity.com>. These devices are currently widely used in research on human daily physical activity. Further, the package can handle accelerometer data file from any other sensor brand providing that the data is stored in csv format. Also the package allows for external function embedding.
 URL: https://github.com/wadpac/GGIR/, https://groups.google.com/forum/#!forum/RpackageGGIR, https://wadpac.github.io/GGIR/

NAMESPACE

@@ -34,14 +34,15 @@ export(g.analyse, g.calibrate,
        CalcSleepRegularityIndex, load_params,
        check_params, extract_params,
        g.imputeTimegaps, g.part5.classifyNaps,
-       tidyup_df, cosinorAnalyses,
-       detect_nonwear_clipping, applyCosinorAnalyses,
+       tidyup_df, cosinor_IS_IV_Analyses,
+       detect_nonwear_clipping, apply_cosinor_IS_IV_Analyses,
        ShellDoc2Vignette, parametersVignette,
        correctOlderMilestoneData,
        convertEpochData, appendRecords, extractID,
        g.part5_analyseSegment, g.part5_initialise_ts,
        g.part5.analyseRest, part6AlignIndividuals,
-       part6PairwiseAggregation, g.part6, g.report.part6, check_log, g.report.part5_dictionary)
+       part6PairwiseAggregation, g.part6, g.report.part6, 
+       check_log, g.report.part5_dictionary, DFA, ABI, SSP)
 
 
 importFrom("grDevices", "colors", "dev.off", "pdf",
@@ -55,8 +56,8 @@ importFrom("utils", "read.csv", "sessionInfo", "write.csv",
            "memory.limit", "memory.size", "tail", "unzip", "lsf.str")
 importFrom("stats", "aggregate.data.frame", "weighted.mean",
            "rnorm","median","aggregate","C", "lm.wfit",
-           "quantile", "sd","coef", "lm", "embed", "na.pass",
-            "residuals", "fitted", "reshape", "cor")
+           "quantile", "sd","coef", "lm", "embed", "na.pass", "na.omit",
+            "residuals", "fitted", "reshape", "cor", "arima", "lm.fit")
 import(data.table)
 importFrom("methods", "is")
 importFrom("utils", "available.packages", "packageVersion", "help", "read.table")

NEWS.md

@@ -1,3 +1,30 @@
+# CHANGES IN GGIR VERSION 3.1-??
+
+- Part 6:
+
+  - Add DFA functionality #839 (credits: Ian Meneghel Danilevicz and Victor Barreto Mesquita)
+
+  - Add fragmentation metrics, same function as in part 5 but now applied at recording level #839
+
+  - Circadian rhythm analysis now ignore invalid nights controlled with new parameter "includecrit.part6".
+
+  - Circadian rhythm analysis now provides consistent results for both regardless of whether time series are stored with or without valid windows with parameter save_ms5raw_without_invalid.
+
+- Part 2 + 6: Revised and simplified IV and IS calculation which now ignores invalid timestamps and also comes with phi statistic (credits: Ian Meneghel Danilevicz). In part 2 we used to have 2 calculation, which is now replaced by just one and applied to all valid data points in the recordings. In part 6 this is repeated by for the time window as specified with parameter part6Window. Further, IVIS now uses argument threshold.lig as threshold to distinguish inactivity from active.
+
+- Part 2: Increase sensitivity clipping detection from 80% to 30% of a window. If the long epoch has more than 30% of its values close to the dynamic range of the sensor then it will be labelled as clipping.
+
+- Part 4: Now also has copy of part 5 variable ACC_spt_mg.
+
+- Part 5:
+
+  - Moving LXMX from part 5 to part 6 as these are circadian rhythm analysis.
+
+  - Omitted fragmentation metrics for spt window as not meaningful given that part 5 only focusses on valid daytime behaviours.
+  - includenightcrit.part5 added to allow for controlling inclusion of windows in part 5 based on their amount of valid data during spt window. Default remains the same.
+
+  - Fix incorrect usage of part 5 inclusion criteria testing, which used fraction as percentage.
+
 # CHANGES IN GGIR VERSION 3.1-5
 
 - Part 5: 

R/ABI.R

@@ -0,0 +1,33 @@
+# NOTE ABOUT DOCUMENTATION:
+# GGIR does not use Roxygen. The documentation below is not used. 
+# All function documentation can be found in the man/*.Rd files.
+# Please edit documentation there.
+#
+#- @title Activity balance index (ABI)
+#- 
+#- @description This function estimates the Activity balance index (ABI), which is a transformation of the self-similarity parameter (SSP), also known as scaling exponent or alpha.  
+#- @param x the estimated self-similarity parameter (SSP)
+#-
+#- @details ABI = exp(-abs(SSP-1)/exp(-2))
+#-
+#- @return The estimated Activity balance index (ABI) is a real number between zero and one.
+#-
+#- @author Ian Meneghel Danilevicz  
+#- 
+#- @references C.-K. Peng, S.V. Buldyrev, S. Havlin, M. Simons, H.E. Stanley, A.L. Goldberger Phys. Rev. E, 49 (1994), p. 1685
+#- Mesquita, Victor & Filho, Florencio & Rodrigues, Paulo. (2020). Detection of crossover points in detrended fluctuation analysis: An application to EEG signals of patients with epilepsy. Bioinformatics. 10.1093/bioinformatics/btaa955. 
+#-
+#- @examples 
+#- # Estimate Activity balance index of a very known time series available on R base: the sunspot.year.
+#- 
+#- ssp = SSP(sunspot.year, scale = 1.2)
+#- abi = ABI(ssp)
+#-
+ABI = function(x){
+  if (is.na(x)) {
+    y = NA
+  } else {
+    y = exp(-abs(x - 1) * exp(2))  
+  }
+  return(y)  
+}

R/DFA.R

@@ -0,0 +1,86 @@
+# NOTE ABOUT DOCUMENTATION:
+# GGIR does not use Roxygen. The documentation below is not used. 
+# All function documentation can be found in the man/*.Rd files.
+# Please edit documentation there.
+#
+#- @title Detrended Fluctuation Analysis
+#- 
+#- @description ...
+#- @usage DFA(data, scale = 2^(1/8), box_size = 4, m = 1)
+#- @param data Univariate time series (must be a vector or data frame)
+#- @param scale Specifies the ratio between successive box sizes (by default scale = 2^(1/8))
+#- @param box_size Vector of box sizes (must be used in conjunction with scale = "F")
+#- @param m An integer of the polynomial order for the detrending (by default m=1)
+#- 
+#- @details The DFA fluctuation can be computed in a geometric scale or for different choices of boxes sizes.
+#- 
+#- @return Estimated alpha is a real number between zero and two.
+#- 
+#- @note It is not possible estimating alpha for multiple time series at once.
+#- 
+#- @author Ian Meneghel Danilevicz  and Victor Mesquita
+#- 
+#- @references C.-K. Peng, S.V. Buldyrev, S. Havlin, M. Simons, H.E. Stanley, A.L. Goldberger Phys. Rev. E, 49 (1994), p. 1685
+#- Mesquita, Victor & Filho, Florencio & Rodrigues, Paulo. (2020). Detection of crossover points in detrended fluctuation analysis: An application to EEG signals of patients with epilepsy. Bioinformatics. 10.1093/bioinformatics/btaa955.
+#- 
+#- @examples
+#- # Estimate self-similarity of a very known time series available on R base: the sunspot.year.
+#- # Then the spend time with each method is compared.
+#- 
+
+DFA = function(data, scale = 2^(1/8), box_size = 4, m = 1){
+  if (inherits(x = data, "data.frame")) {
+    data = data[, 1]
+  }
+  N = length(data)
+  if (scale != "F") {
+    box_size <- NULL
+    n = 1
+    n_aux <- 0
+    box_size[1] <- 4
+    for (n in 1:N) {
+      while (n_aux < round(N/4)) {
+        n_aux <- box_size[1]
+        n = n + 1
+        box_size[n] <- ceiling(scale * box_size[n - 1])
+        n_aux <- box_size[n] + 4
+      }
+    }
+  }
+  ninbox2 <- NULL
+  for (j in 1:length(box_size)) {
+    ninbox2[j] <- N %/% box_size[j]
+  }
+  aux_seq = seq_along(ninbox2)  
+  aux_length = aux_seq[length(aux_seq)]  
+  y_k = cumsum(data) - mean(data)
+  rm(data, n_aux, scale)
+  aux_mat = matrix(nrow = aux_length, ncol = 2)
+  for (j in seq_along(ninbox2)) {
+    W = box_size[j] * ninbox2[j]
+    aux_j = numeric(W)
+    fit = y_k
+    for (i in seq_len(W)) {
+      if (i == 1) {
+        aux_j[1] = box_size[j]
+        tmp1 = i:aux_j[i]
+      } else if (i >= 2) {
+        aux_j[i] = aux_j[i - 1] + box_size[j]
+        tmp1 = (aux_j[i - 1] + 1):(aux_j[i])
+      }
+      if (m == 1) {
+        fit[tmp1] = lm.fit(x = cbind(1, tmp1), y = y_k[tmp1])$fitted.values
+      } else {
+        fit[tmp1] = lm(y_k[tmp1] ~ poly(tmp1, m, raw = TRUE))$fitted.values
+      }
+      if (i >= ninbox2[j]) {
+        aux_j[i] <- 0
+      }
+    }
+    aux_mat[j,] = c(round(box_size[j], digits = 0),
+                    sqrt((1 / N) * sum((y_k[1:W] - fit[1:W]) ^ 2)))
+  }
+  colnames(aux_mat) <- c("box", "DFA")
+  aux_list = aux_mat
+  return(aux_list)
+}

R/GGIR.R

@@ -320,7 +320,7 @@ GGIR = function(mode = 1:5, datadir = c(), outputdir = c(),
     }
     g.part6(datadir = datadir, metadatadir = metadatadir, f0 = f0, f1 = f1,
             params_general = params_general, params_phyact = params_phyact,
-            params_247 = params_247,
+            params_247 = params_247, params_cleaning = params_cleaning,
             verbose = verbose)
   }
 

R/SSP.R

@@ -0,0 +1,45 @@
+# NOTE ABOUT DOCUMENTATION:
+# GGIR does not use Roxygen. The documentation below is not used. 
+# All function documentation can be found in the man/*.Rd files.
+# Please edit documentation there.
+#
+#- @title Estimated self-similarity parameter
+#- 
+#- @description This function estimates the self-similarity parameter (SSP), also known as scaling exponent or alpha.  
+#- @usage alpha_hat(data,scale = 2^(1/8),box_size = 4,m=1) 
+#- @param data Univariate time series (must be a vector or data frame)
+#- @param scale Specifies the ratio between successive box sizes (by default scale = 2^(1/8))
+#- @param box_size Vector of box sizes (must be used in conjunction with scale = "F")
+#- @param m An integer of the polynomial order for the detrending (by default m=1)
+#-
+#- @details The DFA fluctuation can be computed in a geometric scale or for different choices of boxes sizes.
+#-
+#- @return Estimated alpha is a real number between zero and two.
+#- 
+#- @note It is not possible estimating alpha for multiple time series at once. 
+#-
+#- @author Ian Meneghel Danilevicz and Victor Mesquita 
+#- 
+#- @references C.-K. Peng, S.V. Buldyrev, S. Havlin, M. Simons, H.E. Stanley, A.L. Goldberger Phys. Rev. E, 49 (1994), p. 1685
+#- Mesquita, Victor & Filho, Florencio & Rodrigues, Paulo. (2020). Detection of crossover points in detrended fluctuation analysis: An application to EEG signals of patients with epilepsy. Bioinformatics. 10.1093/bioinformatics/btaa955. 
+#-
+#- @examples 
+#- # Estimate self-similarity of a very known time series available on R base: the sunspot.year.
+#- # Then the spend time with each method is compared.
+#- 
+#- SSP(sunspot.year, scale = 2)
+#- SSP(sunspot.year, scale = 1.2)
+
+SSP = function(data, scale = 2^(1/8), box_size = 4, m = 1){
+  if (inherits(x = data, "data.frame")) {
+    data = data[, 1]
+  }
+  if (length(data) <= box_size || any(is.na(data))) {
+    alpha_hat = NA
+  } else {
+    dfa_hat = DFA(data, scale = scale, box_size = box_size, m = m)    
+    est_ols = lm(log(dfa_hat[,2]) ~ log(dfa_hat[,1]))
+    alpha_hat = est_ols$coefficients[[2]]    
+  }
+  return(alpha_hat)
+} 

R/applyCosinorAnalyses.R → R/apply_cosinor_IS_IV_Analyses.R

@@ -1,10 +1,10 @@
-applyCosinorAnalyses = function(ts, qcheck, midnightsi, epochsizes) {
+apply_cosinor_IS_IV_Analyses = function(ts, qcheck, midnightsi, epochsizes, threshold = NULL) {
   # qcheck - vector of length ts to indicate invalid values
   ws2 = epochsizes[2]
   ws3 = epochsizes[1]
   # Re-derive Xi but this time include entire time series
   # Here, we ignore duplicated values (when clock moves backward due to DST)
-  handleDST = !duplicated(ts)
+  handleDST = !duplicated(ts$time)
   qcheck = qcheck[handleDST]
   Xi = ts[handleDST, grep(pattern = "time", x = colnames(ts), invert = TRUE)]
   Nlong_epochs_day =  (1440 * 60) / ws2 # this is 96 by default
@@ -59,11 +59,8 @@ applyCosinorAnalyses = function(ts, qcheck, midnightsi, epochsizes) {
     } else {
       epochsize = ws3
     }
-    # log transform of data in millig
-    notna = !is.na(Xi)
-    Xi[notna] = log((Xi[notna]*1000) + 1)
-
-    cosinor_coef = cosinorAnalyses(Xi = Xi, epochsize = epochsize, timeOffsetHours = timeOffsetHours) 
+    cosinor_coef = cosinor_IS_IV_Analyses(Xi = Xi, epochsize = epochsize, 
+                                   timeOffsetHours = timeOffsetHours, threshold = threshold)
     cosinor_coef$timeOffsetHours = timeOffsetHours
   } else {
     cosinor_coef = c()

R/check_params.R

@@ -84,7 +84,7 @@ check_params = function(params_sleep = c(), params_metrics = c(),
                        "IVIS.activity.metric", "IVIS_acc_threshold",
                        "qM5L5", "MX.ig.min.dur", "M5L5res", "winhr", "LUXthresholds", "LUX_cal_constant",
                        "LUX_cal_exponent", "LUX_day_segments", "L5M5window")
-    boolean_params = c("cosinor", "part6CR", "part6HCA")
+    boolean_params = c("cosinor", "part6CR", "part6HCA", "part6DFA")
     character_params = c("qwindow_dateformat", "part6Window")
     check_class("247", params = params_247, parnames = numeric_params, parclass = "numeric")
     check_class("247", params = params_247, parnames = boolean_params, parclass = "boolean")
@@ -104,7 +104,7 @@ check_params = function(params_sleep = c(), params_metrics = c(),
   if (length(params_cleaning) > 0) {
     numeric_params = c("includedaycrit", "ndayswindow", "data_masking_strategy", "maxdur", "hrs.del.start",
                        "hrs.del.end", "includedaycrit.part5", "minimum_MM_length.part5",
-                       "includenightcrit", "max_calendar_days")
+                       "includenightcrit", "max_calendar_days", "includecrit.part6", "includenightcrit.part5")
     boolean_params = c("excludefirstlast.part5", "do.imp", "excludefirstlast",
                        "excludefirst.part4", "excludelast.part4", "nonWearEdgeCorrection")
     character_params = c("data_cleaning_file", "TimeSegments2ZeroFile")
@@ -277,12 +277,22 @@ check_params = function(params_sleep = c(), params_metrics = c(),
                                                      replacement = "/", x = params_cleaning[["data_cleaning_file"]])
     }
     if (params_cleaning[["includedaycrit.part5"]] < 0) {
-      warning("\nNegative value of includedaycrit.part5 is not allowed, please change.")
+      stop("\nNegative value of includedaycrit.part5 is not allowed, please change.")
     } else if (params_cleaning[["includedaycrit.part5"]]  > 24) {
-      warning(paste0("\nIncorrect value of includedaycrit.part5, this should be",
+      stop(paste0("\nIncorrect value of includedaycrit.part5, this should be",
                      " a fraction of the day between zero and one or the number ",
                      "of hours in a day."))
     }
+    if (params_cleaning[["includenightcrit.part5"]] < 0) {
+      stop("\nNegative value of includenightcrit.part5 is not allowed, please change.")
+    } else if (params_cleaning[["includenightcrit.part5"]]  > 24) {
+      stop(paste0("\nIncorrect value of includenightcrit.part5, this should be",
+                     " a fraction of the day between zero and one or the number ",
+                     "of hours in a day."))
+    }
+    if (any(params_cleaning[["includecrit.part6"]] < 0) | any(params_cleaning[["includecrit.part6"]] > 1)) {
+      stop("Values of includecrit.part6 are not in the range [0, 1]. Please fix.")
+    }
   }
   if (length(params_phyact) > 0) {
     if (length(params_phyact[["bout.metric"]]) > 0 |
@@ -345,7 +355,7 @@ check_params = function(params_sleep = c(), params_metrics = c(),
       }
     }
     # params 247 & params output
-    if (params_247[["part6HCA"]] == TRUE) {
+    if (params_247[["part6HCA"]] == TRUE || params_247[["part6CR"]] == TRUE) {
       # Add RData because part 6 will need it
       params_output[["save_ms5raw_format"]] = unique(c(params_output[["save_ms5raw_format"]], "RData"))
       params_output[["save_ms5rawlevels"]] = TRUE