IVIS now use revised version throughout GGIR, remove log transformati…

…on, update documentation #955
wadpac · Feb 14, 2024 · 2a601c4 · 2a601c4
1 parent 1d368ec
commit 2a601c4
Show file tree

Hide file tree

Showing 17 changed files with 93 additions and 123 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -2,6 +2,9 @@
 
 - Part 6: Add DFA functionality #839 (credits: Ian Meneghel Danilevicz and Victor Barreto Mesquita)
 
+- 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.
+
+
 # CHANGES IN GGIR VERSION 3.0-5
 
 - Part 5: Fix bug in functionality for Sensewear data (externally derived epoch data) #1030

diff --git a/R/applyCosinorAnalyses.R b/R/applyCosinorAnalyses.R
@@ -1,4 +1,4 @@
-applyCosinorAnalyses = function(ts, qcheck, midnightsi, epochsizes) {
+applyCosinorAnalyses = function(ts, qcheck, midnightsi, epochsizes, threshold = NULL) {
   # qcheck - vector of length ts to indicate invalid values
   ws2 = epochsizes[2]
   ws3 = epochsizes[1]
@@ -59,11 +59,13 @@ 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) 
+    # transform data to millig if data is stored in g-units
+    if (max(Xi, na.rm = TRUE) < 8) Xi[notna] = Xi[notna] * 1000
+    # log transform, add 1 to avoid log(0)
+    Xi[notna] = log(Xi[notna] + 1)
+    cosinor_coef = cosinorAnalyses(Xi = Xi, epochsize = epochsize, 
+                                   timeOffsetHours = timeOffsetHours, threshold = threshold) 
     cosinor_coef$timeOffsetHours = timeOffsetHours
   } else {
     cosinor_coef = c()

diff --git a/R/cosinorAnalyses.R b/R/cosinorAnalyses.R
@@ -1,4 +1,4 @@
-cosinorAnalyses = function(Xi, epochsize = 60, timeOffsetHours = 0) {
+cosinorAnalyses = function(Xi, epochsize = 60, timeOffsetHours = 0, threshold = NULL) {
   # Apply Cosinor function from ActRC
   N = 1440 * (60 / epochsize) # Number of epochs per day
   Xi = Xi[1:(N * floor(length(Xi) / N))] # ActCR expects integer number of days
@@ -23,9 +23,9 @@ cosinorAnalyses = function(Xi, epochsize = 60, timeOffsetHours = 0) {
   k = ceiling(abs(coef$params$acr) / (pi * 2))
   if (coef$params$acr < 0) coef$params$acr = coef$params$acr + (k * 2 * pi)
   # Perform IVIS on the same input signal to allow for direct comparison
-  IVIS = g.IVIS(Xi = Xi / 1000, # divide by 1000 because function g.IVIS internally multiplies by 1000 when IVIS.activity.metric = 2
+  IVIS = g.IVIS(Xi = (exp(Xi) - 1), # undo log transformation for IV IS
                 epochSize = epochsize, 
-                threshold = log(20 + 1)) # take log, because Xi is logtransformed with offset of 1
+                threshold = threshold) # take log, because Xi is logtransformed with offset of 1
 
   coefext$params$R2 = cor(coefext$cosinor_ts$original, coefext$cosinor_ts$fittedYext)^2
   coef$params$R2 = cor(coefext$cosinor_ts$original, coefext$cosinor_ts$fittedY)^2

diff --git a/R/g.IVIS.R b/R/g.IVIS.R
@@ -1,10 +1,9 @@
 g.IVIS = function(Xi, epochSize = 60, threshold = NULL) {
   if (!is.null(threshold)) {
-    if (max(Xi, na.rm = TRUE) < 8) Xi = Xi * 1000
     Xi = ifelse(Xi > threshold, 0, 1)  
   }
   if (epochSize < 3600) {
-    # Always aggregate Xi to 1 hour resolution
+    # Always aggregate Xi to 1 hour resolution if it is not provided at 1 hour resolution
     df = data.frame(Xi = Xi, time = numeric(length(Xi)))
     time = seq(0, length(Xi) * epochSize, by = epochSize)
     df$time = time[1:nrow(df)]
@@ -24,7 +23,6 @@ g.IVIS = function(Xi, epochSize = 60, threshold = NULL) {
       aveDay = aggregate(. ~ hour, data = ts, mean, na.action = na.omit)
       Xh = aveDay$Xi
       Xm = suppressWarnings(mean(Xh,na.rm = TRUE)) # Average acceleration per day
-      N = length(Xi[!is.na(Xi)])
       deltaXi = diff(Xi)^2
       N = length(Xi[!is.na(Xi)])
 

diff --git a/R/g.analyse.R b/R/g.analyse.R
@@ -207,7 +207,7 @@ g.analyse =  function(I, C, M, IMP, params_247 = c(), params_phyact = c(),
                                  quantiletype = quantiletype, ws3 = ws3,
                                  doiglevels = doiglevels, firstmidnighti = firstmidnighti, ws2 = ws2,
                                  midnightsi = midnightsi, params_247 = params_247, qcheck = qcheck,
-                                 acc.metric = acc.metric)
+                                 acc.metric = acc.metric, params_phyact = params_phyact)
   cosinor_coef = output_avday$cosinor_coef
 
   #--------------------------------------------------------------

diff --git a/R/g.analyse.avday.R b/R/g.analyse.avday.R
@@ -1,6 +1,6 @@
 g.analyse.avday = function(doquan, averageday, M, IMP, t_TWDI, quantiletype,
                            ws3, doiglevels, firstmidnighti, ws2, midnightsi, params_247 = c(), 
-                           qcheck = c(), acc.metric = c(), ...) {
+                           qcheck = c(), acc.metric = c(), params_phyact, ...) {
   #get input variables
   input = list(...)
   if (length(input) > 0 || length(params_247) == 0) {
@@ -99,32 +99,21 @@ g.analyse.avday = function(doquan, averageday, M, IMP, t_TWDI, quantiletype,
       }
     }
   }
-  # IS and IV variables
-  fmn = midnightsi[1] * (ws2/ws3) # select data from first midnight to last midnight because we need full calendar days to compare
-  lmn = (midnightsi[length(midnightsi)] * (ws2/ws3)) - 1
-  # By using the metashort from the IMP we do not need to ignore segments, because data imputed
-  Xi = IMP$metashort[fmn:lmn, acc.metric]
-  # IV IS
-  IVISout = g.IVIS(Xi, epochSize = ws3,
-                   threshold = params_247[["IVIS_acc_threshold"]])
-  InterdailyStability = IVISout$InterdailyStability
-  IntradailyVariability = IVISout$IntradailyVariability
-  rm(Xi)
-
   #----------------------------------
-  # (Extended) Cosinor analysis
+  # Cosinor analysis, (Extended) Cosinor analysis, including IV, IS, and phi
+  # based on time series where invalid data points are set to NA
   if (params_247[["cosinor"]] == TRUE) {
     cosinor_coef = applyCosinorAnalyses(ts = IMP$metashort[, c("timestamp", acc.metric)],
                                         qcheck = qcheck,
-                                        midnightsi, epochsizes = c(ws3, ws2))
+                                        midnightsi, epochsizes = c(ws3, ws2),
+                                        threshold = params_phyact[["threshold.lig"]])
   } else {
     cosinor_coef = c()
   }
-  invisible(list(InterdailyStability = InterdailyStability,
-                 IntradailyVariability = IntradailyVariability, 
-                 igfullr_names = igfullr_names,
+  invisible(list(igfullr_names = igfullr_names,
                  igfullr = igfullr, QUAN = QUAN,
                  qlevels_names = qlevels_names,
                  ML5AD = ML5AD,
-                 ML5AD_names = ML5AD_names, cosinor_coef = cosinor_coef))
+                 ML5AD_names = ML5AD_names,
+                 cosinor_coef = cosinor_coef))
 }
diff --git a/R/g.analyse.perfile.R b/R/g.analyse.perfile.R
@@ -154,14 +154,7 @@ g.analyse.perfile = function(I, C, metrics_nav,
     if (length(iNA) > 0) filesummary[(vi:(vi + 1))[iNA]] = 0
     s_names[vi:(vi + 1)] = c("N valid WEdays","N valid WKdays")
     vi = vi + 2
-    # Add ISIV to filesummary
-    filesummary[vi:(vi + 2)] = c(output_avday$InterdailyStability, output_avday$IntradailyVariability,
-                                 params_247[["IVIS_windowsize_minutes"]])
-    iNA = which(is.na(filesummary[vi:(vi + 3)]) == TRUE)
-    if (length(iNA) > 0) filesummary[(vi:(vi + 3))[iNA]] = " "
-    s_names[vi:(vi + 2)] = c("IS_interdailystability", "IV_intradailyvariability", "IVIS_windowsize_minutes")
-    vi = vi + 4
-    # Cosinor analysis
+    # Cosinor analysis + IV + IS + phi
     if (length(cosinor_coef) > 0) {
       filesummary[vi]  = c(cosinor_coef$timeOffsetHours)
       s_names[vi] = c("cosinor_timeOffsetHours")
@@ -191,12 +184,13 @@ g.analyse.perfile = function(I, C, metrics_nav,
                                 "cosinorExt_DownMesor", "cosinorExt_MESOR",
                                 "cosinorExt_ndays", "cosinorExt_F_pseudo", "cosinorExt_R2")
       vi = vi + 11
-      filesummary[vi:(vi + 1)]  = c(cosinor_coef$IVIS$InterdailyStability,
-                                    cosinor_coef$IVIS$IntradailyVariability)
-      s_names[vi:(vi + 1)] = c("cosinorIS", "cosinorIV")
-      vi = vi + 2
+      filesummary[vi:(vi + 2)]  = c(cosinor_coef$IVIS$InterdailyStability,
+                                    cosinor_coef$IVIS$IntradailyVariability,
+                                    cosinor_coef$IVIS$phi)
+      s_names[vi:(vi + 2)] = c("IS", "IV", "phi")
+      vi = vi + 3
     } else {
-      vi = vi + 20
+      vi = vi + 21
     }
 
     # Variables per metric - summarise with stratification to weekdays and weekend days

diff --git a/R/g.part6.R b/R/g.part6.R
@@ -201,21 +201,20 @@ g.part6 = function(datadir = c(), metadatadir = c(), f0 = c(), f1 = c(),
                            no = length(which(ts$invalidepoch == 0)) / ((3600 * 24) / epochSize))
       s_names[fi] = "N_valid_days"
       fi = fi + 1
+      #=================================================================
+      # Cosinor analysis which comes with IV IS estimtes
       if (do.cr == TRUE) {
-        # Cosinor analysis (which includes IVIS)
+        # Note: applyCosinorAnalyses below uses column invalidepoch to turn
+        # imputed values to NA.
         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))
+                                            epochsizes = rep(epochSize, 2),
+                                            threshold = params_phyact[["threshold.lig"]])
         rm(acc4cos)
-      } else {
-        cosinor_coef = NULL
-      }
-
-      if (length(cosinor_coef) > 0) {
         summary[fi] = cosinor_coef$timeOffsetHours
         s_names[fi] = "cosinor_timeOffsetHours"
         fi = fi + 1
@@ -253,22 +252,24 @@ g.part6 = function(datadir = c(), metadatadir = c(), f0 = c(), f1 = c(),
         s_names[fi:(fi + 2)] = c("IS", "IV", "phi")
         fi = fi + 3
       } else {
-        cosinor_coef = c()
-        s_names[fi:(fi + 19)] = c("cosinor_timeOffsetHours", "cosinor_mes", 
+        cosinor_coef = NULL
+        s_names[fi:(fi + 20)] = c("cosinor_timeOffsetHours", "cosinor_mes", 
                                   "cosinor_amp", "cosinor_acrophase",
                                   "cosinor_acrotime", "cosinor_ndays", "cosinor_R2", 
                                   "cosinorExt_minimum", "cosinorExt_amp", 
                                   "cosinorExt_alpha", "cosinorExt_beta", 
                                   "cosinorExt_acrotime", "cosinorExt_UpMesor",
                                   "cosinorExt_DownMesor", "cosinorExt_MESOR",
                                   "cosinorExt_ndays", "cosinorExt_F_pseudo", 
-                                  "cosinorExt_R2", "cosinorIS", "cosinorIV")
-        fi = fi + 20
+                                  "cosinorExt_R2", "IS", "IV", "phi")
+        fi = fi + 21
       }
-      if (params_247[["part6DFA"]] == TRUE) {
+
+      #=======================================================================
+      # DFA analyses is used independent of do.cr because it is time consuming
+      if (params_247[["part6DFA"]] == TRUE && do.cr == TRUE) {
         ssp = SSP(ts$ACC)
         abi = ABI(ssp)
-
         summary[fi:(fi + 1)] = c(ssp, abi)
         s_names[fi:(fi + 1)] = c("SSP", "ABI")
         fi = fi + 2

diff --git a/man/GGIR.Rd b/man/GGIR.Rd
@@ -1320,10 +1320,14 @@ GGIR(mode = 1:5,
         over which L5M5 needs to be calculated. Now this is done with argument qwindow.}
 
       \item{cosinor}{
-        Boolean (default = FALSE). Whether to apply the cosinor analysis from the ActCR package.}
+        Boolean (default = FALSE). Whether to apply the cosinor analysis from the ActCR package in part 2.
+        In part 6 cosinor analysis is applied by default and cannot be turned off.}
         
       \item{part6CR}{
-        Boolean (default = FALSE) to indicate whether circadian rhythm analysis should be run by part 6.
+        Boolean (default = FALSE) to indicate whether circadian rhythm analysis should
+        be run by part 6, this includes: cosinor analysis, extended cosinor analysis,
+        IS, IV, and phi. Optionally this can be expanded with detrended fluctutation analysis
+        which is controlled by parameter `part6DFA`.
       }
       \item{part6HCA}{
         Boolean (default = FALSE) to indicate whether Household Co Analysis should 
@@ -1342,7 +1346,7 @@ GGIR(mode = 1:5,
       }
       \item{part6DFA}{
         Boolean (default = FALSE) to indicate whether to perform Detrended Fluctuation
-        Analysis. Turned off by default because it can be very time consuming.
+        Analysis. Turned off by default because it can be time consuming.
       }
     }
   }

diff --git a/man/applyCosinorAnalyses.Rd b/man/applyCosinorAnalyses.Rd
@@ -8,7 +8,7 @@
   before applying the cosinorAnlayses
 }
 \usage{
-  applyCosinorAnalyses(ts, qcheck, midnightsi, epochsizes)
+  applyCosinorAnalyses(ts, qcheck, midnightsi, epochsizes, threshold = NULL)
 }
 \arguments{
   \item{ts}{
@@ -24,6 +24,9 @@
   \item{epochsizes}{
     Epoch size for ts and qcheck respectively
   }
+  \item{threshold}{
+    See \link{cosinorAnalyses}
+  }
 }
 \author{
   Vincent T van Hees <v.vanhees@accelting.com>

diff --git a/man/cosinorAnalyses.Rd b/man/cosinorAnalyses.Rd
@@ -4,10 +4,11 @@
   Apply cosinor anlaysis and extended cosinor analysis
 }
 \description{
-  Applies cosinor anlaysis from the ActCR package to the time series
+  Applies cosinor anlaysis from the ActCR package to the time series, as well
+  as IV, IS and phi estimates.
 }
 \usage{
-  cosinorAnalyses(Xi, epochsize = 60, timeOffsetHours = 0)
+  cosinorAnalyses(Xi, epochsize = 60, timeOffsetHours = 0, threshold = NULL)
 }
 \arguments{
   \item{Xi}{
@@ -19,6 +20,10 @@
   \item{timeOffsetHours}{
     Numeric time in hours relative to next midnight
   }
+  \item{threshold}{
+    Numeric value to use as threshold to distinguish inactivity from active behaviour
+    for the IV and IS analysis. GGIR uses parameter threshold.lig to set this threshold.
+  }
 }
 \keyword{internal}
 \author{

diff --git a/man/g.IVIS.Rd b/man/g.IVIS.Rd
@@ -31,6 +31,7 @@
 }
 \keyword{internal}
 \author{
+  Ian Meneghel Danilevicz <ian.meneghel-danilevicz@inserm.fr>
   Vincent T van Hees <v.vanhees@accelting.com>
 }
 \references{

diff --git a/man/g.analyse.avday.Rd b/man/g.analyse.avday.Rd
@@ -10,7 +10,8 @@
 \usage{
   g.analyse.avday(doquan, averageday, M, IMP, t_TWDI, quantiletype,
                    ws3, doiglevels, firstmidnighti, ws2, midnightsi,
-                   params_247 = c(), qcheck = c(), acc.metric = c(), ...)
+                   params_247 = c(), qcheck = c(), acc.metric = c(),
+                   params_phyact = NULL, ...)
 }
 \arguments{
   \item{doquan}{Boolean whether quantile analysis should be done}
@@ -25,27 +26,28 @@
   \item{ws2}{see \link{g.weardec}}
   \item{midnightsi}{see \link{g.detecmidnight}}
   \item{params_247}{
-    See \link{g.part2}
+    See \link{GGIR}
   }
   \item{qcheck}{
     Vector with indicators of when data is valid (value=0) or invalid (value=1).
   }
   \item{acc.metric}{
-    Character, see \link{g.part1}. Here, it is used to decided which acceleration metric
+    Character, see \link{GGIR}. Here, it is used to decided which acceleration metric
     to use for IVIS and cosinor analyses.
   }
+  \item{params_phyact}{
+    See \link{GGIR}
+  }
   \item{...}{
    Any argument used in the previous version of g.analyse.avday, which will now
    be used to overrule the arguments specified with the parameter objects.
   }
 }
 \value{
-  \item{\code{InterdailyStability}}{}
-  \item{\code{IntradailyVariability}}{}
-  \item{\code{igfullr_names}}{}
-  \item{\code{igfullr}}{}
-  \item{\code{QUAN}}{}
-  \item{\code{qlevels_names}}{}
+  \item{\code{igfullr_names}}{Intensity gradient variable names}
+  \item{\code{igfullr}}{Intensity gradient values}
+  \item{\code{QUAN}}{Quantiles}
+  \item{\code{qlevels_names}}{Quantile level names}
   \item{\code{ML5AD}}{}
   \item{\code{ML5AD_names}}{}
 }

diff --git a/tests/testthat/test_cosinor.R b/tests/testthat/test_cosinor.R
@@ -13,6 +13,8 @@ test_that("cosinorAnalyses provides expected output", {
     if (missingdata == TRUE) {
       is.na(counts[N:floor(N * 2)]) = TRUE
     }
+    # Plot code below commented out, but useful for debugging to 
+    # understand how dummy time series looks like.
     # x11()
     # plot(log((counts * 1000) + 1), type = "l")
     mod = cosinorAnalyses(Xi = log((counts * 1000) + 1), epochsize = epochSizeSeconds, timeOffsetHours = timeOffsetHours)
@@ -51,8 +53,8 @@ test_that("cosinorAnalyses provides expected output", {
   expect_equal(coef60$coefext$params$R2, 0.8976208, tolerance  = 0.01)
 
   # IV IS
-  expect_equal(coef60$IVIS$InterdailyStability, 0.9945789, tolerance  = 0.01)
-  expect_equal(coef60$IVIS$IntradailyVariability, 1.10, tolerance  = 0.01)
+  expect_equal(coef60$IVIS$InterdailyStability, 1, tolerance  = 0.01)
+  expect_equal(coef60$IVIS$IntradailyVariability, 0.06774028, tolerance  = 0.01)
 
   fields_to_compare = c("minimum", "amp", "alpha", "beta", "acrotime", "DownMesor", "MESOR")
   expect_equal(coef60$coefext[fields_to_compare], coef300$coefext[fields_to_compare], tolerance = 0.1)

diff --git a/tests/testthat/test_g.IVIS.R b/tests/testthat/test_g.IVIS.R
@@ -2,7 +2,7 @@ library(GGIR)
 context("g.IVIS")
 test_that("g.IVIS returns expected output value without missing data", {
   set.seed(1234)
-  xx = rnorm(n = 1440)
+  xx = rnorm(n = 1440) * 1000
   Xi1 = rep(xx, 7)
   Xi2 = rep(0:7, each = 1440)
   T1 = g.IVIS(Xi1, epochSize = 60)

diff --git a/tests/testthat/test_part6.R b/tests/testthat/test_part6.R
@@ -87,7 +87,7 @@ test_that("Part 6 with household co-analysis", {
   expect_equal(output_part6$cosinor_mes, 2.451769, tolerance = 0.00001)
   expect_equal(output_part6$cosinorExt_minimum, 1.288636, tolerance = 0.00001)
   expect_equal(output_part6$cosinorExt_MESOR, 2.164644, tolerance = 0.00001)
-  expect_equal(sum(output_part6[5:26]), 327.9476, tolerance = 0.0001)
+  expect_equal(sum(output_part6[5:26]), 327.9062, tolerance = 0.0001)
 
 
   # Run Circadian rhythm analysis with non-default window