Skip to content

Commit

Permalink
Merge branch 'issue708_centralise_cosinor' into issue906_create_part6
Browse files Browse the repository at this point in the history
  • Loading branch information
vincentvanhees committed Sep 26, 2023
2 parents 10e5f35 + 90a3946 commit 0ec4efa
Show file tree
Hide file tree
Showing 10 changed files with 177 additions and 113 deletions.
2 changes: 1 addition & 1 deletion NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ export(g.analyse, g.calibrate,
check_params, extract_params,
g.imputeTimegaps, g.part5.classifyNaps,
tidyup_df, cosinorAnalyses,
detect_nonwear_clipping,
detect_nonwear_clipping, applyCosinorAnalyses,
ShellDoc2Vignette, parametersVignette,
correctOlderMilestoneData,
convertEpochData, appendRecords, extractID,
Expand Down
4 changes: 4 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
72 changes: 72 additions & 0 deletions R/applyCosinorAnalyses.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,72 @@
applyCosinorAnalyses = function(ts, qcheck, midnightsi, epochsizes) {
# qcheck - vector of length ts to indicate invalid values
ws2 = epochsizes[2]
ws3 = epochsizes[1]
# Re-derive Xi but this time include entire time series
# Here, we ignore duplicated values (when clock moves backward due to DST)
handleDST = !duplicated(ts)
qcheck = qcheck[handleDST]
Xi = ts[handleDST, grep(pattern = "time", x = colnames(ts), invert = TRUE)]
Nlong_epochs_day = (1440 * 60) / ws2 # this is 96 by default
dstgap = which(diff(midnightsi) != Nlong_epochs_day)
if (length(dstgap) > 0) {
# Time moved forward due to DST
gaplocation = ((midnightsi[dstgap[1]] * ws2) / ws3) + (2 * (3600/ws3))
# Insert NA values
Xi = c(Xi[1:gaplocation], rep(NA, 3600/ws3), Xi[(gaplocation + 1):length(Xi)])
qcheck = c(qcheck[1:gaplocation], rep(NA, 3600/ws3), qcheck[(gaplocation + 1):length(qcheck)])
}

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

cosinor_coef = cosinorAnalyses(Xi = Xi, epochsize = epochsize, timeOffsetHours = timeOffsetHours)
cosinor_coef$timeOffsetHours = timeOffsetHours
} else {
cosinor_coef = c()
}
return(cosinor_coef)
}
11 changes: 11 additions & 0 deletions R/cosinorAnalyses.R
Original file line number Diff line number Diff line change
Expand Up @@ -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))
}
72 changes: 3 additions & 69 deletions R/g.analyse.avday.R
Original file line number Diff line number Diff line change
Expand Up @@ -120,75 +120,9 @@ g.analyse.avday = function(doquan, averageday, M, IMP, t_TWDI, quantiletype,
#----------------------------------
# (Extended) Cosinor analysis
if (params_247[["cosinor"]] == TRUE) {
# Re-derive Xi but this time include entire time series
# Here, we ignore duplicated values (when clock moves backward due to DST)
handleDST = !duplicated(IMP$metashort)
qcheck = qcheck[handleDST]
Xi = IMP$metashort[handleDST, acc.metric]
Nlong_epochs_day = (1440 * 60) / ws2 # this is 96 by default
dstgap = which(diff(midnightsi) != Nlong_epochs_day)
if (length(dstgap) > 0) {
# Time moved forward due to DST
gaplocation = ((midnightsi[dstgap[1]] * ws2) / ws3) + (2 * (3600/ws3))
# Insert NA values
Xi = c(Xi[1:gaplocation], rep(NA, 3600/ws3), Xi[(gaplocation + 1):length(Xi)])
qcheck = c(qcheck[1:gaplocation], rep(NA, 3600/ws3), qcheck[(gaplocation + 1):length(qcheck)])
}

# Xi = log((Xi * 1000) + 1) # log transformed to be more robust against peaks in the data
# set non-wear to missing values, because for Cosinor fit
# it seems more logical to only fit with real data
# this comes at the price of not being able to extract F_pseudo
firstvalid = 1
if (length(which(qcheck == 1)) > 0) {
is.na(Xi[which(qcheck == 1)]) = TRUE
# ignore invalid start of recording (if applicable)
# such that 24 hour blocks start from first valid value
firstvalid = which(qcheck == 0)[1]
if (is.na(firstvalid) == FALSE) {
if (firstvalid != 1) {
Xi = Xi[firstvalid:length(Xi)]
}
}
}
if (length(which(is.na(Xi) == FALSE)) > (1440 * (60/ws3))) { # Only attempt cosinor analyses if there is more than 24 hours of data
midnightsi_ws3 = (midnightsi - 1) * (ws2 / ws3)
timeOffsetHours = (midnightsi_ws3[which(midnightsi_ws3 >= firstvalid - 1)[1]] - (firstvalid - 1)) / (3600 / ws3)
if (ws3 < 60) {
# If epochsize < 1 minute then aggregate to 1 minute by taking maximum value
# but keep NA values
XTtime = rep(1:length(Xi), each = 60 / ws3)
XT = data.frame(Xi = Xi, time = XTtime[1:length(Xi)])
custommean = function(x) {
y = NA
if (length(x) > 0) {
if (length(which(is.na(x) == FALSE) ) > 0) {
y = mean(x, na.rm = TRUE)
}
}
return(y)
}
XT = aggregate(x = XT, by = list(XT$time), FUN = custommean)
if (length(which(is.nan(XT$Xi) == TRUE)) > 0) {
is.na(XT$Xi[which(is.nan(XT$Xi) == TRUE)]) = TRUE
}
# experimental: clip all peaks above Xth percentile?
# Q9 = quantile(x = XT$Xi, probs = 0.75, na.rm = TRUE)
# XT$Xi[which(XT$Xi >= Q9)] = Q9

# log transform of data in millig
notna = !is.na(XT$Xi)
XT$Xi[notna] = log((XT$Xi[notna]*1000) + 1)
Xi = XT$Xi
epochsize = 60
} else {
epochsize = ws3
}
cosinor_coef = cosinorAnalyses(Xi = Xi, epochsize = epochsize, timeOffsetHours = timeOffsetHours)
cosinor_coef$timeOffsetHours = timeOffsetHours
} else {
cosinor_coef = c()
}
cosinor_coef = applyCosinorAnalyses(ts = IMP$metashort[, c("timestamp", acc.metric)],
qcheck = qcheck,
midnightsi, epochsizes = c(ws3, ws2))
} else {
cosinor_coef = c()
}
Expand Down
24 changes: 14 additions & 10 deletions R/g.analyse.perfile.R
Original file line number Diff line number Diff line change
Expand Up @@ -156,15 +156,16 @@ g.analyse.perfile = function(I, C, metrics_nav,
filesummary[vi] = c(cosinor_coef$timeOffsetHours)
s_names[vi] = c("cosinor_timeOffsetHours")
vi = vi + 1
try(expr = {filesummary[vi:(vi + 4)] = as.numeric(c(cosinor_coef$coef$params$mes,
try(expr = {filesummary[vi:(vi + 5)] = as.numeric(c(cosinor_coef$coef$params$mes,
cosinor_coef$coef$params$amp,
cosinor_coef$coef$params$acr,
cosinor_coef$coef$params$acrotime,
cosinor_coef$coef$params$ndays))}, silent = TRUE)
s_names[vi:(vi + 4)] = c("cosinor_mes", "cosinor_amp", "cosinor_acrophase",
"cosinor_acrotime", "cosinor_ndays")
vi = vi + 5
try(expr = {filesummary[vi:(vi + 9)] = c(cosinor_coef$coefext$params$minimum,
cosinor_coef$coef$params$ndays,
cosinor_coef$coef$params$R2))}, silent = TRUE)
s_names[vi:(vi + 5)] = c("cosinor_mes", "cosinor_amp", "cosinor_acrophase",
"cosinor_acrotime", "cosinor_ndays", "cosinor_R2")
vi = vi + 6
try(expr = {filesummary[vi:(vi + 10)] = c(cosinor_coef$coefext$params$minimum,
cosinor_coef$coefext$params$amp,
cosinor_coef$coefext$params$alpha,
cosinor_coef$coefext$params$beta,
Expand All @@ -173,16 +174,19 @@ g.analyse.perfile = function(I, C, metrics_nav,
cosinor_coef$coefext$params$DownMesor,
cosinor_coef$coefext$params$MESOR,
cosinor_coef$coefext$params$ndays,
cosinor_coef$coefext$params$F_pseudo)}, silent = TRUE)
s_names[vi:(vi + 9)] = c("cosinorExt_minimum", "cosinorExt_amp", "cosinorExt_alpha",
cosinor_coef$coefext$params$F_pseudo,
cosinor_coef$coefext$params$R2)}, silent = TRUE)
s_names[vi:(vi + 10)] = c("cosinorExt_minimum", "cosinorExt_amp", "cosinorExt_alpha",
"cosinorExt_beta", "cosinorExt_acrotime", "cosinorExt_UpMesor",
"cosinorExt_DownMesor", "cosinorExt_MESOR",
"cosinorExt_ndays", "cosinorExt_F_pseudo")
vi = vi + 10
"cosinorExt_ndays", "cosinorExt_F_pseudo", "cosinorExt_R2")
vi = vi + 11
filesummary[vi:(vi + 1)] = c(cosinor_coef$IVIS$InterdailyStability,
cosinor_coef$IVIS$IntradailyVariability)
s_names[vi:(vi + 1)] = c("cosinorIS", "cosinorIV")
vi = vi + 2
} else {
vi = vi + 20
}

# Variables per metric - summarise with stratification to weekdays and weekend days
Expand Down
27 changes: 16 additions & 11 deletions R/g.part4.R
Original file line number Diff line number Diff line change
Expand Up @@ -62,11 +62,15 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1,
# initialize output variable names
colnamesnightsummary = c("ID", "night", "sleeponset", "wakeup", "SptDuration", "sleepparam", "guider_onset",
"guider_wakeup", "guider_SptDuration", "error_onset", "error_wake", "error_dur", "fraction_night_invalid",
"SleepDurationInSpt", "WASO", "duration_sib_wakinghours", "number_sib_sleepperiod", "number_of_awakenings",
"number_sib_wakinghours", "duration_sib_wakinghours_atleast15min", "sleeponset_ts", "wakeup_ts", "guider_onset_ts",
"guider_wakeup_ts", "sleeplatency", "sleepefficiency", "page", "daysleeper", "weekday", "calendar_date",
"SleepDurationInSpt", "WASO", "duration_sib_wakinghours",
"number_sib_sleepperiod", "number_of_awakenings",
"number_sib_wakinghours", "duration_sib_wakinghours_atleast15min",
"sleeponset_ts", "wakeup_ts", "guider_onset_ts", "guider_wakeup_ts",
"sleeplatency", "sleepefficiency", "page", "daysleeper", "weekday", "calendar_date",
"filename", "cleaningcode", "sleeplog_used", "sleeplog_ID", "acc_available", "guider", "SleepRegularityIndex", "SriFractionValid",
"longitudinal_axis")
"longitudinal_axis") #


if (params_output[["storefolderstructure"]] == TRUE) {
colnamesnightsummary = c(colnamesnightsummary, "filename_dir", "foldername")
}
Expand Down Expand Up @@ -344,10 +348,10 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1,
# initialize dataframe to hold sleep period overview:
spocum = data.frame(nb = numeric(0), start = numeric(0), end = numeric(0),
dur = numeric(0), def = character(0))

spocumi = 1 # counter for sleep periods
# continue now with the specific data of the night

guider.df2 = guider.df[which(guider.df$night == j), ]
# ================================================================================ get
# sleeplog (or HDCZA or L5+/-6hr algorithm) onset and waking time and assess whether it
Expand Down Expand Up @@ -457,6 +461,7 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1,
} else {
acc_available = TRUE
}

