Move time series preparation for cosinor analysis to new function

NAMESPACE

@@ -35,7 +35,7 @@ export(g.analyse, g.calibrate,
        check_params, extract_params,
        g.imputeTimegaps, g.part5.classifyNaps,
        tidyup_df, cosinorAnalyses,
-       detect_nonwear_clipping,
+       detect_nonwear_clipping, applyCosinorAnalyses,
        ShellDoc2Vignette, parametersVignette,
        correctOlderMilestoneData,
        convertEpochData, appendRecords, extractID,
@@ -53,7 +53,7 @@ importFrom("utils", "read.csv", "sessionInfo", "write.csv",
 importFrom("stats", "aggregate.data.frame", "weighted.mean",
            "rnorm","median","aggregate","C", "lm.wfit",
            "quantile", "sd","coef", "lm", "embed", "na.pass",
-            "residuals", "fitted")
+            "residuals", "fitted", "cor")
 import(data.table)
 importFrom("methods", "is")
 importFrom("utils", "available.packages", "packageVersion", "help", "read.table")

NEWS.md

@@ -1,9 +1,12 @@
 # CHANGES IN GGIR VERSION 2.10-5
 
+- Part 2: Move cosinor analysis code to its own function in order to ease re-using it in both part 2 and part 6.
+
+- Part2: Expand cosinor analysis with R2
+
 - Part 2: Arguments hrs.del.start and hrs.del.end when combined with strategy = 3 and strategy = 5 now count
 relative to start and end of the most active time window as identified.  #905
 
-
 # CHANGES IN GGIR VERSION 2.10-4
 
 - Part 4: Now better able to handle nights without sustained inactivity bouts (rest) #911

R/applyCosinorAnalyses.R

@@ -0,0 +1,72 @@
+applyCosinorAnalyses = function(ts, qcheck, midnightsi, epochsizes) {
+  # 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)
+  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
+  dstgap = which(diff(midnightsi) != Nlong_epochs_day)
+  if (length(dstgap) > 0) {
+    # Time moved forward due to DST
+    gaplocation = ((midnightsi[dstgap[1]] * ws2) / ws3) + (2 * (3600/ws3))
+    # Insert NA values
+    Xi = c(Xi[1:gaplocation], rep(NA, 3600/ws3), Xi[(gaplocation + 1):length(Xi)])
+    qcheck = c(qcheck[1:gaplocation], rep(NA, 3600/ws3), qcheck[(gaplocation + 1):length(qcheck)])
+  }
+
+  # Xi = log((Xi * 1000) + 1)  # log transformed to be more robust against peaks in the data
+  # set non-wear to missing values, because for Cosinor fit
+  # it seems more logical to only fit with real data
+  # this comes at the price of not being able to extract F_pseudo
+  firstvalid = 1
+  if (length(which(qcheck == 1)) > 0) {
+    is.na(Xi[which(qcheck == 1)]) = TRUE
+    # ignore invalid start of recording (if applicable)
+    # such that 24 hour blocks start from first valid value
+    firstvalid = which(qcheck == 0)[1]
+    if (is.na(firstvalid) == FALSE) {
+      if (firstvalid != 1) {
+        Xi = Xi[firstvalid:length(Xi)]
+      }
+    }
+  }
+  if (length(which(is.na(Xi) == FALSE)) > (1440 * (60/ws3))) { # Only attempt cosinor analyses if there is more than 24 hours of data
+    midnightsi_ws3 = (midnightsi - 1) * (ws2 / ws3)
+    timeOffsetHours = (midnightsi_ws3[which(midnightsi_ws3 >= firstvalid - 1)[1]] - (firstvalid - 1)) / (3600 / ws3)
+    if (ws3 < 60) {
+      # If epochsize < 1 minute then aggregate to 1 minute by taking average value
+      # but keep NA values
+      XTtime = rep(1:length(Xi), each = 60 / ws3)
+      XT = data.frame(Xi = Xi, time = XTtime[1:length(Xi)])
+      custommean = function(x) {
+        y = NA
+        if (length(x) > 0) {
+          if (length(which(is.na(x) == FALSE) ) > 0) {
+            y = mean(x, na.rm = TRUE)
+          }
+        }
+        return(y)
+      }
+      XT = aggregate(x = XT, by = list(XT$time), FUN = custommean)
+      if (length(which(is.nan(XT$Xi) == TRUE)) > 0) {
+        is.na(XT$Xi[which(is.nan(XT$Xi) == TRUE)]) = TRUE
+      }
+      Xi = XT$Xi
+      epochsize = 60
+    } 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$timeOffsetHours = timeOffsetHours
+  } else {
+    cosinor_coef = c()
+  }
+  return(cosinor_coef)
+}

R/cosinorAnalyses.R

@@ -30,5 +30,16 @@ cosinorAnalyses = function(Xi, epochsize = 60, timeOffsetHours = 0) {
                 IVIS_acc_threshold = log(20 + 1),
                 IVIS_per_daypair = TRUE) # 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
+
+  # # this should equal: https://en.wikipedia.org/wiki/Coefficient_of_determination
+  # yi = coefext$cosinor_ts$original
+  # fi = coefext$cosinor_ts$fittedY
+  # meanY = mean(coefext$cosinor_ts$original)
+  # SSres = sum((yi - fi)^2)
+  # SStot = sum((y - meanY)^2)
+  # R2 = 1 - (SSres / SStot)
+
   invisible(list(coef = coef, coefext = coefext, IVIS = IVIS))
 }

