Issue1499 speed up DFA calculations

NEWS.md

@@ -2,8 +2,12 @@
 
 - Documentation: Remove suggestion that configfile can be used to store annotations inside the configuration file to keep focus on the main purpose of the configfile, which is to facilitate reproducing analysis. #1473
 
-- Part 6: Enable extraction of participant ID from filename consistent with the other parts #1478
+- Part 6: 
+
+  - Enable extraction of participant ID from filename consistent with the other parts #1478
 
+  - Improve speed calculation of DFA analysis using a vectorized approach #1499
+
 - Part 5: Variable dictionary, now also documents ACC_spt_mg, ACC_spt_mg_median, and ACC_spt_mg_stdev. #1490
 
 # CHANGES IN GGIR VERSION 3.3-5

R/DFA.R

@@ -29,58 +29,61 @@
 #- 
 
 DFA = function(data, scale = 2^(1/8), box_size = 4, m = 1){
+
+  # Helper to replicate the fluctuation analysis using sapply later
+  calc_fluctuation = function(n, y, N_total) {
+    num_boxes = N_total %/% n
+    W = n * num_boxes
+    y_truncated = y[1:W]
+
+    # Reshape y into a matrix: each column is one box
+    y_matrix = matrix(y_truncated, nrow = n)
+
+    # Create the x-vector for a single box: 1:n
+    x = 1:n
+    # Slope (beta) for each box: beta = cov(x,y) / var(x)
+    x_detrend = x - mean(x)
+    sum_x2 = sum(x_detrend^2)
+
+    # betas is a vector of slopes (one for each box)
+    betas = colSums(y_matrix * x_detrend) / sum_x2
+
+    # intercepts is a vector of intercepts (one for each box)
+    y_means = colMeans(y_matrix)
+    intercepts = y_means - (betas * mean(x))
+
+    # Fits for every point
+    # fit_matrix[i, j] = intercepts[j] + betas[j] * x[i]
+    betas_mtx = matrix(betas, nrow = n, ncol = num_boxes, byrow = TRUE)
+    intercepts_mtx = matrix(intercepts, nrow = n, ncol = num_boxes, byrow = TRUE)
+    fit_matrix = intercepts_mtx + betas_mtx * x
+
+    # RMS calculation as in DFA (average fluctuation)
+    return(sqrt(sum((y_matrix - fit_matrix)^2) / N_total))
+  }
+
+  # MAIN SCRIPT ------
   if (inherits(x = data, "data.frame")) {
     data = data[, 1]
   }
+
+  # Box sizes calculation
   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
-      }
+    box_sizes = 4
+    while (tail(box_sizes, 1) + 4 < round(N/4)) {
+      box_sizes <- c(box_sizes, ceiling(scale * tail(box_sizes, 1)))
     }
+  } else {
+    box_sizes = box_size
   }
-  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: integrated acceleration - global mean
   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)
+
+  # Calculate the detrended fluctuation amplitude across all chosen time scales
+  dfa_values = sapply(box_sizes, calc_fluctuation, y = y_k, N_total = N)
+
+  aux_mat = cbind(box = box_sizes, DFA = dfa_values)
+  return(aux_mat)
 }

R/SSP.R

@@ -11,10 +11,11 @@
 #- @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)
+#- @param epochSize An integer with epoch length in seconds
 #-
 #- @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.
+#- @return Estimated alpha, alpha_1, and alpha_2: real number between zero and two.
 #- 
 #- @note It is not possible estimating alpha for multiple time series at once. 
 #-
@@ -30,16 +31,43 @@
 #- SSP(sunspot.year, scale = 2)
 #- SSP(sunspot.year, scale = 1.2)
 