if (nrow(sleepdet) == 0) next
ki = which(sleepdet$definition == k)
if (length(ki) == 0) next
Expand Down Expand Up @@ -681,13 +686,13 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1,
}
delta_t1 = diff(as.numeric(spocum.t$end))
spocum.t$dur = correct01010pattern(spocum.t$dur)

#----------------------------
nightsummary[sumi, 1] = accid
nightsummary[sumi, 2] = j #night
# remove double rows
spocum.t = spocum.t[!duplicated(spocum.t), ]

#------------------------------------
# ACCELEROMETER
if (length(which(as.numeric(spocum.t$dur) == 1)) > 0) {
Expand Down Expand Up @@ -868,6 +873,7 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1,
spocum.t.dur_sibd = 0
spocum.t.dur_sibd_atleast15min = 0
}

nightsummary[sumi, 14] = spocum.t.dur.noc #total nocturnalsleep /accumulated sleep duration
nightsummary[sumi, 15] = nightsummary[sumi, 5] - spocum.t.dur.noc #WASO
nightsummary[sumi, 16] = spocum.t.dur_sibd #total sib (sustained inactivty bout) duration during wakinghours
Expand Down Expand Up @@ -1043,11 +1049,10 @@ g.part4 = function(datadir = c(), metadatadir = c(), f0 = f0, f1 = f1,
if (length(nnights.list) == 0) {
# if there were no nights to analyse
nightsummary[sumi, 1:2] = c(accid, 0)
nightsummary[sumi, 3:30] = NA
nightsummary[sumi, c(3:30, 34, 36:39)] = NA
nightsummary[sumi, 31] = fnames[i]
nightsummary[sumi, 32] = 4 #cleaningcode = 4 (no nights of accelerometer available)
nightsummary[sumi, c(33,35)] = c(FALSE, TRUE) #sleeplog_used acc_available
nightsummary[sumi, 36:39] = NA
nightsummary[sumi, c(33, 35)] = c(FALSE, TRUE) #sleeplog_used acc_available
if (params_output[["storefolderstructure"]] == TRUE) {
nightsummary[sumi, 40:41] = c(ffd[i], ffp[i]) #full filename structure and use the lowest foldername as foldername name
}
Expand Down
20 changes: 12 additions & 8 deletions R/g.report.part4.R
Original file line number Diff line number Diff line change
Expand Up @@ -295,13 +295,18 @@ g.report.part4 = function(datadir = c(), metadatadir = c(), loglocation = c(),
if (dotwice == 1) {
nightsummary.tmp = turn_numeric(x = nightsummary.tmp, varnames = gdn)
}
varnames_tmp = c("SptDuration", "sleeponset",
"wakeup", "WASO", "SleepDurationInSpt",
"number_sib_sleepperiod", "duration_sib_wakinghours",
"number_of_awakenings", "number_sib_wakinghours",
"duration_sib_wakinghours_atleast15min",
"sleeplatency", "sleepefficiency", "number_of_awakenings",
"guider_inbedDuration", "guider_inbedStart",
"guider_inbedEnd", "guider_SptDuration", "guider_onset",
"guider_wakeup", "SleepRegularityIndex",
"SriFractionValid")
nightsummary.tmp = turn_numeric(x = nightsummary.tmp,
varnames = c("SptDuration", "sleeponset",
"wakeup", "WASO", "SleepDurationInSpt", "number_sib_sleepperiod", "duration_sib_wakinghours",
"number_of_awakenings", "number_sib_wakinghours", "duration_sib_wakinghours_atleast15min",
"sleeplatency", "sleepefficiency", "number_of_awakenings", "guider_inbedDuration", "guider_inbedStart",
"guider_inbedEnd", "guider_SptDuration", "guider_onset", "guider_wakeup", "SleepRegularityIndex",
"SriFractionValid"))
varnames = varnames_tmp)
weekday = nightsummary.tmp$weekday[this_sleepparam]
if (dotwice == 1) {
for (k in 1:3) {
Expand Down Expand Up @@ -428,8 +433,7 @@ g.report.part4 = function(datadir = c(), metadatadir = c(), loglocation = c(),
udefn[j], "_mn", sep = ""), paste("average_dur_sib_wakinghours_", TW, "_", udefn[j],
"_sd", sep = ""))
NDAYsibd = length(which(nightsummary.tmp$number_sib_wakinghours[indexUdef] > 0))
if (length(NDAYsibd) == 0)
NDAYsibd = 0
if (length(NDAYsibd) == 0) NDAYsibd = 0
personSummary[i, (cnt + 19)] = NDAYsibd
personSummarynames = c(personSummarynames, paste("n_days_w_sib_wakinghours_", TW, "_",
udefn[j], sep = ""))
Expand Down
Loading

0 comments on commit 0ec4efa

Please sign in to comment.