R/g.analyse.avday.R

@@ -117,75 +117,9 @@ g.analyse.avday = function(doquan, averageday, M, IMP, t_TWDI, quantiletype,
   #----------------------------------
   # (Extended) Cosinor analysis
   if (params_247[["cosinor"]] == TRUE) {
-    # Re-derive Xi but this time include entire time series
-    # Here, we ignore duplicated values (when clock moves backward due to DST)
-    handleDST = !duplicated(IMP$metashort)
-    qcheck = qcheck[handleDST]
-    Xi = IMP$metashort[handleDST, acc.metric]
-    Nlong_epochs_day =  (1440 * 60) / ws2 # this is 96 by default
-    dstgap = which(diff(midnightsi) != Nlong_epochs_day)
-    if (length(dstgap) > 0) {
-      # Time moved forward due to DST
-      gaplocation = ((midnightsi[dstgap[1]] * ws2) / ws3) + (2 * (3600/ws3))
-      # Insert NA values
-      Xi = c(Xi[1:gaplocation], rep(NA, 3600/ws3), Xi[(gaplocation + 1):length(Xi)])
-      qcheck = c(qcheck[1:gaplocation], rep(NA, 3600/ws3), qcheck[(gaplocation + 1):length(qcheck)])
-    }
-
-    # Xi = log((Xi * 1000) + 1)  # log transformed to be more robust against peaks in the data
-    # set non-wear to missing values, because for Cosinor fit
-    # it seems more logical to only fit with real data
-    # this comes at the price of not being able to extract F_pseudo
-    firstvalid = 1
-    if (length(which(qcheck == 1)) > 0) {
-      is.na(Xi[which(qcheck == 1)]) = TRUE
-      # ignore invalid start of recording (if applicable)
-      # such that 24 hour blocks start from first valid value
-      firstvalid = which(qcheck == 0)[1]
-      if (is.na(firstvalid) == FALSE) {
-        if (firstvalid != 1) {
-          Xi = Xi[firstvalid:length(Xi)]
-        }
-      }
-    }
-    if (length(which(is.na(Xi) == FALSE)) > (1440 * (60/ws3))) { # Only attempt cosinor analyses if there is more than 24 hours of data
-      midnightsi_ws3 = (midnightsi - 1) * (ws2 / ws3)
-      timeOffsetHours = (midnightsi_ws3[which(midnightsi_ws3 >= firstvalid - 1)[1]] - (firstvalid - 1)) / (3600 / ws3)
-      if (ws3 < 60) {
-        # If epochsize < 1 minute then aggregate to 1 minute by taking maximum value
-        # but keep NA values
-        XTtime = rep(1:length(Xi), each = 60 / ws3)
-        XT = data.frame(Xi = Xi, time = XTtime[1:length(Xi)])
-        custommean = function(x) {
-          y = NA
-          if (length(x) > 0) {
-            if (length(which(is.na(x) == FALSE) ) > 0) {
-              y = mean(x, na.rm = TRUE)
-            }
-          }
-          return(y)
-        }
-        XT = aggregate(x = XT, by = list(XT$time), FUN = custommean)
-        if (length(which(is.nan(XT$Xi) == TRUE)) > 0) {
-          is.na(XT$Xi[which(is.nan(XT$Xi) == TRUE)]) = TRUE
-        }
-        # experimental: clip all peaks above Xth percentile?
-        # Q9 = quantile(x = XT$Xi, probs = 0.75, na.rm = TRUE)
-        # XT$Xi[which(XT$Xi >= Q9)] = Q9
-
-        # log transform of data in millig
-        notna = !is.na(XT$Xi)
-        XT$Xi[notna] = log((XT$Xi[notna]*1000) + 1) 
-        Xi = XT$Xi
-        epochsize = 60
-      } else {
-        epochsize = ws3
-      }
-      cosinor_coef = cosinorAnalyses(Xi = Xi, epochsize = epochsize, timeOffsetHours = timeOffsetHours) 
-      cosinor_coef$timeOffsetHours = timeOffsetHours
-    } else {
-      cosinor_coef = c()
-    }
+    cosinor_coef = applyCosinorAnalyses(ts = IMP$metashort[, c("timestamp", acc.metric)],
+                                        qcheck = qcheck,
+                                        midnightsi, epochsizes = c(ws3, ws2))
   } else {
     cosinor_coef = c()
   }

R/g.analyse.perfile.R

@@ -156,15 +156,16 @@ g.analyse.perfile = function(I, C, metrics_nav,
       filesummary[vi]  = c(cosinor_coef$timeOffsetHours)
       s_names[vi] = c("cosinor_timeOffsetHours")
       vi = vi + 1
-      try(expr = {filesummary[vi:(vi + 4)]  = as.numeric(c(cosinor_coef$coef$params$mes,
+      try(expr = {filesummary[vi:(vi + 5)]  = as.numeric(c(cosinor_coef$coef$params$mes,
                                                            cosinor_coef$coef$params$amp,
                                                            cosinor_coef$coef$params$acr,
                                                            cosinor_coef$coef$params$acrotime,
-                                                           cosinor_coef$coef$params$ndays))}, silent = TRUE)
-      s_names[vi:(vi + 4)] = c("cosinor_mes", "cosinor_amp", "cosinor_acrophase",
-                               "cosinor_acrotime", "cosinor_ndays")
-      vi = vi + 5
-      try(expr = {filesummary[vi:(vi + 9)]  = c(cosinor_coef$coefext$params$minimum,
+                                                           cosinor_coef$coef$params$ndays,
+                                                           cosinor_coef$coef$params$R2))}, silent = TRUE)
+      s_names[vi:(vi + 5)] = c("cosinor_mes", "cosinor_amp", "cosinor_acrophase",
+                               "cosinor_acrotime", "cosinor_ndays", "cosinor_R2")
+      vi = vi + 6
+      try(expr = {filesummary[vi:(vi + 10)]  = c(cosinor_coef$coefext$params$minimum,
                                                 cosinor_coef$coefext$params$amp,
                                                 cosinor_coef$coefext$params$alpha,
                                                 cosinor_coef$coefext$params$beta,
@@ -173,16 +174,19 @@ g.analyse.perfile = function(I, C, metrics_nav,
                                                 cosinor_coef$coefext$params$DownMesor,
                                                 cosinor_coef$coefext$params$MESOR,
                                                 cosinor_coef$coefext$params$ndays,
-                                                cosinor_coef$coefext$params$F_pseudo)}, silent = TRUE)
-      s_names[vi:(vi + 9)] = c("cosinorExt_minimum", "cosinorExt_amp", "cosinorExt_alpha",
+                                                cosinor_coef$coefext$params$F_pseudo,
+                                                cosinor_coef$coefext$params$R2)}, silent = TRUE)
+      s_names[vi:(vi + 10)] = c("cosinorExt_minimum", "cosinorExt_amp", "cosinorExt_alpha",
                                "cosinorExt_beta", "cosinorExt_acrotime", "cosinorExt_UpMesor",
                                "cosinorExt_DownMesor", "cosinorExt_MESOR",
-                               "cosinorExt_ndays", "cosinorExt_F_pseudo")
-      vi = vi + 10
+                               "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
+    } else {
+      vi = vi + 20
     }
 
     # Variables per metric - summarise with stratification to weekdays and weekend days

R/g.part4.R

@@ -62,11 +62,15 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1,
   # initialize output variable names
   colnamesnightsummary = c("ID", "night", "sleeponset", "wakeup", "SptDuration", "sleepparam", "guider_onset",
                            "guider_wakeup", "guider_SptDuration", "error_onset", "error_wake", "error_dur", "fraction_night_invalid",
-                           "SleepDurationInSpt", "WASO", "duration_sib_wakinghours", "number_sib_sleepperiod", "number_of_awakenings",
-                           "number_sib_wakinghours", "duration_sib_wakinghours_atleast15min", "sleeponset_ts", "wakeup_ts", "guider_onset_ts",
-                           "guider_wakeup_ts", "sleeplatency", "sleepefficiency", "page", "daysleeper", "weekday", "calendar_date",
+                           "SleepDurationInSpt", "WASO", "duration_sib_wakinghours",
+                           "number_sib_sleepperiod", "number_of_awakenings",
+                           "number_sib_wakinghours", "duration_sib_wakinghours_atleast15min",
+                           "sleeponset_ts", "wakeup_ts", "guider_onset_ts",  "guider_wakeup_ts",
+                           "sleeplatency", "sleepefficiency", "page", "daysleeper", "weekday", "calendar_date",
                            "filename", "cleaningcode", "sleeplog_used", "sleeplog_ID", "acc_available", "guider", "SleepRegularityIndex", "SriFractionValid",
-                           "longitudinal_axis")
+                           "longitudinal_axis") #
+
+
   if (params_output[["storefolderstructure"]] == TRUE) {
     colnamesnightsummary = c(colnamesnightsummary, "filename_dir", "foldername")
   }
@@ -344,10 +348,10 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1,
           # initialize dataframe to hold sleep period overview:
           spocum = data.frame(nb = numeric(0), start = numeric(0),  end = numeric(0),
                               dur = numeric(0), def = character(0))
-
+          
           spocumi = 1  # counter for sleep periods
           # continue now with the specific data of the night
-
+          
           guider.df2 = guider.df[which(guider.df$night == j), ]
           # ================================================================================ get
           # sleeplog (or HDCZA or L5+/-6hr algorithm) onset and waking time and assess whether it
@@ -457,6 +461,7 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1,
               } else {
                 acc_available = TRUE
               }
+
               if (nrow(sleepdet) == 0) next
               ki = which(sleepdet$definition == k)
               if (length(ki) == 0) next
@@ -681,13 +686,13 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1,
                   }
                   delta_t1 = diff(as.numeric(spocum.t$end))
                   spocum.t$dur = correct01010pattern(spocum.t$dur)
