diff --git a/NAMESPACE b/NAMESPACE index 9a71c1811..8a4d2ca78 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -39,7 +39,9 @@ export(g.analyse, g.calibrate, ShellDoc2Vignette, parametersVignette, correctOlderMilestoneData, convertEpochData, appendRecords, extractID, - g.part5_analyseSegment, g.part5_initialise_ts) + g.part5_analyseSegment, g.part5_initialise_ts, + g.part5.analyseRest) + importFrom("grDevices", "colors", "dev.off", "pdf","rainbow","rgb") importFrom("graphics", "abline", "axis", "par", "plot", "plot.new", "rect","axis.POSIXct", "barplot", "box", "legend", diff --git a/NEWS.md b/NEWS.md index eacd19d42..70cb3de99 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,8 @@ # CHANGES IN GGIR VERSION 2.10-3 +- Part 5: Now able to assess overlap between self-reported naps and nonwear +with sustained inactivity bouts in order to facilitate research on nap detection. Fixes #687 + - Part 5: Time series now also exported if recording only includes one night, even though this is not sufficient for the main part 5 analyses. #894 Further, the time series now also come with lightpeak (LUX). - Visualreport: Bug fixed extracting the numeric value of the days to exclude in g.plot5 #879 diff --git a/R/check_params.R b/R/check_params.R index af67978f3..04a9f3e25 100644 --- a/R/check_params.R +++ b/R/check_params.R @@ -30,7 +30,7 @@ check_params = function(params_sleep = c(), params_metrics = c(), #----------------------------------------------------------------------------------------- if (length(params_sleep) > 0) { # Check class of sleep parameters numeric_params = c("anglethreshold", "timethreshold", "longitudinal_axis", "possible_nap_window", "possible_nap_dur", - "colid", "coln1", "def.noc.sleep", "nnights", "sleepefficiency.metric") + "colid", "coln1", "def.noc.sleep", "nnights", "sleepefficiency.metric", "possible_nap_edge_acc") boolean_params = c("ignorenonwear", "constrain2range", "HASPT.ignore.invalid", "relyonguider", "sleeplogidnum") character_params = c("HASPT.algo", "HASIB.algo", "Sadeh_axis", "nap_model", diff --git a/R/g.part5.R b/R/g.part5.R index 78b1fb8d0..2f1a906a9 100644 --- a/R/g.part5.R +++ b/R/g.part5.R @@ -228,6 +228,7 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), } for (sibDef 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) { ts = ts_backup @@ -315,7 +316,7 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), #=============================================== # Use sib.report to classify naps, non-wear and integrate these in time series # does not depend on bout detection criteria or window definitions. - if (params_output[["do.sibreport"]] == TRUE & length(params_sleep[["nap_model"]]) > 0) { + if (params_output[["do.sibreport"]] == TRUE) { IDtmp = as.character(ID) sibreport = g.sibreport(ts, ID = IDtmp, epochlength = ws3new, logs_diaries, desiredtz = params_general[["desiredtz"]]) @@ -325,66 +326,66 @@ 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, "_",sibDef,".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") { - warning("\nNap classification currently assumes acc.metric = ENMO and HASIB.algo = vanHees2015, so output may not be meaningful") - } - naps_nonwear = g.part5.classifyNaps(sibreport = sibreport, - desiredtz = params_general[["desiredtz"]], - possible_nap_window = params_sleep[["possible_nap_window"]], - possible_nap_dur = params_sleep[["possible_nap_dur"]], - nap_model = params_sleep[["nap_model"]], - HASIB.algo = params_sleep[["HASIB.algo"]]) - # store in ts object, such that it is exported in as time series - ts$nap1_nonwear2 = 0 - # napsindices = which(naps_nonwear$probability_nap == 1) - # if (length(napsindices) > 0) { - if (length(naps_nonwear) > 0) { - for (nni in 1:nrow(naps_nonwear)) { - nnc_window = which(time_POSIX >= naps_nonwear$start[nni] & time_POSIX <= naps_nonwear$end[nni] & ts$diur == 0) - if (length(nnc_window) > 0) { - if (naps_nonwear$probability_nap[nni] == 1) { - ts$nap1_nonwear2[nnc_window] = 1 # nap - } else if (naps_nonwear$probability_nap[nni] == 0) { - ts$nap1_nonwear2[nnc_window] = 2 # nonwear + if (length(params_sleep[["nap_model"]]) > 0) { + # nap detection + if (params_general[["acc.metric"]] != "ENMO" | + params_sleep[["HASIB.algo"]] != "vanHees2015") { + warning("\nNap classification currently assumes acc.metric = ENMO and HASIB.algo = vanHees2015, so output may not be meaningful") + } + naps_nonwear = g.part5.classifyNaps(sibreport = sibreport, + desiredtz = params_general[["desiredtz"]], + possible_nap_window = params_sleep[["possible_nap_window"]], + possible_nap_dur = params_sleep[["possible_nap_dur"]], + nap_model = params_sleep[["nap_model"]], + HASIB.algo = params_sleep[["HASIB.algo"]]) + # store in ts object, such that it is exported in as time series + ts$nap1_nonwear2 = 0 + # napsindices = which(naps_nonwear$probability_nap == 1) + # if (length(napsindices) > 0) { + if (length(naps_nonwear) > 0) { + for (nni in 1:nrow(naps_nonwear)) { + nnc_window = which(time_POSIX >= naps_nonwear$start[nni] & time_POSIX <= naps_nonwear$end[nni] & ts$diur == 0) + if (length(nnc_window) > 0) { + if (naps_nonwear$probability_nap[nni] == 1) { + ts$nap1_nonwear2[nnc_window] = 1 # nap + } else if (naps_nonwear$probability_nap[nni] == 0) { + ts$nap1_nonwear2[nnc_window] = 2 # nonwear + } } } } - } - # impute non-naps episodes (non-wear) - nonwearindices = which(naps_nonwear$probability_nap == 0) - if (length(nonwearindices) > 0) { - for (nni in nonwearindices) { - nwwindow_start = which(time_POSIX >= naps_nonwear$start[nni] & time_POSIX <= naps_nonwear$end[nni] & ts$diur == 0) - if (length(nwwindow_start) > 0) { - Nepochsin24Hours = (60/ws3new) * 60 * 24 - if (nwwindow_start[1] > Nepochsin24Hours) { - nwwindow = nwwindow_start - Nepochsin24Hours # impute time series with preceding day - if (length(which(ts$nap1_nonwear2[nwwindow] == 2)) / length(nwwindow) > 0.5) { - # if there is also a lot of overlap with non-wear there then do next day - nwwindow = nwwindow_start + Nepochsin24Hours + # impute non-naps episodes (non-wear) + nonwearindices = which(naps_nonwear$probability_nap == 0) + if (length(nonwearindices) > 0) { + for (nni in nonwearindices) { + nwwindow_start = which(time_POSIX >= naps_nonwear$start[nni] & time_POSIX <= naps_nonwear$end[nni] & ts$diur == 0) + if (length(nwwindow_start) > 0) { + Nepochsin24Hours = (60/ws3new) * 60 * 24 + if (nwwindow_start[1] > Nepochsin24Hours) { + nwwindow = nwwindow_start - Nepochsin24Hours # impute time series with preceding day + if (length(which(ts$nap1_nonwear2[nwwindow] == 2)) / length(nwwindow) > 0.5) { + # if there is also a lot of overlap with non-wear there then do next day + nwwindow = nwwindow_start + Nepochsin24Hours + } + } else { + nwwindow = nwwindow_start + Nepochsin24Hours # if there is not preceding day use next day } - } else { - nwwindow = nwwindow_start + Nepochsin24Hours # if there is not preceding day use next day - } - if (max(nwwindow) <= nrow(ts)) { # only attempt imputation if possible - # check again that there is not a lot of overlap with non-wear - if (length(which(ts$nap1_nonwear2[nwwindow] == 2)) / length(nwwindow) > 0.5) { - ts$ACC[nwwindow_start] = ts$ACC[nwwindow] # impute + if (max(nwwindow) <= nrow(ts)) { # only attempt imputation if possible + # check again that there is not a lot of overlap with non-wear + if (length(which(ts$nap1_nonwear2[nwwindow] == 2)) / length(nwwindow) > 0.5) { + ts$ACC[nwwindow_start] = ts$ACC[nwwindow] # impute + } } } } } } + } else { + sibreport = NULL } ts$window = 0 # backup of nightsi outside threshold defintions to avoid @@ -395,7 +396,7 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), for (TRVi in params_phyact[["threshold.vig"]]) { # derive behavioral levels (class), e.g. MVPA, inactivity bouts, etc. levelList = identify_levels(ts = ts, TRLi = TRLi, TRMi = TRMi, TRVi = TRVi, - ws3 = ws3new, params_phyact = params_phyact) + ws3 = ws3new, params_phyact = params_phyact) LEVELS = levelList$LEVELS OLEVELS = levelList$OLEVELS Lnames = levelList$Lnames @@ -503,7 +504,8 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), fullFilename = fullfilenames[i], add_one_day_to_next_date, lightpeak_available, tail_expansion_log, - foldernamei = foldername[i]) + foldernamei = foldername[i], + sibreport = sibreport) # Extract essential object to be used as input for the next # segment indexlog = gas$indexlog @@ -523,6 +525,17 @@ g.part5 = function(datadir = c(), metadatadir = c(), f0=c(), f1=c(), ws3new = timeList$epochSize if (doNext == TRUE) next } + #=============================================== + # FOLDER STRUCTURE + if (params_output[["storefolderstructure"]] == TRUE) { + if ("filename_dir" %in% ds_names) fi = which( ds_names == "filename_dir") + dsummary[di,fi] = fullfilenames[i] #full filename structure + ds_names[fi] = "filename_dir"; fi = fi + 1 + dsummary[di,fi] = foldername[i] #store the lowest foldername + if ("foldername" %in% ds_names) fi = which( ds_names == "foldername") + ds_names[fi] = "foldername"; fi = fi + 1 + } + di = di + 1 } } di = di + 1 diff --git a/R/g.part5.analyseRest.R b/R/g.part5.analyseRest.R new file mode 100644 index 000000000..fefedf6c8 --- /dev/null +++ b/R/g.part5.analyseRest.R @@ -0,0 +1,196 @@ +g.part5.analyseRest = function(sibreport = NULL, dsummary = NULL, + ds_names = NULL, fi = NULL, di = NULL, + time = NULL, tz = NULL, possible_nap_dur = 0, + possible_nap_edge_acc = Inf) { + # transform time to POSIX + if (is.ISO8601(as.character(time[1]))) { + time = iso8601chartime2POSIX(time, tz = tz) + } + + # Only consider sib episodes with minimum duration + if (length(grep(pattern = "mean_acc_1min", x = colnames(sibreport))) > 0) { + sibreport$acc_edge = pmax(sibreport$mean_acc_1min_before, sibreport$mean_acc_1min_after) + } else { + sibreport$acc_edge = 0 + } + longboutsi = which((sibreport$type == "sib" & + sibreport$duration >= possible_nap_dur[1] & + sibreport$acc_edge <= possible_nap_edge_acc) | + (sibreport$type != "sib" & sibreport$duration >= 1)) + # for qc purposes: + dsummary[di,fi] = length(longboutsi) + ds_names[fi] = "sibreport_n_items" + fi = fi + 1 + if (length(longboutsi) > 0) { + sibreport = sibreport[longboutsi,] + srep_tmp = sibreport[which(sibreport$start > min(time) & + sibreport$end < max(time)),] + + # account for possibility that some of these categories do not exis + # identify overlapping and non-overlapping, (nap-sib, non-wear-sib, sib, nap, nonwear) + # calculate for all five categories number, total duration, mean duration + # for qc purposes: + dsummary[di,fi] = nrow(srep_tmp) + ds_names[fi] = "sibreport_n_items_day" + fi = fi + 1 + + if (nrow(srep_tmp) > 0) { + sibs = which(srep_tmp$type == "sib") + srep_tmp$SIBoverlapNonwear = 0 + srep_tmp$SIBoverlapNap = 0 + srep_tmp$NonwearOverlapSIB = 0 + srep_tmp$NapOverlapSIB = 0 + srep_tmp$start = as.POSIXlt(srep_tmp$start, tz = tz) + srep_tmp$end = as.POSIXlt(srep_tmp$end, tz = tz) + # # for qc purposes: + # dsummary[di,fi] = nrow(srep_tmp) + # ds_names[fi] = "n_sibs_sibreport" + # fi = fi + 1 + if (length(sibs) > 0) { + classes = unique(srep_tmp$type) + selfreport = which(srep_tmp$type == "nonwear" | srep_tmp$type == "nap") + if (length(selfreport) > 0) { + for (si in sibs) { + for (sr in selfreport) { + if (srep_tmp$start[si] <= srep_tmp$end[sr] & + srep_tmp$end[si] >= srep_tmp$start[sr]) { + end_overlap = as.numeric(pmin(srep_tmp$end[si], srep_tmp$end[sr])) + start_overlap = as.numeric(pmax(srep_tmp$start[si], srep_tmp$start[sr])) + duration_overlap = end_overlap - start_overlap + duration_sib = as.numeric(srep_tmp$end[si]) - as.numeric(srep_tmp$start[si]) + perc_overlap = round(100 * (duration_overlap / duration_sib), digits = 1) + if (srep_tmp$type[sr] == "nonwear") { + srep_tmp$SIBoverlapNonwear[si] = perc_overlap + } else if (srep_tmp$type[sr] == "nap") { + srep_tmp$SIBoverlapNap[si] = perc_overlap + } + } + + if (srep_tmp$start[sr] <= srep_tmp$end[si] & + srep_tmp$end[sr] >= srep_tmp$start[si]) { + end_overlap = as.numeric(pmin(srep_tmp$end[si], srep_tmp$end[sr])) + start_overlap = as.numeric(pmax(srep_tmp$start[si], srep_tmp$start[sr])) + duration_overlap = end_overlap - start_overlap + duration_sr = as.numeric(srep_tmp$end[sr]) - as.numeric(srep_tmp$start[sr]) + perc_overlap = round(100 * (duration_overlap / duration_sr), digits = 1) + if (srep_tmp$type[sr] == "nonwear") { + srep_tmp$NonwearOverlapSIB[sr] = perc_overlap + } else if (srep_tmp$type[sr] == "nap") { + srep_tmp$NapOverlapSIB[sr] = perc_overlap + } + } + + } + } + } + } + + sibs_indices = which(srep_tmp$type == "sib") + nap_indices = which(srep_tmp$type == "nap") + nonwear_indices = which(srep_tmp$type == "nonwear") + SIBoverlapNap_indices = which(srep_tmp$SIBoverlapNap != 0) + SIBoverlapNonwear_indices = which(srep_tmp$SIBoverlapNonwear != 0) + NapOverlapSIB_indices = which(srep_tmp$NapOverlapSIB != 0) + NonwearOverlapSIB_indices = which(srep_tmp$NonwearOverlapSIB != 0) + dsummary[di,fi:(fi + 6)] = c(length(sibs_indices), + length(nap_indices), + length(nonwear_indices), + length(SIBoverlapNap_indices), + length(SIBoverlapNonwear_indices), + length(NapOverlapSIB_indices), + length(NonwearOverlapSIB_indices)) + + ds_names[fi:(fi + 6)] = c("nbouts_day_sib", "nbouts_day_srnap", "nbouts_day_srnonw", + "noverl_sib_srnap", "noverl_sib_srnonw", + "noverl_srnap_sib", "noverl_srnonw_sib") + fi = fi + 7 + if (length(sibs_indices) > 0) { + dsummary[di,fi:(fi + 1)] = c(mean(srep_tmp$duration[sibs_indices]), + sum(srep_tmp$duration[sibs_indices])) + } else { + dsummary[di,fi:(fi + 1)] = c(0, 0) + } + ds_names[fi:(fi + 1)] = c("frag_mean_dur_sib_day", "dur_day_sib_min") + fi = fi + 2 + if (length(nap_indices) > 0) { + srep_tmp$duration[nap_indices] = (as.numeric(srep_tmp$end[nap_indices]) - + as.numeric(srep_tmp$start[nap_indices])) / 60 + dsummary[di,fi:(fi + 1)] = c(mean(srep_tmp$duration[nap_indices]), + sum(srep_tmp$duration[nap_indices])) + } else { + dsummary[di,fi:(fi + 1)] = c(0, 0) + } + + ds_names[fi:(fi + 1)] = c("frag_mean_dur_srnap_day", "dur_day_srnap_min") + fi = fi + 2 + + if (length(nonwear_indices) > 0) { + dsummary[di,fi:(fi + 1)] = c(mean(srep_tmp$duration[nonwear_indices]), + sum(srep_tmp$duration[nonwear_indices])) + } else { + dsummary[di,fi:(fi + 1)] = c(0, 0) + } + ds_names[fi:(fi + 1)] = c("frag_mean_dur_srnonw_day", "dur_day_srnonw_min") + fi = fi + 2 + # Overlap sib with srnap + calcOverlapPercentage = function(overlap, duration) { + return(sum(overlap * duration) / sum(duration)) + } + if (length(SIBoverlapNap_indices) > 0) { + overlap_perc = calcOverlapPercentage(overlap = srep_tmp$SIBoverlapNap[SIBoverlapNap_indices], + duration = srep_tmp$duration[SIBoverlapNap_indices]) + dsummary[di,fi:(fi + 2)] = c(mean(srep_tmp$duration[SIBoverlapNap_indices]), + sum(srep_tmp$duration[SIBoverlapNap_indices]), + overlap_perc) + } else { + dsummary[di,fi:(fi + 2)] = c(0, 0, 0) + } + ds_names[fi:(fi + 2)] = c("mdur_sib_overl_srnap", "tdur_sib_overl_srnap", "perc_sib_overl_srnap") + fi = fi + 3 + + # Overlap srnap with sib + calcOverlapPercentage = function(overlap, duration) { + return(sum(overlap * duration) / sum(duration)) + } + if (length(SIBoverlapNap_indices) > 0) { + overlap_perc = calcOverlapPercentage(overlap = srep_tmp$NapOverlapSIB[NapOverlapSIB_indices], + duration = srep_tmp$duration[NapOverlapSIB_indices]) + dsummary[di,fi:(fi + 2)] = c(mean(srep_tmp$duration[NapOverlapSIB_indices]), + sum(srep_tmp$duration[NapOverlapSIB_indices]), + overlap_perc) + } else { + dsummary[di,fi:(fi + 2)] = c(0, 0, 0) + } + ds_names[fi:(fi + 2)] = c("mdur_srnap_overl_sib", "tdur_srnap_overl_sib", "perc_srnap_overl_sib") + fi = fi + 3 + # Overlap sib with srnonw + if (length(SIBoverlapNonwear_indices) > 0) { + overlap_perc = calcOverlapPercentage(overlap = srep_tmp$SIBoverlapNonwear[SIBoverlapNonwear_indices], + duration = srep_tmp$duration[SIBoverlapNonwear_indices]) + dsummary[di,fi:(fi + 2)] = c(mean(srep_tmp$duration[SIBoverlapNonwear_indices]), + sum(srep_tmp$duration[SIBoverlapNonwear_indices]), + overlap_perc) + } else { + dsummary[di,fi:(fi + 2)] = c(0, 0, 0) + } + ds_names[fi:(fi + 2)] = c("mdur_sib_overl_srnonw", "tdur_sib_overl_srnonw", "perc_sib_overl_srnonw") + fi = fi + 3 + # Overlap srnonw with sib + if (length(SIBoverlapNonwear_indices) > 0) { + overlap_perc = calcOverlapPercentage(overlap = srep_tmp$NonwearOverlapSIB[NonwearOverlapSIB_indices], + duration = srep_tmp$duration[NonwearOverlapSIB_indices]) + dsummary[di,fi:(fi + 2)] = c(mean(srep_tmp$duration[NonwearOverlapSIB_indices]), + sum(srep_tmp$duration[NonwearOverlapSIB_indices]), + overlap_perc) + } else { + dsummary[di,fi:(fi + 2)] = c(0, 0, 0) + } + ds_names[fi:(fi + 2)] = c("mdur_srnonw_overl_sib", "tdur_srnonw_overl_sib", "perc_srnonw_overl_sib") + fi = fi + 3 + rm(srep_tmp) + } else { + fi = fi + 26 + } + } + invisible(list(fi = fi, di = di, ds_names = ds_names, dsummary = dsummary)) +} diff --git a/R/g.part5_analyseSegment.R b/R/g.part5_analyseSegment.R index faac20ade..cf3cce9a0 100644 --- a/R/g.part5_analyseSegment.R +++ b/R/g.part5_analyseSegment.R @@ -10,7 +10,7 @@ g.part5_analyseSegment = function(indexlog, timeList, levelList, add_one_day_to_next_date, lightpeak_available, tail_expansion_log, - foldernamei) { + foldernamei, sibreport = NULL) { # unpack list objects: # indexlog fileIndex = indexlog$fileIndex @@ -449,6 +449,20 @@ g.part5_analyseSegment = function(indexlog, timeList, levelList, fi = fi + length(luxperseg$values) } } + #======================================================= + # nap/sib/nonwear overlap analysis + #======================================================= + if (params_output[["do.sibreport"]] == TRUE & !is.null(sibreport)) { + restAnalyses = g.part5.analyseRest(sibreport = sibreport, dsummary = dsummary, + ds_names = ds_names, fi = fi, di = si, + time = ts$time[sse[ts$diur[sse] == 0]], + tz = params_general[["desiredtz"]], + possible_nap_dur = params_sleep[["possible_nap_dur"]]) + fi = restAnalyses$fi + si = restAnalyses$di + dsummary = restAnalyses$dsummary + ds_names = restAnalyses$ds_names + } #=============================================== # FOLDER STRUCTURE if (params_output[["storefolderstructure"]] == TRUE) { diff --git a/R/g.sibreport.R b/R/g.sibreport.R index 87848db80..ff86ea4d8 100644 --- a/R/g.sibreport.R +++ b/R/g.sibreport.R @@ -5,7 +5,7 @@ g.sibreport = function(ts, ID, epochlength, logs_diaries=c(), desiredtz="") { Nsibs = length(sib_starts) sibreport = c() if (Nsibs > 0) { - sibreport = data.frame(ID= rep(ID, Nsibs), + sibreport = data.frame(ID = rep(ID, Nsibs), type = rep("sib", Nsibs), start = character(Nsibs), end = character(Nsibs), @@ -16,13 +16,13 @@ g.sibreport = function(ts, ID, epochlength, logs_diaries=c(), desiredtz="") { sibreport$start[sibi] = format(ts$time[dayind[sib_starts[sibi]]]) sibreport$end[sibi] = format(ts$time[dayind[sib_ends[sibi]]]) if (is.ISO8601(sibreport$start[sibi])) { - sibreport$start[sibi] = format(iso8601chartime2POSIX(sibreport$start[sibi],tz=desiredtz)) - sibreport$end[sibi] = format(iso8601chartime2POSIX(sibreport$end[sibi],tz=desiredtz)) + sibreport$start[sibi] = format(iso8601chartime2POSIX(sibreport$start[sibi], tz = desiredtz)) + sibreport$end[sibi] = format(iso8601chartime2POSIX(sibreport$end[sibi], tz = desiredtz)) } sibreport$duration[sibi] = ((sib_ends[sibi] - sib_starts[sibi]) + 1) / (60/epochlength) boutind = sib_starts[sibi]:sib_ends[sibi] - minute_before = (sib_starts[sibi]-(60/epochlength)):(sib_starts[sibi]-1) - minute_after = (sib_ends[sibi]+1):(sib_ends[sibi]+(60/epochlength)) + minute_before = (sib_starts[sibi] - (60/epochlength)):(sib_starts[sibi] - 1) + minute_after = (sib_ends[sibi] + 1):(sib_ends[sibi] + (60/epochlength)) if (min(minute_before) > 1) { sibreport$mean_acc_1min_before[sibi] = round(mean(ts$ACC[dayind[minute_before]]), digits = 3) } @@ -41,15 +41,19 @@ g.sibreport = function(ts, ID, epochlength, logs_diaries=c(), desiredtz="") { logreport = c() if (length(log) > 0) { relevant_rows = which(log$ID == ID) + if (length(relevant_rows) > 0) { log = log[relevant_rows,] # extract ID for (i in 1:nrow(log)) { # loop over lines (days) + tmp = log[i,] # convert into timestamps # only attempt if there are at least 2 timestamps to process nonempty = which(tmp[3:ncol(tmp)] != "" & tmp[3:ncol(tmp)] != "NA") if (length(nonempty) > 1) { date = as.Date(as.character(tmp[1,2]), format = dateformat) times = format(unlist(tmp[1,3:ncol(tmp)])) + times = grep(pattern = "NA", value = TRUE, invert = TRUE, x = times) + times = gsub(pattern = " ", replacement = "", x = times) times = times[which(times %in% c("", "NA") == FALSE)] # ignore entries without start and/or end time t_to_remove = c() @@ -66,13 +70,17 @@ g.sibreport = function(ts, ID, epochlength, logs_diaries=c(), desiredtz="") { if (length(times) > 1) { Nevents = floor(length(times) / 2) timestamps = sort(as.POSIXlt(paste0(date, " ", times), tz = desiredtz)) - logreport_tmp = data.frame(ID= rep(ID, Nevents), + logreport_tmp = data.frame(ID = rep(ID, Nevents), type = rep(logname, Nevents), start = rep("", Nevents), - end = rep("", Nevents), stringsAsFactors = FALSE) + end = rep("", Nevents), + duration = rep(0, Nevents), stringsAsFactors = FALSE) for (bi in 1:Nevents) { - logreport_tmp$start[bi] = format(timestamps[(bi*2)-1]) - logreport_tmp$end[bi] = format(timestamps[(bi*2)]) + tt1 = as.POSIXlt(timestamps[(bi * 2) - 1], tz = desiredtz) + tt2 = as.POSIXlt(timestamps[(bi * 2)], tz = desiredtz) + logreport_tmp$start[bi] = format(tt1) + logreport_tmp$end[bi] = format(tt2) + logreport_tmp$duration[bi] = abs(as.numeric(difftime(time1 = tt1, time2 = tt2, units = "mins"))) } } if (length(logreport) == 0) { @@ -86,17 +94,17 @@ g.sibreport = function(ts, ID, epochlength, logs_diaries=c(), desiredtz="") { } return(logreport) } - naplogreport = extract_logs(naplog, ID, logname="nap") - nonwearlogreport = extract_logs(nonwearlog, ID, logname="nonwear") + naplogreport = extract_logs(naplog, ID, logname = "nap") + nonwearlogreport = extract_logs(nonwearlog, ID, logname = "nonwear") logreport = sibreport # append all together in one output data.frame if (length(logreport) > 0 & length(naplogreport) > 0) { - logreport = merge(logreport, naplogreport, by=c("ID", "type", "start", "end"), all=TRUE) + logreport = merge(logreport, naplogreport, by = c("ID", "type", "start", "end", "duration"), all = TRUE) } else if (length(logreport) == 0 & length(naplogreport) > 0) { logreport = naplogreport } if (length(logreport) > 0 & length(nonwearlogreport) > 0) { - logreport = merge(logreport, nonwearlogreport, by=c("ID", "type", "start", "end"), all=TRUE) + logreport = merge(logreport, nonwearlogreport, by = c("ID", "type", "start", "end", "duration"), all = TRUE) } else if (length(logreport) == 0 & length(nonwearlogreport) > 0) { logreport = nonwearlogreport } diff --git a/R/load_params.R b/R/load_params.R index 5dfa437c1..c653d918f 100644 --- a/R/load_params.R +++ b/R/load_params.R @@ -26,7 +26,8 @@ load_params = function(group = c("sleep", "metrics", "rawdata", sleeplogsep = NULL, sleepwindowType = "SPT", possible_nap_window = c(9, 18), possible_nap_dur = c(15, 240), - nap_model = c(), sleepefficiency.metric = 1) + nap_model = c(), sleepefficiency.metric = 1, + possible_nap_edge_acc = Inf) } if ("metrics" %in% group) { params_metrics = list(do.anglex = FALSE, do.angley = FALSE, do.anglez = TRUE, diff --git a/man/GGIR.Rd b/man/GGIR.Rd index 8ee4a5f05..9c7107c0d 100755 --- a/man/GGIR.Rd +++ b/man/GGIR.Rd @@ -1104,7 +1104,10 @@ GGIR(mode = 1:5, wake times). If 2, sleep efficiency is calculated as detected nocturnal sleep time divided by detected duration in sleep period time plus sleep latency (where sleep latency refers to the difference between log-derived sleep onset and - accelerometer-detected sleep onset). + accelerometer-detected sleep onset).} + + \item{possible_nap_edge_acc}{ + Minimum acceleration before or after the SIB for the nap to be considered. } } } @@ -1291,7 +1294,11 @@ GGIR(mode = 1:5, \item{do.sibreport}{ Boolean (default = FALSE). In \link{g.part4}: To indicate whether to generate report for the sustained - inactivity bouts (SIB).} + inactivity bouts (SIB). If set to TRUE and when an advanced sleep diary is + available in part 4 then part 5 will use this to generate summary statistics + on the overlap between self-reported nonwear and napping with SIB. Here, + SIB cna be filter based on argument possible_nap_edge_acc and the first value + of possible_nap_dur} \item{do.visual}{ Boolean (default = TRUE). diff --git a/man/g.part5.analyseRest.Rd b/man/g.part5.analyseRest.Rd new file mode 100644 index 000000000..d4809f421 --- /dev/null +++ b/man/g.part5.analyseRest.Rd @@ -0,0 +1,53 @@ +\name{g.part5.analyseRest} +\alias{g.part5.analyseRest} +\title{ + Analyse rest (internal function) +} +\description{ + Analyses overlap self-reported napping, non-wear and sib. + Internal function not intended for direct use by GGIR end-user. +} +\usage{ + g.part5.analyseRest(sibreport = NULL, dsummary = NULL, + ds_names = NULL, fi = NULL, + di = NULL, time = NULL, tz = NULL, + possible_nap_dur = 0, + possible_nap_edge_acc = Inf) +} +\arguments{ + \item{sibreport}{ + sibreport data.frame produced by \link{g.sibreport} + } + \item{dsummary}{ + matrix created internally by \link{g.part5} + } + \item{ds_names}{ + character vector with variable names corresponding to dsummary + created internally by \link{g.part5} + } + \item{fi}{ + Numeric scalar to indicate variable index, created internally by \link{g.part5} + } + \item{di}{ + Numeric scalar to indicate recording index, created internally by \link{g.part5} + } + \item{time}{ + Daytime section of time column from the ts object, created internally by \link{g.part5}, + } + \item{tz}{ + Timezone database name + } + \item{possible_nap_dur}{ + Minimum sib duration to be considered. For self-reported naps/nonwear + there is a minimum duration of 1 epoch. + } + \item{possible_nap_edge_acc}{ + Maximum acceleration before or after the SIB for it to be considered. + } +} +\value{ + List with updated objects dsummary, ds_names, fi, and di +} +\author{ + Vincent T van Hees +} \ No newline at end of file diff --git a/man/g.part5_analyseSegment.Rd b/man/g.part5_analyseSegment.Rd index cef196d71..6a24a92af 100644 --- a/man/g.part5_analyseSegment.Rd +++ b/man/g.part5_analyseSegment.Rd @@ -19,7 +19,7 @@ add_one_day_to_next_date, lightpeak_available, tail_expansion_log, - foldernamei) + foldernamei, sibreport = NULL) } \arguments{ \item{indexlog}{ @@ -83,6 +83,9 @@ \item{foldernamei}{ Character with folder name in which the data file is stored. } + \item{sibreport}{ + Sibreport object as passed on from \link{g.part5} + } } \value{ List with objects: indexlog, timeList, and the matrix with the prelimenary results diff --git a/tests/testthat/test_chainof5parts.R b/tests/testthat/test_chainof5parts.R index 9cd0ac68a..8ee2a9937 100644 --- a/tests/testthat/test_chainof5parts.R +++ b/tests/testthat/test_chainof5parts.R @@ -162,7 +162,7 @@ test_that("chainof5parts", { expect_true(dir.exists(dirname)) expect_true(file.exists(rn[1])) expect_that(nrow(output),equals(4)) # 2023-09-04: changed because part5 now gives also first and last day - expect_that(ncol(output),equals(160)) # changed because intensity gradient now included in the test + expect_that(ncol(output),equals(187)) # changed because intensity gradient now included in the test expect_that(round(as.numeric(output$wakeup[2]), digits = 4), equals(36)) dirname_raw = "output_test/meta/ms5.outraw/40_100_400" rn2 = dir(dirname_raw,full.names = TRUE, recursive = T) diff --git a/tests/testthat/test_load_check_params.R b/tests/testthat/test_load_check_params.R index 4ee135ad1..4f5bffac4 100644 --- a/tests/testthat/test_load_check_params.R +++ b/tests/testthat/test_load_check_params.R @@ -10,7 +10,7 @@ test_that("load_params can load parameters", { # Test length of objects expect_equal(length(params), 8) - expect_equal(length(params$params_sleep), 21) + expect_equal(length(params$params_sleep), 22) expect_equal(length(params$params_metrics), 41) expect_equal(length(params$params_rawdata), 37) expect_equal(length(params$params_247), 20) diff --git a/tests/testthat/test_part5_analyseRest.R b/tests/testthat/test_part5_analyseRest.R new file mode 100644 index 000000000..1115e75d8 --- /dev/null +++ b/tests/testthat/test_part5_analyseRest.R @@ -0,0 +1,100 @@ +library(GGIR) +context("g.part5.analyseRest") +tz = "Europe/Amsterdam" +ds_names = rep("", 40) +dsummary = matrix("", 1, 40) +startday = as.POSIXlt(x = "2022-06-02 08:00:00", tz = tz) +time = seq(startday, startday + (16 * 3600), by = 60) + + + +test_that("Overlap 1 nap and 1 sib", { + fi = 1 + di = 1 + sibreport = data.frame(ID = rep("test123", 2), type = c("nap", "sib"), + start = c("2022-06-02 14:00:00", "2022-06-02 14:05:00"), + end = c("2022-06-02 14:20:00", "2022-06-02 14:20:00")) + sibreport$duration = as.numeric(difftime(time1 = sibreport$end, time2 = sibreport$start, units = "mins", tz = tz)) + restAnalyses = g.part5.analyseRest(sibreport = sibreport, dsummary = dsummary, + ds_names = ds_names, fi = fi, di = di, + time = time, + tz = tz) + fi = restAnalyses$fi + di = restAnalyses$di + dsummary = restAnalyses$dsummary + ds_names = restAnalyses$ds_names + dsummary = as.numeric(dsummary[, which(ds_names != "")]) + ds_names = ds_names[ds_names != ""] + names(dsummary) = ds_names + dsummary = as.data.frame(t(dsummary)) + + expect_equal(dsummary$nbouts_day_sib, 1) + expect_equal(dsummary$nbouts_day_srnap, 1) + expect_equal(dsummary$frag_mean_dur_sib_day, 15) + expect_equal(dsummary$frag_mean_dur_srnap_day, 20) + expect_equal(dsummary$perc_sib_overl_srnap, 100) + expect_equal(dsummary$perc_srnap_overl_sib, 75) + expect_equal(sum(dsummary), 323) +}) + +test_that("Overlap 1 nonwear and 1 sib", { + fi = 1 + di = 1 + sibreport = data.frame(ID = rep("test123", 2), type = c("nonwear", "sib"), + start = c("2022-06-02 14:00:00", "2022-06-02 14:05:00"), + end = c("2022-06-02 14:20:00", "2022-06-02 14:20:00")) + sibreport$duration = as.numeric(difftime(time1 = sibreport$end, time2 = sibreport$start, units = "mins", tz = tz)) + restAnalyses = g.part5.analyseRest(sibreport = sibreport, dsummary = dsummary, + ds_names = ds_names, fi = fi, di = di, + time = time, + tz = tz) + fi = restAnalyses$fi + di = restAnalyses$di + dsummary = restAnalyses$dsummary + ds_names = restAnalyses$ds_names + dsummary = as.numeric(dsummary[, which(ds_names != "")]) + ds_names = ds_names[ds_names != ""] + names(dsummary) = ds_names + dsummary = as.data.frame(t(dsummary)) + + expect_equal(dsummary$nbouts_day_sib, 1) + expect_equal(dsummary$nbouts_day_srnonw, 1) + expect_equal(dsummary$frag_mean_dur_sib_day, 15) + expect_equal(dsummary$frag_mean_dur_srnonw_day, 20) + expect_equal(dsummary$perc_sib_overl_srnonw, 100) + expect_equal(dsummary$perc_srnonw_overl_sib, 75) + expect_equal(sum(dsummary), 323) +}) + + +test_that("No overlap 1 nonwear, 1 nap, and 1 sib", { + fi = 1 + di = 1 + sibreport = data.frame(ID = rep("test123", 3), type = c("nonwear", "nap", "sib"), + start = c("2022-06-02 12:00:00", "2022-06-02 13:00:00", "2022-06-02 15:00:00"), + end = c("2022-06-02 12:20:00", "2022-06-02 13:20:00", "2022-06-02 15:20:00")) + sibreport$duration = as.numeric(difftime(time1 = sibreport$end, time2 = sibreport$start, units = "mins", tz = tz)) + restAnalyses = g.part5.analyseRest(sibreport = sibreport, dsummary = dsummary, + ds_names = ds_names, fi = fi, di = di, + time = time, + tz = tz) + fi = restAnalyses$fi + di = restAnalyses$di + dsummary = restAnalyses$dsummary + ds_names = restAnalyses$ds_names + dsummary = as.numeric(dsummary[, which(ds_names != "")]) + ds_names = ds_names[ds_names != ""] + names(dsummary) = ds_names + dsummary = as.data.frame(t(dsummary)) + + expect_equal(dsummary$nbouts_day_sib, 1) + expect_equal(dsummary$nbouts_day_srnap, 1) + expect_equal(dsummary$nbouts_day_srnonw, 1) + expect_equal(dsummary$frag_mean_dur_sib_day, 20) + expect_equal(dsummary$frag_mean_dur_srnap_day, 20) + expect_equal(dsummary$frag_mean_dur_srnonw_day, 20) + expect_equal(dsummary$perc_sib_overl_srnap, 0) + expect_equal(dsummary$perc_sib_overl_srnonw, 0) + expect_equal(dsummary$perc_srnonw_overl_sib, 0) + expect_equal(sum(dsummary), 129) +}) diff --git a/tests/testthat/test_sibreport.R b/tests/testthat/test_sibreport.R index dfa0262f6..7b2680479 100644 --- a/tests/testthat/test_sibreport.R +++ b/tests/testthat/test_sibreport.R @@ -2,22 +2,29 @@ library(GGIR) context("g.sibreportc") test_that("g.sibreport creates expected object", { ts = data.frame(sibdetection = c(0, 0, 1, 1, 1, 0, 0), - diur = rep(0, 7), - time = seq(as.POSIXlt(x = "2021-3-3 15:00:00", tz="Europe/Amsterdam"), - as.POSIXlt(x = "2021-3-3 15:06:00", tz="Europe/Amsterdam"), by = 60), - ACC = 1:7, angle = 7:1) + diur = rep(0, 7), + time = seq(as.POSIXlt(x = "2021-3-3 15:00:00", tz = "Europe/Amsterdam"), + as.POSIXlt(x = "2021-3-3 15:06:00", tz = "Europe/Amsterdam"), + by = 60), + ACC = 1:7, angle = 7:1) ID = 12345 epochlength = 60 - nonwearlog = data.frame(ID=rep(ID, 5), date=seq(as.Date("2021-3-3"),as.Date("2021-3-7"), by=1), - nonwear1=rep("15:00:00", 5), nonwear2=rep("15:30:00", 5)) - naplog = data.frame(ID=rep(ID, 2), date=seq(as.Date("2021-3-3"),as.Date("2021-3-4"), by=1), - nonwear1=rep("13:00:00", 2), nonwear2=rep("13:30:00", 2)) + nonwearlog = data.frame(ID = rep(ID, 5), + date = seq(as.Date("2021-3-3"), as.Date("2021-3-7"), by = 1), + nonwear1 = rep("15:00:00", 5), + nonwear2 = rep("15:30:00", 5)) + naplog = data.frame(ID = rep(ID, 2), + date = seq(as.Date("2021-3-3"),as.Date("2021-3-4"), by = 1), + nonwear1 = rep("13:00:00", 2), + nonwear2 = rep("13:30:00", 2)) + + logs_diaries = list(sleeplog = c(), nonwearlog = nonwearlog, + naplog = naplog, dateformat = "%Y-%m-%d") - logs_diaries = list(sleeplog=c(), nonwearlog=nonwearlog, - naplog=naplog, dateformat="%Y-%m-%d") SIBREPORT = g.sibreport(ts, ID, epochlength, logs_diaries) + expect_equal(nrow(SIBREPORT), 8) - expect_equal(mean(SIBREPORT$duration, na.rm = TRUE), 3) + expect_equal(mean(SIBREPORT$duration, na.rm = TRUE), 26.625) expect_equal(mean(SIBREPORT$mean_acc_1min_before, na.rm = TRUE), 2) expect_equal(mean(SIBREPORT$mean_acc_1min_after, na.rm = TRUE), 6) expect_equal(SIBREPORT$start[8], "2021-03-03 15:02:00") diff --git a/vignettes/GGIR.Rmd b/vignettes/GGIR.Rmd index 15a213d22..c6712fc60 100644 --- a/vignettes/GGIR.Rmd +++ b/vignettes/GGIR.Rmd @@ -438,7 +438,8 @@ comes with the following changes: - You can add columns relating to self-reported napping time and nonwear time. These are not used for the sleep analysis in g.part3 and g.part4, but used in g.part5 to facilitate napping analysis, see argument - `do.sibreport`. Multiple naps and multiple nonwear periods can be + `do.sibreport` and the paragraph on naps. Multiple naps and multiple + nonwear periods can be entered per day. - Leave cells for missing values blank. - Column names are critical for the advanced sleeplog format: Date @@ -975,6 +976,16 @@ vigorous activity: | `W0.5` | Derived from power law exponent alpha, see [Chastin et al. 201 0](https://www.sciencedirect.com/science/article/abs/pii/S096663620900602X?via%3Dihub) | | nap_count | Total number of naps, only calculated when argument do.sibreport = TRUE, currently optimised for 3.5-year olds. See function documentation for function `g.part5.classifyNaps` in the [GGIR function documentation (pdf)](https:/CRAN.R-project.org/package=GGIR/GGIR.pdf). | | nap_totalduration | Total nap duration, only calculated when argument do.sibreport = TRUE, currently optimised for 3.5-year old. See function documentation for function `g.part5.classifyNaps` in the [GGIR function documentation (pdf)](https:/CRAN.R-project.org/package=GGIR/GGIR.pdf). | +| sibreport_n_items | Only created if `do.sibreport = TRUE`. Number of items in the sibreport | +| sibreport_n_items_day | Only created if `do.sibreport = TRUE`. Number of items in the sibreport for this specific day | +| nbouts_day_X | Only created if `do.sibreport = TRUE`. Number of bouts in a day of X where X can be sib (sustained inactivity bout), srnap (self-reported nap) or srnonw (self-reported nonwear) | +| noverl_X | Only created if `do.sibreport = TRUE`. Number of overlapping bouts in a day of X where X can be sib_srnap, sib_srnonw, srnap_sib, or srnonw_sib | +| frag_mean_dur_X_day | Only created if `do.sibreport = TRUE`. Mean duration of X per day, where X can be sib, srnap or srnonw | +| dur_day_X_min | Only created if `do.sibreport = TRUE`. Total duration in minutes of X per day, where X can be sib, srnap or srnonw | +| mdur_X_overl_Y | Only created if `do.sibreport = TRUE`. Mean duration of the overlap between X and Y, which are combinations of sib, srnap or srnonw | +| tdur_X_overl_Y | Only created if `do.sibreport = TRUE`. Total duration in minutes of the overlap between X and Y, which are combinations of sib, srnap or srnonw | +| perc_X_overl_Y | Only created if `do.sibreport = TRUE`. Percentage of overlap between X and Y, which are combinations of sib, srnap or srnonw | + **Special note if you are working on compositional data analysis:**