-SSP = function(data, scale = 2^(1/8), box_size = 4, m = 1){
+SSP = function(data, scale = 2^(1/8), box_size = 4, m = 1, epochSize){
   if (inherits(x = data, "data.frame")) {
     data = data[, 1]
   }
+
   if (length(data) <= box_size || any(is.na(data))) {
-    alpha_hat = NA
+    return(list(alpha_overall = NA, alpha_1 = NA, alpha_2 = 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]]    
+    # Calculate fluctuation values for all box sizes
+    dfa_hat = DFA(data, scale = scale, box_size = box_size, m = m)
+
+    box_sizes = dfa_hat[, "box"]
+    dfa_values = dfa_hat[, "DFA"]
+    box_sizes_min = box_sizes * epochSize / 60
+
+    # Overall Alpha
+    model_overall = lm(log(dfa_values) ~ log(box_sizes))
+    alpha_overall = as.numeric(coef(model_overall)[2])
+
+    # Alpha_1 (<= 90 minutes)
+    idx_1 = which(box_sizes_min <= 90)
+    if (length(idx_1) >= 3) {
+      model_1 = lm(log(dfa_values[idx_1]) ~ log(box_sizes[idx_1]))
+      alpha_1 = as.numeric(coef(model_1)[2])
+    } else {
+      alpha_1 = NA
+    }
+
+    # Alpha_2 (Between 120 and 600 minutes)
+    idx_2 = which(box_sizes_min >= 120 & box_sizes_min <= 600)
+    if (length(idx_2) >= 3) {
+      model_2 = lm(log(dfa_values[idx_2]) ~ log(box_sizes[idx_2]))
+      alpha_2 = as.numeric(coef(model_2)[2])
+    } else {
+      alpha_2 = NA
+    }
   }
-  return(alpha_hat)
+
+  return(list(alpha_overall = alpha_overall, alpha_1 = alpha_1, alpha_2 = alpha_2))
 } 

R/g.part6.R

@@ -506,11 +506,16 @@ g.part6 = function(datadir = c(), metadatadir = c(), f0 = c(), f1 = c(),
         #=======================================================================
         # DFA analyses is used independent of do.cr because it is time consuming
         if (params_247[["part6DFA"]] == TRUE) {
-          ssp = SSP(ts$ACC[!is.na(ts$ACC)])
-          abi = ABI(ssp)
-          summary[fi:(fi + 1)] = c(ssp, abi)
-          s_names[fi:(fi + 1)] = c("SSP", "ABI")
-          fi = fi + 2
+          SSP = SSP(ts$ACC[!is.na(ts$ACC)], epochSize = epochSize)
+          ssp = SSP$alpha_overall; ssp_short = SSP$alpha_1; ssp_long = SSP$alpha_2
+          abi = ABI(ssp); abi_short = ABI(ssp_short); abi_long = ABI(ssp_long)
+          ssp_diff = ssp_short - ssp_long
+          summary[fi:(fi + 6)] = c(ssp, abi, ssp_short, abi_short, ssp_long, abi_long, ssp_diff)
+          s_names[fi:(fi + 6)] = c("SSP", "ABI", 
+                                   "SSP_short", "ABI_short", 
+                                   "SSP_long", "ABI_long", 
+                                   "SSP_diff")
+          fi = fi + 7
         }
         #------------------------------------------------------------
         # Sleep Regularity Index

man/DFA.Rd

@@ -22,7 +22,9 @@
   }
 }
 \value{
-  Estimated alpha is a real number between zero and two.
+A matrix with two columns:
+  \item{box}{The box sizes (time scales) used for the analysis.}
+  \item{DFA}{The corresponding detrended fluctuation amplitudes (Root Mean Square).}
 }
 \details{
   The DFA fluctuation can be computed in a geometric scale or for different choices of boxes sizes.
@@ -45,4 +47,5 @@
 \author{
   Ian Meneghel Danilevicz <ian.meneghel-danilevicz@inserm.fr>
   Victor Barreto Mesquita <victormesquita40@hotmail.com>
+  Jairo H Migueles <jairo@jhmigueles.com>
 }

man/SSP.Rd

@@ -2,10 +2,12 @@
 \alias{SSP}
 \title{Estimated self-similarity parameter}
 \description{
-  This function estimates the self-similarity parameter (SSP), also known as scaling exponent or alpha.
+  This function estimates the overall self-similarity parameter (SSP), 
+  also known as the scaling exponent or alpha, as well as the piecewise 
+  short-term (alpha_1) and long-term (alpha_2) scaling exponents.
 }
 \usage{
-  SSP(data,scale = 2^(1/8),box_size = 4,m=1)
+  SSP(data, scale = 2^(1/8), box_size = 4, m = 1, epochSize)
 }
 \arguments{
   \item{data}{
@@ -18,14 +20,26 @@
     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)
+    An integer of the polynomial order for the detrending (by default m = 1)
+  }
+  \item{epochSize}{
+    The epoch size of the data in seconds. Used to convert box sizes to minutes 
+    to apply time boundaries for alpha_1 and alpha_2.
   }
 }
 \value{
-  Estimated alpha is a real number between zero and two.
+  A list containing three real numbers (all of them between 0 and 2):
+    \item{alpha_overall}{The overall estimated scaling exponent across all box sizes.}
+    \item{alpha_1}{The short-term scaling exponent, calculated for time scales smaller than or equal to 90 minutes.}
+    \item{alpha_2}{The long-term scaling exponent, calculated for time scales between 120 and 600 minutes.}
 }
 \details{
-  The DFA fluctuation can be computed in a geometric scale or for different choices of boxes sizes.
+  The DFA fluctuation can be computed in a geometric scale or for different 
+  choices of boxes sizes. In human motor activity, the scaling behavior exhibits 
+  a crossover point around 1.5 to 2 hours. Therefore, alpha_1 captures the 
+  short-term temporal correlations regulated by multiple physiological controls, 
+  while alpha_2 captures long-term fluctuations that are heavily reliant on 
+  the central circadian pacemaker.
 }
 \note{
   It is not possible estimating alpha for multiple time series at once.
@@ -34,7 +48,10 @@
   # 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{
-    ssp = SSP(sunspot.year)
+    ssp_results = SSP(sunspot.year)
+    ssp_results$alpha_overall
+    ssp_results$alpha_1
+    ssp_results$alpha_2
   }
 }
 \references{
@@ -44,4 +61,5 @@
 \author{
   Ian Meneghel Danilevicz <ian.meneghel-danilevicz@inserm.fr>
   Victor Barreto Mesquita <victormesquita40@hotmail.com>
+  Jairo H Migueles <jairo@jhmigueles.com>
 }

tests/testthat/test_dfa.R

@@ -2,26 +2,28 @@ library(GGIR)
 context("DFA, SSP and ABI")
 test_that("DFA produces expected abi and ssp from sunspot.year dataset", {
   skip_on_cran()
-  ssp = SSP(sunspot.year)
-  abi = ABI(ssp)
+  ssp = SSP(sunspot.year, epochSize = 60)
+  abi = ABI(ssp$alpha_overall)
   expect_equal(abi, 0.145756, tolerance = 0.0001)
-  expect_equal(ssp, 0.7393684, tolerance = 0.0001)
+  expect_equal(ssp$alpha_overall, 0.7393684, tolerance = 0.0001)
+  expect_equal(length(ssp), 3)
+  expect_all_true(c("alpha_overall", "alpha_1", "alpha_2") %in% names(ssp))
 })
 
 test_that("DFA skips calculation if there are NA values", {
   skip_on_cran()
   sy_tmp = sunspot.year
   sy_tmp[5] = NA
-  ssp = SSP(sy_tmp)
-  abi = ABI(ssp)
+  ssp = SSP(sy_tmp, epochSize = 60)
+  abi = ABI(ssp$alpha_overall)
   expect_equal(abi, NA)
-  expect_equal(ssp, NA)
+  expect_equal(ssp$alpha_overall, NA)
 })
 
 test_that("DFA can handle too short time series", {
   skip_on_cran()
-  ssp = SSP(sunspot.year[1:2])
-  abi = ABI(ssp)
+  ssp = SSP(sunspot.year[1:2], epochSize = 60)
+  abi = ABI(ssp$alpha_overall)
   expect_equal(abi, NA)
-  expect_equal(ssp, NA)
+  expect_equal(ssp$alpha_overall, NA)
 })

vignettes/GGIRoutput.Rmd

@@ -432,6 +432,11 @@ These variables names are consistent with the variable names [as discussed for G
 | `FRAG_` | Fragmentation variables, as also [discussed for part 5 output](https://wadpac.github.io/GGIR/articles/GGIRoutput.html#fragmentation). Only difference now is that fragmentation variables name ending with `_day` are specific to the waking hours of a day, while variable names ending iwht `_spt` are specific to the SPT window.  |
 | SSP | See [section on SSP in chapter 13](https://wadpac.github.io/GGIR/articles/chapter13_CircadianRhythm.html#self-similarity-paramerter-ssp) |
 | ABI | See [section on ABI in chapter 13](https://wadpac.github.io/GGIR/articles/chapter13_CircadianRhythm.html#activity-balance-index-abi) |
+| SSP_short | See [section on SSP in chapter 13](https://wadpac.github.io/GGIR/articles/chapter13_CircadianRhythm.html#self-similarity-paramerter-ssp) |
+| ABI_short | See [section on ABI in chapter 13](https://wadpac.github.io/GGIR/articles/chapter13_CircadianRhythm.html#activity-balance-index-abi) |
+| SSP_long | See [section on SSP in chapter 13](https://wadpac.github.io/GGIR/articles/chapter13_CircadianRhythm.html#self-similarity-paramerter-ssp) |
+| ABI_long | See [section on ABI in chapter 13](https://wadpac.github.io/GGIR/articles/chapter13_CircadianRhythm.html#activity-balance-index-abi) |
+| SSP_diff | See [section on SSP in chapter 13](https://wadpac.github.io/GGIR/articles/chapter13_CircadianRhythm.html#self-similarity-paramerter-ssp) |
 | SleepRegularityIndex2 | The Sleep Regularity Index inspired by [Phillips et al. 2017](https://www.nature.com/articles/s41598-017-03171-4), but calculated per day-pair to enable user to study patterns across days, and calculated based on classified naps and sleep. See [Chapter 10](https://wadpac.github.io/GGIR/articles/chapter10_SleepAnalysis.html#sleep-regularity-index-sri) for details. |
 | SleepRegularityIndex2_Ndaypairs | Number of day pairs used for calculating SleepRegularityIndex2. |
 

vignettes/chapter13_CircadianRhythm.Rmd

@@ -18,8 +18,8 @@ Circadian rhythm are the physical, mental, and behavioral changes an organism ex
 
 Most of the techniques are applied in both GGIR part 2 and 6. The main differences between these implementations are:
 
-- In GGIR part 2 the circadian rhythm analysis tries to use every valid data point in the recording.
-- In GGIR part 6 the circadian rhythm analysis uses a specific time window as defined by the user, e.g. from first wake-up time till before last wake-up time, and omits entire days that are considered invalid or incomplete.
+-   In GGIR part 2 the circadian rhythm analysis tries to use every valid data point in the recording.
+-   In GGIR part 6 the circadian rhythm analysis uses a specific time window as defined by the user, e.g. from first wake-up time till before last wake-up time, and omits entire days that are considered invalid or incomplete.
 
 As a result, when your intention is to compare physical activity or sleep estimates derived in GGIR, part 6 estimates ensure that the same days and time points are used as in the other parts of GGIR. Whereas, part 2 does not attempt to do so.
 
@@ -78,16 +78,27 @@ Phi indicates how correlated the multi-day acceleration time series is with itse
 
 ### Self-similarity paramerter (SSP)
 
-The self-similarity paramter (SSP) is also known as scaling exponent or alpha. SSP is a real number between zero and two. Values in the range (0, 1) indicate stationary motion behaviour. Values int he range (1, 2 indicate nonstationary motion behaviour. For details see [Mesquita et al 2020](https://doi.org/10.1093/bioinformatics/btaa955) and [Danilevicz et al. 2024](https://doi.org/10.1186/s12874-024-02255-w).
+The self-similarity parameter (SSP) is also known as scaling exponent or alpha. The overall SSP is a real number between zero and two. Values in the range (0, 1) indicate stationary motion behaviour. Values int he range (1, 2 indicate non-stationary motion behaviour. For details see [Mesquita et al 2020](https://doi.org/10.1093/bioinformatics/btaa955) and [Danilevicz et al. 2024](https://doi.org/10.1186/s12874-024-02255-w).
+
+In human motor activity, researchers often observe that the scaling behaviour is not uniform across all time scales but exhibits a "crossover point" around 1.5 to 2 hours. To capture distinct biological mechanisms, GGIR splits the DFA into two distinct time-scale regions:
+
+-   **Short-term scaling exponent (**α1**or SSP_short):** Measures the fractal temporal correlations at time scales smaller than or equal to 90 minutes. Fluctuations at this scale are regulated by multiple physiological controls, including general mobility, mood, and cognition.
+
+-   **Long-term scaling exponent (**α2**or SSP_long):** Measures the temporal correlations at time scales between 120 and 600 minutes (2 to 10 hours). Fluctuations at these larger time scales are heavily reliant on the central circadian pattern.
+
+Evaluating these regions separately allows researchers to disentangle physical frailty from central circadian degradation. For example, the difference between the two regions (Δα=α1−α2) quantifies the specific breakdown of multi-scale fractal patterns and can serve as a biomarker for SCN neurodegeneration and dementia severity. For more details on the clinical application of these piecewise metrics, see [Li et al. 2019](https://www.google.com/url?sa=E&q=https%3A%2F%2Fdoi.org%2F10.1126%2Fscitranslmed.aax1977).
+
+Additionally, the output includes the difference between the short-term and long-term scaling exponents (Δα=α1 − α2) as SSP_diff in the part 6 report. In healthy individuals, fractal regulation is uniform across all time scales, yielding a difference close to zero. A larger difference indicates a disruption in multiscale fractal regulation
 
 ### Activity Balance Index (ABI)
 
 The Activity Balance Index (ABI) was introduced by [Danilevicz et al. 2024](https://doi.org/10.1186/s12874-024-02255-w) and is a transformation of SSP. ABI measures how the activity over the observed period is balanced, higher values reflect a more balanced pattern of activity. ABI is a real number between zero and one and calculated from the acceleration metric time series directly without the need for cut-points.
 
+In addition to the overall ABI (**ABI_overall**), GGIR calculates **ABI_short** and **ABI_long**, which are the Activity Balance Indices transformed from the short-term (SSP_short) and long-term (SSP2) scaling exponents, respectively. This allows for a direct evaluation of behavioural balance at both the general mobility and circadian levels.
+
 ### Sleep Regularity Index (SRI)
 
-A discussion of the Sleep Regularity Index can be found in
-<https://wadpac.github.io/GGIR/articles/SleepRegularityIndex.html>
+A discussion of the Sleep Regularity Index can be found in <https://wadpac.github.io/GGIR/articles/SleepRegularityIndex.html>
 
 ## Related output