From e0f5b0e4423749e09fbd43e6780aff88d2548c1d Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Tue, 13 Dec 2022 15:48:50 +0100 Subject: [PATCH 01/34] moving the application of the cosinor anlyses inside g.analyse.avday.R to a seperate function to ease re-using it in g.part5 #708 --- R/applyCosinorAnalyses.R | 75 +++++++++++++++++++ R/g.analyse.avday.R | 141 ++++++++++++++++++------------------ R/g.part5.R | 8 ++ man/applyCosinorAnalyses.Rd | 34 +++++++++ 4 files changed, 189 insertions(+), 69 deletions(-) create mode 100644 R/applyCosinorAnalyses.R create mode 100644 man/applyCosinorAnalyses.Rd diff --git a/R/applyCosinorAnalyses.R b/R/applyCosinorAnalyses.R new file mode 100644 index 000000000..6d1b7b441 --- /dev/null +++ b/R/applyCosinorAnalyses.R @@ -0,0 +1,75 @@ +applyCosinorAnalyses = function(ts, qcheck, midnightsi, shortEpoch, longEpoch) { + # qcheck - vector of length ts to indicate invalid values + ws2 = longEpoch + ws3 = shortEpoch + # 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, which(colnames(ts) %in% "timestamp" == FALSE)] + 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() + } + return(cosinor_coef) +} \ No newline at end of file diff --git a/R/g.analyse.avday.R b/R/g.analyse.avday.R index 47017fb0a..a2df794cf 100644 --- a/R/g.analyse.avday.R +++ b/R/g.analyse.avday.R @@ -120,75 +120,78 @@ 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, shortEpoch = ws3, longEpoch = ws2) + # # 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() + # } } else { cosinor_coef = c() } diff --git a/R/g.part5.R b/R/g.part5.R index cbba4be76..4b0a4b4ad 100644 --- a/R/g.part5.R +++ b/R/g.part5.R @@ -412,6 +412,14 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), } } ts$window = 0 + #=============================================================== + # Cosinor analyses based on only the data used for GGIR part5 + + + + + + #=============================================================== for (TRLi in params_phyact[["threshold.lig"]]) { for (TRMi in params_phyact[["threshold.mod"]]) { for (TRVi in params_phyact[["threshold.vig"]]) { diff --git a/man/applyCosinorAnalyses.Rd b/man/applyCosinorAnalyses.Rd new file mode 100644 index 000000000..1d733a635 --- /dev/null +++ b/man/applyCosinorAnalyses.Rd @@ -0,0 +1,34 @@ +\name{applyCosinorAnalyses} +\alias{applyCosinorAnalyses} +\title{ + Apply Cosinor Analyses to time series +} +\description{ + Wrapper function around \link{cosinorAnalyses} that first prepares the time series + before applying the cosinorAnlayses +} +\usage{ + applyCosinorAnalyses(ts, qcheck, midnightsi, shortEpoch, longEpoch) +} +\arguments{ + \item{ts}{ + Data.frame with timestamps and acceleration metric. + } + \item{qcheck}{ + Vector of equal length as number of rows in ts with value 1 for invalid + timestamps, 0 otherwise. + } + \item{midnightsi}{ + Indices for midnights in the time series + } + \item{shortEpoch}{ + First element of argument windowsizes, the resolution of the ts and qcheck values. + } + \item{longEpoch}{ + Second element of argument windowsizes, the resolution of the midnightsi + values in seconds. + } +} +\author{ + Vincent T van Hees +} \ No newline at end of file From 7eed1098aea112b0d303ce0b46e38a5d6033fb09 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Tue, 13 Dec 2022 17:07:47 +0100 Subject: [PATCH 02/34] #708 implementing applyCosinorAnlayses in g.part5 and making sure it works with both 5 second and 60second aggregated epochs --- R/applyCosinorAnalyses.R | 8 ++--- R/g.analyse.avday.R | 71 +------------------------------------ R/g.part5.R | 39 +++++++++++++++----- man/applyCosinorAnalyses.Rd | 10 ++---- 4 files changed, 39 insertions(+), 89 deletions(-) diff --git a/R/applyCosinorAnalyses.R b/R/applyCosinorAnalyses.R index 6d1b7b441..1d2235ff1 100644 --- a/R/applyCosinorAnalyses.R +++ b/R/applyCosinorAnalyses.R @@ -1,12 +1,12 @@ -applyCosinorAnalyses = function(ts, qcheck, midnightsi, shortEpoch, longEpoch) { +applyCosinorAnalyses = function(ts, qcheck, midnightsi, epochsizes) { # qcheck - vector of length ts to indicate invalid values - ws2 = longEpoch - ws3 = shortEpoch + 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, which(colnames(ts) %in% "timestamp" == FALSE)] + 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) { diff --git a/R/g.analyse.avday.R b/R/g.analyse.avday.R index a2df794cf..c4fab41a2 100644 --- a/R/g.analyse.avday.R +++ b/R/g.analyse.avday.R @@ -122,76 +122,7 @@ g.analyse.avday = function(doquan, averageday, M, IMP, t_TWDI, quantiletype, if (params_247[["cosinor"]] == TRUE) { cosinor_coef = applyCosinorAnalyses(ts = IMP$metashort[, c("timestamp", acc.metric)], qcheck = qcheck, - midnightsi, shortEpoch = ws3, longEpoch = ws2) - # # 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() - # } + midnightsi, epochsizes = c(ws3, ws2)) } else { cosinor_coef = c() } diff --git a/R/g.part5.R b/R/g.part5.R index 4b0a4b4ad..71188880f 100644 --- a/R/g.part5.R +++ b/R/g.part5.R @@ -332,6 +332,24 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), } if (length(tail_expansion_log) != 0 & nrow(ts) > max(nightsi)) nightsi[length(nightsi) + 1] = nrow(ts) # include last window Nts = nrow(ts) + } else { + ts$time = iso8601chartime2POSIX(ts$time,tz = params_general[["desiredtz"]]) + ws3new = ws3 # change because below it is used to decide how many epochs are there in + # extract nightsi again + tempp = unclass(ts$time) + if (is.na(tempp$sec[1]) == TRUE) { + tempp = unclass(as.POSIXlt(ts$time, tz = params_general[["desiredtz"]])) + } + sec = tempp$sec + min = tempp$min + hour = tempp$hour + if (params_general[["dayborder"]] == 0) { + nightsi = which(sec == 0 & min == 0 & hour == 0) + } else { + nightsi = which(sec == 0 & min == (params_general[["dayborder"]] - floor(params_general[["dayborder"]])) * 60 & hour == floor(params_general[["dayborder"]])) #shift the definition of midnight if required + } + if (length(tail_expansion_log) != 0 & nrow(ts) > max(nightsi)) nightsi[length(nightsi) + 1] = nrow(ts) # include last window + Nts = nrow(ts) } # if ("angle" %in% colnames(ts)) { # ts = ts[, -which(colnames(ts) == "angle")] @@ -412,14 +430,7 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), } } ts$window = 0 - #=============================================================== - # Cosinor analyses based on only the data used for GGIR part5 - - - - - - #=============================================================== + for (TRLi in params_phyact[["threshold.lig"]]) { for (TRMi in params_phyact[["threshold.mod"]]) { for (TRVi in params_phyact[["threshold.vig"]]) { @@ -930,6 +941,18 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), } } } + #=============================================================== + # Cosinor analyses based on only the data used for GGIR part5 + if (params_247[["cosinor"]] == TRUE) { + cosinor_coef = applyCosinorAnalyses(ts = ts[, c("time", "ACC")], + qcheck = ts$nonwear, + midnightsi = nightsi, + epochsizes = c(ws3, ws3)) + } else { + cosinor_coef = c() + } + + #=============================================================== } } if ("angle" %in% colnames(ts)) { diff --git a/man/applyCosinorAnalyses.Rd b/man/applyCosinorAnalyses.Rd index 1d733a635..77e3198e8 100644 --- a/man/applyCosinorAnalyses.Rd +++ b/man/applyCosinorAnalyses.Rd @@ -8,7 +8,7 @@ before applying the cosinorAnlayses } \usage{ - applyCosinorAnalyses(ts, qcheck, midnightsi, shortEpoch, longEpoch) + applyCosinorAnalyses(ts, qcheck, midnightsi, epochsizes) } \arguments{ \item{ts}{ @@ -21,12 +21,8 @@ \item{midnightsi}{ Indices for midnights in the time series } - \item{shortEpoch}{ - First element of argument windowsizes, the resolution of the ts and qcheck values. - } - \item{longEpoch}{ - Second element of argument windowsizes, the resolution of the midnightsi - values in seconds. + \item{epochsizes}{ + Epoch seiz for ts and qcheck respectively } } \author{ From 5ef3924bac701f17e53f426640fa0d6bbdd548f9 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Mon, 19 Dec 2022 11:27:59 +0100 Subject: [PATCH 03/34] adding cosinor variables to part 5 output --- R/g.part5.R | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/R/g.part5.R b/R/g.part5.R index 71188880f..a36309705 100644 --- a/R/g.part5.R +++ b/R/g.part5.R @@ -948,6 +948,39 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), qcheck = ts$nonwear, midnightsi = nightsi, epochsizes = c(ws3, ws3)) + if (length(cosinor_coef) > 0) { + # assign same value to all rows to ease creating reports + dsummary[,fi] = c(cosinor_coef$timeOffsetHours) + ds_names[fi] = c("cosinor_timeOffsetHours") + fi = fi + 1 + try(expr = {dsummary[, fi:(fi + 4)] = 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) + ds_names[fi:(fi + 4)] = c("cosinor_mes", "cosinor_amp", "cosinor_acrophase", + "cosinor_acrotime", "cosinor_ndays") + fi = fi + 5 + try(expr = {dsummary[, fi:(fi + 9)] = c(cosinor_coef$coefext$params$minimum, + cosinor_coef$coefext$params$amp, + cosinor_coef$coefext$params$alpha, + cosinor_coef$coefext$params$beta, + cosinor_coef$coefext$params$acrotime, + cosinor_coef$coefext$params$UpMesor, + cosinor_coef$coefext$params$DownMesor, + cosinor_coef$coefext$params$MESOR, + cosinor_coef$coefext$params$ndays, + cosinor_coef$coefext$params$F_pseudo)}, silent = TRUE) + ds_names[fi:(fi + 9)] = c("cosinorExt_minimum", "cosinorExt_amp", "cosinorExt_alpha", + "cosinorExt_beta", "cosinorExt_acrotime", "cosinorExt_UpMesor", + "cosinorExt_DownMesor", "cosinorExt_MESOR", + "cosinorExt_ndays", "cosinorExt_F_pseudo") + fi = fi + 10 + dsummary[fi:(fi + 1)] = c(cosinor_coef$IVIS$InterdailyStability, + cosinor_coef$IVIS$IntradailyVariability) + ds_names[fi:(fi + 1)] = c("cosinorIS", "cosinorIV") + fi = fi + 2 + } } else { cosinor_coef = c() } From 07ed472d4dc2ca229b2681907906b0b4edbf3609 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Mon, 19 Dec 2022 11:28:18 +0100 Subject: [PATCH 04/34] tidying up syntax --- R/g.report.part5.R | 266 ++++++++++++++++++++++++++++++++------------- 1 file changed, 189 insertions(+), 77 deletions(-) diff --git a/R/g.report.part5.R b/R/g.report.part5.R index 1d5bb3cc3..b1c035026 100644 --- a/R/g.report.part5.R +++ b/R/g.report.part5.R @@ -63,8 +63,8 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c return(indices) } ms5.out = "/meta/ms5.out" - if (file.exists(paste(metadatadir,ms5.out,sep=""))) { - if (length(dir(paste(metadatadir,ms5.out,sep=""))) == 0) { + if (file.exists(paste0(metadatadir, ms5.out))) { + if (length(dir(paste0(metadatadir, ms5.out))) == 0) { try.generate.report = FALSE #do not run this function if there is no milestone data from g.part5 } else { try.generate.report = TRUE @@ -76,15 +76,15 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c if (try.generate.report == TRUE) { #====================================================================== # loop through meta-files - fnames.ms5 = list.files(paste0(metadatadir,ms5.out),full.names=TRUE) - if(f1 > length(fnames.ms5)) f1 = length(fnames.ms5) + fnames.ms5 = list.files(paste0(metadatadir,ms5.out), full.names = TRUE) + if (f1 > length(fnames.ms5)) f1 = length(fnames.ms5) cat(" loading all the milestone data from part 5 this can take a few minutes\n") myfun = function(x, expectedCols = c()) { tail_expansion_log = NULL load(file = x) cut = which(output[, 1] == "") if (length(cut) > 0 & length(cut) < nrow(output)) { - output = output[-cut,which(colnames(output) != "")] + output = output[-cut, which(colnames(output) != "")] } out = as.matrix(output) if (length(expectedCols) > 0) { @@ -109,7 +109,7 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c expectedCols = colnames(out_try) # print(fnames.ms5[f0:f1]) outputfinal = as.data.frame(do.call(rbind,lapply(fnames.ms5[f0:f1],myfun, expectedCols)), stringsAsFactors = FALSE) - cut = which(sapply(outputfinal, function(x) all(x=="")) == TRUE) # Find columns filled with missing values which(output[1,] == "" & output[2,] == "") + cut = which(sapply(outputfinal, function(x) all(x == "")) == TRUE) # Find columns filled with missing values which(output[1,] == "" & output[2,] == "") if (length(cut) > 0) { outputfinal = outputfinal[,-cut] } @@ -160,57 +160,61 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c CN = colnames(outputfinal) outputfinal2 = outputfinal colnames(outputfinal2) = CN - delcol = which(colnames(outputfinal2) == "window" | colnames(outputfinal2) == "TRLi" | - colnames(outputfinal2) == "TRMi" | colnames(outputfinal2) == "TRVi" | - colnames(outputfinal2) == "sleepparam") + delcol = grep(pattern = "window|TRLi|TRMi|TRVi|sleepparam", + x = colnames(outputfinal2)) outputfinal2 = outputfinal2[,-delcol] OF3 = outputfinal2[seluwi,] OF3 = as.data.frame(OF3, stringsAsFactors = TRUE) #------------------------------------------------------------- # store all summaries in csv files without cleaning criteria OF3_clean = tidyup_df(OF3) - data.table::fwrite(OF3_clean,paste(metadatadir,"/results/QC/part5_daysummary_full_", - uwi[j],"_L",uTRLi[h1],"M",uTRMi[h2],"V",uTRVi[h3], - "_",usleepparam[h4],".csv",sep=""),row.names=FALSE, na = "") + colsWithoutCosinor = grep(pattern = "cosinor", x = colnames(OF3_clean), invert = TRUE) + data.table::fwrite(x = OF3_clean[, colsWithoutCosinor], + file = paste(metadatadir,"/results/QC/part5_daysummary_full_", + uwi[j], "_L", uTRLi[h1], "M", uTRMi[h2], "V", uTRVi[h3], + "_", usleepparam[h4], ".csv", sep = ""), row.names = FALSE, na = "") # store all summaries in csv files with cleaning criteria validdaysi = getValidDayIndices(OF3,includedaycrit.part5, excludefirstlast.part5, window = uwi[j]) - data.table::fwrite(OF3_clean[validdaysi,],paste(metadatadir,"/results/part5_daysummary_", - uwi[j],"_L",uTRLi[h1],"M",uTRMi[h2],"V", - uTRVi[h3],"_",usleepparam[h4],".csv",sep=""), row.names=FALSE, na = "") + data.table::fwrite(x = OF3_clean[validdaysi, colsWithoutCosinor], + file = paste(metadatadir, "/results/part5_daysummary_", + uwi[j], "_L", uTRLi[h1], "M", uTRMi[h2], "V", + uTRVi[h3], "_", usleepparam[h4], ".csv", sep = ""), + row.names = FALSE, na = "") #------------------------------------------------------------------------------------ #also compute summary per person - agg_plainNweighted = function(df,filename="filename",daytype="daytype") { + agg_plainNweighted = function(df, filename = "filename", daytype = "daytype") { # function to take both the weighted (by weekday/weekendday) and plain average of all numeric variables # df: input data.frame (OF3 outside this function) - ignorevar = c("daysleeper","cleaningcode","night_number","sleeplog_used","ID","acc_available","window_number", + ignorevar = c("daysleeper", "cleaningcode", "night_number", "sleeplog_used", + "ID", "acc_available", "window_number", "boutcriter.mvpa", "boutcriter.lig", "boutcriter.in", "bout.metric") for (ee in 1:ncol(df)) { # make sure that numeric columns have class numeric nr = nrow(df) if (nr > 30) nr = 30 - options(warn=-1) - trynum = as.numeric(as.character(df[1:nr,ee])) - options(warn=0) + options(warn = -1) + trynum = as.numeric(as.character(df[1:nr, ee])) + options(warn = 0) if (length(which(is.na(trynum) == TRUE)) != nr & length(which(ignorevar == names(df)[ee])) == 0) { - options(warn=-1) - class(df[,ee]) = "numeric" - options(warn=0) + options(warn = -1) + class(df[, ee]) = "numeric" + options(warn = 0) } } plain_mean = function(x) { - options(warn=-1) - plain_mean = mean(x,na.rm=TRUE) - options(warn=0) + options(warn = -1) + plain_mean = mean(x, na.rm = TRUE) + options(warn = 0) if (is.na(plain_mean) == TRUE) { plain_mean = x[1] } return(plain_mean) } # aggregate across all days - PlainAggregate = aggregate.data.frame(df,by=list(df$filename),FUN=plain_mean) + PlainAggregate = aggregate.data.frame(df, by = list(df$filename), FUN = plain_mean) PlainAggregate = PlainAggregate[,-1] # aggregate per day type (weekday or weekenddays) - AggregateWDWE = aggregate.data.frame(df,by=list(df$filename,df$daytype),plain_mean) + AggregateWDWE = aggregate.data.frame(df, by = list(df$filename, df$daytype), plain_mean) AggregateWDWE = AggregateWDWE[,-c(1:2)] # Add counted number of days for Gini, Cov, alpha Fragmentation variables, because # days are dropped if there are not enough fragments: @@ -221,8 +225,8 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c if (length(vars_with_mininum_Nfrag_i) > 0) { varname_minfrag = vars_with_mininum_Nfrag[vars_with_mininum_Nfrag_i[1]] DAYCOUNT_Frag_Multiclass = aggregate.data.frame(df[,varname_minfrag], - by=list(df$filename,df$daytype), - FUN=function(x) length(which(is.na(x) == FALSE))) + by = list(df$filename,df$daytype), + FUN = function(x) length(which(is.na(x) == FALSE))) colnames(DAYCOUNT_Frag_Multiclass)[1:2] = c("filename","daytype") colnames(DAYCOUNT_Frag_Multiclass)[3] = "Nvaliddays_AL10F" # AL10F, abbreviation for: at least 10 fragments AggregateWDWE = merge(AggregateWDWE, DAYCOUNT_Frag_Multiclass, by.x = c("filename","daytype")) @@ -231,22 +235,40 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c AggregateWDWE$len <- 0 AggregateWDWE$len[which(as.character(AggregateWDWE$daytype) == "WD")] = 5 #weighting of weekdays AggregateWDWE$len[which(as.character(AggregateWDWE$daytype) == "WE")] = 2 #weighting of weekend days - dt <- data.table::as.data.table(AggregateWDWE[,which(lapply(AggregateWDWE, class)=="numeric" | - names(AggregateWDWE) == filename)]) - options(warn=-1) + dt <- + data.table::as.data.table(AggregateWDWE[, which(lapply(AggregateWDWE, class) == + "numeric" | + names(AggregateWDWE) == filename)]) + options(warn = -1) .SD <- .N <- count <- a <- NULL - WeightedAggregate <- dt[,lapply(.SD,weighted.mean,w=len,na.rm=TRUE),by=list(filename)] - options(warn=0) + WeightedAggregate <- + dt[, lapply(.SD, weighted.mean, w = len, na.rm = TRUE), by = list(filename)] + options(warn = 0) LUXmetrics = c("above1000", "timeawake", "mean", "imputed", "ignored") add_missing_LUX = function(x, LUX_day_segments, weeksegment=c(), LUXmetrics) { # missing columns, add these: NLUXseg = length(LUX_day_segments) if (length(weeksegment) > 0) { - LUX_segment_vars_expected = paste0("LUX_",LUXmetrics,"_",LUX_day_segments[1:(NLUXseg-1)],"-",LUX_day_segments[2:(NLUXseg)],"hr_day_",weeksegment) + LUX_segment_vars_expected = paste0( + "LUX_", + LUXmetrics, + "_", + LUX_day_segments[1:(NLUXseg - 1)], + "-", + LUX_day_segments[2:(NLUXseg)], + "hr_day_", + weeksegment + ) } else { - LUX_segment_vars_expected = paste0("LUX_",LUXmetrics,"_",LUX_day_segments[1:(NLUXseg-1)],"-",LUX_day_segments[2:(NLUXseg)],"hr_day") + LUX_segment_vars_expected = paste0("LUX_", + LUXmetrics, + "_", + LUX_day_segments[1:(NLUXseg - 1)], + "-", + LUX_day_segments[2:(NLUXseg)], + "hr_day") } - dummy_df = as.data.frame(matrix(NaN,1, (NLUXseg-1))) + dummy_df = as.data.frame(matrix(NaN, 1, (NLUXseg - 1))) colnames(dummy_df) = LUX_segment_vars_expected if (length(which(LUX_segment_vars_expected %in% colnames(x))) > 0) { x = as.data.frame(merge(x, dummy_df, all.x = T)) @@ -259,48 +281,73 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c } LUX_segment_vars = c() for (li in 1:length(LUXmetrics)) { - LUX_segment_vars = c(LUX_segment_vars, grep(pattern = paste0("LUX_",LUXmetrics[li]),x = colnames(WeightedAggregate), value=TRUE)) + LUX_segment_vars = c(LUX_segment_vars, + grep( + pattern = paste0("LUX_", LUXmetrics[li]), + x = colnames(WeightedAggregate), + value = TRUE + )) } if (length(LUX_segment_vars) > 0 & length(LUX_segment_vars) < 24 & length(LUX_day_segments) > 0) { - WeightedAggregate = add_missing_LUX(WeightedAggregate, LUX_day_segments, weeksegment=c(), LUXmetrics = LUXmetrics) + WeightedAggregate = add_missing_LUX( + WeightedAggregate, + LUX_day_segments, + weeksegment = c(), + LUXmetrics = LUXmetrics + ) } # merge them into one output data.frame (G) LUX_segment_vars = c() for (li in 1:length(LUXmetrics)) { - LUX_segment_vars = colnames(PlainAggregate) %in% grep(x = colnames(PlainAggregate), pattern=paste0("LUX_",LUXmetrics[li]), value=TRUE) + LUX_segment_vars = colnames(PlainAggregate) %in% grep( + x = colnames(PlainAggregate), + pattern = paste0("LUX_", LUXmetrics[li]), + value = TRUE + ) } - charcol = which(lapply(PlainAggregate, class) != "numeric" & names(PlainAggregate) != filename & !(LUX_segment_vars)) - numcol = which(lapply(PlainAggregate, class) == "numeric" | LUX_segment_vars) + charcol = which( + lapply(PlainAggregate, class) != "numeric" & + names(PlainAggregate) != filename & !(LUX_segment_vars) + ) + numcol = which(lapply(PlainAggregate, class) == "numeric" | + LUX_segment_vars) WeightedAggregate = as.data.frame(WeightedAggregate, stringsAsFactors = TRUE) - G = base::merge(PlainAggregate,WeightedAggregate,by="filename",all.x=TRUE) - p0b = paste0(names(PlainAggregate[,charcol]),".x") - p1 = paste0(names(PlainAggregate[,numcol]),".x") - p2 = paste0(names(PlainAggregate[,numcol]),".y") + G = base::merge(PlainAggregate, + WeightedAggregate, + by = "filename", + all.x = TRUE) + p0b = paste0(names(PlainAggregate[, charcol]), ".x") + p1 = paste0(names(PlainAggregate[, numcol]), ".x") + p2 = paste0(names(PlainAggregate[, numcol]), ".y") for (i in 1:length(p0b)) { - names(G)[which(names(G)==p0b[i])] = paste0(names(PlainAggregate[,charcol])[i]) + names(G)[which(names(G) == p0b[i])] = paste0(names(PlainAggregate[, charcol])[i]) } for (i in 1:length(p1)) { - names(G)[which(names(G)==p1[i])] = paste0(names(PlainAggregate[,numcol])[i],"_pla") + names(G)[which(names(G) == p1[i])] = paste0(names(PlainAggregate[, numcol])[i], "_pla") } for (i in 1:length(p2)) { - names(G)[which(names(G)==p2[i])] = paste0(names(PlainAggregate[,numcol])[i],"_wei") + names(G)[which(names(G) == p2[i])] = paste0(names(PlainAggregate[, numcol])[i], "_wei") } # expand output with weekday (WD) and weekend (WE) day aggregates for (weeksegment in c("WD", "WE")) { - temp_aggregate = AggregateWDWE[which(AggregateWDWE$daytype==weeksegment),] + temp_aggregate = AggregateWDWE[which(AggregateWDWE$daytype == weeksegment), ] charcol = which(lapply(temp_aggregate, class) != "numeric" & names(temp_aggregate) != filename) numcol = which(lapply(temp_aggregate, class) %in% c("numeric", "integer") == TRUE) names(temp_aggregate)[numcol] = paste0(names(temp_aggregate)[numcol], "_", weeksegment) temp_aggregate = temp_aggregate[,c(which(colnames(temp_aggregate) == "filename"), numcol)] LUX_segment_vars = c() for (li in 1:length(LUXmetrics)) { - LUX_segment_vars = grep(pattern = paste0("LUX_",LUXmetrics[li]),x = colnames(temp_aggregate), value=TRUE) + LUX_segment_vars = grep( + pattern = paste0("LUX_", LUXmetrics[li]), + x = colnames(temp_aggregate), + value = TRUE + ) } if (length(LUX_segment_vars) > 0 & length(LUX_segment_vars) < 24 & length(LUX_day_segments) > 0) { temp_aggregate = add_missing_LUX(temp_aggregate, LUX_day_segments, weeksegment, LUXmetrics) } G = base::merge(G, temp_aggregate, - by="filename", all.x=TRUE) + by = "filename", all.x = TRUE) } G = G[,-which(names(G) %in% c("len", "daytype", "len_WE", "len_WD"))] return(G) @@ -310,13 +357,26 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c # add column to define what are weekenddays and weekdays as needed for function agg_plainNweighted # before processing OF3, first identify which days have enough monitor wear time validdaysi = getValidDayIndices(OF3,includedaycrit.part5, excludefirstlast.part5, window = uwi[j]) - if (length(validdaysi) >0) { # do not attempt to aggregate if there are no valid days + if (length(validdaysi) > 0) { + # do not attempt to aggregate if there are no valid days # aggregate OF3 (days) to person summaries in OF4 - OF4 = agg_plainNweighted(OF3[validdaysi,],filename="filename",day="daytype") + OF4 = agg_plainNweighted(OF3[validdaysi, ], + filename = "filename", daytype = "daytype") # calculate additional variables - OF3tmp = OF3[,c("filename","night_number","daysleeper","cleaningcode","sleeplog_used","guider", - "acc_available","nonwear_perc_day","nonwear_perc_spt","daytype","dur_day_min", - "dur_spt_min")] + OF3tmp = OF3[, c( + "filename", + "night_number", + "daysleeper", + "cleaningcode", + "sleeplog_used", + "guider", + "acc_available", + "nonwear_perc_day", + "nonwear_perc_spt", + "daytype", + "dur_day_min", + "dur_spt_min" + )] foo34 = function(df,aggPerIndividual,nameold,namenew,cval) { # function to help with calculating additinal variables # related to counting how many days of measurement there are @@ -326,12 +386,20 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c # df is the non-aggregated data (days across individuals # we want to extra the number of days per individuals that meet the # criteria in df, and make it allign with aggPerIndividual. - df2 = function(x) df2 = length(which(x==cval)) # check which values meets criterion - mmm = as.data.frame(aggregate.data.frame(df,by=list(df$filename),FUN = df2), + df2 = function(x) { + df2 = length(which(x == cval)) # check which values meets criterion + } + mmm = as.data.frame(aggregate.data.frame(df, by = list(df$filename), FUN = df2), stringsAsFactors = TRUE) - mmm2 = data.frame(filename=mmm$Group.1, cc=mmm[,nameold], stringsAsFactors = TRUE) - aggPerIndividual = merge(aggPerIndividual, mmm2,by="filename") - names(aggPerIndividual)[which(names(aggPerIndividual)=="cc")] = namenew + mmm2 = data.frame( + filename = mmm$Group.1, + cc = mmm[, nameold], + stringsAsFactors = TRUE + ) + aggPerIndividual = merge(aggPerIndividual, mmm2, by = + "filename") + names(aggPerIndividual)[which(names(aggPerIndividual) == + "cc")] = namenew foo34 = aggPerIndividual } # # calculate number of valid days (both night and day criteria met) @@ -340,33 +408,77 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c OF3tmp$validdays[validdaysi] = 1 # now we have a label for the valid days, we can create a new variable # in OF4 that is a count of the number of valid days: - OF4 = foo34(df=OF3tmp,aggPerIndividual=OF4,nameold="validdays",namenew="Nvaliddays",cval=1) + OF4 = foo34(df = OF3tmp, aggPerIndividual = OF4, + nameold = "validdays", + namenew = "Nvaliddays", cval = 1) # do the same for WE (weekend days): OF3tmp$validdays = 0 OF3tmp$validdays[validdaysi[which(OF3tmp$daytype[validdaysi] == "WE")]] = 1 - OF4 = foo34(df=OF3tmp,aggPerIndividual=OF4,nameold="validdays",namenew="Nvaliddays_WE",cval=1) + OF4 = foo34( + df = OF3tmp, + aggPerIndividual = OF4, + nameold = "validdays", + namenew = "Nvaliddays_WE", + cval = 1 + ) # do the same for WD (weekdays): OF3tmp$validdays = 0 OF3tmp$validdays[validdaysi[which(OF3tmp$daytype[validdaysi] == "WD")]] = 1 - OF4 = foo34(df=OF3tmp,aggPerIndividual=OF4,nameold="validdays",namenew="Nvaliddays_WD",cval=1) # create variable from it + OF4 = foo34( + df = OF3tmp, + aggPerIndividual = OF4, + nameold = "validdays", + namenew = "Nvaliddays_WD", + cval = 1 + ) # create variable from it # do the same for daysleeper,cleaningcode, sleeplog_used, acc_available: OF3tmp$validdays = 1 - OF4 = foo34(df=OF3tmp[validdaysi,],aggPerIndividual=OF4,nameold="daysleeper",namenew="Ndaysleeper",cval=1) - OF4 = foo34(df=OF3tmp[validdaysi,],aggPerIndividual=OF4,nameold="cleaningcode",namenew="Ncleaningcodezero",cval=0) + OF4 = foo34( + df = OF3tmp[validdaysi, ], + aggPerIndividual = OF4, + nameold = "daysleeper", + namenew = "Ndaysleeper", + cval = 1 + ) + OF4 = foo34( + df = OF3tmp[validdaysi, ], + aggPerIndividual = OF4, + nameold = "cleaningcode", + namenew = "Ncleaningcodezero", + cval = 0 + ) for (ccode in 1:6) { - OF4 = foo34(df=OF3tmp[validdaysi,], aggPerIndividual=OF4, nameold="cleaningcode", - namenew=paste0("Ncleaningcode", ccode), cval=ccode) + OF4 = foo34( + df = OF3tmp[validdaysi, ], + aggPerIndividual = OF4, + nameold = "cleaningcode", + namenew = paste0("Ncleaningcode", ccode), + cval = ccode + ) } - OF4 = foo34(df=OF3tmp[validdaysi,],aggPerIndividual=OF4,nameold="sleeplog_used",namenew="Nsleeplog_used",cval=TRUE) - OF4 = foo34(df=OF3tmp[validdaysi,],aggPerIndividual=OF4,nameold="acc_available",namenew="Nacc_available",cval=1) + OF4 = foo34( + df = OF3tmp[validdaysi, ], + aggPerIndividual = OF4, + nameold = "sleeplog_used", + namenew = "Nsleeplog_used", + cval = TRUE + ) + OF4 = foo34( + df = OF3tmp[validdaysi, ], + aggPerIndividual = OF4, + nameold = "acc_available", + namenew = "Nacc_available", + cval = 1 + ) # Move valid day count variables to beginning of dataframe - OF4 = cbind(OF4[,1:5],OF4[,(ncol(OF4)-10):ncol(OF4)],OF4[,6:(ncol(OF4)-11)]) + OF4 = cbind(OF4[, 1:5], OF4[, (ncol(OF4) - 10):ncol(OF4)], OF4[, 6:(ncol(OF4) - + 11)]) nom = names(OF4) cut = which(nom == "sleeponset_ts" | nom == "wakeup_ts" | nom == "night_number" | nom == "window_number" | nom == "daysleeper" | nom == "cleaningcode" | nom == "acc_available" | nom == "guider" | nom == "L5TIME" | nom == "M5TIME" | nom == "L10TIME" | nom == "M10TIME" | nom == "acc_available" | nom == "daytype") - names(OF4)[which(names(OF4)=="weekday")] = "startday" + names(OF4)[which(names(OF4) == "weekday")] = "startday" OF4 = OF4[,-cut] for (col4 in 1:ncol(OF4)) { navalues = which(is.na(OF4[,col4]) == TRUE) @@ -376,9 +488,9 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c } #Move Nvaliddays variables to the front of the spreadsheet Nvaliddays_variables = grep(x = colnames(OF4), pattern = "Nvaliddays", value = FALSE) - Nvaliddays_variables = unique(c(which(colnames(OF4) =="Nvaliddays"), - which(colnames(OF4) =="Nvaliddays_WD"), - which(colnames(OF4) =="Nvaliddays_WE"), Nvaliddays_variables)) + Nvaliddays_variables = unique(c(which(colnames(OF4) == "Nvaliddays"), + which(colnames(OF4) == "Nvaliddays_WD"), + which(colnames(OF4) == "Nvaliddays_WE"), Nvaliddays_variables)) OF4 = OF4[,unique(c(1:4, Nvaliddays_variables, 5:ncol(OF4)))] #------------------------------------------------------------- # store all summaries in csv files From 56af54a4c80d649006d16de0c18efb7743b36867 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Mon, 19 Dec 2022 11:57:23 +0100 Subject: [PATCH 05/34] add cosinor analyses to namespace --- NAMESPACE | 2 +- R/g.part5.R | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 5a7b9da8e..043bd798d 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -35,7 +35,7 @@ export(g.analyse, g.calibrate, CalcSleepRegularityIndex, load_params, check_params, extract_params, g.imputeTimegaps, g.part5.classifyNaps, - tidyup_df, cosinorAnalyses, + tidyup_df, cosinorAnalyses, applyCosinorAnalyses, ShellDoc2Vignette, parametersVignette) importFrom("grDevices", "colors", "dev.off", "pdf","rainbow","rgb") importFrom("graphics", "abline", "axis", "par", "plot", "plot.new", diff --git a/R/g.part5.R b/R/g.part5.R index a36309705..3be91040a 100644 --- a/R/g.part5.R +++ b/R/g.part5.R @@ -947,7 +947,7 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), cosinor_coef = applyCosinorAnalyses(ts = ts[, c("time", "ACC")], qcheck = ts$nonwear, midnightsi = nightsi, - epochsizes = c(ws3, ws3)) + epochsizes = c(ws3new, ws3new)) if (length(cosinor_coef) > 0) { # assign same value to all rows to ease creating reports dsummary[,fi] = c(cosinor_coef$timeOffsetHours) @@ -1084,7 +1084,8 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), "g.part5.onsetwaketiming", "g.part5.wakesleepwindows", "g.part5.savetimeseries", "g.fragmentation", "g.intensitygradient", "g.part5.handle_lux_extremes", "g.part5.lux_persegment", "g.sibreport", - "extract_params", "load_params", "check_params") + "extract_params", "load_params", "check_params", "cosinorAnalyses", + "applyCosinorAnalyses") errhand = 'stop' } i = 0 # declare i because foreach uses it, without declaring it From a94cc013b41e39d020a07148d01b46d6567d3d76 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Mon, 19 Dec 2022 14:44:23 +0100 Subject: [PATCH 06/34] adding mean and number of daytime sibs that last at least 15 minutes --- R/g.part4.R | 78 +++++++++++++++++++++++++--------------------- R/g.report.part4.R | 45 +++++++++++++++++++------- 2 files changed, 75 insertions(+), 48 deletions(-) diff --git a/R/g.part4.R b/R/g.part4.R index d8cf3d184..123553009 100644 --- a/R/g.part4.R +++ b/R/g.part4.R @@ -65,7 +65,10 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1, 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", + "number_sib_wakinghours", "duration_sib_wakinghours_atleast15min", + "meandur_sib_wakinghours_atleast15min", + "number_sib_wakinghours_atleast15min", + "sleeponset_ts", "wakeup_ts", "guider_onset_ts", "guider_wakeup_ts", "sleeplatency", "sleepefficiency", "page", "daysleeper", "weekday", "calendar_date", "filename", "cleaningcode", "sleeplog_used", "acc_available", "guider", "SleepRegularityIndex", "SriFractionValid", "longitudinal_axis") @@ -676,7 +679,7 @@ 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 @@ -850,18 +853,19 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1, } # ====================================================================================================== sibds_atleast15min = 0 + spocum.t.dur_sibd = 0 + spocum.t.dur_sibd_atleast15min = 0 + spocum.m.dur_sibd_atleast15min = 0 + spocum.n.dur_sibd_atleast15min = 0 if (length(sibds) > 0) { spocum.t.dur_sibd = sum(sibds) atleast15min = which(sibds >= 1/4) if (length(atleast15min) > 0) { sibds_atleast15min = sibds[atleast15min] spocum.t.dur_sibd_atleast15min = sum(sibds_atleast15min) - } else { - spocum.t.dur_sibd_atleast15min = 0 - } - } else { - spocum.t.dur_sibd = 0 - spocum.t.dur_sibd_atleast15min = 0 + spocum.m.dur_sibd_atleast15min = mean(sibds_atleast15min, na.rm = TRUE) + spocum.n.dur_sibd_atleast15min = length(sibds_atleast15min) + } } nightsummary[sumi, 14] = spocum.t.dur.noc #total nocturnalsleep /accumulated sleep duration nightsummary[sumi, 15] = nightsummary[sumi, 5] - spocum.t.dur.noc #WASO @@ -870,6 +874,8 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1, nightsummary[sumi, 18] = nightsummary[sumi, 17] - 1 #number of awakenings nightsummary[sumi, 19] = length(which(spocum.t$dur == 0)) #number of sib (sustained inactivty bout) during wakinghours nightsummary[sumi, 20] = as.numeric(spocum.t.dur_sibd_atleast15min) #total sib (sustained inactivty bout) duration during wakinghours of at least 5 minutes + nightsummary[sumi, 21] = as.numeric(spocum.m.dur_sibd_atleast15min) + nightsummary[sumi, 22] = as.numeric(spocum.n.dur_sibd_atleast15min) #------------------------------------------------------- # Also report timestamps in non-numeric format: acc_onset = nightsummary[sumi, 3] @@ -880,24 +886,24 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1, # convert into clocktime acc_onsetTS = convertHRsinceprevMN2Clocktime(acc_onset) acc_wakeTS = convertHRsinceprevMN2Clocktime(acc_wake) - nightsummary[sumi, 21] = acc_onsetTS - nightsummary[sumi, 22] = acc_wakeTS + nightsummary[sumi, 23] = acc_onsetTS + nightsummary[sumi, 24] = acc_wakeTS #---------------------------------------------- - nightsummary[sumi, 23] = tmp1 #guider_onset_ts - nightsummary[sumi, 24] = tmp4 #guider_onset_ts + nightsummary[sumi, 25] = tmp1 #guider_onset_ts + nightsummary[sumi, 26] = tmp4 #guider_onset_ts if (params_sleep[["sleepwindowType"]] == "TimeInBed") { # If guider isa sleeplog and if the sleeplog recorded time in bed then # calculate: sleep latency: - nightsummary[sumi, 25] = round(nightsummary[sumi, 3] - nightsummary[sumi, 7], + nightsummary[sumi, 27] = round(nightsummary[sumi, 3] - nightsummary[sumi, 7], digits = 7) #sleeponset - guider_onset # sleep efficiency: - nightsummary[sumi, 26] = round(nightsummary[sumi, 14]/nightsummary[sumi, 9], digits = 5) #accumulated nocturnal sleep / guider + nightsummary[sumi, 28] = round(nightsummary[sumi, 14]/nightsummary[sumi, 9], digits = 5) #accumulated nocturnal sleep / guider } - nightsummary[sumi, 27] = pagei - nightsummary[sumi, 28] = daysleeper[j] - nightsummary[sumi, 29] = wdayname[j] - nightsummary[sumi, 30] = calendar_date[j] - nightsummary[sumi, 31] = fnames[i] + nightsummary[sumi, 29] = pagei + nightsummary[sumi, 30] = daysleeper[j] + nightsummary[sumi, 31] = wdayname[j] + nightsummary[sumi, 32] = calendar_date[j] + nightsummary[sumi, 33] = fnames[i] # nightsummary #------------------------------------------------------------------------ # PLOT @@ -974,12 +980,12 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1, } # PLOT #------------------------------------------------------------------------ - nightsummary[sumi, 32] = cleaningcode - nightsummary[sumi, 33] = sleeplog_used - nightsummary[sumi, 34] = acc_available - nightsummary[sumi, 35] = guider + nightsummary[sumi, 34] = cleaningcode + nightsummary[sumi, 35] = sleeplog_used + nightsummary[sumi, 36] = acc_available + nightsummary[sumi, 37] = guider # Extract SRI for this night - nightsummary[sumi, 36:37] = NA + nightsummary[sumi, 38:39] = NA if (!exists("SleepRegularityIndex")) { SleepRegularityIndex = NA } @@ -989,18 +995,18 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1, calendar_date_reformat = as.character(format(x = calendar_date_asDate, format = "%d/%m/%Y")) SRIindex = which(SRI$date == calendar_date_reformat & SRI$frac_valid > (params_cleaning[["includenightcrit"]]/24)) if (length(SRIindex) > 0) { - nightsummary[sumi, 36] = SRI$SleepRegularityIndex[SRIindex[1]] - nightsummary[sumi, 37] = SRI$frac_valid[SRIindex[1]] + nightsummary[sumi, 38] = SRI$SleepRegularityIndex[SRIindex[1]] + nightsummary[sumi, 39] = SRI$frac_valid[SRIindex[1]] } } if (length(longitudinal_axis) == 0) { - nightsummary[sumi, 38] = NA + nightsummary[sumi, 40] = NA } else { - nightsummary[sumi, 38] = longitudinal_axis + nightsummary[sumi, 40] = longitudinal_axis } if (params_output[["storefolderstructure"]] == TRUE) { - nightsummary[sumi, 39] = ffd[i] #full filename structure - nightsummary[sumi, 40] = ffp[i] #use the lowest foldername as foldername name + nightsummary[sumi, 41] = ffd[i] #full filename structure + nightsummary[sumi, 42] = ffp[i] #use the lowest foldername as foldername name } sumi = sumi + 1 } #run through definitions @@ -1033,13 +1039,13 @@ 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, 31] = fnames[i] - nightsummary[sumi, 32] = 4 #cleaningcode = 4 (no nights of accelerometer available) - nightsummary[sumi, 33:34] = c(FALSE, TRUE) #sleeplog_used acc_available - nightsummary[sumi, 35:38] = NA + nightsummary[sumi, 3:32] = NA + nightsummary[sumi, 33] = fnames[i] + nightsummary[sumi, 34] = 4 #cleaningcode = 4 (no nights of accelerometer available) + nightsummary[sumi, 35:36] = c(FALSE, TRUE) #sleeplog_used acc_available + nightsummary[sumi, 37:40] = NA if (params_output[["storefolderstructure"]] == TRUE) { - nightsummary[sumi, 39:40] = c(ffd[i], ffp[i]) #full filename structure and use the lowest foldername as foldername name + nightsummary[sumi, 41:42] = c(ffd[i], ffp[i]) #full filename structure and use the lowest foldername as foldername name } sumi = sumi + 1 } diff --git a/R/g.report.part4.R b/R/g.report.part4.R index 57ec1a919..150aec844 100644 --- a/R/g.report.part4.R +++ b/R/g.report.part4.R @@ -292,7 +292,10 @@ g.report.part4 = function(datadir = c(), metadatadir = c(), loglocation = c(), f 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", + "number_of_awakenings", "number_sib_wakinghours", + "duration_sib_wakinghours_atleast15min", + "meandur_sib_wakinghours_atleast15min", + "number_sib_wakinghours_atleast15min", "sleeplatency", "sleepefficiency", "number_of_awakenings", "guider_inbedDuration", "guider_inbedStart", "guider_inbedEnd", "guider_SptDuration", "guider_onset", "guider_wakeup", "SleepRegularityIndex", "SriFractionValid")) @@ -404,6 +407,7 @@ g.report.part4 = function(datadir = c(), metadatadir = c(), loglocation = c(), f na.rm = TRUE) personSummarynames = c(personSummarynames, paste("number_sib_wakinghours_", TW, "_", udefn[j], "_mn", sep = ""), paste("number_sib_wakinghours_", TW, "_", udefn[j], "_sd", sep = "")) + personSummary[i, (cnt + 15)] = mean(nightsummary.tmp$duration_sib_wakinghours_atleast15min[indexUdef], na.rm = TRUE) personSummary[i, (cnt + 16)] = sd(nightsummary.tmp$duration_sib_wakinghours_atleast15min[indexUdef], @@ -411,37 +415,54 @@ g.report.part4 = function(datadir = c(), metadatadir = c(), loglocation = c(), f personSummarynames = c(personSummarynames, paste("duration_sib_wakinghours_atleast15min_", TW, "_", udefn[j], "_mn", sep = ""), paste("duration_sib_wakinghours_atleast15min_", TW, "_", udefn[j], "_sd", sep = "")) + personSummary[i, (cnt + 17)] = mean(nightsummary.tmp$meandur_sib_wakinghours_atleast15min[indexUdef], + na.rm = TRUE) + personSummary[i, (cnt + 18)] = sd(nightsummary.tmp$meandur_sib_wakinghours_atleast15min[indexUdef], + na.rm = TRUE) + personSummarynames = c(personSummarynames, paste("meanduration_sib_wakinghours_atleast15min_", + TW, "_", udefn[j], "_mn", sep = ""), + paste("meanduration_sib_wakinghours_atleast15min_", + TW, "_", udefn[j], "_sd", sep = "")) + personSummary[i, (cnt + 19)] = mean(nightsummary.tmp$number_sib_wakinghours_atleast15min[indexUdef], + na.rm = TRUE) + personSummary[i, (cnt + 20)] = sd(nightsummary.tmp$number_sib_wakinghours_atleast15min[indexUdef], + na.rm = TRUE) + personSummarynames = c(personSummarynames, paste("number_sib_wakinghours_atleast15min_", + TW, "_", udefn[j], "_mn", sep = ""), + paste("number_sib_wakinghours_atleast15min_", + TW, "_", udefn[j], "_sd", sep = "")) + # average sibd during the day AVEsibdDUR = c(nightsummary.tmp$duration_sib_wakinghours[indexUdef]/nightsummary.tmp$number_sib_wakinghours[indexUdef]) if (length(which(nightsummary.tmp$number_sib_wakinghours[indexUdef] == 0))) { AVEsibdDUR[which(nightsummary.tmp$number_sib_wakinghours[indexUdef] == 0)] = 0 } - personSummary[i, (cnt + 17)] = mean(AVEsibdDUR, na.rm = TRUE) - personSummary[i, (cnt + 18)] = sd(AVEsibdDUR, na.rm = TRUE) + personSummary[i, (cnt + 21)] = mean(AVEsibdDUR, na.rm = TRUE) + personSummary[i, (cnt + 22)] = sd(AVEsibdDUR, na.rm = TRUE) personSummarynames = c(personSummarynames, paste("average_dur_sib_wakinghours_", TW, "_", 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 - personSummary[i, (cnt + 19)] = NDAYsibd + personSummary[i, (cnt + 23)] = NDAYsibd personSummarynames = c(personSummarynames, paste("n_days_w_sib_wakinghours_", TW, "_", udefn[j], sep = "")) - personSummary[i, (cnt + 20)] = mean(nightsummary.tmp$sleeponset[indexUdef], na.rm = TRUE) - personSummary[i, (cnt + 21)] = sd(nightsummary.tmp$sleeponset[indexUdef], na.rm = TRUE) + personSummary[i, (cnt + 24)] = mean(nightsummary.tmp$sleeponset[indexUdef], na.rm = TRUE) + personSummary[i, (cnt + 25)] = sd(nightsummary.tmp$sleeponset[indexUdef], na.rm = TRUE) personSummarynames = c(personSummarynames, paste("sleeponset_", TW, "_", udefn[j], "_mn", sep = ""), paste("sleeponset_", TW, "_", udefn[j], "_sd", sep = "")) - personSummary[i, (cnt + 22)] = mean(nightsummary.tmp$wakeup[indexUdef], na.rm = TRUE) - personSummary[i, (cnt + 23)] = sd(nightsummary.tmp$wakeup[indexUdef], na.rm = TRUE) + personSummary[i, (cnt + 26)] = mean(nightsummary.tmp$wakeup[indexUdef], na.rm = TRUE) + personSummary[i, (cnt + 27)] = sd(nightsummary.tmp$wakeup[indexUdef], na.rm = TRUE) personSummarynames = c(personSummarynames, paste("wakeup_", TW, "_", udefn[j], "_mn", sep = ""), paste("wakeup_", TW, "_", udefn[j], "_sd", sep = "")) - personSummary[i, (cnt + 24)] = mean(nightsummary.tmp$SleepRegularityIndex[indexUdef], + personSummary[i, (cnt + 28)] = mean(nightsummary.tmp$SleepRegularityIndex[indexUdef], na.rm = TRUE) - personSummary[i, (cnt + 25)] = sd(nightsummary.tmp$SleepRegularityIndex[indexUdef], na.rm = TRUE) + personSummary[i, (cnt + 29)] = sd(nightsummary.tmp$SleepRegularityIndex[indexUdef], na.rm = TRUE) personSummarynames = c(personSummarynames, paste("SleepRegularityIndex_", TW, "_", udefn[j], "_mn", sep = ""), paste("SleepRegularityIndex_", TW, "_", udefn[j], "_sd", sep = "")) - personSummary[i, (cnt + 26)] = mean(nightsummary.tmp$SriFractionValid[indexUdef], na.rm = TRUE) - personSummary[i, (cnt + 27)] = sd(nightsummary.tmp$SriFractionValid[indexUdef], na.rm = TRUE) + personSummary[i, (cnt + 30)] = mean(nightsummary.tmp$SriFractionValid[indexUdef], na.rm = TRUE) + personSummary[i, (cnt + 31)] = sd(nightsummary.tmp$SriFractionValid[indexUdef], na.rm = TRUE) personSummarynames = c(personSummarynames, paste("SriFractionValid_", TW, "_", udefn[j], "_mn", sep = ""), paste("SriFractionValid_", TW, "_", udefn[j], "_sd", sep = "")) cnt = cnt + 27 From d3838ea83cd9e3acf0dbdf9f4990d62af79235ef Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Mon, 19 Dec 2022 18:10:02 +0100 Subject: [PATCH 07/34] adding R2 --- R/cosinorAnalyses.R | 11 ++++++++ R/g.analyse.perfile.R | 24 +++++++++------- R/g.part5.R | 64 ++++++++++++++++++++++++------------------- R/g.report.part5.R | 2 +- 4 files changed, 62 insertions(+), 39 deletions(-) diff --git a/R/cosinorAnalyses.R b/R/cosinorAnalyses.R index eb7d93702..c66c54972 100644 --- a/R/cosinorAnalyses.R +++ b/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)) } \ No newline at end of file diff --git a/R/g.analyse.perfile.R b/R/g.analyse.perfile.R index 00459fa04..b1384e528 100644 --- a/R/g.analyse.perfile.R +++ b/R/g.analyse.perfile.R @@ -111,15 +111,16 @@ g.analyse.perfile = function(ID, fname, deviceSerialNumber, sensor.location, sta 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, @@ -128,16 +129,19 @@ g.analyse.perfile = function(ID, fname, deviceSerialNumber, sensor.location, sta 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 diff --git a/R/g.part5.R b/R/g.part5.R index 3be91040a..75a92a19d 100644 --- a/R/g.part5.R +++ b/R/g.part5.R @@ -950,43 +950,52 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), epochsizes = c(ws3new, ws3new)) if (length(cosinor_coef) > 0) { # assign same value to all rows to ease creating reports - dsummary[,fi] = c(cosinor_coef$timeOffsetHours) + dsummary[1:di, fi] = c(cosinor_coef$timeOffsetHours) ds_names[fi] = c("cosinor_timeOffsetHours") fi = fi + 1 - try(expr = {dsummary[, fi:(fi + 4)] = 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) - ds_names[fi:(fi + 4)] = c("cosinor_mes", "cosinor_amp", "cosinor_acrophase", - "cosinor_acrotime", "cosinor_ndays") - fi = fi + 5 - try(expr = {dsummary[, fi:(fi + 9)] = c(cosinor_coef$coefext$params$minimum, - cosinor_coef$coefext$params$amp, - cosinor_coef$coefext$params$alpha, - cosinor_coef$coefext$params$beta, - cosinor_coef$coefext$params$acrotime, - cosinor_coef$coefext$params$UpMesor, - cosinor_coef$coefext$params$DownMesor, - cosinor_coef$coefext$params$MESOR, - cosinor_coef$coefext$params$ndays, - cosinor_coef$coefext$params$F_pseudo)}, silent = TRUE) - ds_names[fi:(fi + 9)] = c("cosinorExt_minimum", "cosinorExt_amp", "cosinorExt_alpha", - "cosinorExt_beta", "cosinorExt_acrotime", "cosinorExt_UpMesor", - "cosinorExt_DownMesor", "cosinorExt_MESOR", - "cosinorExt_ndays", "cosinorExt_F_pseudo") - fi = fi + 10 - dsummary[fi:(fi + 1)] = c(cosinor_coef$IVIS$InterdailyStability, - cosinor_coef$IVIS$IntradailyVariability) + try(expr = {dsummary[1:di, fi:(fi + 5)] = matrix(rep(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, + cosinor_coef$coef$params$R2)), times = di), nrow = di, byrow = TRUE)}, + silent = TRUE) + + ds_names[fi:(fi + 5)] = c("cosinor_mes", "cosinor_amp", "cosinor_acrophase", + "cosinor_acrotime", "cosinor_ndays", "cosinor_R2") + fi = fi + 6 + try(expr = {dsummary[1:di, fi:(fi + 10)] = matrix(rep(c(cosinor_coef$coefext$params$minimum, + cosinor_coef$coefext$params$amp, + cosinor_coef$coefext$params$alpha, + cosinor_coef$coefext$params$beta, + cosinor_coef$coefext$params$acrotime, + cosinor_coef$coefext$params$UpMesor, + cosinor_coef$coefext$params$DownMesor, + cosinor_coef$coefext$params$MESOR, + cosinor_coef$coefext$params$ndays, + cosinor_coef$coefext$params$F_pseudo, + cosinor_coef$coefext$params$R2), times = di), nrow = di, byrow = TRUE)}, + silent = TRUE) + ds_names[fi:(fi + 10)] = c("cosinorExt_minimum", "cosinorExt_amp", "cosinorExt_alpha", + "cosinorExt_beta", "cosinorExt_acrotime", "cosinorExt_UpMesor", + "cosinorExt_DownMesor", "cosinorExt_MESOR", + "cosinorExt_ndays", "cosinorExt_F_pseudo", "cosinorExt_R2") + fi = fi + 11 + dsummary[1:di, fi:(fi + 1)] = matrix(rep(c(cosinor_coef$IVIS$InterdailyStability, + cosinor_coef$IVIS$IntradailyVariability), times = di), nrow = di) ds_names[fi:(fi + 1)] = c("cosinorIS", "cosinorIV") fi = fi + 2 + } else { + cosinor_coef = c() + fi = fi + 20 } } else { cosinor_coef = c() + fi = fi + 20 } - #=============================================================== } + } if ("angle" %in% colnames(ts)) { ts = ts[, -which(colnames(ts) == "angle")] @@ -998,7 +1007,6 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), } output = data.frame(dsummary,stringsAsFactors = FALSE) names(output) = ds_names - # correct definition of sleep log availability for window = WW, because now it # also relies on sleep log from previous night whoareWW = which(output$window == "WW") # look up WW diff --git a/R/g.report.part5.R b/R/g.report.part5.R index b1c035026..d8bedc7d5 100644 --- a/R/g.report.part5.R +++ b/R/g.report.part5.R @@ -174,7 +174,7 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c uwi[j], "_L", uTRLi[h1], "M", uTRMi[h2], "V", uTRVi[h3], "_", usleepparam[h4], ".csv", sep = ""), row.names = FALSE, na = "") # store all summaries in csv files with cleaning criteria - validdaysi = getValidDayIndices(OF3,includedaycrit.part5, excludefirstlast.part5, window = uwi[j]) + validdaysi = getValidDayIndices(OF3, includedaycrit.part5, excludefirstlast.part5, window = uwi[j]) data.table::fwrite(x = OF3_clean[validdaysi, colsWithoutCosinor], file = paste(metadatadir, "/results/part5_daysummary_", uwi[j], "_L", uTRLi[h1], "M", uTRMi[h2], "V", From 5bcb882eccf195f1c80a3ea0effa3c1b2e332bf7 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Thu, 22 Dec 2022 10:31:09 +0100 Subject: [PATCH 08/34] making sure that cosinor analyses in part5 only uses windows and is calculated seperately for MM and WW window --- R/g.part5.R | 128 ++++++++++++++++++++++++++++++---------------------- 1 file changed, 74 insertions(+), 54 deletions(-) diff --git a/R/g.part5.R b/R/g.part5.R index 75a92a19d..b7c1767df 100644 --- a/R/g.part5.R +++ b/R/g.part5.R @@ -430,7 +430,7 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), } } ts$window = 0 - + backup_cosinor_MM = backup_cosinor_WW = NULL for (TRLi in params_phyact[["threshold.lig"]]) { for (TRMi in params_phyact[["threshold.mod"]]) { for (TRVi in params_phyact[["threshold.vig"]]) { @@ -911,7 +911,78 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), } } } + + #=============================================================== + # Cosinor analyses based on only the data used for GGIR part5 + if (params_247[["cosinor"]] == TRUE) { + if ((timewindowi == "MM" & length(backup_cosinor_MM) == 0) | + (timewindowi == "WW" & length(backup_cosinor_WW) == 0)) { + # avoid computing same parameter twice because this part of the code is + # not dependent on the lig, mod, vig thresholds + cosinor_coef = applyCosinorAnalyses(ts = ts[which(ts$window != 0), c("time", "ACC")], + qcheck = ts$nonwear[which(ts$window != 0)], + midnightsi = nightsi, + epochsizes = c(ws3new, ws3new)) + if (timewindowi == "MM") { + backup_cosinor_MM = cosinor_coef + } else if (timewindowi == "WW") { + backup_cosinor_WW = cosinor_coef + } + } else { + if (timewindowi == "MM") { + cosinor_coef = backup_cosinor_MM + } else if (timewindowi == "WW") { + cosinor_coef = backup_cosinor_WW + } + } + if (length(cosinor_coef) > 0) { + # assign same value to all rows to ease creating reports + dsummary[1:di, fi] = c(cosinor_coef$timeOffsetHours) + ds_names[fi] = c("cosinor_timeOffsetHours") + fi = fi + 1 + try(expr = {dsummary[1:di, fi:(fi + 5)] = matrix(rep(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, + cosinor_coef$coef$params$R2)), times = di), nrow = di, byrow = TRUE)}, + silent = TRUE) + + ds_names[fi:(fi + 5)] = c("cosinor_mes", "cosinor_amp", "cosinor_acrophase", + "cosinor_acrotime", "cosinor_ndays", "cosinor_R2") + fi = fi + 6 + try(expr = {dsummary[1:di, fi:(fi + 10)] = matrix(rep(c(cosinor_coef$coefext$params$minimum, + cosinor_coef$coefext$params$amp, + cosinor_coef$coefext$params$alpha, + cosinor_coef$coefext$params$beta, + cosinor_coef$coefext$params$acrotime, + cosinor_coef$coefext$params$UpMesor, + cosinor_coef$coefext$params$DownMesor, + cosinor_coef$coefext$params$MESOR, + cosinor_coef$coefext$params$ndays, + cosinor_coef$coefext$params$F_pseudo, + cosinor_coef$coefext$params$R2), times = di), nrow = di, byrow = TRUE)}, + silent = TRUE) + ds_names[fi:(fi + 10)] = c("cosinorExt_minimum", "cosinorExt_amp", "cosinorExt_alpha", + "cosinorExt_beta", "cosinorExt_acrotime", "cosinorExt_UpMesor", + "cosinorExt_DownMesor", "cosinorExt_MESOR", + "cosinorExt_ndays", "cosinorExt_F_pseudo", "cosinorExt_R2") + fi = fi + 11 + dsummary[1:di, fi:(fi + 1)] = matrix(rep(c(cosinor_coef$IVIS$InterdailyStability, + cosinor_coef$IVIS$IntradailyVariability), times = di), nrow = di) + ds_names[fi:(fi + 1)] = c("cosinorIS", "cosinorIV") + fi = fi + 2 + } else { + cosinor_coef = c() + fi = fi + 20 + } + } else { + cosinor_coef = c() + fi = fi + 20 + } + } + if (params_output[["save_ms5rawlevels"]] == TRUE) { legendfile = paste0(metadatadir,ms5.outraw,"/behavioralcodes",as.Date(Sys.time()),".csv") if (file.exists(legendfile) == FALSE) { @@ -941,58 +1012,7 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), } } } - #=============================================================== - # Cosinor analyses based on only the data used for GGIR part5 - if (params_247[["cosinor"]] == TRUE) { - cosinor_coef = applyCosinorAnalyses(ts = ts[, c("time", "ACC")], - qcheck = ts$nonwear, - midnightsi = nightsi, - epochsizes = c(ws3new, ws3new)) - if (length(cosinor_coef) > 0) { - # assign same value to all rows to ease creating reports - dsummary[1:di, fi] = c(cosinor_coef$timeOffsetHours) - ds_names[fi] = c("cosinor_timeOffsetHours") - fi = fi + 1 - try(expr = {dsummary[1:di, fi:(fi + 5)] = matrix(rep(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, - cosinor_coef$coef$params$R2)), times = di), nrow = di, byrow = TRUE)}, - silent = TRUE) - - ds_names[fi:(fi + 5)] = c("cosinor_mes", "cosinor_amp", "cosinor_acrophase", - "cosinor_acrotime", "cosinor_ndays", "cosinor_R2") - fi = fi + 6 - try(expr = {dsummary[1:di, fi:(fi + 10)] = matrix(rep(c(cosinor_coef$coefext$params$minimum, - cosinor_coef$coefext$params$amp, - cosinor_coef$coefext$params$alpha, - cosinor_coef$coefext$params$beta, - cosinor_coef$coefext$params$acrotime, - cosinor_coef$coefext$params$UpMesor, - cosinor_coef$coefext$params$DownMesor, - cosinor_coef$coefext$params$MESOR, - cosinor_coef$coefext$params$ndays, - cosinor_coef$coefext$params$F_pseudo, - cosinor_coef$coefext$params$R2), times = di), nrow = di, byrow = TRUE)}, - silent = TRUE) - ds_names[fi:(fi + 10)] = c("cosinorExt_minimum", "cosinorExt_amp", "cosinorExt_alpha", - "cosinorExt_beta", "cosinorExt_acrotime", "cosinorExt_UpMesor", - "cosinorExt_DownMesor", "cosinorExt_MESOR", - "cosinorExt_ndays", "cosinorExt_F_pseudo", "cosinorExt_R2") - fi = fi + 11 - dsummary[1:di, fi:(fi + 1)] = matrix(rep(c(cosinor_coef$IVIS$InterdailyStability, - cosinor_coef$IVIS$IntradailyVariability), times = di), nrow = di) - ds_names[fi:(fi + 1)] = c("cosinorIS", "cosinorIV") - fi = fi + 2 - } else { - cosinor_coef = c() - fi = fi + 20 - } - } else { - cosinor_coef = c() - fi = fi + 20 - } + #=============================================================== } @@ -1093,7 +1113,7 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), "g.part5.savetimeseries", "g.fragmentation", "g.intensitygradient", "g.part5.handle_lux_extremes", "g.part5.lux_persegment", "g.sibreport", "extract_params", "load_params", "check_params", "cosinorAnalyses", - "applyCosinorAnalyses") + "applyCosinorAnalyses", "g.IVIS") errhand = 'stop' } i = 0 # declare i because foreach uses it, without declaring it From fc45400047b81c6837d6a311b73e7eee17e3e058 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Mon, 6 Feb 2023 13:27:52 +0100 Subject: [PATCH 09/34] remove parseGT3Xggir as function no longer used --- R/g.part1.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/g.part1.R b/R/g.part1.R index dc063a2c4..8bffcaab2 100644 --- a/R/g.part1.R +++ b/R/g.part1.R @@ -462,7 +462,7 @@ g.part1 = function(datadir = c(), outputdir = c(), f0 = 1, f1 = c(), "g.getstarttime", "POSIXtime2iso8601", "iso8601chartime2POSIX", "datadir2fnames", "read.myacc.csv", "get_nw_clip_block_params", "get_starttime_weekday_meantemp_truncdata", "ismovisens", - "g.extractheadervars", "g.imputeTimegaps", "parseGT3Xggir") + "g.extractheadervars", "g.imputeTimegaps") errhand = 'stop' # Note: This will not work for cwa files, because those also need Rcpp functions. # So, it is probably best to turn off parallel when debugging cwa data. From eee7aca9bdf503c4beda674471a599ddeed6da9f Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Mon, 6 Feb 2023 14:45:49 +0100 Subject: [PATCH 10/34] moving log transform outside ifelse statement as this always need to be applied, and adding division by 1000 because analyses expect g-units #708 --- R/applyCosinorAnalyses.R | 13 +++++-------- R/g.part5.R | 5 ++++- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/R/applyCosinorAnalyses.R b/R/applyCosinorAnalyses.R index 1d2235ff1..24dbb4a64 100644 --- a/R/applyCosinorAnalyses.R +++ b/R/applyCosinorAnalyses.R @@ -37,7 +37,7 @@ applyCosinorAnalyses = function(ts, qcheck, midnightsi, epochsizes) { 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 + # 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)]) @@ -54,18 +54,15 @@ applyCosinorAnalyses = function(ts, qcheck, midnightsi, epochsizes) { 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 } + # 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 { diff --git a/R/g.part5.R b/R/g.part5.R index b7c1767df..9e2a33f5e 100644 --- a/R/g.part5.R +++ b/R/g.part5.R @@ -919,10 +919,13 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), (timewindowi == "WW" & length(backup_cosinor_WW) == 0)) { # avoid computing same parameter twice because this part of the code is # not dependent on the lig, mod, vig thresholds - cosinor_coef = applyCosinorAnalyses(ts = ts[which(ts$window != 0), c("time", "ACC")], + acc4cos = ts[which(ts$window != 0), c("time", "ACC")] + acc4cos$ACC = acc4cos$ACC / 1000 # convert to mg because that is what applyCosinorAnalyses expects + cosinor_coef = applyCosinorAnalyses(ts = acc4cos, qcheck = ts$nonwear[which(ts$window != 0)], midnightsi = nightsi, epochsizes = c(ws3new, ws3new)) + rm(acc4cos) if (timewindowi == "MM") { backup_cosinor_MM = cosinor_coef } else if (timewindowi == "WW") { From 030a22a403acb9a898a93b4c100bef773d728d5f Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Mon, 6 Feb 2023 15:00:55 +0100 Subject: [PATCH 11/34] skip calculating weighted, WD and WE values for cosinor as estimates apply to the entire recording #708 --- NAMESPACE | 2 +- R/g.report.part5.R | 4 +++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 043bd798d..00b4b34c4 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -48,7 +48,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") diff --git a/R/g.report.part5.R b/R/g.report.part5.R index d8bedc7d5..3f97bc58a 100644 --- a/R/g.report.part5.R +++ b/R/g.report.part5.R @@ -187,7 +187,9 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c # df: input data.frame (OF3 outside this function) ignorevar = c("daysleeper", "cleaningcode", "night_number", "sleeplog_used", "ID", "acc_available", "window_number", - "boutcriter.mvpa", "boutcriter.lig", "boutcriter.in", "bout.metric") + "boutcriter.mvpa", "boutcriter.lig", "boutcriter.in", "bout.metric", + grep(pattern = "cosinor", x = names(df), ignore.case = TRUE, value = TRUE)) # skip cosinor variables + print(ignorevar) for (ee in 1:ncol(df)) { # make sure that numeric columns have class numeric nr = nrow(df) if (nr > 30) nr = 30 From caafd63712df558b78a038a2616953f69867281d Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Mon, 6 Feb 2023 16:11:20 +0100 Subject: [PATCH 12/34] fix misallignment in output columns caused by changes in this branch #708 --- R/g.report.part4.R | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/R/g.report.part4.R b/R/g.report.part4.R index 150aec844..9620552d0 100644 --- a/R/g.report.part4.R +++ b/R/g.report.part4.R @@ -407,7 +407,6 @@ g.report.part4 = function(datadir = c(), metadatadir = c(), loglocation = c(), f na.rm = TRUE) personSummarynames = c(personSummarynames, paste("number_sib_wakinghours_", TW, "_", udefn[j], "_mn", sep = ""), paste("number_sib_wakinghours_", TW, "_", udefn[j], "_sd", sep = "")) - personSummary[i, (cnt + 15)] = mean(nightsummary.tmp$duration_sib_wakinghours_atleast15min[indexUdef], na.rm = TRUE) personSummary[i, (cnt + 16)] = sd(nightsummary.tmp$duration_sib_wakinghours_atleast15min[indexUdef], @@ -443,8 +442,7 @@ g.report.part4 = function(datadir = c(), metadatadir = c(), loglocation = c(), f 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 + 23)] = NDAYsibd personSummarynames = c(personSummarynames, paste("n_days_w_sib_wakinghours_", TW, "_", udefn[j], sep = "")) @@ -465,7 +463,7 @@ g.report.part4 = function(datadir = c(), metadatadir = c(), loglocation = c(), f personSummary[i, (cnt + 31)] = sd(nightsummary.tmp$SriFractionValid[indexUdef], na.rm = TRUE) personSummarynames = c(personSummarynames, paste("SriFractionValid_", TW, "_", udefn[j], "_mn", sep = ""), paste("SriFractionValid_", TW, "_", udefn[j], "_sd", sep = "")) - cnt = cnt + 27 + cnt = cnt + 31 if (sleepwindowType == "TimeInBed") { personSummary[i, (cnt + 1)] = mean(nightsummary.tmp$sleepefficiency[indexUdef], na.rm = TRUE) personSummary[i, (cnt + 2)] = sd(nightsummary.tmp$sleepefficiency[indexUdef], na.rm = TRUE) From 91b592d7333ca62ea9e464c73ada09e1845bb810 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Tue, 13 Dec 2022 15:48:50 +0100 Subject: [PATCH 13/34] moving the application of the cosinor anlyses inside g.analyse.avday.R to a seperate function to ease re-using it in g.part5 #708 --- R/applyCosinorAnalyses.R | 75 +++++++++++++++++++ R/g.analyse.avday.R | 141 ++++++++++++++++++------------------ R/g.part5.R | 8 ++ man/applyCosinorAnalyses.Rd | 34 +++++++++ 4 files changed, 189 insertions(+), 69 deletions(-) create mode 100644 R/applyCosinorAnalyses.R create mode 100644 man/applyCosinorAnalyses.Rd diff --git a/R/applyCosinorAnalyses.R b/R/applyCosinorAnalyses.R new file mode 100644 index 000000000..6d1b7b441 --- /dev/null +++ b/R/applyCosinorAnalyses.R @@ -0,0 +1,75 @@ +applyCosinorAnalyses = function(ts, qcheck, midnightsi, shortEpoch, longEpoch) { + # qcheck - vector of length ts to indicate invalid values + ws2 = longEpoch + ws3 = shortEpoch + # 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, which(colnames(ts) %in% "timestamp" == FALSE)] + 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() + } + return(cosinor_coef) +} \ No newline at end of file diff --git a/R/g.analyse.avday.R b/R/g.analyse.avday.R index 47017fb0a..a2df794cf 100644 --- a/R/g.analyse.avday.R +++ b/R/g.analyse.avday.R @@ -120,75 +120,78 @@ 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, shortEpoch = ws3, longEpoch = ws2) + # # 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() + # } } else { cosinor_coef = c() } diff --git a/R/g.part5.R b/R/g.part5.R index cbba4be76..4b0a4b4ad 100644 --- a/R/g.part5.R +++ b/R/g.part5.R @@ -412,6 +412,14 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), } } ts$window = 0 + #=============================================================== + # Cosinor analyses based on only the data used for GGIR part5 + + + + + + #=============================================================== for (TRLi in params_phyact[["threshold.lig"]]) { for (TRMi in params_phyact[["threshold.mod"]]) { for (TRVi in params_phyact[["threshold.vig"]]) { diff --git a/man/applyCosinorAnalyses.Rd b/man/applyCosinorAnalyses.Rd new file mode 100644 index 000000000..1d733a635 --- /dev/null +++ b/man/applyCosinorAnalyses.Rd @@ -0,0 +1,34 @@ +\name{applyCosinorAnalyses} +\alias{applyCosinorAnalyses} +\title{ + Apply Cosinor Analyses to time series +} +\description{ + Wrapper function around \link{cosinorAnalyses} that first prepares the time series + before applying the cosinorAnlayses +} +\usage{ + applyCosinorAnalyses(ts, qcheck, midnightsi, shortEpoch, longEpoch) +} +\arguments{ + \item{ts}{ + Data.frame with timestamps and acceleration metric. + } + \item{qcheck}{ + Vector of equal length as number of rows in ts with value 1 for invalid + timestamps, 0 otherwise. + } + \item{midnightsi}{ + Indices for midnights in the time series + } + \item{shortEpoch}{ + First element of argument windowsizes, the resolution of the ts and qcheck values. + } + \item{longEpoch}{ + Second element of argument windowsizes, the resolution of the midnightsi + values in seconds. + } +} +\author{ + Vincent T van Hees +} \ No newline at end of file From 64323bdda41e00f2be6c24be188915687c9a01f0 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Tue, 13 Dec 2022 17:07:47 +0100 Subject: [PATCH 14/34] #708 implementing applyCosinorAnlayses in g.part5 and making sure it works with both 5 second and 60second aggregated epochs --- R/applyCosinorAnalyses.R | 8 ++--- R/g.analyse.avday.R | 71 +------------------------------------ R/g.part5.R | 39 +++++++++++++++----- man/applyCosinorAnalyses.Rd | 10 ++---- 4 files changed, 39 insertions(+), 89 deletions(-) diff --git a/R/applyCosinorAnalyses.R b/R/applyCosinorAnalyses.R index 6d1b7b441..1d2235ff1 100644 --- a/R/applyCosinorAnalyses.R +++ b/R/applyCosinorAnalyses.R @@ -1,12 +1,12 @@ -applyCosinorAnalyses = function(ts, qcheck, midnightsi, shortEpoch, longEpoch) { +applyCosinorAnalyses = function(ts, qcheck, midnightsi, epochsizes) { # qcheck - vector of length ts to indicate invalid values - ws2 = longEpoch - ws3 = shortEpoch + 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, which(colnames(ts) %in% "timestamp" == FALSE)] + 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) { diff --git a/R/g.analyse.avday.R b/R/g.analyse.avday.R index a2df794cf..c4fab41a2 100644 --- a/R/g.analyse.avday.R +++ b/R/g.analyse.avday.R @@ -122,76 +122,7 @@ g.analyse.avday = function(doquan, averageday, M, IMP, t_TWDI, quantiletype, if (params_247[["cosinor"]] == TRUE) { cosinor_coef = applyCosinorAnalyses(ts = IMP$metashort[, c("timestamp", acc.metric)], qcheck = qcheck, - midnightsi, shortEpoch = ws3, longEpoch = ws2) - # # 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() - # } + midnightsi, epochsizes = c(ws3, ws2)) } else { cosinor_coef = c() } diff --git a/R/g.part5.R b/R/g.part5.R index 4b0a4b4ad..71188880f 100644 --- a/R/g.part5.R +++ b/R/g.part5.R @@ -332,6 +332,24 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), } if (length(tail_expansion_log) != 0 & nrow(ts) > max(nightsi)) nightsi[length(nightsi) + 1] = nrow(ts) # include last window Nts = nrow(ts) + } else { + ts$time = iso8601chartime2POSIX(ts$time,tz = params_general[["desiredtz"]]) + ws3new = ws3 # change because below it is used to decide how many epochs are there in + # extract nightsi again + tempp = unclass(ts$time) + if (is.na(tempp$sec[1]) == TRUE) { + tempp = unclass(as.POSIXlt(ts$time, tz = params_general[["desiredtz"]])) + } + sec = tempp$sec + min = tempp$min + hour = tempp$hour + if (params_general[["dayborder"]] == 0) { + nightsi = which(sec == 0 & min == 0 & hour == 0) + } else { + nightsi = which(sec == 0 & min == (params_general[["dayborder"]] - floor(params_general[["dayborder"]])) * 60 & hour == floor(params_general[["dayborder"]])) #shift the definition of midnight if required + } + if (length(tail_expansion_log) != 0 & nrow(ts) > max(nightsi)) nightsi[length(nightsi) + 1] = nrow(ts) # include last window + Nts = nrow(ts) } # if ("angle" %in% colnames(ts)) { # ts = ts[, -which(colnames(ts) == "angle")] @@ -412,14 +430,7 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), } } ts$window = 0 - #=============================================================== - # Cosinor analyses based on only the data used for GGIR part5 - - - - - - #=============================================================== + for (TRLi in params_phyact[["threshold.lig"]]) { for (TRMi in params_phyact[["threshold.mod"]]) { for (TRVi in params_phyact[["threshold.vig"]]) { @@ -930,6 +941,18 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), } } } + #=============================================================== + # Cosinor analyses based on only the data used for GGIR part5 + if (params_247[["cosinor"]] == TRUE) { + cosinor_coef = applyCosinorAnalyses(ts = ts[, c("time", "ACC")], + qcheck = ts$nonwear, + midnightsi = nightsi, + epochsizes = c(ws3, ws3)) + } else { + cosinor_coef = c() + } + + #=============================================================== } } if ("angle" %in% colnames(ts)) { diff --git a/man/applyCosinorAnalyses.Rd b/man/applyCosinorAnalyses.Rd index 1d733a635..77e3198e8 100644 --- a/man/applyCosinorAnalyses.Rd +++ b/man/applyCosinorAnalyses.Rd @@ -8,7 +8,7 @@ before applying the cosinorAnlayses } \usage{ - applyCosinorAnalyses(ts, qcheck, midnightsi, shortEpoch, longEpoch) + applyCosinorAnalyses(ts, qcheck, midnightsi, epochsizes) } \arguments{ \item{ts}{ @@ -21,12 +21,8 @@ \item{midnightsi}{ Indices for midnights in the time series } - \item{shortEpoch}{ - First element of argument windowsizes, the resolution of the ts and qcheck values. - } - \item{longEpoch}{ - Second element of argument windowsizes, the resolution of the midnightsi - values in seconds. + \item{epochsizes}{ + Epoch seiz for ts and qcheck respectively } } \author{ From 003e561f060dd5c7ba12eec308d8c6e72daf4412 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Mon, 19 Dec 2022 11:27:59 +0100 Subject: [PATCH 15/34] adding cosinor variables to part 5 output --- R/g.part5.R | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/R/g.part5.R b/R/g.part5.R index 71188880f..a36309705 100644 --- a/R/g.part5.R +++ b/R/g.part5.R @@ -948,6 +948,39 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), qcheck = ts$nonwear, midnightsi = nightsi, epochsizes = c(ws3, ws3)) + if (length(cosinor_coef) > 0) { + # assign same value to all rows to ease creating reports + dsummary[,fi] = c(cosinor_coef$timeOffsetHours) + ds_names[fi] = c("cosinor_timeOffsetHours") + fi = fi + 1 + try(expr = {dsummary[, fi:(fi + 4)] = 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) + ds_names[fi:(fi + 4)] = c("cosinor_mes", "cosinor_amp", "cosinor_acrophase", + "cosinor_acrotime", "cosinor_ndays") + fi = fi + 5 + try(expr = {dsummary[, fi:(fi + 9)] = c(cosinor_coef$coefext$params$minimum, + cosinor_coef$coefext$params$amp, + cosinor_coef$coefext$params$alpha, + cosinor_coef$coefext$params$beta, + cosinor_coef$coefext$params$acrotime, + cosinor_coef$coefext$params$UpMesor, + cosinor_coef$coefext$params$DownMesor, + cosinor_coef$coefext$params$MESOR, + cosinor_coef$coefext$params$ndays, + cosinor_coef$coefext$params$F_pseudo)}, silent = TRUE) + ds_names[fi:(fi + 9)] = c("cosinorExt_minimum", "cosinorExt_amp", "cosinorExt_alpha", + "cosinorExt_beta", "cosinorExt_acrotime", "cosinorExt_UpMesor", + "cosinorExt_DownMesor", "cosinorExt_MESOR", + "cosinorExt_ndays", "cosinorExt_F_pseudo") + fi = fi + 10 + dsummary[fi:(fi + 1)] = c(cosinor_coef$IVIS$InterdailyStability, + cosinor_coef$IVIS$IntradailyVariability) + ds_names[fi:(fi + 1)] = c("cosinorIS", "cosinorIV") + fi = fi + 2 + } } else { cosinor_coef = c() } From a817df07b2659bca6bfb8743ecdcff29cce1302a Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Mon, 19 Dec 2022 11:28:18 +0100 Subject: [PATCH 16/34] tidying up syntax --- R/g.report.part5.R | 266 ++++++++++++++++++++++++++++++++------------- 1 file changed, 189 insertions(+), 77 deletions(-) diff --git a/R/g.report.part5.R b/R/g.report.part5.R index 1d5bb3cc3..b1c035026 100644 --- a/R/g.report.part5.R +++ b/R/g.report.part5.R @@ -63,8 +63,8 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c return(indices) } ms5.out = "/meta/ms5.out" - if (file.exists(paste(metadatadir,ms5.out,sep=""))) { - if (length(dir(paste(metadatadir,ms5.out,sep=""))) == 0) { + if (file.exists(paste0(metadatadir, ms5.out))) { + if (length(dir(paste0(metadatadir, ms5.out))) == 0) { try.generate.report = FALSE #do not run this function if there is no milestone data from g.part5 } else { try.generate.report = TRUE @@ -76,15 +76,15 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c if (try.generate.report == TRUE) { #====================================================================== # loop through meta-files - fnames.ms5 = list.files(paste0(metadatadir,ms5.out),full.names=TRUE) - if(f1 > length(fnames.ms5)) f1 = length(fnames.ms5) + fnames.ms5 = list.files(paste0(metadatadir,ms5.out), full.names = TRUE) + if (f1 > length(fnames.ms5)) f1 = length(fnames.ms5) cat(" loading all the milestone data from part 5 this can take a few minutes\n") myfun = function(x, expectedCols = c()) { tail_expansion_log = NULL load(file = x) cut = which(output[, 1] == "") if (length(cut) > 0 & length(cut) < nrow(output)) { - output = output[-cut,which(colnames(output) != "")] + output = output[-cut, which(colnames(output) != "")] } out = as.matrix(output) if (length(expectedCols) > 0) { @@ -109,7 +109,7 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c expectedCols = colnames(out_try) # print(fnames.ms5[f0:f1]) outputfinal = as.data.frame(do.call(rbind,lapply(fnames.ms5[f0:f1],myfun, expectedCols)), stringsAsFactors = FALSE) - cut = which(sapply(outputfinal, function(x) all(x=="")) == TRUE) # Find columns filled with missing values which(output[1,] == "" & output[2,] == "") + cut = which(sapply(outputfinal, function(x) all(x == "")) == TRUE) # Find columns filled with missing values which(output[1,] == "" & output[2,] == "") if (length(cut) > 0) { outputfinal = outputfinal[,-cut] } @@ -160,57 +160,61 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c CN = colnames(outputfinal) outputfinal2 = outputfinal colnames(outputfinal2) = CN - delcol = which(colnames(outputfinal2) == "window" | colnames(outputfinal2) == "TRLi" | - colnames(outputfinal2) == "TRMi" | colnames(outputfinal2) == "TRVi" | - colnames(outputfinal2) == "sleepparam") + delcol = grep(pattern = "window|TRLi|TRMi|TRVi|sleepparam", + x = colnames(outputfinal2)) outputfinal2 = outputfinal2[,-delcol] OF3 = outputfinal2[seluwi,] OF3 = as.data.frame(OF3, stringsAsFactors = TRUE) #------------------------------------------------------------- # store all summaries in csv files without cleaning criteria OF3_clean = tidyup_df(OF3) - data.table::fwrite(OF3_clean,paste(metadatadir,"/results/QC/part5_daysummary_full_", - uwi[j],"_L",uTRLi[h1],"M",uTRMi[h2],"V",uTRVi[h3], - "_",usleepparam[h4],".csv",sep=""),row.names=FALSE, na = "") + colsWithoutCosinor = grep(pattern = "cosinor", x = colnames(OF3_clean), invert = TRUE) + data.table::fwrite(x = OF3_clean[, colsWithoutCosinor], + file = paste(metadatadir,"/results/QC/part5_daysummary_full_", + uwi[j], "_L", uTRLi[h1], "M", uTRMi[h2], "V", uTRVi[h3], + "_", usleepparam[h4], ".csv", sep = ""), row.names = FALSE, na = "") # store all summaries in csv files with cleaning criteria validdaysi = getValidDayIndices(OF3,includedaycrit.part5, excludefirstlast.part5, window = uwi[j]) - data.table::fwrite(OF3_clean[validdaysi,],paste(metadatadir,"/results/part5_daysummary_", - uwi[j],"_L",uTRLi[h1],"M",uTRMi[h2],"V", - uTRVi[h3],"_",usleepparam[h4],".csv",sep=""), row.names=FALSE, na = "") + data.table::fwrite(x = OF3_clean[validdaysi, colsWithoutCosinor], + file = paste(metadatadir, "/results/part5_daysummary_", + uwi[j], "_L", uTRLi[h1], "M", uTRMi[h2], "V", + uTRVi[h3], "_", usleepparam[h4], ".csv", sep = ""), + row.names = FALSE, na = "") #------------------------------------------------------------------------------------ #also compute summary per person - agg_plainNweighted = function(df,filename="filename",daytype="daytype") { + agg_plainNweighted = function(df, filename = "filename", daytype = "daytype") { # function to take both the weighted (by weekday/weekendday) and plain average of all numeric variables # df: input data.frame (OF3 outside this function) - ignorevar = c("daysleeper","cleaningcode","night_number","sleeplog_used","ID","acc_available","window_number", + ignorevar = c("daysleeper", "cleaningcode", "night_number", "sleeplog_used", + "ID", "acc_available", "window_number", "boutcriter.mvpa", "boutcriter.lig", "boutcriter.in", "bout.metric") for (ee in 1:ncol(df)) { # make sure that numeric columns have class numeric nr = nrow(df) if (nr > 30) nr = 30 - options(warn=-1) - trynum = as.numeric(as.character(df[1:nr,ee])) - options(warn=0) + options(warn = -1) + trynum = as.numeric(as.character(df[1:nr, ee])) + options(warn = 0) if (length(which(is.na(trynum) == TRUE)) != nr & length(which(ignorevar == names(df)[ee])) == 0) { - options(warn=-1) - class(df[,ee]) = "numeric" - options(warn=0) + options(warn = -1) + class(df[, ee]) = "numeric" + options(warn = 0) } } plain_mean = function(x) { - options(warn=-1) - plain_mean = mean(x,na.rm=TRUE) - options(warn=0) + options(warn = -1) + plain_mean = mean(x, na.rm = TRUE) + options(warn = 0) if (is.na(plain_mean) == TRUE) { plain_mean = x[1] } return(plain_mean) } # aggregate across all days - PlainAggregate = aggregate.data.frame(df,by=list(df$filename),FUN=plain_mean) + PlainAggregate = aggregate.data.frame(df, by = list(df$filename), FUN = plain_mean) PlainAggregate = PlainAggregate[,-1] # aggregate per day type (weekday or weekenddays) - AggregateWDWE = aggregate.data.frame(df,by=list(df$filename,df$daytype),plain_mean) + AggregateWDWE = aggregate.data.frame(df, by = list(df$filename, df$daytype), plain_mean) AggregateWDWE = AggregateWDWE[,-c(1:2)] # Add counted number of days for Gini, Cov, alpha Fragmentation variables, because # days are dropped if there are not enough fragments: @@ -221,8 +225,8 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c if (length(vars_with_mininum_Nfrag_i) > 0) { varname_minfrag = vars_with_mininum_Nfrag[vars_with_mininum_Nfrag_i[1]] DAYCOUNT_Frag_Multiclass = aggregate.data.frame(df[,varname_minfrag], - by=list(df$filename,df$daytype), - FUN=function(x) length(which(is.na(x) == FALSE))) + by = list(df$filename,df$daytype), + FUN = function(x) length(which(is.na(x) == FALSE))) colnames(DAYCOUNT_Frag_Multiclass)[1:2] = c("filename","daytype") colnames(DAYCOUNT_Frag_Multiclass)[3] = "Nvaliddays_AL10F" # AL10F, abbreviation for: at least 10 fragments AggregateWDWE = merge(AggregateWDWE, DAYCOUNT_Frag_Multiclass, by.x = c("filename","daytype")) @@ -231,22 +235,40 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c AggregateWDWE$len <- 0 AggregateWDWE$len[which(as.character(AggregateWDWE$daytype) == "WD")] = 5 #weighting of weekdays AggregateWDWE$len[which(as.character(AggregateWDWE$daytype) == "WE")] = 2 #weighting of weekend days - dt <- data.table::as.data.table(AggregateWDWE[,which(lapply(AggregateWDWE, class)=="numeric" | - names(AggregateWDWE) == filename)]) - options(warn=-1) + dt <- + data.table::as.data.table(AggregateWDWE[, which(lapply(AggregateWDWE, class) == + "numeric" | + names(AggregateWDWE) == filename)]) + options(warn = -1) .SD <- .N <- count <- a <- NULL - WeightedAggregate <- dt[,lapply(.SD,weighted.mean,w=len,na.rm=TRUE),by=list(filename)] - options(warn=0) + WeightedAggregate <- + dt[, lapply(.SD, weighted.mean, w = len, na.rm = TRUE), by = list(filename)] + options(warn = 0) LUXmetrics = c("above1000", "timeawake", "mean", "imputed", "ignored") add_missing_LUX = function(x, LUX_day_segments, weeksegment=c(), LUXmetrics) { # missing columns, add these: NLUXseg = length(LUX_day_segments) if (length(weeksegment) > 0) { - LUX_segment_vars_expected = paste0("LUX_",LUXmetrics,"_",LUX_day_segments[1:(NLUXseg-1)],"-",LUX_day_segments[2:(NLUXseg)],"hr_day_",weeksegment) + LUX_segment_vars_expected = paste0( + "LUX_", + LUXmetrics, + "_", + LUX_day_segments[1:(NLUXseg - 1)], + "-", + LUX_day_segments[2:(NLUXseg)], + "hr_day_", + weeksegment + ) } else { - LUX_segment_vars_expected = paste0("LUX_",LUXmetrics,"_",LUX_day_segments[1:(NLUXseg-1)],"-",LUX_day_segments[2:(NLUXseg)],"hr_day") + LUX_segment_vars_expected = paste0("LUX_", + LUXmetrics, + "_", + LUX_day_segments[1:(NLUXseg - 1)], + "-", + LUX_day_segments[2:(NLUXseg)], + "hr_day") } - dummy_df = as.data.frame(matrix(NaN,1, (NLUXseg-1))) + dummy_df = as.data.frame(matrix(NaN, 1, (NLUXseg - 1))) colnames(dummy_df) = LUX_segment_vars_expected if (length(which(LUX_segment_vars_expected %in% colnames(x))) > 0) { x = as.data.frame(merge(x, dummy_df, all.x = T)) @@ -259,48 +281,73 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c } LUX_segment_vars = c() for (li in 1:length(LUXmetrics)) { - LUX_segment_vars = c(LUX_segment_vars, grep(pattern = paste0("LUX_",LUXmetrics[li]),x = colnames(WeightedAggregate), value=TRUE)) + LUX_segment_vars = c(LUX_segment_vars, + grep( + pattern = paste0("LUX_", LUXmetrics[li]), + x = colnames(WeightedAggregate), + value = TRUE + )) } if (length(LUX_segment_vars) > 0 & length(LUX_segment_vars) < 24 & length(LUX_day_segments) > 0) { - WeightedAggregate = add_missing_LUX(WeightedAggregate, LUX_day_segments, weeksegment=c(), LUXmetrics = LUXmetrics) + WeightedAggregate = add_missing_LUX( + WeightedAggregate, + LUX_day_segments, + weeksegment = c(), + LUXmetrics = LUXmetrics + ) } # merge them into one output data.frame (G) LUX_segment_vars = c() for (li in 1:length(LUXmetrics)) { - LUX_segment_vars = colnames(PlainAggregate) %in% grep(x = colnames(PlainAggregate), pattern=paste0("LUX_",LUXmetrics[li]), value=TRUE) + LUX_segment_vars = colnames(PlainAggregate) %in% grep( + x = colnames(PlainAggregate), + pattern = paste0("LUX_", LUXmetrics[li]), + value = TRUE + ) } - charcol = which(lapply(PlainAggregate, class) != "numeric" & names(PlainAggregate) != filename & !(LUX_segment_vars)) - numcol = which(lapply(PlainAggregate, class) == "numeric" | LUX_segment_vars) + charcol = which( + lapply(PlainAggregate, class) != "numeric" & + names(PlainAggregate) != filename & !(LUX_segment_vars) + ) + numcol = which(lapply(PlainAggregate, class) == "numeric" | + LUX_segment_vars) WeightedAggregate = as.data.frame(WeightedAggregate, stringsAsFactors = TRUE) - G = base::merge(PlainAggregate,WeightedAggregate,by="filename",all.x=TRUE) - p0b = paste0(names(PlainAggregate[,charcol]),".x") - p1 = paste0(names(PlainAggregate[,numcol]),".x") - p2 = paste0(names(PlainAggregate[,numcol]),".y") + G = base::merge(PlainAggregate, + WeightedAggregate, + by = "filename", + all.x = TRUE) + p0b = paste0(names(PlainAggregate[, charcol]), ".x") + p1 = paste0(names(PlainAggregate[, numcol]), ".x") + p2 = paste0(names(PlainAggregate[, numcol]), ".y") for (i in 1:length(p0b)) { - names(G)[which(names(G)==p0b[i])] = paste0(names(PlainAggregate[,charcol])[i]) + names(G)[which(names(G) == p0b[i])] = paste0(names(PlainAggregate[, charcol])[i]) } for (i in 1:length(p1)) { - names(G)[which(names(G)==p1[i])] = paste0(names(PlainAggregate[,numcol])[i],"_pla") + names(G)[which(names(G) == p1[i])] = paste0(names(PlainAggregate[, numcol])[i], "_pla") } for (i in 1:length(p2)) { - names(G)[which(names(G)==p2[i])] = paste0(names(PlainAggregate[,numcol])[i],"_wei") + names(G)[which(names(G) == p2[i])] = paste0(names(PlainAggregate[, numcol])[i], "_wei") } # expand output with weekday (WD) and weekend (WE) day aggregates for (weeksegment in c("WD", "WE")) { - temp_aggregate = AggregateWDWE[which(AggregateWDWE$daytype==weeksegment),] + temp_aggregate = AggregateWDWE[which(AggregateWDWE$daytype == weeksegment), ] charcol = which(lapply(temp_aggregate, class) != "numeric" & names(temp_aggregate) != filename) numcol = which(lapply(temp_aggregate, class) %in% c("numeric", "integer") == TRUE) names(temp_aggregate)[numcol] = paste0(names(temp_aggregate)[numcol], "_", weeksegment) temp_aggregate = temp_aggregate[,c(which(colnames(temp_aggregate) == "filename"), numcol)] LUX_segment_vars = c() for (li in 1:length(LUXmetrics)) { - LUX_segment_vars = grep(pattern = paste0("LUX_",LUXmetrics[li]),x = colnames(temp_aggregate), value=TRUE) + LUX_segment_vars = grep( + pattern = paste0("LUX_", LUXmetrics[li]), + x = colnames(temp_aggregate), + value = TRUE + ) } if (length(LUX_segment_vars) > 0 & length(LUX_segment_vars) < 24 & length(LUX_day_segments) > 0) { temp_aggregate = add_missing_LUX(temp_aggregate, LUX_day_segments, weeksegment, LUXmetrics) } G = base::merge(G, temp_aggregate, - by="filename", all.x=TRUE) + by = "filename", all.x = TRUE) } G = G[,-which(names(G) %in% c("len", "daytype", "len_WE", "len_WD"))] return(G) @@ -310,13 +357,26 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c # add column to define what are weekenddays and weekdays as needed for function agg_plainNweighted # before processing OF3, first identify which days have enough monitor wear time validdaysi = getValidDayIndices(OF3,includedaycrit.part5, excludefirstlast.part5, window = uwi[j]) - if (length(validdaysi) >0) { # do not attempt to aggregate if there are no valid days + if (length(validdaysi) > 0) { + # do not attempt to aggregate if there are no valid days # aggregate OF3 (days) to person summaries in OF4 - OF4 = agg_plainNweighted(OF3[validdaysi,],filename="filename",day="daytype") + OF4 = agg_plainNweighted(OF3[validdaysi, ], + filename = "filename", daytype = "daytype") # calculate additional variables - OF3tmp = OF3[,c("filename","night_number","daysleeper","cleaningcode","sleeplog_used","guider", - "acc_available","nonwear_perc_day","nonwear_perc_spt","daytype","dur_day_min", - "dur_spt_min")] + OF3tmp = OF3[, c( + "filename", + "night_number", + "daysleeper", + "cleaningcode", + "sleeplog_used", + "guider", + "acc_available", + "nonwear_perc_day", + "nonwear_perc_spt", + "daytype", + "dur_day_min", + "dur_spt_min" + )] foo34 = function(df,aggPerIndividual,nameold,namenew,cval) { # function to help with calculating additinal variables # related to counting how many days of measurement there are @@ -326,12 +386,20 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c # df is the non-aggregated data (days across individuals # we want to extra the number of days per individuals that meet the # criteria in df, and make it allign with aggPerIndividual. - df2 = function(x) df2 = length(which(x==cval)) # check which values meets criterion - mmm = as.data.frame(aggregate.data.frame(df,by=list(df$filename),FUN = df2), + df2 = function(x) { + df2 = length(which(x == cval)) # check which values meets criterion + } + mmm = as.data.frame(aggregate.data.frame(df, by = list(df$filename), FUN = df2), stringsAsFactors = TRUE) - mmm2 = data.frame(filename=mmm$Group.1, cc=mmm[,nameold], stringsAsFactors = TRUE) - aggPerIndividual = merge(aggPerIndividual, mmm2,by="filename") - names(aggPerIndividual)[which(names(aggPerIndividual)=="cc")] = namenew + mmm2 = data.frame( + filename = mmm$Group.1, + cc = mmm[, nameold], + stringsAsFactors = TRUE + ) + aggPerIndividual = merge(aggPerIndividual, mmm2, by = + "filename") + names(aggPerIndividual)[which(names(aggPerIndividual) == + "cc")] = namenew foo34 = aggPerIndividual } # # calculate number of valid days (both night and day criteria met) @@ -340,33 +408,77 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c OF3tmp$validdays[validdaysi] = 1 # now we have a label for the valid days, we can create a new variable # in OF4 that is a count of the number of valid days: - OF4 = foo34(df=OF3tmp,aggPerIndividual=OF4,nameold="validdays",namenew="Nvaliddays",cval=1) + OF4 = foo34(df = OF3tmp, aggPerIndividual = OF4, + nameold = "validdays", + namenew = "Nvaliddays", cval = 1) # do the same for WE (weekend days): OF3tmp$validdays = 0 OF3tmp$validdays[validdaysi[which(OF3tmp$daytype[validdaysi] == "WE")]] = 1 - OF4 = foo34(df=OF3tmp,aggPerIndividual=OF4,nameold="validdays",namenew="Nvaliddays_WE",cval=1) + OF4 = foo34( + df = OF3tmp, + aggPerIndividual = OF4, + nameold = "validdays", + namenew = "Nvaliddays_WE", + cval = 1 + ) # do the same for WD (weekdays): OF3tmp$validdays = 0 OF3tmp$validdays[validdaysi[which(OF3tmp$daytype[validdaysi] == "WD")]] = 1 - OF4 = foo34(df=OF3tmp,aggPerIndividual=OF4,nameold="validdays",namenew="Nvaliddays_WD",cval=1) # create variable from it + OF4 = foo34( + df = OF3tmp, + aggPerIndividual = OF4, + nameold = "validdays", + namenew = "Nvaliddays_WD", + cval = 1 + ) # create variable from it # do the same for daysleeper,cleaningcode, sleeplog_used, acc_available: OF3tmp$validdays = 1 - OF4 = foo34(df=OF3tmp[validdaysi,],aggPerIndividual=OF4,nameold="daysleeper",namenew="Ndaysleeper",cval=1) - OF4 = foo34(df=OF3tmp[validdaysi,],aggPerIndividual=OF4,nameold="cleaningcode",namenew="Ncleaningcodezero",cval=0) + OF4 = foo34( + df = OF3tmp[validdaysi, ], + aggPerIndividual = OF4, + nameold = "daysleeper", + namenew = "Ndaysleeper", + cval = 1 + ) + OF4 = foo34( + df = OF3tmp[validdaysi, ], + aggPerIndividual = OF4, + nameold = "cleaningcode", + namenew = "Ncleaningcodezero", + cval = 0 + ) for (ccode in 1:6) { - OF4 = foo34(df=OF3tmp[validdaysi,], aggPerIndividual=OF4, nameold="cleaningcode", - namenew=paste0("Ncleaningcode", ccode), cval=ccode) + OF4 = foo34( + df = OF3tmp[validdaysi, ], + aggPerIndividual = OF4, + nameold = "cleaningcode", + namenew = paste0("Ncleaningcode", ccode), + cval = ccode + ) } - OF4 = foo34(df=OF3tmp[validdaysi,],aggPerIndividual=OF4,nameold="sleeplog_used",namenew="Nsleeplog_used",cval=TRUE) - OF4 = foo34(df=OF3tmp[validdaysi,],aggPerIndividual=OF4,nameold="acc_available",namenew="Nacc_available",cval=1) + OF4 = foo34( + df = OF3tmp[validdaysi, ], + aggPerIndividual = OF4, + nameold = "sleeplog_used", + namenew = "Nsleeplog_used", + cval = TRUE + ) + OF4 = foo34( + df = OF3tmp[validdaysi, ], + aggPerIndividual = OF4, + nameold = "acc_available", + namenew = "Nacc_available", + cval = 1 + ) # Move valid day count variables to beginning of dataframe - OF4 = cbind(OF4[,1:5],OF4[,(ncol(OF4)-10):ncol(OF4)],OF4[,6:(ncol(OF4)-11)]) + OF4 = cbind(OF4[, 1:5], OF4[, (ncol(OF4) - 10):ncol(OF4)], OF4[, 6:(ncol(OF4) - + 11)]) nom = names(OF4) cut = which(nom == "sleeponset_ts" | nom == "wakeup_ts" | nom == "night_number" | nom == "window_number" | nom == "daysleeper" | nom == "cleaningcode" | nom == "acc_available" | nom == "guider" | nom == "L5TIME" | nom == "M5TIME" | nom == "L10TIME" | nom == "M10TIME" | nom == "acc_available" | nom == "daytype") - names(OF4)[which(names(OF4)=="weekday")] = "startday" + names(OF4)[which(names(OF4) == "weekday")] = "startday" OF4 = OF4[,-cut] for (col4 in 1:ncol(OF4)) { navalues = which(is.na(OF4[,col4]) == TRUE) @@ -376,9 +488,9 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c } #Move Nvaliddays variables to the front of the spreadsheet Nvaliddays_variables = grep(x = colnames(OF4), pattern = "Nvaliddays", value = FALSE) - Nvaliddays_variables = unique(c(which(colnames(OF4) =="Nvaliddays"), - which(colnames(OF4) =="Nvaliddays_WD"), - which(colnames(OF4) =="Nvaliddays_WE"), Nvaliddays_variables)) + Nvaliddays_variables = unique(c(which(colnames(OF4) == "Nvaliddays"), + which(colnames(OF4) == "Nvaliddays_WD"), + which(colnames(OF4) == "Nvaliddays_WE"), Nvaliddays_variables)) OF4 = OF4[,unique(c(1:4, Nvaliddays_variables, 5:ncol(OF4)))] #------------------------------------------------------------- # store all summaries in csv files From 6e6050ce0dafdee9255c565e92d3d51b07df1185 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Mon, 19 Dec 2022 11:57:23 +0100 Subject: [PATCH 17/34] add cosinor analyses to namespace --- NAMESPACE | 2 +- R/g.part5.R | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 5a7b9da8e..043bd798d 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -35,7 +35,7 @@ export(g.analyse, g.calibrate, CalcSleepRegularityIndex, load_params, check_params, extract_params, g.imputeTimegaps, g.part5.classifyNaps, - tidyup_df, cosinorAnalyses, + tidyup_df, cosinorAnalyses, applyCosinorAnalyses, ShellDoc2Vignette, parametersVignette) importFrom("grDevices", "colors", "dev.off", "pdf","rainbow","rgb") importFrom("graphics", "abline", "axis", "par", "plot", "plot.new", diff --git a/R/g.part5.R b/R/g.part5.R index a36309705..3be91040a 100644 --- a/R/g.part5.R +++ b/R/g.part5.R @@ -947,7 +947,7 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), cosinor_coef = applyCosinorAnalyses(ts = ts[, c("time", "ACC")], qcheck = ts$nonwear, midnightsi = nightsi, - epochsizes = c(ws3, ws3)) + epochsizes = c(ws3new, ws3new)) if (length(cosinor_coef) > 0) { # assign same value to all rows to ease creating reports dsummary[,fi] = c(cosinor_coef$timeOffsetHours) @@ -1084,7 +1084,8 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), "g.part5.onsetwaketiming", "g.part5.wakesleepwindows", "g.part5.savetimeseries", "g.fragmentation", "g.intensitygradient", "g.part5.handle_lux_extremes", "g.part5.lux_persegment", "g.sibreport", - "extract_params", "load_params", "check_params") + "extract_params", "load_params", "check_params", "cosinorAnalyses", + "applyCosinorAnalyses") errhand = 'stop' } i = 0 # declare i because foreach uses it, without declaring it From 32638bce00d1be8b4cab868fe778e18924b3bca4 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Mon, 19 Dec 2022 14:44:23 +0100 Subject: [PATCH 18/34] adding mean and number of daytime sibs that last at least 15 minutes --- R/g.part4.R | 78 +++++++++++++++++++++++++--------------------- R/g.report.part4.R | 45 +++++++++++++++++++------- 2 files changed, 75 insertions(+), 48 deletions(-) diff --git a/R/g.part4.R b/R/g.part4.R index d8cf3d184..123553009 100644 --- a/R/g.part4.R +++ b/R/g.part4.R @@ -65,7 +65,10 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1, 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", + "number_sib_wakinghours", "duration_sib_wakinghours_atleast15min", + "meandur_sib_wakinghours_atleast15min", + "number_sib_wakinghours_atleast15min", + "sleeponset_ts", "wakeup_ts", "guider_onset_ts", "guider_wakeup_ts", "sleeplatency", "sleepefficiency", "page", "daysleeper", "weekday", "calendar_date", "filename", "cleaningcode", "sleeplog_used", "acc_available", "guider", "SleepRegularityIndex", "SriFractionValid", "longitudinal_axis") @@ -676,7 +679,7 @@ 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 @@ -850,18 +853,19 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1, } # ====================================================================================================== sibds_atleast15min = 0 + spocum.t.dur_sibd = 0 + spocum.t.dur_sibd_atleast15min = 0 + spocum.m.dur_sibd_atleast15min = 0 + spocum.n.dur_sibd_atleast15min = 0 if (length(sibds) > 0) { spocum.t.dur_sibd = sum(sibds) atleast15min = which(sibds >= 1/4) if (length(atleast15min) > 0) { sibds_atleast15min = sibds[atleast15min] spocum.t.dur_sibd_atleast15min = sum(sibds_atleast15min) - } else { - spocum.t.dur_sibd_atleast15min = 0 - } - } else { - spocum.t.dur_sibd = 0 - spocum.t.dur_sibd_atleast15min = 0 + spocum.m.dur_sibd_atleast15min = mean(sibds_atleast15min, na.rm = TRUE) + spocum.n.dur_sibd_atleast15min = length(sibds_atleast15min) + } } nightsummary[sumi, 14] = spocum.t.dur.noc #total nocturnalsleep /accumulated sleep duration nightsummary[sumi, 15] = nightsummary[sumi, 5] - spocum.t.dur.noc #WASO @@ -870,6 +874,8 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1, nightsummary[sumi, 18] = nightsummary[sumi, 17] - 1 #number of awakenings nightsummary[sumi, 19] = length(which(spocum.t$dur == 0)) #number of sib (sustained inactivty bout) during wakinghours nightsummary[sumi, 20] = as.numeric(spocum.t.dur_sibd_atleast15min) #total sib (sustained inactivty bout) duration during wakinghours of at least 5 minutes + nightsummary[sumi, 21] = as.numeric(spocum.m.dur_sibd_atleast15min) + nightsummary[sumi, 22] = as.numeric(spocum.n.dur_sibd_atleast15min) #------------------------------------------------------- # Also report timestamps in non-numeric format: acc_onset = nightsummary[sumi, 3] @@ -880,24 +886,24 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1, # convert into clocktime acc_onsetTS = convertHRsinceprevMN2Clocktime(acc_onset) acc_wakeTS = convertHRsinceprevMN2Clocktime(acc_wake) - nightsummary[sumi, 21] = acc_onsetTS - nightsummary[sumi, 22] = acc_wakeTS + nightsummary[sumi, 23] = acc_onsetTS + nightsummary[sumi, 24] = acc_wakeTS #---------------------------------------------- - nightsummary[sumi, 23] = tmp1 #guider_onset_ts - nightsummary[sumi, 24] = tmp4 #guider_onset_ts + nightsummary[sumi, 25] = tmp1 #guider_onset_ts + nightsummary[sumi, 26] = tmp4 #guider_onset_ts if (params_sleep[["sleepwindowType"]] == "TimeInBed") { # If guider isa sleeplog and if the sleeplog recorded time in bed then # calculate: sleep latency: - nightsummary[sumi, 25] = round(nightsummary[sumi, 3] - nightsummary[sumi, 7], + nightsummary[sumi, 27] = round(nightsummary[sumi, 3] - nightsummary[sumi, 7], digits = 7) #sleeponset - guider_onset # sleep efficiency: - nightsummary[sumi, 26] = round(nightsummary[sumi, 14]/nightsummary[sumi, 9], digits = 5) #accumulated nocturnal sleep / guider + nightsummary[sumi, 28] = round(nightsummary[sumi, 14]/nightsummary[sumi, 9], digits = 5) #accumulated nocturnal sleep / guider } - nightsummary[sumi, 27] = pagei - nightsummary[sumi, 28] = daysleeper[j] - nightsummary[sumi, 29] = wdayname[j] - nightsummary[sumi, 30] = calendar_date[j] - nightsummary[sumi, 31] = fnames[i] + nightsummary[sumi, 29] = pagei + nightsummary[sumi, 30] = daysleeper[j] + nightsummary[sumi, 31] = wdayname[j] + nightsummary[sumi, 32] = calendar_date[j] + nightsummary[sumi, 33] = fnames[i] # nightsummary #------------------------------------------------------------------------ # PLOT @@ -974,12 +980,12 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1, } # PLOT #------------------------------------------------------------------------ - nightsummary[sumi, 32] = cleaningcode - nightsummary[sumi, 33] = sleeplog_used - nightsummary[sumi, 34] = acc_available - nightsummary[sumi, 35] = guider + nightsummary[sumi, 34] = cleaningcode + nightsummary[sumi, 35] = sleeplog_used + nightsummary[sumi, 36] = acc_available + nightsummary[sumi, 37] = guider # Extract SRI for this night - nightsummary[sumi, 36:37] = NA + nightsummary[sumi, 38:39] = NA if (!exists("SleepRegularityIndex")) { SleepRegularityIndex = NA } @@ -989,18 +995,18 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1, calendar_date_reformat = as.character(format(x = calendar_date_asDate, format = "%d/%m/%Y")) SRIindex = which(SRI$date == calendar_date_reformat & SRI$frac_valid > (params_cleaning[["includenightcrit"]]/24)) if (length(SRIindex) > 0) { - nightsummary[sumi, 36] = SRI$SleepRegularityIndex[SRIindex[1]] - nightsummary[sumi, 37] = SRI$frac_valid[SRIindex[1]] + nightsummary[sumi, 38] = SRI$SleepRegularityIndex[SRIindex[1]] + nightsummary[sumi, 39] = SRI$frac_valid[SRIindex[1]] } } if (length(longitudinal_axis) == 0) { - nightsummary[sumi, 38] = NA + nightsummary[sumi, 40] = NA } else { - nightsummary[sumi, 38] = longitudinal_axis + nightsummary[sumi, 40] = longitudinal_axis } if (params_output[["storefolderstructure"]] == TRUE) { - nightsummary[sumi, 39] = ffd[i] #full filename structure - nightsummary[sumi, 40] = ffp[i] #use the lowest foldername as foldername name + nightsummary[sumi, 41] = ffd[i] #full filename structure + nightsummary[sumi, 42] = ffp[i] #use the lowest foldername as foldername name } sumi = sumi + 1 } #run through definitions @@ -1033,13 +1039,13 @@ 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, 31] = fnames[i] - nightsummary[sumi, 32] = 4 #cleaningcode = 4 (no nights of accelerometer available) - nightsummary[sumi, 33:34] = c(FALSE, TRUE) #sleeplog_used acc_available - nightsummary[sumi, 35:38] = NA + nightsummary[sumi, 3:32] = NA + nightsummary[sumi, 33] = fnames[i] + nightsummary[sumi, 34] = 4 #cleaningcode = 4 (no nights of accelerometer available) + nightsummary[sumi, 35:36] = c(FALSE, TRUE) #sleeplog_used acc_available + nightsummary[sumi, 37:40] = NA if (params_output[["storefolderstructure"]] == TRUE) { - nightsummary[sumi, 39:40] = c(ffd[i], ffp[i]) #full filename structure and use the lowest foldername as foldername name + nightsummary[sumi, 41:42] = c(ffd[i], ffp[i]) #full filename structure and use the lowest foldername as foldername name } sumi = sumi + 1 } diff --git a/R/g.report.part4.R b/R/g.report.part4.R index 57ec1a919..150aec844 100644 --- a/R/g.report.part4.R +++ b/R/g.report.part4.R @@ -292,7 +292,10 @@ g.report.part4 = function(datadir = c(), metadatadir = c(), loglocation = c(), f 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", + "number_of_awakenings", "number_sib_wakinghours", + "duration_sib_wakinghours_atleast15min", + "meandur_sib_wakinghours_atleast15min", + "number_sib_wakinghours_atleast15min", "sleeplatency", "sleepefficiency", "number_of_awakenings", "guider_inbedDuration", "guider_inbedStart", "guider_inbedEnd", "guider_SptDuration", "guider_onset", "guider_wakeup", "SleepRegularityIndex", "SriFractionValid")) @@ -404,6 +407,7 @@ g.report.part4 = function(datadir = c(), metadatadir = c(), loglocation = c(), f na.rm = TRUE) personSummarynames = c(personSummarynames, paste("number_sib_wakinghours_", TW, "_", udefn[j], "_mn", sep = ""), paste("number_sib_wakinghours_", TW, "_", udefn[j], "_sd", sep = "")) + personSummary[i, (cnt + 15)] = mean(nightsummary.tmp$duration_sib_wakinghours_atleast15min[indexUdef], na.rm = TRUE) personSummary[i, (cnt + 16)] = sd(nightsummary.tmp$duration_sib_wakinghours_atleast15min[indexUdef], @@ -411,37 +415,54 @@ g.report.part4 = function(datadir = c(), metadatadir = c(), loglocation = c(), f personSummarynames = c(personSummarynames, paste("duration_sib_wakinghours_atleast15min_", TW, "_", udefn[j], "_mn", sep = ""), paste("duration_sib_wakinghours_atleast15min_", TW, "_", udefn[j], "_sd", sep = "")) + personSummary[i, (cnt + 17)] = mean(nightsummary.tmp$meandur_sib_wakinghours_atleast15min[indexUdef], + na.rm = TRUE) + personSummary[i, (cnt + 18)] = sd(nightsummary.tmp$meandur_sib_wakinghours_atleast15min[indexUdef], + na.rm = TRUE) + personSummarynames = c(personSummarynames, paste("meanduration_sib_wakinghours_atleast15min_", + TW, "_", udefn[j], "_mn", sep = ""), + paste("meanduration_sib_wakinghours_atleast15min_", + TW, "_", udefn[j], "_sd", sep = "")) + personSummary[i, (cnt + 19)] = mean(nightsummary.tmp$number_sib_wakinghours_atleast15min[indexUdef], + na.rm = TRUE) + personSummary[i, (cnt + 20)] = sd(nightsummary.tmp$number_sib_wakinghours_atleast15min[indexUdef], + na.rm = TRUE) + personSummarynames = c(personSummarynames, paste("number_sib_wakinghours_atleast15min_", + TW, "_", udefn[j], "_mn", sep = ""), + paste("number_sib_wakinghours_atleast15min_", + TW, "_", udefn[j], "_sd", sep = "")) + # average sibd during the day AVEsibdDUR = c(nightsummary.tmp$duration_sib_wakinghours[indexUdef]/nightsummary.tmp$number_sib_wakinghours[indexUdef]) if (length(which(nightsummary.tmp$number_sib_wakinghours[indexUdef] == 0))) { AVEsibdDUR[which(nightsummary.tmp$number_sib_wakinghours[indexUdef] == 0)] = 0 } - personSummary[i, (cnt + 17)] = mean(AVEsibdDUR, na.rm = TRUE) - personSummary[i, (cnt + 18)] = sd(AVEsibdDUR, na.rm = TRUE) + personSummary[i, (cnt + 21)] = mean(AVEsibdDUR, na.rm = TRUE) + personSummary[i, (cnt + 22)] = sd(AVEsibdDUR, na.rm = TRUE) personSummarynames = c(personSummarynames, paste("average_dur_sib_wakinghours_", TW, "_", 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 - personSummary[i, (cnt + 19)] = NDAYsibd + personSummary[i, (cnt + 23)] = NDAYsibd personSummarynames = c(personSummarynames, paste("n_days_w_sib_wakinghours_", TW, "_", udefn[j], sep = "")) - personSummary[i, (cnt + 20)] = mean(nightsummary.tmp$sleeponset[indexUdef], na.rm = TRUE) - personSummary[i, (cnt + 21)] = sd(nightsummary.tmp$sleeponset[indexUdef], na.rm = TRUE) + personSummary[i, (cnt + 24)] = mean(nightsummary.tmp$sleeponset[indexUdef], na.rm = TRUE) + personSummary[i, (cnt + 25)] = sd(nightsummary.tmp$sleeponset[indexUdef], na.rm = TRUE) personSummarynames = c(personSummarynames, paste("sleeponset_", TW, "_", udefn[j], "_mn", sep = ""), paste("sleeponset_", TW, "_", udefn[j], "_sd", sep = "")) - personSummary[i, (cnt + 22)] = mean(nightsummary.tmp$wakeup[indexUdef], na.rm = TRUE) - personSummary[i, (cnt + 23)] = sd(nightsummary.tmp$wakeup[indexUdef], na.rm = TRUE) + personSummary[i, (cnt + 26)] = mean(nightsummary.tmp$wakeup[indexUdef], na.rm = TRUE) + personSummary[i, (cnt + 27)] = sd(nightsummary.tmp$wakeup[indexUdef], na.rm = TRUE) personSummarynames = c(personSummarynames, paste("wakeup_", TW, "_", udefn[j], "_mn", sep = ""), paste("wakeup_", TW, "_", udefn[j], "_sd", sep = "")) - personSummary[i, (cnt + 24)] = mean(nightsummary.tmp$SleepRegularityIndex[indexUdef], + personSummary[i, (cnt + 28)] = mean(nightsummary.tmp$SleepRegularityIndex[indexUdef], na.rm = TRUE) - personSummary[i, (cnt + 25)] = sd(nightsummary.tmp$SleepRegularityIndex[indexUdef], na.rm = TRUE) + personSummary[i, (cnt + 29)] = sd(nightsummary.tmp$SleepRegularityIndex[indexUdef], na.rm = TRUE) personSummarynames = c(personSummarynames, paste("SleepRegularityIndex_", TW, "_", udefn[j], "_mn", sep = ""), paste("SleepRegularityIndex_", TW, "_", udefn[j], "_sd", sep = "")) - personSummary[i, (cnt + 26)] = mean(nightsummary.tmp$SriFractionValid[indexUdef], na.rm = TRUE) - personSummary[i, (cnt + 27)] = sd(nightsummary.tmp$SriFractionValid[indexUdef], na.rm = TRUE) + personSummary[i, (cnt + 30)] = mean(nightsummary.tmp$SriFractionValid[indexUdef], na.rm = TRUE) + personSummary[i, (cnt + 31)] = sd(nightsummary.tmp$SriFractionValid[indexUdef], na.rm = TRUE) personSummarynames = c(personSummarynames, paste("SriFractionValid_", TW, "_", udefn[j], "_mn", sep = ""), paste("SriFractionValid_", TW, "_", udefn[j], "_sd", sep = "")) cnt = cnt + 27 From 3308d4eeec6c50813b6936e05a9d870dafd0e663 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Mon, 19 Dec 2022 18:10:02 +0100 Subject: [PATCH 19/34] adding R2 --- R/cosinorAnalyses.R | 11 ++++++++ R/g.analyse.perfile.R | 24 +++++++++------- R/g.part5.R | 64 ++++++++++++++++++++++++------------------- R/g.report.part5.R | 2 +- 4 files changed, 62 insertions(+), 39 deletions(-) diff --git a/R/cosinorAnalyses.R b/R/cosinorAnalyses.R index eb7d93702..c66c54972 100644 --- a/R/cosinorAnalyses.R +++ b/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)) } \ No newline at end of file diff --git a/R/g.analyse.perfile.R b/R/g.analyse.perfile.R index 00459fa04..b1384e528 100644 --- a/R/g.analyse.perfile.R +++ b/R/g.analyse.perfile.R @@ -111,15 +111,16 @@ g.analyse.perfile = function(ID, fname, deviceSerialNumber, sensor.location, sta 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, @@ -128,16 +129,19 @@ g.analyse.perfile = function(ID, fname, deviceSerialNumber, sensor.location, sta 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 diff --git a/R/g.part5.R b/R/g.part5.R index 3be91040a..75a92a19d 100644 --- a/R/g.part5.R +++ b/R/g.part5.R @@ -950,43 +950,52 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), epochsizes = c(ws3new, ws3new)) if (length(cosinor_coef) > 0) { # assign same value to all rows to ease creating reports - dsummary[,fi] = c(cosinor_coef$timeOffsetHours) + dsummary[1:di, fi] = c(cosinor_coef$timeOffsetHours) ds_names[fi] = c("cosinor_timeOffsetHours") fi = fi + 1 - try(expr = {dsummary[, fi:(fi + 4)] = 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) - ds_names[fi:(fi + 4)] = c("cosinor_mes", "cosinor_amp", "cosinor_acrophase", - "cosinor_acrotime", "cosinor_ndays") - fi = fi + 5 - try(expr = {dsummary[, fi:(fi + 9)] = c(cosinor_coef$coefext$params$minimum, - cosinor_coef$coefext$params$amp, - cosinor_coef$coefext$params$alpha, - cosinor_coef$coefext$params$beta, - cosinor_coef$coefext$params$acrotime, - cosinor_coef$coefext$params$UpMesor, - cosinor_coef$coefext$params$DownMesor, - cosinor_coef$coefext$params$MESOR, - cosinor_coef$coefext$params$ndays, - cosinor_coef$coefext$params$F_pseudo)}, silent = TRUE) - ds_names[fi:(fi + 9)] = c("cosinorExt_minimum", "cosinorExt_amp", "cosinorExt_alpha", - "cosinorExt_beta", "cosinorExt_acrotime", "cosinorExt_UpMesor", - "cosinorExt_DownMesor", "cosinorExt_MESOR", - "cosinorExt_ndays", "cosinorExt_F_pseudo") - fi = fi + 10 - dsummary[fi:(fi + 1)] = c(cosinor_coef$IVIS$InterdailyStability, - cosinor_coef$IVIS$IntradailyVariability) + try(expr = {dsummary[1:di, fi:(fi + 5)] = matrix(rep(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, + cosinor_coef$coef$params$R2)), times = di), nrow = di, byrow = TRUE)}, + silent = TRUE) + + ds_names[fi:(fi + 5)] = c("cosinor_mes", "cosinor_amp", "cosinor_acrophase", + "cosinor_acrotime", "cosinor_ndays", "cosinor_R2") + fi = fi + 6 + try(expr = {dsummary[1:di, fi:(fi + 10)] = matrix(rep(c(cosinor_coef$coefext$params$minimum, + cosinor_coef$coefext$params$amp, + cosinor_coef$coefext$params$alpha, + cosinor_coef$coefext$params$beta, + cosinor_coef$coefext$params$acrotime, + cosinor_coef$coefext$params$UpMesor, + cosinor_coef$coefext$params$DownMesor, + cosinor_coef$coefext$params$MESOR, + cosinor_coef$coefext$params$ndays, + cosinor_coef$coefext$params$F_pseudo, + cosinor_coef$coefext$params$R2), times = di), nrow = di, byrow = TRUE)}, + silent = TRUE) + ds_names[fi:(fi + 10)] = c("cosinorExt_minimum", "cosinorExt_amp", "cosinorExt_alpha", + "cosinorExt_beta", "cosinorExt_acrotime", "cosinorExt_UpMesor", + "cosinorExt_DownMesor", "cosinorExt_MESOR", + "cosinorExt_ndays", "cosinorExt_F_pseudo", "cosinorExt_R2") + fi = fi + 11 + dsummary[1:di, fi:(fi + 1)] = matrix(rep(c(cosinor_coef$IVIS$InterdailyStability, + cosinor_coef$IVIS$IntradailyVariability), times = di), nrow = di) ds_names[fi:(fi + 1)] = c("cosinorIS", "cosinorIV") fi = fi + 2 + } else { + cosinor_coef = c() + fi = fi + 20 } } else { cosinor_coef = c() + fi = fi + 20 } - #=============================================================== } + } if ("angle" %in% colnames(ts)) { ts = ts[, -which(colnames(ts) == "angle")] @@ -998,7 +1007,6 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), } output = data.frame(dsummary,stringsAsFactors = FALSE) names(output) = ds_names - # correct definition of sleep log availability for window = WW, because now it # also relies on sleep log from previous night whoareWW = which(output$window == "WW") # look up WW diff --git a/R/g.report.part5.R b/R/g.report.part5.R index b1c035026..d8bedc7d5 100644 --- a/R/g.report.part5.R +++ b/R/g.report.part5.R @@ -174,7 +174,7 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c uwi[j], "_L", uTRLi[h1], "M", uTRMi[h2], "V", uTRVi[h3], "_", usleepparam[h4], ".csv", sep = ""), row.names = FALSE, na = "") # store all summaries in csv files with cleaning criteria - validdaysi = getValidDayIndices(OF3,includedaycrit.part5, excludefirstlast.part5, window = uwi[j]) + validdaysi = getValidDayIndices(OF3, includedaycrit.part5, excludefirstlast.part5, window = uwi[j]) data.table::fwrite(x = OF3_clean[validdaysi, colsWithoutCosinor], file = paste(metadatadir, "/results/part5_daysummary_", uwi[j], "_L", uTRLi[h1], "M", uTRMi[h2], "V", From 25cb509282e017bad5033c92b50deee86c0ccfe8 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Thu, 22 Dec 2022 10:31:09 +0100 Subject: [PATCH 20/34] making sure that cosinor analyses in part5 only uses windows and is calculated seperately for MM and WW window --- R/g.part5.R | 128 ++++++++++++++++++++++++++++++---------------------- 1 file changed, 74 insertions(+), 54 deletions(-) diff --git a/R/g.part5.R b/R/g.part5.R index 75a92a19d..b7c1767df 100644 --- a/R/g.part5.R +++ b/R/g.part5.R @@ -430,7 +430,7 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), } } ts$window = 0 - + backup_cosinor_MM = backup_cosinor_WW = NULL for (TRLi in params_phyact[["threshold.lig"]]) { for (TRMi in params_phyact[["threshold.mod"]]) { for (TRVi in params_phyact[["threshold.vig"]]) { @@ -911,7 +911,78 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), } } } + + #=============================================================== + # Cosinor analyses based on only the data used for GGIR part5 + if (params_247[["cosinor"]] == TRUE) { + if ((timewindowi == "MM" & length(backup_cosinor_MM) == 0) | + (timewindowi == "WW" & length(backup_cosinor_WW) == 0)) { + # avoid computing same parameter twice because this part of the code is + # not dependent on the lig, mod, vig thresholds + cosinor_coef = applyCosinorAnalyses(ts = ts[which(ts$window != 0), c("time", "ACC")], + qcheck = ts$nonwear[which(ts$window != 0)], + midnightsi = nightsi, + epochsizes = c(ws3new, ws3new)) + if (timewindowi == "MM") { + backup_cosinor_MM = cosinor_coef + } else if (timewindowi == "WW") { + backup_cosinor_WW = cosinor_coef + } + } else { + if (timewindowi == "MM") { + cosinor_coef = backup_cosinor_MM + } else if (timewindowi == "WW") { + cosinor_coef = backup_cosinor_WW + } + } + if (length(cosinor_coef) > 0) { + # assign same value to all rows to ease creating reports + dsummary[1:di, fi] = c(cosinor_coef$timeOffsetHours) + ds_names[fi] = c("cosinor_timeOffsetHours") + fi = fi + 1 + try(expr = {dsummary[1:di, fi:(fi + 5)] = matrix(rep(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, + cosinor_coef$coef$params$R2)), times = di), nrow = di, byrow = TRUE)}, + silent = TRUE) + + ds_names[fi:(fi + 5)] = c("cosinor_mes", "cosinor_amp", "cosinor_acrophase", + "cosinor_acrotime", "cosinor_ndays", "cosinor_R2") + fi = fi + 6 + try(expr = {dsummary[1:di, fi:(fi + 10)] = matrix(rep(c(cosinor_coef$coefext$params$minimum, + cosinor_coef$coefext$params$amp, + cosinor_coef$coefext$params$alpha, + cosinor_coef$coefext$params$beta, + cosinor_coef$coefext$params$acrotime, + cosinor_coef$coefext$params$UpMesor, + cosinor_coef$coefext$params$DownMesor, + cosinor_coef$coefext$params$MESOR, + cosinor_coef$coefext$params$ndays, + cosinor_coef$coefext$params$F_pseudo, + cosinor_coef$coefext$params$R2), times = di), nrow = di, byrow = TRUE)}, + silent = TRUE) + ds_names[fi:(fi + 10)] = c("cosinorExt_minimum", "cosinorExt_amp", "cosinorExt_alpha", + "cosinorExt_beta", "cosinorExt_acrotime", "cosinorExt_UpMesor", + "cosinorExt_DownMesor", "cosinorExt_MESOR", + "cosinorExt_ndays", "cosinorExt_F_pseudo", "cosinorExt_R2") + fi = fi + 11 + dsummary[1:di, fi:(fi + 1)] = matrix(rep(c(cosinor_coef$IVIS$InterdailyStability, + cosinor_coef$IVIS$IntradailyVariability), times = di), nrow = di) + ds_names[fi:(fi + 1)] = c("cosinorIS", "cosinorIV") + fi = fi + 2 + } else { + cosinor_coef = c() + fi = fi + 20 + } + } else { + cosinor_coef = c() + fi = fi + 20 + } + } + if (params_output[["save_ms5rawlevels"]] == TRUE) { legendfile = paste0(metadatadir,ms5.outraw,"/behavioralcodes",as.Date(Sys.time()),".csv") if (file.exists(legendfile) == FALSE) { @@ -941,58 +1012,7 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), } } } - #=============================================================== - # Cosinor analyses based on only the data used for GGIR part5 - if (params_247[["cosinor"]] == TRUE) { - cosinor_coef = applyCosinorAnalyses(ts = ts[, c("time", "ACC")], - qcheck = ts$nonwear, - midnightsi = nightsi, - epochsizes = c(ws3new, ws3new)) - if (length(cosinor_coef) > 0) { - # assign same value to all rows to ease creating reports - dsummary[1:di, fi] = c(cosinor_coef$timeOffsetHours) - ds_names[fi] = c("cosinor_timeOffsetHours") - fi = fi + 1 - try(expr = {dsummary[1:di, fi:(fi + 5)] = matrix(rep(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, - cosinor_coef$coef$params$R2)), times = di), nrow = di, byrow = TRUE)}, - silent = TRUE) - - ds_names[fi:(fi + 5)] = c("cosinor_mes", "cosinor_amp", "cosinor_acrophase", - "cosinor_acrotime", "cosinor_ndays", "cosinor_R2") - fi = fi + 6 - try(expr = {dsummary[1:di, fi:(fi + 10)] = matrix(rep(c(cosinor_coef$coefext$params$minimum, - cosinor_coef$coefext$params$amp, - cosinor_coef$coefext$params$alpha, - cosinor_coef$coefext$params$beta, - cosinor_coef$coefext$params$acrotime, - cosinor_coef$coefext$params$UpMesor, - cosinor_coef$coefext$params$DownMesor, - cosinor_coef$coefext$params$MESOR, - cosinor_coef$coefext$params$ndays, - cosinor_coef$coefext$params$F_pseudo, - cosinor_coef$coefext$params$R2), times = di), nrow = di, byrow = TRUE)}, - silent = TRUE) - ds_names[fi:(fi + 10)] = c("cosinorExt_minimum", "cosinorExt_amp", "cosinorExt_alpha", - "cosinorExt_beta", "cosinorExt_acrotime", "cosinorExt_UpMesor", - "cosinorExt_DownMesor", "cosinorExt_MESOR", - "cosinorExt_ndays", "cosinorExt_F_pseudo", "cosinorExt_R2") - fi = fi + 11 - dsummary[1:di, fi:(fi + 1)] = matrix(rep(c(cosinor_coef$IVIS$InterdailyStability, - cosinor_coef$IVIS$IntradailyVariability), times = di), nrow = di) - ds_names[fi:(fi + 1)] = c("cosinorIS", "cosinorIV") - fi = fi + 2 - } else { - cosinor_coef = c() - fi = fi + 20 - } - } else { - cosinor_coef = c() - fi = fi + 20 - } + #=============================================================== } @@ -1093,7 +1113,7 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), "g.part5.savetimeseries", "g.fragmentation", "g.intensitygradient", "g.part5.handle_lux_extremes", "g.part5.lux_persegment", "g.sibreport", "extract_params", "load_params", "check_params", "cosinorAnalyses", - "applyCosinorAnalyses") + "applyCosinorAnalyses", "g.IVIS") errhand = 'stop' } i = 0 # declare i because foreach uses it, without declaring it From e18a5c15c632b825b30838af78416aea9bc8e9e1 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Mon, 6 Feb 2023 13:27:52 +0100 Subject: [PATCH 21/34] remove parseGT3Xggir as function no longer used --- R/g.part1.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/g.part1.R b/R/g.part1.R index dc063a2c4..8bffcaab2 100644 --- a/R/g.part1.R +++ b/R/g.part1.R @@ -462,7 +462,7 @@ g.part1 = function(datadir = c(), outputdir = c(), f0 = 1, f1 = c(), "g.getstarttime", "POSIXtime2iso8601", "iso8601chartime2POSIX", "datadir2fnames", "read.myacc.csv", "get_nw_clip_block_params", "get_starttime_weekday_meantemp_truncdata", "ismovisens", - "g.extractheadervars", "g.imputeTimegaps", "parseGT3Xggir") + "g.extractheadervars", "g.imputeTimegaps") errhand = 'stop' # Note: This will not work for cwa files, because those also need Rcpp functions. # So, it is probably best to turn off parallel when debugging cwa data. From be3f73a662e88b3fa7c640e066d12d413513a2cf Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Mon, 6 Feb 2023 14:45:49 +0100 Subject: [PATCH 22/34] moving log transform outside ifelse statement as this always need to be applied, and adding division by 1000 because analyses expect g-units #708 --- R/applyCosinorAnalyses.R | 13 +++++-------- R/g.part5.R | 5 ++++- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/R/applyCosinorAnalyses.R b/R/applyCosinorAnalyses.R index 1d2235ff1..24dbb4a64 100644 --- a/R/applyCosinorAnalyses.R +++ b/R/applyCosinorAnalyses.R @@ -37,7 +37,7 @@ applyCosinorAnalyses = function(ts, qcheck, midnightsi, epochsizes) { 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 + # 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)]) @@ -54,18 +54,15 @@ applyCosinorAnalyses = function(ts, qcheck, midnightsi, epochsizes) { 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 } + # 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 { diff --git a/R/g.part5.R b/R/g.part5.R index b7c1767df..9e2a33f5e 100644 --- a/R/g.part5.R +++ b/R/g.part5.R @@ -919,10 +919,13 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), (timewindowi == "WW" & length(backup_cosinor_WW) == 0)) { # avoid computing same parameter twice because this part of the code is # not dependent on the lig, mod, vig thresholds - cosinor_coef = applyCosinorAnalyses(ts = ts[which(ts$window != 0), c("time", "ACC")], + acc4cos = ts[which(ts$window != 0), c("time", "ACC")] + acc4cos$ACC = acc4cos$ACC / 1000 # convert to mg because that is what applyCosinorAnalyses expects + cosinor_coef = applyCosinorAnalyses(ts = acc4cos, qcheck = ts$nonwear[which(ts$window != 0)], midnightsi = nightsi, epochsizes = c(ws3new, ws3new)) + rm(acc4cos) if (timewindowi == "MM") { backup_cosinor_MM = cosinor_coef } else if (timewindowi == "WW") { From 3a5da8f26e260d3ab39dc1b6dd82497ce5406e83 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Mon, 6 Feb 2023 15:00:55 +0100 Subject: [PATCH 23/34] skip calculating weighted, WD and WE values for cosinor as estimates apply to the entire recording #708 --- NAMESPACE | 2 +- R/g.report.part5.R | 4 +++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 043bd798d..00b4b34c4 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -48,7 +48,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") diff --git a/R/g.report.part5.R b/R/g.report.part5.R index d8bedc7d5..3f97bc58a 100644 --- a/R/g.report.part5.R +++ b/R/g.report.part5.R @@ -187,7 +187,9 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c # df: input data.frame (OF3 outside this function) ignorevar = c("daysleeper", "cleaningcode", "night_number", "sleeplog_used", "ID", "acc_available", "window_number", - "boutcriter.mvpa", "boutcriter.lig", "boutcriter.in", "bout.metric") + "boutcriter.mvpa", "boutcriter.lig", "boutcriter.in", "bout.metric", + grep(pattern = "cosinor", x = names(df), ignore.case = TRUE, value = TRUE)) # skip cosinor variables + print(ignorevar) for (ee in 1:ncol(df)) { # make sure that numeric columns have class numeric nr = nrow(df) if (nr > 30) nr = 30 From f5ea91c01fdbf38b013793f003d13bf88236e285 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Mon, 6 Feb 2023 16:11:20 +0100 Subject: [PATCH 24/34] fix misallignment in output columns caused by changes in this branch #708 --- R/g.report.part4.R | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/R/g.report.part4.R b/R/g.report.part4.R index 150aec844..9620552d0 100644 --- a/R/g.report.part4.R +++ b/R/g.report.part4.R @@ -407,7 +407,6 @@ g.report.part4 = function(datadir = c(), metadatadir = c(), loglocation = c(), f na.rm = TRUE) personSummarynames = c(personSummarynames, paste("number_sib_wakinghours_", TW, "_", udefn[j], "_mn", sep = ""), paste("number_sib_wakinghours_", TW, "_", udefn[j], "_sd", sep = "")) - personSummary[i, (cnt + 15)] = mean(nightsummary.tmp$duration_sib_wakinghours_atleast15min[indexUdef], na.rm = TRUE) personSummary[i, (cnt + 16)] = sd(nightsummary.tmp$duration_sib_wakinghours_atleast15min[indexUdef], @@ -443,8 +442,7 @@ g.report.part4 = function(datadir = c(), metadatadir = c(), loglocation = c(), f 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 + 23)] = NDAYsibd personSummarynames = c(personSummarynames, paste("n_days_w_sib_wakinghours_", TW, "_", udefn[j], sep = "")) @@ -465,7 +463,7 @@ g.report.part4 = function(datadir = c(), metadatadir = c(), loglocation = c(), f personSummary[i, (cnt + 31)] = sd(nightsummary.tmp$SriFractionValid[indexUdef], na.rm = TRUE) personSummarynames = c(personSummarynames, paste("SriFractionValid_", TW, "_", udefn[j], "_mn", sep = ""), paste("SriFractionValid_", TW, "_", udefn[j], "_sd", sep = "")) - cnt = cnt + 27 + cnt = cnt + 31 if (sleepwindowType == "TimeInBed") { personSummary[i, (cnt + 1)] = mean(nightsummary.tmp$sleepefficiency[indexUdef], na.rm = TRUE) personSummary[i, (cnt + 2)] = sd(nightsummary.tmp$sleepefficiency[indexUdef], na.rm = TRUE) From ce7267621ebd4755d0f4a694850020afdaee660b Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Tue, 14 Feb 2023 18:07:43 +0100 Subject: [PATCH 25/34] comment out print statement --- R/g.report.part5.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/g.report.part5.R b/R/g.report.part5.R index 3f97bc58a..0038c9232 100644 --- a/R/g.report.part5.R +++ b/R/g.report.part5.R @@ -189,7 +189,7 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c "ID", "acc_available", "window_number", "boutcriter.mvpa", "boutcriter.lig", "boutcriter.in", "bout.metric", grep(pattern = "cosinor", x = names(df), ignore.case = TRUE, value = TRUE)) # skip cosinor variables - print(ignorevar) + # print(ignorevar) for (ee in 1:ncol(df)) { # make sure that numeric columns have class numeric nr = nrow(df) if (nr > 30) nr = 30 From 98f00f58eff51f5481925963657a15092001820b Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Sun, 18 Jun 2023 14:08:03 +0200 Subject: [PATCH 26/34] fix merge problem with previous commit --- R/g.part4.R | 22 +++++++++++----------- R/g.part5.R | 26 +++++++++++++------------- 2 files changed, 24 insertions(+), 24 deletions(-) diff --git a/R/g.part4.R b/R/g.part4.R index 5a14d11cb..708df267a 100644 --- a/R/g.part4.R +++ b/R/g.part4.R @@ -72,8 +72,8 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1, "sleeplatency", "sleepefficiency", "page", "daysleeper", "weekday", "calendar_date", "filename", "cleaningcode", "sleeplog_used", "sleeplog_ID", "acc_available", "guider", "SleepRegularityIndex", "SriFractionValid", "longitudinal_axis") # - - + + if (params_output[["storefolderstructure"]] == TRUE) { colnamesnightsummary = c(colnamesnightsummary, "filename_dir", "foldername") } @@ -349,10 +349,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 @@ -462,7 +462,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) sleepdet.t = sleepdet[ki, ] @@ -684,13 +684,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) { @@ -872,7 +872,7 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1, spocum.n.dur_sibd_atleast15min = length(sibds_atleast15min) } } - + 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 @@ -895,12 +895,12 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1, nightsummary[sumi, 23] = acc_onsetTS nightsummary[sumi, 24] = acc_wakeTS #---------------------------------------------- - nightsummary[sumi, 23] = tmp1 #guider_onset_ts - nightsummary[sumi, 24] = tmp4 #guider_onset_ts + nightsummary[sumi, 25] = tmp1 #guider_onset_ts + nightsummary[sumi, 26] = tmp4 #guider_wake_ts if (params_sleep[["sleepwindowType"]] == "TimeInBed") { # If guider isa sleeplog and if the sleeplog recorded time in bed then # calculate: sleep latency: - nightsummary[sumi, 25] = round(nightsummary[sumi, 3] - nightsummary[sumi, 7], + nightsummary[sumi, 27] = round(nightsummary[sumi, 3] - nightsummary[sumi, 7], digits = 7) #sleeponset - guider_onset # sleep efficiency: nightsummary[sumi, 28] = round(nightsummary[sumi, 14]/nightsummary[sumi, 9], digits = 5) #accumulated nocturnal sleep / guider diff --git a/R/g.part5.R b/R/g.part5.R index 31f5c9bc3..2289cda19 100644 --- a/R/g.part5.R +++ b/R/g.part5.R @@ -62,7 +62,7 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), #====================================================================== # compile lists of milestone data filenames fnames.ms3 = dir(paste(metadatadir, "/meta/ms3.out", sep = "")) - + fnames.ms5 = dir(paste(metadatadir, "/meta/ms5.out", sep = "")) # path to sleeplog milestonedata, if it exists: sleeplogRDA = paste(metadatadir, "/meta/sleeplog.RData", sep = "") @@ -269,7 +269,7 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), S$sib.end.time[replaceLastWakeup] = ts$time[nrow(ts)] } } - + for (j in def) { # loop through sleep definitions (defined by angle and time threshold in g.part3) ws3new = ws3 # reset wse3new, because if part5_agg2_60seconds is TRUE then this will have been change in the previous iteration of the loop if (params_general[["part5_agg2_60seconds"]] == TRUE) { @@ -311,7 +311,7 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), Nepochsinhour, Nts, SPTE_end, ws3) if (params_general[["part5_agg2_60seconds"]] == TRUE) { # Optionally aggregate to 1 minute epoch: ts$time_num = floor(as.numeric(iso8601chartime2POSIX(ts$time,tz = params_general[["desiredtz"]])) / 60) * 60 - + # only include angle if angle is present angleColName = ifelse("angle" %in% names(ts), yes = "angle", no = NULL) if (lightpeak_available == TRUE) { @@ -379,14 +379,14 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), dir.create(file.path(metadatadir, ms5.sibreport)) } shortendFname = gsub(pattern = "[.]|RData|csv|cwa|bin", replacement = "", x = fnames.ms3[i], ignore.case = TRUE) - + sibreport_fname = paste0(metadatadir,ms5.sibreport,"/sib_report_", shortendFname, "_",j,".csv") data.table::fwrite(x = sibreport, file = sibreport_fname, row.names = FALSE, sep = params_output[["sep_reports"]]) # nap/sib/nonwear overlap analysis - + # TO DO - + # nap detection if (params_general[["acc.metric"]] != "ENMO" | params_sleep[["HASIB.algo"]] != "vanHees2015") { @@ -441,7 +441,7 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), } } ts$window = 0 - + backup_cosinor_MM = backup_cosinor_WW = NULL # 2023-04-23 - backup of nightsi outside threshold look to avoid # overwriting the backup after the first iteration @@ -936,7 +936,7 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), } } } - + #=============================================================== # Cosinor analyses based on only the data used for GGIR part5 if (params_247[["cosinor"]] == TRUE) { @@ -975,7 +975,7 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), cosinor_coef$coef$params$ndays, cosinor_coef$coef$params$R2)), times = di), nrow = di, byrow = TRUE)}, silent = TRUE) - + ds_names[fi:(fi + 5)] = c("cosinor_mes", "cosinor_amp", "cosinor_acrophase", "cosinor_acrotime", "cosinor_ndays", "cosinor_R2") fi = fi + 6 @@ -1008,9 +1008,9 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), cosinor_coef = c() fi = fi + 20 } - + } - + if (params_output[["save_ms5rawlevels"]] == TRUE) { legendfile = paste0(metadatadir,ms5.outraw,"/behavioralcodes",as.Date(Sys.time()),".csv") if (file.exists(legendfile) == FALSE) { @@ -1042,10 +1042,10 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), } } } - + #=============================================================== } - + } if ("angle" %in% colnames(ts)) { ts = ts[, -which(colnames(ts) == "angle")] From 7d474389a223aa1ee868caa08cc1d6aadb95983e Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Tue, 11 Jul 2023 12:01:50 +0200 Subject: [PATCH 27/34] undo indentation changes --- R/g.report.part5.R | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/R/g.report.part5.R b/R/g.report.part5.R index 70a72390a..58b9d3761 100644 --- a/R/g.report.part5.R +++ b/R/g.report.part5.R @@ -52,11 +52,11 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c # identify first and last day per file first_days = aggregate(window_number ~ filename, data = x, FUN = min, na.rm = TRUE) last_days = aggregate(window_number ~ filename, data = x, FUN = max, na.rm = TRUE) - + # match first and last days with the output dataframe exclude_firsts = which(paste(x$filename, x$window_number) %in% paste(first_days$filename, first_days$window_number)) exclude_lasts = which(paste(x$filename, x$window_number) %in% paste(last_days$filename, last_days$window_number)) - + # keep only indices that do not match with first and last days indices2exclude = which(indices %in% c(exclude_firsts, exclude_lasts)) if (length(indices2exclude) > 0) indices = indices[-indices2exclude] @@ -116,7 +116,7 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c # order data.frame outputfinal$window_number = as.numeric(gsub(" ", "", outputfinal$window_number)) outputfinal = outputfinal[order(outputfinal$filename, outputfinal$window, outputfinal$window_number),] - + # Find columns filled with missing values cut = which(sapply(outputfinal, function(x) all(x == "")) == TRUE) if (length(cut) > 0) { @@ -199,7 +199,7 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c agg_plainNweighted = function(df, filename = "filename", daytype = "daytype") { - + # function to take both the weighted (by weekday/weekendday) and plain average of all numeric variables # df: input data.frame (OF3 outside this function) ignorevar = c("daysleeper", "cleaningcode", "night_number", "sleeplog_used", @@ -512,12 +512,12 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c # store all summaries in csv files OF4_clean = tidyup_df(OF4) data.table::fwrite(OF4_clean,paste(metadatadir,"/results/part5_personsummary_", - uwi[j],"_L",uTRLi[h1],"M", uTRMi[h2], "V", uTRVi[h3], "_", usleepparam[h4], ".csv", sep = ""), + uwi[j],"_L",uTRLi[h1],"M", uTRMi[h2], "V", uTRVi[h3], "_", usleepparam[h4], ".csv", sep = ""), row.names = FALSE, na = "", sep = sep_reports) } } } - + } } } From f26e5069b43e0bd46ea9bbb15cdf3a109b433287 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Tue, 11 Jul 2023 14:57:10 +0200 Subject: [PATCH 28/34] #708 omit addition of extra daytime sib variables to part 4, the plan is to add them to part 5 --- R/g.part4.R | 72 +++++++++++++++++++++------------------------- R/g.report.part4.R | 45 ++++++++--------------------- 2 files changed, 45 insertions(+), 72 deletions(-) diff --git a/R/g.part4.R b/R/g.part4.R index 22b86f8b3..c226d6dd1 100644 --- a/R/g.part4.R +++ b/R/g.part4.R @@ -66,8 +66,6 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1, "SleepDurationInSpt", "WASO", "duration_sib_wakinghours", "number_sib_sleepperiod", "number_of_awakenings", "number_sib_wakinghours", "duration_sib_wakinghours_atleast15min", - "meandur_sib_wakinghours_atleast15min", - "number_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", @@ -860,19 +858,18 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1, } # ====================================================================================================== sibds_atleast15min = 0 - spocum.t.dur_sibd = 0 - spocum.t.dur_sibd_atleast15min = 0 - spocum.m.dur_sibd_atleast15min = 0 - spocum.n.dur_sibd_atleast15min = 0 if (length(sibds) > 0) { spocum.t.dur_sibd = sum(sibds) atleast15min = which(sibds >= 1/4) if (length(atleast15min) > 0) { sibds_atleast15min = sibds[atleast15min] spocum.t.dur_sibd_atleast15min = sum(sibds_atleast15min) - spocum.m.dur_sibd_atleast15min = mean(sibds_atleast15min, na.rm = TRUE) - spocum.n.dur_sibd_atleast15min = length(sibds_atleast15min) + } else { + spocum.t.dur_sibd_atleast15min = 0 } + } else { + spocum.t.dur_sibd = 0 + spocum.t.dur_sibd_atleast15min = 0 } nightsummary[sumi, 14] = spocum.t.dur.noc #total nocturnalsleep /accumulated sleep duration @@ -882,8 +879,6 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1, nightsummary[sumi, 18] = nightsummary[sumi, 17] - 1 #number of awakenings nightsummary[sumi, 19] = length(which(spocum.t$dur == 0)) #number of sib (sustained inactivty bout) during wakinghours nightsummary[sumi, 20] = as.numeric(spocum.t.dur_sibd_atleast15min) #total sib (sustained inactivty bout) duration during wakinghours of at least 5 minutes - nightsummary[sumi, 21] = as.numeric(spocum.m.dur_sibd_atleast15min) - nightsummary[sumi, 22] = as.numeric(spocum.n.dur_sibd_atleast15min) #------------------------------------------------------- # Also report timestamps in non-numeric format: acc_onset = nightsummary[sumi, 3] @@ -894,24 +889,24 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1, # convert into clocktime acc_onsetTS = convertHRsinceprevMN2Clocktime(acc_onset) acc_wakeTS = convertHRsinceprevMN2Clocktime(acc_wake) - nightsummary[sumi, 23] = acc_onsetTS - nightsummary[sumi, 24] = acc_wakeTS + nightsummary[sumi, 21] = acc_onsetTS + nightsummary[sumi, 22] = acc_wakeTS #---------------------------------------------- - nightsummary[sumi, 25] = tmp1 #guider_onset_ts - nightsummary[sumi, 26] = tmp4 #guider_wake_ts + nightsummary[sumi, 23] = tmp1 #guider_onset_ts + nightsummary[sumi, 24] = tmp4 #guider_wake_ts if (params_sleep[["sleepwindowType"]] == "TimeInBed") { # If guider isa sleeplog and if the sleeplog recorded time in bed then # calculate: sleep latency: - nightsummary[sumi, 27] = round(nightsummary[sumi, 3] - nightsummary[sumi, 7], + nightsummary[sumi, 25] = round(nightsummary[sumi, 3] - nightsummary[sumi, 7], digits = 7) #sleeponset - guider_onset # sleep efficiency: - nightsummary[sumi, 28] = round(nightsummary[sumi, 14]/nightsummary[sumi, 9], digits = 5) #accumulated nocturnal sleep / guider + nightsummary[sumi, 26] = round(nightsummary[sumi, 14]/nightsummary[sumi, 9], digits = 5) #accumulated nocturnal sleep / guider } - nightsummary[sumi, 29] = pagei - nightsummary[sumi, 30] = daysleeper[j] - nightsummary[sumi, 31] = wdayname[j] - nightsummary[sumi, 32] = calendar_date[j] - nightsummary[sumi, 33] = fnames[i] + nightsummary[sumi, 27] = pagei + nightsummary[sumi, 28] = daysleeper[j] + nightsummary[sumi, 29] = wdayname[j] + nightsummary[sumi, 30] = calendar_date[j] + nightsummary[sumi, 31] = fnames[i] # nightsummary #------------------------------------------------------------------------ # PLOT @@ -988,13 +983,13 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1, } # PLOT #------------------------------------------------------------------------ - nightsummary[sumi, 34] = cleaningcode - nightsummary[sumi, 35] = sleeplog_used - nightsummary[sumi, 36] = logid - nightsummary[sumi, 37] = acc_available - nightsummary[sumi, 38] = guider + nightsummary[sumi, 32] = cleaningcode + nightsummary[sumi, 33] = sleeplog_used + nightsummary[sumi, 34] = logid + nightsummary[sumi, 35] = acc_available + nightsummary[sumi, 36] = guider # Extract SRI for this night - nightsummary[sumi, 39:40] = NA + nightsummary[sumi, 37:38] = NA if (!exists("SleepRegularityIndex")) { SleepRegularityIndex = NA } @@ -1004,18 +999,18 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1, calendar_date_reformat = format(x = calendar_date_asDate, format = "%d/%m/%Y") SRIindex = which(SRI$date == calendar_date_reformat & SRI$frac_valid > (params_cleaning[["includenightcrit"]]/24)) if (length(SRIindex) > 0) { - nightsummary[sumi, 39] = SRI$SleepRegularityIndex[SRIindex[1]] - nightsummary[sumi, 40] = SRI$frac_valid[SRIindex[1]] + nightsummary[sumi, 37] = SRI$SleepRegularityIndex[SRIindex[1]] + nightsummary[sumi, 38] = SRI$frac_valid[SRIindex[1]] } } if (length(longitudinal_axis) == 0) { - nightsummary[sumi, 41] = NA + nightsummary[sumi, 39] = NA } else { - nightsummary[sumi, 41] = longitudinal_axis + nightsummary[sumi, 39] = longitudinal_axis } if (params_output[["storefolderstructure"]] == TRUE) { - nightsummary[sumi, 42] = ffd[i] #full filename structure - nightsummary[sumi, 43] = ffp[i] #use the lowest foldername as foldername name + nightsummary[sumi, 40] = ffd[i] #full filename structure + nightsummary[sumi, 41] = ffp[i] #use the lowest foldername as foldername name } sumi = sumi + 1 } #run through definitions @@ -1048,13 +1043,12 @@ 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:32] = NA - nightsummary[sumi, 33] = fnames[i] - nightsummary[sumi, 34] = 4 #cleaningcode = 4 (no nights of accelerometer available) - nightsummary[sumi, 35:36] = c(FALSE, TRUE) #sleeplog_used acc_available - nightsummary[sumi, 37:40] = 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 if (params_output[["storefolderstructure"]] == TRUE) { - nightsummary[sumi, 41:42] = c(ffd[i], ffp[i]) #full filename structure and use the lowest foldername as foldername name + nightsummary[sumi, 40:41] = c(ffd[i], ffp[i]) #full filename structure and use the lowest foldername as foldername name } sumi = sumi + 1 } diff --git a/R/g.report.part4.R b/R/g.report.part4.R index 89e078e45..fc97ec3e0 100644 --- a/R/g.report.part4.R +++ b/R/g.report.part4.R @@ -298,10 +298,6 @@ g.report.part4 = function(datadir = c(), metadatadir = c(), loglocation = c(), 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", - "meandur_sib_wakinghours_atleast15min", - "number_sib_wakinghours_atleast15min", "sleeplatency", "sleepefficiency", "number_of_awakenings", "guider_inbedDuration", "guider_inbedStart", "guider_inbedEnd", "guider_SptDuration", "guider_onset", "guider_wakeup", "SleepRegularityIndex", "SriFractionValid")) @@ -420,56 +416,39 @@ g.report.part4 = function(datadir = c(), metadatadir = c(), loglocation = c(), personSummarynames = c(personSummarynames, paste("duration_sib_wakinghours_atleast15min_", TW, "_", udefn[j], "_mn", sep = ""), paste("duration_sib_wakinghours_atleast15min_", TW, "_", udefn[j], "_sd", sep = "")) - personSummary[i, (cnt + 17)] = mean(nightsummary.tmp$meandur_sib_wakinghours_atleast15min[indexUdef], - na.rm = TRUE) - personSummary[i, (cnt + 18)] = sd(nightsummary.tmp$meandur_sib_wakinghours_atleast15min[indexUdef], - na.rm = TRUE) - personSummarynames = c(personSummarynames, paste("meanduration_sib_wakinghours_atleast15min_", - TW, "_", udefn[j], "_mn", sep = ""), - paste("meanduration_sib_wakinghours_atleast15min_", - TW, "_", udefn[j], "_sd", sep = "")) - personSummary[i, (cnt + 19)] = mean(nightsummary.tmp$number_sib_wakinghours_atleast15min[indexUdef], - na.rm = TRUE) - personSummary[i, (cnt + 20)] = sd(nightsummary.tmp$number_sib_wakinghours_atleast15min[indexUdef], - na.rm = TRUE) - personSummarynames = c(personSummarynames, paste("number_sib_wakinghours_atleast15min_", - TW, "_", udefn[j], "_mn", sep = ""), - paste("number_sib_wakinghours_atleast15min_", - TW, "_", udefn[j], "_sd", sep = "")) - # average sibd during the day AVEsibdDUR = c(nightsummary.tmp$duration_sib_wakinghours[indexUdef]/nightsummary.tmp$number_sib_wakinghours[indexUdef]) if (length(which(nightsummary.tmp$number_sib_wakinghours[indexUdef] == 0))) { AVEsibdDUR[which(nightsummary.tmp$number_sib_wakinghours[indexUdef] == 0)] = 0 } - personSummary[i, (cnt + 21)] = mean(AVEsibdDUR, na.rm = TRUE) - personSummary[i, (cnt + 22)] = sd(AVEsibdDUR, na.rm = TRUE) + personSummary[i, (cnt + 17)] = mean(AVEsibdDUR, na.rm = TRUE) + personSummary[i, (cnt + 18)] = sd(AVEsibdDUR, na.rm = TRUE) personSummarynames = c(personSummarynames, paste("average_dur_sib_wakinghours_", TW, "_", 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 - personSummary[i, (cnt + 23)] = NDAYsibd + personSummary[i, (cnt + 19)] = NDAYsibd personSummarynames = c(personSummarynames, paste("n_days_w_sib_wakinghours_", TW, "_", udefn[j], sep = "")) - personSummary[i, (cnt + 24)] = mean(nightsummary.tmp$sleeponset[indexUdef], na.rm = TRUE) - personSummary[i, (cnt + 25)] = sd(nightsummary.tmp$sleeponset[indexUdef], na.rm = TRUE) + personSummary[i, (cnt + 20)] = mean(nightsummary.tmp$sleeponset[indexUdef], na.rm = TRUE) + personSummary[i, (cnt + 21)] = sd(nightsummary.tmp$sleeponset[indexUdef], na.rm = TRUE) personSummarynames = c(personSummarynames, paste("sleeponset_", TW, "_", udefn[j], "_mn", sep = ""), paste("sleeponset_", TW, "_", udefn[j], "_sd", sep = "")) - personSummary[i, (cnt + 26)] = mean(nightsummary.tmp$wakeup[indexUdef], na.rm = TRUE) - personSummary[i, (cnt + 27)] = sd(nightsummary.tmp$wakeup[indexUdef], na.rm = TRUE) + personSummary[i, (cnt + 22)] = mean(nightsummary.tmp$wakeup[indexUdef], na.rm = TRUE) + personSummary[i, (cnt + 23)] = sd(nightsummary.tmp$wakeup[indexUdef], na.rm = TRUE) personSummarynames = c(personSummarynames, paste("wakeup_", TW, "_", udefn[j], "_mn", sep = ""), paste("wakeup_", TW, "_", udefn[j], "_sd", sep = "")) - personSummary[i, (cnt + 28)] = mean(nightsummary.tmp$SleepRegularityIndex[indexUdef], + personSummary[i, (cnt + 24)] = mean(nightsummary.tmp$SleepRegularityIndex[indexUdef], na.rm = TRUE) - personSummary[i, (cnt + 29)] = sd(nightsummary.tmp$SleepRegularityIndex[indexUdef], na.rm = TRUE) + personSummary[i, (cnt + 25)] = sd(nightsummary.tmp$SleepRegularityIndex[indexUdef], na.rm = TRUE) personSummarynames = c(personSummarynames, paste("SleepRegularityIndex_", TW, "_", udefn[j], "_mn", sep = ""), paste("SleepRegularityIndex_", TW, "_", udefn[j], "_sd", sep = "")) - personSummary[i, (cnt + 30)] = mean(nightsummary.tmp$SriFractionValid[indexUdef], na.rm = TRUE) - personSummary[i, (cnt + 31)] = sd(nightsummary.tmp$SriFractionValid[indexUdef], na.rm = TRUE) + personSummary[i, (cnt + 26)] = mean(nightsummary.tmp$SriFractionValid[indexUdef], na.rm = TRUE) + personSummary[i, (cnt + 27)] = sd(nightsummary.tmp$SriFractionValid[indexUdef], na.rm = TRUE) personSummarynames = c(personSummarynames, paste("SriFractionValid_", TW, "_", udefn[j], "_mn", sep = ""), paste("SriFractionValid_", TW, "_", udefn[j], "_sd", sep = "")) - cnt = cnt + 31 + cnt = cnt + 27 if (sleepwindowType == "TimeInBed") { personSummary[i, (cnt + 1)] = mean(nightsummary.tmp$sleepefficiency[indexUdef], na.rm = TRUE) personSummary[i, (cnt + 2)] = sd(nightsummary.tmp$sleepefficiency[indexUdef], na.rm = TRUE) From 7150cd35e4f2d5f12c3ce9dd90151430efc8aa35 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Tue, 11 Jul 2023 15:12:38 +0200 Subject: [PATCH 29/34] minor improvement to syntax --- R/g.report.part4.R | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/R/g.report.part4.R b/R/g.report.part4.R index fc97ec3e0..6482d6458 100644 --- a/R/g.report.part4.R +++ b/R/g.report.part4.R @@ -295,12 +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", - "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) { From 8e5a4f6a8a78f6f6ce37dd789fe3beb9f1c6e572 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Tue, 26 Sep 2023 15:11:58 +0200 Subject: [PATCH 30/34] undo addition of cosinor to part 5 as this will be moved to part 6 --- R/g.part5.R | 94 ----------------------------------------------------- 1 file changed, 94 deletions(-) diff --git a/R/g.part5.R b/R/g.part5.R index 2050559f8..ad075e80a 100644 --- a/R/g.part5.R +++ b/R/g.part5.R @@ -312,24 +312,6 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), } if (length(tail_expansion_log) != 0 & nrow(ts) > max(nightsi)) nightsi[length(nightsi) + 1] = nrow(ts) # include last window Nts = nrow(ts) - } else { - ts$time = iso8601chartime2POSIX(ts$time,tz = params_general[["desiredtz"]]) - ws3new = ws3 # change because below it is used to decide how many epochs are there in - # extract nightsi again - tempp = unclass(ts$time) - if (is.na(tempp$sec[1]) == TRUE) { - tempp = unclass(as.POSIXlt(ts$time, tz = params_general[["desiredtz"]])) - } - sec = tempp$sec - min = tempp$min - hour = tempp$hour - if (params_general[["dayborder"]] == 0) { - nightsi = which(sec == 0 & min == 0 & hour == 0) - } else { - nightsi = which(sec == 0 & min == (params_general[["dayborder"]] - floor(params_general[["dayborder"]])) * 60 & hour == floor(params_general[["dayborder"]])) #shift the definition of midnight if required - } - if (length(tail_expansion_log) != 0 & nrow(ts) > max(nightsi)) nightsi[length(nightsi) + 1] = nrow(ts) # include last window - Nts = nrow(ts) } #=============================================== # Use sib.report to classify naps, non-wear and integrate these in time series @@ -558,79 +540,6 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), } di = di + 1 } - - #=============================================================== - # Cosinor analyses based on only the data used for GGIR part5 - if (params_247[["cosinor"]] == TRUE) { - if ((timewindowi == "MM" & length(backup_cosinor_MM) == 0) | - (timewindowi == "WW" & length(backup_cosinor_WW) == 0)) { - # avoid computing same parameter twice because this part of the code is - # not dependent on the lig, mod, vig thresholds - acc4cos = ts[which(ts$window != 0), c("time", "ACC")] - acc4cos$ACC = acc4cos$ACC / 1000 # convert to mg because that is what applyCosinorAnalyses expects - cosinor_coef = applyCosinorAnalyses(ts = acc4cos, - qcheck = ts$nonwear[which(ts$window != 0)], - midnightsi = nightsi, - epochsizes = c(ws3new, ws3new)) - rm(acc4cos) - if (timewindowi == "MM") { - backup_cosinor_MM = cosinor_coef - } else if (timewindowi == "WW") { - backup_cosinor_WW = cosinor_coef - } - } else { - if (timewindowi == "MM") { - cosinor_coef = backup_cosinor_MM - } else if (timewindowi == "WW") { - cosinor_coef = backup_cosinor_WW - } - } - if (length(cosinor_coef) > 0) { - # assign same value to all rows to ease creating reports - dsummary[1:di, fi] = c(cosinor_coef$timeOffsetHours) - ds_names[fi] = c("cosinor_timeOffsetHours") - fi = fi + 1 - try(expr = {dsummary[1:di, fi:(fi + 5)] = matrix(rep(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, - cosinor_coef$coef$params$R2)), times = di), nrow = di, byrow = TRUE)}, - silent = TRUE) - - ds_names[fi:(fi + 5)] = c("cosinor_mes", "cosinor_amp", "cosinor_acrophase", - "cosinor_acrotime", "cosinor_ndays", "cosinor_R2") - fi = fi + 6 - try(expr = {dsummary[1:di, fi:(fi + 10)] = matrix(rep(c(cosinor_coef$coefext$params$minimum, - cosinor_coef$coefext$params$amp, - cosinor_coef$coefext$params$alpha, - cosinor_coef$coefext$params$beta, - cosinor_coef$coefext$params$acrotime, - cosinor_coef$coefext$params$UpMesor, - cosinor_coef$coefext$params$DownMesor, - cosinor_coef$coefext$params$MESOR, - cosinor_coef$coefext$params$ndays, - cosinor_coef$coefext$params$F_pseudo, - cosinor_coef$coefext$params$R2), times = di), nrow = di, byrow = TRUE)}, - silent = TRUE) - ds_names[fi:(fi + 10)] = c("cosinorExt_minimum", "cosinorExt_amp", "cosinorExt_alpha", - "cosinorExt_beta", "cosinorExt_acrotime", "cosinorExt_UpMesor", - "cosinorExt_DownMesor", "cosinorExt_MESOR", - "cosinorExt_ndays", "cosinorExt_F_pseudo", "cosinorExt_R2") - fi = fi + 11 - dsummary[1:di, fi:(fi + 1)] = matrix(rep(c(cosinor_coef$IVIS$InterdailyStability, - cosinor_coef$IVIS$IntradailyVariability), times = di), nrow = di) - ds_names[fi:(fi + 1)] = c("cosinorIS", "cosinorIV") - fi = fi + 2 - } else { - cosinor_coef = c() - fi = fi + 20 - } - } else { - cosinor_coef = c() - fi = fi + 20 - } - } if (params_output[["save_ms5rawlevels"]] == TRUE) { @@ -671,10 +580,7 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), } } } - - #=============================================================== } - } if ("angle" %in% colnames(ts)) { ts = ts[, -which(colnames(ts) == "angle")] From 16e838851b09ee3be86984544643143b09187f26 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Tue, 26 Sep 2023 15:17:31 +0200 Subject: [PATCH 31/34] remove cosinor from part5 --- R/g.part5.R | 2 +- R/g.report.part5.R | 12 +++--------- 2 files changed, 4 insertions(+), 10 deletions(-) diff --git a/R/g.part5.R b/R/g.part5.R index ad075e80a..c33dac70b 100644 --- a/R/g.part5.R +++ b/R/g.part5.R @@ -541,7 +541,6 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), di = di + 1 } } - if (params_output[["save_ms5rawlevels"]] == TRUE) { legendfile = paste0(metadatadir,ms5.outraw,"/behavioralcodes",as.Date(Sys.time()),".csv") if (file.exists(legendfile) == FALSE) { @@ -592,6 +591,7 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), } output = data.frame(dsummary,stringsAsFactors = FALSE) names(output) = ds_names + # correct definition of sleep log availability for window = WW, because now it # also relies on sleep log from previous night whoareWW = which(output$window == "WW") # look up WW diff --git a/R/g.report.part5.R b/R/g.report.part5.R index 4a89971d4..1469b4efd 100644 --- a/R/g.report.part5.R +++ b/R/g.report.part5.R @@ -215,16 +215,14 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c if (uwi[j] != "Segments") { delcol = c(delcol, which(colnames(outputfinal2) == "window")) } - outputfinal2 = outputfinal2[,-delcol] OF3 = outputfinal2[seluwi,] OF3 = as.data.frame(OF3, stringsAsFactors = TRUE) #------------------------------------------------------------- # store all summaries in csv files without cleaning criteria OF3_clean = tidyup_df(OF3) - colsWithoutCosinor = grep(pattern = "cosinor", x = colnames(OF3_clean), invert = TRUE) data.table::fwrite( - OF3_clean[, colsWithoutCosinor], + OF3_clean, paste( metadatadir, "/results/QC/part5_daysummary_full_", @@ -253,9 +251,7 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c ignorevar = c("daysleeper", "cleaningcode", "night_number", "sleeplog_used", "ID", "acc_available", "window_number", - "boutcriter.mvpa", "boutcriter.lig", "boutcriter.in", "bout.metric", - grep(pattern = "cosinor", x = names(df), ignore.case = TRUE, value = TRUE)) # skip cosinor variables - # print(ignorevar) + "boutcriter.mvpa", "boutcriter.lig", "boutcriter.in", "bout.metric") # skip cosinor variables for (ee in 1:ncol(df)) { # make sure that numeric columns have class numeric nr = nrow(df) if (nr > 30) nr = 30 @@ -295,7 +291,6 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c AggregateWDWE = aggregate.data.frame(df, by = by, plain_mean) AggregateWDWE = AggregateWDWE[, -grep("^Group", colnames(AggregateWDWE))] # Add counted number of days for Gini, Cov, alpha Fragmentation variables, because - # days are dropped if there are not enough fragments: vars_with_mininum_Nfrag = c("FRAG_Gini_dur_PA_day", "FRAG_CoV_dur_PA_day", "FRAG_alpha_dur_PA_day", "FRAG_Gini_dur_IN_day", @@ -333,6 +328,7 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c dt <- data.table::as.data.table(AggregateWDWE[,which(lapply(AggregateWDWE, class) == "numeric" | names(AggregateWDWE) == filename)]) } + options(warn = -1) .SD <- .N <- count <- a <- NULL if (window == "Segments") { @@ -547,7 +543,6 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c namenew = namenew, cval = 1, window = uwi[j]) # create variable from it - # do the same for daysleeper,cleaningcode, sleeplog_used, acc_available: OF3tmp$validdays = 1 # redefine by considering only valid days @@ -626,7 +621,6 @@ g.report.part5 = function(metadatadir = c(), f0 = c(), f1 = c(), loglocation = c # store all summaries in csv files OF4_clean = tidyup_df(OF4) data.table::fwrite(OF4_clean,paste(metadatadir,"/results/part5_personsummary_", - uwi[j],"_L",uTRLi[h1],"M", uTRMi[h2], "V", uTRVi[h3], "_", usleepparam[h4], ".csv", sep = ""), From 9efe0168484fabf5a1f5431f32dc3bd1a146f13d Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Tue, 26 Sep 2023 15:20:55 +0200 Subject: [PATCH 32/34] Update NEWS.md --- NEWS.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/NEWS.md b/NEWS.md index b3e6aa27a..c4ccfd655 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,10 @@ - Part 5: Function g.fragmentation code cleaned up, TP formulas revised following feedback Ian Danilevicz, and code now also produces sleep fragmentation #817 +- 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 + # CHANGES IN GGIR VERSION 2.10-3 - Part 1: Fixed minor bug in ismovisens that failed when datadir started with "./" #897 From 90a39462d4922531a3974c891156cfb031e6f173 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Tue, 26 Sep 2023 15:32:44 +0200 Subject: [PATCH 33/34] fix typo in documentation --- man/applyCosinorAnalyses.Rd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/man/applyCosinorAnalyses.Rd b/man/applyCosinorAnalyses.Rd index 77e3198e8..90109af20 100644 --- a/man/applyCosinorAnalyses.Rd +++ b/man/applyCosinorAnalyses.Rd @@ -22,7 +22,7 @@ Indices for midnights in the time series } \item{epochsizes}{ - Epoch seiz for ts and qcheck respectively + Epoch size for ts and qcheck respectively } } \author{ From 0fb6f924e700cb4f8fede5111d127738a35c7213 Mon Sep 17 00:00:00 2001 From: Vincent van Hees Date: Wed, 11 Oct 2023 09:42:42 +0200 Subject: [PATCH 34/34] add unit test for cosinor R2, remove incorrect statement about deprecation in GGIR.Rd, and move NEWS item to new 2.10-5 section #708 --- NEWS.md | 10 ++++++---- man/GGIR.Rd | 1 - tests/testthat/test_cosinor.R | 1 + 3 files changed, 7 insertions(+), 5 deletions(-) diff --git a/NEWS.md b/NEWS.md index 9b2f04717..f26c91b47 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,13 +1,15 @@ +# 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 + # CHANGES IN GGIR VERSION 2.10-4 - Part 4: Now better able to handle nights without sustained inactivity bouts (rest) #911 - Part 5: Function g.fragmentation code cleaned up, TP formulas revised following feedback Ian Danilevicz, and code now also produces sleep fragmentation #817 -- 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 1: Fixed issue #914 relating to specifying timezone for processing ad-hoc csv format raw data. - Part 1: Improved recognition of ActiGraph csv that occassionally confused for Axivity csv #918 diff --git a/man/GGIR.Rd b/man/GGIR.Rd index 670d23aa1..d3278919a 100755 --- a/man/GGIR.Rd +++ b/man/GGIR.Rd @@ -1249,7 +1249,6 @@ GGIR(mode = 1:5, over which L5M5 needs to be calculated. Now this is done with argument qwindow.} \item{cosinor}{ - Argument depricated after version 1.5-24. Boolean (default = FALSE). Whether to apply the cosinor analysis from the ActCR package.} } diff --git a/tests/testthat/test_cosinor.R b/tests/testthat/test_cosinor.R index a4fce97d0..11537498d 100644 --- a/tests/testthat/test_cosinor.R +++ b/tests/testthat/test_cosinor.R @@ -48,6 +48,7 @@ test_that("cosinorAnalyses provides expected output", { expect_equal(coef60$coefext$params$UpMesor, 19.52358, tolerance = 0.01) expect_equal(coef60$coefext$params$DownMesor, 16.47642, tolerance = 0.01) expect_equal(coef60$coefext$params$MESOR, 3.622534, tolerance = 0.01) + expect_equal(coef60$coefext$params$R2, 0.8976208, tolerance = 0.01) # IV IS expect_equal(coef60$IVIS$InterdailyStability, 0.9945789, tolerance = 0.01)