-
+                  
                   #----------------------------
                   nightsummary[sumi, 1] = accid
                   nightsummary[sumi, 2] = j  #night
                   # remove double rows
                   spocum.t = spocum.t[!duplicated(spocum.t), ]
-
+                  
                   #------------------------------------
                   # ACCELEROMETER
                   if (length(which(as.numeric(spocum.t$dur) == 1)) > 0) {
@@ -868,6 +873,7 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1,
                     spocum.t.dur_sibd = 0
                     spocum.t.dur_sibd_atleast15min = 0
                   }
+
                   nightsummary[sumi, 14] = spocum.t.dur.noc  #total nocturnalsleep /accumulated sleep duration
                   nightsummary[sumi, 15] = nightsummary[sumi, 5] - spocum.t.dur.noc  #WASO
                   nightsummary[sumi, 16] = spocum.t.dur_sibd  #total sib (sustained inactivty bout) duration during wakinghours
@@ -1043,11 +1049,10 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1,
         if (length(nnights.list) == 0) {
           # if there were no nights to analyse
           nightsummary[sumi, 1:2] = c(accid, 0)
-          nightsummary[sumi, 3:30] = NA
+          nightsummary[sumi, c(3:30, 34, 36:39)] = NA
           nightsummary[sumi, 31] = fnames[i]
           nightsummary[sumi, 32] = 4  #cleaningcode = 4 (no nights of accelerometer available)
-          nightsummary[sumi, c(33,35)] = c(FALSE, TRUE)  #sleeplog_used acc_available
-          nightsummary[sumi, 36:39] = NA
+          nightsummary[sumi, c(33, 35)] = c(FALSE, TRUE)  #sleeplog_used acc_available
           if (params_output[["storefolderstructure"]] == TRUE) {
             nightsummary[sumi, 40:41] = c(ffd[i], ffp[i])  #full filename structure and use the lowest foldername as foldername name
           }

R/g.report.part4.R

@@ -295,13 +295,18 @@ g.report.part4 = function(datadir = c(), metadatadir = c(), loglocation = c(),
             if (dotwice == 1) {
               nightsummary.tmp = turn_numeric(x = nightsummary.tmp, varnames = gdn)
             }
+            varnames_tmp = c("SptDuration", "sleeponset",
+                             "wakeup", "WASO", "SleepDurationInSpt",
+                             "number_sib_sleepperiod", "duration_sib_wakinghours",
+                             "number_of_awakenings", "number_sib_wakinghours", 
+                             "duration_sib_wakinghours_atleast15min",
+                             "sleeplatency", "sleepefficiency", "number_of_awakenings",
+                             "guider_inbedDuration", "guider_inbedStart",
+                             "guider_inbedEnd", "guider_SptDuration", "guider_onset",
+                             "guider_wakeup", "SleepRegularityIndex",
+                             "SriFractionValid")
             nightsummary.tmp = turn_numeric(x = nightsummary.tmp, 
-                                            varnames = c("SptDuration", "sleeponset",
-                                                         "wakeup", "WASO", "SleepDurationInSpt", "number_sib_sleepperiod", "duration_sib_wakinghours",
-                                                         "number_of_awakenings", "number_sib_wakinghours", "duration_sib_wakinghours_atleast15min",
-                                                         "sleeplatency", "sleepefficiency", "number_of_awakenings", "guider_inbedDuration", "guider_inbedStart",
-                                                         "guider_inbedEnd", "guider_SptDuration", "guider_onset", "guider_wakeup", "SleepRegularityIndex",
-                                                         "SriFractionValid"))
+                                            varnames = varnames_tmp)
             weekday = nightsummary.tmp$weekday[this_sleepparam]
             if (dotwice == 1) {
               for (k in 1:3) {
@@ -428,8 +433,7 @@ g.report.part4 = function(datadir = c(), metadatadir = c(), loglocation = c(),
                                                                    udefn[j], "_mn", sep = ""), paste("average_dur_sib_wakinghours_", TW, "_", udefn[j],
                                                                                                      "_sd", sep = ""))
                   NDAYsibd = length(which(nightsummary.tmp$number_sib_wakinghours[indexUdef] > 0))
-                  if (length(NDAYsibd) == 0)
-                    NDAYsibd = 0
+                  if (length(NDAYsibd) == 0) NDAYsibd = 0
                   personSummary[i, (cnt + 19)] = NDAYsibd
                   personSummarynames = c(personSummarynames, paste("n_days_w_sib_wakinghours_", TW, "_",
                                                                    udefn[j], sep = ""))