add DFA functionality, fixes #839

Co-Authored-By: Victor Barreto Mesquita <54644227+victormesquita40@users.noreply.github.com>
Co-Authored-By: Ian Danilevicz <43249776+iandanilevicz@users.noreply.github.com>

DESCRIPTION

@@ -20,10 +20,16 @@ Authors@R: c(person("Vincent T","van Hees",role=c("aut","cre"),
                     role = "ctb", email = "egpbos@gmail.com",
                     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="victormesquita40@hotmail.com"),
              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("Chenxuan","Zhao",role="ctb"))
+             person("French National Research Agency",  role = c("cph", "fnd"))
+             )
 Maintainer: Vincent T van Hees <v.vanhees@accelting.com>
 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/> and GENEA devices (not for sale), .csv-export data from  'Actigraph' <https://theactigraph.com> devices, and .cwa and .wav-format 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 and has either no header or a two column header. 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

@@ -41,7 +41,8 @@ export(g.analyse, g.calibrate,
        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",

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,90 @@
+# 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){
+  DFA_aux = function(j, box_size, ninbox2, data, y_k, m, N){
+    aux_j = numeric(box_size[j] * ninbox2[j])
+    fit = y_k
+    for (i in seq_len(box_size[j] * ninbox2[j])) {
+      t1 = Sys.time()
+      if (i == 1) {
+        aux_j[1] = box_size[j]
+        mod_i = lm(y_k[i:aux_j[i]] ~ poly(c(i:aux_j[i]), m, raw = TRUE))
+        fit[i:(aux_j[i])] = mod_i$fitted.values
+      }
+      if (i >= 2) {
+        aux_j[i] = aux_j[i - 1] + box_size[j]
+        mod_i = lm(y_k[(aux_j[i - 1] + 1):(aux_j[i])] ~ poly(c((aux_j[i - 1] +1):(aux_j[i])), m, raw = TRUE))
+        fit[(aux_j[i - 1] + 1):(aux_j[i])] = mod_i$fitted.values
+      }
+      if (i >= ninbox2[j]) {
+        aux_j[i] <- 0
+      }
+      t2 = Sys.time()
+    }
+    DFA = sqrt((1 / N) * sum((y_k[1:(box_size[j] * ninbox2[j])] -
+                                fit[1:(box_size[j] * ninbox2[j])]) ^ 2))
+    Results = c(round(box_size[j], digits = 0), DFA)
+    return(Results)
+  }
+
+  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)
+  aux_mat = matrix(nrow = aux_length, ncol = 2)
+  for (j in seq_along(ninbox2)) {
+    aux_mat[j,] = DFA_aux(j, box_size, ninbox2, data, y_k, m, N)
+  }
+  colnames(aux_mat) <- c("box", "DFA")
+  aux_list = aux_mat
+  return(aux_list)
+}

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 (length(data) <= box_size && 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/check_params.R

@@ -79,7 +79,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")

R/g.part6.R

@@ -202,18 +202,25 @@ g.part6 = function(datadir = c(), metadatadir = c(), f0 = c(), f1 = c(),
       s_names[fi] = "N_valid_days"
       fi = fi + 1
 
-      # Cosinor analysis
+
       if (do.cr == TRUE) {
+        # Cosinor analysis
         colnames(ts)[which(colnames(ts) == "timenum")] = "time"
         acc4cos = ts[, c("time", "ACC")]
         acc4cos$ACC  = acc4cos$ACC / 1000 # convert to mg because that is what applyCosinorAnalyses expects
         cosinor_coef = applyCosinorAnalyses(ts = acc4cos,
                                             qcheck = ts$invalidepoch,
                                             midnightsi = nightsi,
                                             epochsizes = rep(epochSize, 2))
+        # DFA
+        if (params_247[["part6DFA"]] == TRUE) { 
+          ssp = SSP(ts$ACC)
+          abi = ABI(ssp)
+        }
         rm(acc4cos)
       } else {
         cosinor_coef = NULL
+        ssp = abi = NULL
       }
       if (length(cosinor_coef) > 0) {
         summary[fi] = cosinor_coef$timeOffsetHours
@@ -264,6 +271,11 @@ g.part6 = function(datadir = c(), metadatadir = c(), f0 = c(), f1 = c(),
                                   "cosinorExt_R2", "cosinorIS", "cosinorIV")
         fi = fi + 20
       }
+      if (params_247[["part6DFA"]] == TRUE) {
+        summary[fi:(fi + 1)] = c(ssp, abi)
+        s_names[fi:(fi + 1)] = c("SSP", "ABI")
+        fi = fi + 2
+      }
       #=============================================
       # Store results in milestone data
       summary = summary[1:(fi - 1),]

R/load_params.R

@@ -75,7 +75,8 @@ load_params = function(topic = c("sleep", "metrics", "rawdata",
                       LUX_cal_constant = c(), LUX_cal_exponent = c(), LUX_day_segments = c(),
                       L5M5window = c(0, 24), cosinor = FALSE,
                       part6CR = FALSE, part6HCA = FALSE,
-                      part6Window = c("start", "end"))
+                      part6Window = c("start", "end"),
+                      part6DFA = FALSE)
   }
   if ("phyact" %in% topic) {
     params_phyact = list(mvpathreshold = 100, boutcriter = 0.8,

man/ABI.Rd

@@ -0,0 +1,33 @@
+\name{ABI}
+\alias{ABI}
+\title{Activity balance index (ABI)}
+\usage{
+  ABI(x)
+}
+\arguments{
+  \item{x}{the estimated self-similarity parameter (SSP)}
+}
+\value{
+  The estimated Activity balance index (ABI) is a real number between zero and one.
+}
+\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.
+}
+\details{
+  ABI = exp(-abs(SSP-1)/exp(-2))
+}
+\examples{
+  # Estimate Activity balance index of a very known time series
+  # available on R base: the sunspot.year.
+  \dontrun{  
+    ssp = SSP(sunspot.year)
+    abi = ABI(ssp)
+  }
+}
+\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.
+}
+\author{
+  Ian Meneghel Danilevicz <ian.meneghel-danilevicz@inserm.fr>
+}

man/DFA.Rd

@@ -0,0 +1,48 @@
+\name{DFA}
+\alias{DFA}
+\title{Detrended Fluctuation Analysis}
+\description{
+  Detrended Fluctuation Analysis (DFA)
+}
+\usage{
+  DFA(data, scale = 2^(1/8), box_size = 4, m = 1)
+}
+\arguments{
+  \item{data}{
+    Univariate time series (must be a vector or data frame)
+  }
+  \item{scale}{
+    Specifies the ratio between successive box sizes (by default scale = 2^(1/8))
+  }
+  \item{box_size}{
+    Vector of box sizes (must be used in conjunction with scale = "F")
+  }
+  \item{m}{
+    An integer of the polynomial order for the detrending (by default m=1)
+  }
+}
+\value{
+  Estimated alpha is a real number between zero and two.
+}
+\details{
+  The DFA fluctuation can be computed in a geometric scale or for different choices of boxes sizes.
+}
+\note{
+  It is not possible estimating alpha for multiple time series at once.
+}
+\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.
+  \dontrun{
+    dfa = DFA(sunspot.year)
+  }
+}
+\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.
+}
+\author{
+  Ian Meneghel Danilevicz <ian.meneghel-danilevicz@inserm.fr>
+  Victor Barreto Mesquita <victormesquita40@hotmail.com>
+}

man/GGIR.Rd

@@ -1348,6 +1348,10 @@ GGIR(mode = 1:5,
         ignores the first and last 5 hours, and c("O2", "W10") goes from the second
         onset till the 10th wakeup time.
       }
+      \item{part6DFA}{
+        Boolean (default = FALSE) to indicate whether to perform Detrended Fluctuation
+        Analysis. Turned off by default because it can be very time consuming.
+      }
     }
   }