R/summarise.R
In vincentvanhees/ActChronicFatigue: Process accelerometer data from individuals with chronic fatigue
Defines functions
summarise
Documented in
summarise
#' summarise
#'
#' @param outputdir ath to output directory generated by GGIR
#' @param part5_summary Content of GGIR part 5 personsummary file
#' @param model_threshold Acceleration threshold corresponding to the model, which can be derived as negative of the first coefficient divided by the second coefficient.
#' @param referentiewaarden Mean and standard deviation ENMO_fullRecordingMean in a reference population.
#' @param sleepwindowType See, GGIR documentation.
#' @param MVPAdefinition MVPA variable name from the GGIR part2 daysummary output
#' @param threshold_wrist Threshold to classify active or passive
#' @param sep separator used by csv files stored by GGIR
#' @param dec dec separator used by csv files stored by GGIR
#' @return no object is returned, only a summary is printed to the screen
#' @export
#' @importFrom grDevices dev.off pdf adjustcolor
#' @importFrom graphics abline axis barplot par grid
#' @importFrom graphics legend lines text mtext rasterImage
#' @importFrom utils packageVersion
summarise = function(outputdir,
part5_summary,
model_threshold = 75,
referentiewaarden,
sleepwindowType="SPT",
MVPAdefinition = "MVPA_E5S_B10M80%_T100_BFEN_0-24hr",
threshold_wrist = 0.5,
sep = sep, dec = dec) {
part2_summary_file = grep(dir(paste0(outputdir,"/results"), full.names = TRUE),
pattern = "part2_summary", value = T)
part2_summary = read.csv(file = part2_summary_file, sep = sep, dec = dec)
extractid = function(x) {
tmp = unlist(strsplit(x," "))
if (length(tmp) > 1) {
tmp = tmp[1]
} else {
tmp = x
}
return(tmp)
}
if (is.character(part2_summary$ID) == TRUE) {
part2_summary$ID = sapply(part2_summary$ID, FUN = extractid)
}
part2_summary = part2_summary[,c("ID", "ENMO_fullRecordingMean")] # "BFEN_fullRecordingMean" "ENMO_fullRecordingMean"
colnames(part2_summary) = c("ID","Activity_zscore")
part2_summary$Activity_zscore = round((part2_summary$Activity - referentiewaarden[1]) / referentiewaarden[2], digits = 1) # - 23.64) / 15.5 for ENMO
part5_summary = merge(part5_summary, part2_summary, by.x = "ID2", by.y = "ID")
# check for which IDs reports have already been generated:
filesinresults = dir(paste0(outputdir, "/results"))
report_names = grep(pattern = "eweeg_en_slaap_rapport_", x = filesinresults, value = TRUE)
IDdone = c()
if (length(report_names) > 0) {
extractID = function(x) {
return(unlist(strsplit(x, "_rapport_|[.]pd"))[2])
}
IDdone = unlist(lapply(X = report_names, FUN = extractID))
}
if (length(IDdone) > 0) {
most_recent_recordings = which(part5_summary$ID2 %in% IDdone == FALSE)
} else {
most_recent_recordings = 1:nrow(part5_summary)
}
if (length(most_recent_recordings) == 0) {
cat("\nDe bestanden waren al eerder genalyseerd en de rapporten staan in de resultaten. Mocht je de rapporten opnieuw willen creeeren, verwijder dan eerst de oude pdf bestanden.")
} else {
cat(paste0("Klassificaties voor meest recente ", length(most_recent_recordings), " metingen: "))
part5_summary = part5_summary[order(as.Date(part5_summary$calendar_date)),]
if ("sleep_efficiency_after_onset_pla" %in% colnames(part5_summary)) {
colnames(part5_summary)[which(colnames(part5_summary) == "sleep_efficiency_after_onset_pla")] = "sleep_efficiency_pla"
}
recent_recording = part5_summary[most_recent_recordings,
c("ID2", "calendar_date", "prop_perv_passive", "Activity_zscore",
"dur_spt_sleep_min_pla", "dur_spt_min_pla", "sleep_efficiency_pla",
"sleeponset_pla", "wakeup_pla")]
recent_recording$dur_spt_sleep_min_pla = round(recent_recording$dur_spt_sleep_min_pla / 60, digits = 2)
recent_recording$dur_spt_min_pla = round(recent_recording$dur_spt_min_pla / 60, digits = 2)
recent_recording$sleep_efficiency_pla = round(recent_recording$sleep_efficiency_pla * 100, digits = 0)
recent_recording$sleeponset_pla = round(recent_recording$sleeponset_pla - 24, digits = 2)
recent_recording$wakeup_pla = round(recent_recording$wakeup_pla - 24, digits = 2)
recent_recording$prop_perv_passive = round(recent_recording$prop_perv_passive * 100)
recent_recording$classification = "Actief"
pp = which(recent_recording$prop_perv_passive >= threshold_wrist * 100)
if (length(pp) > 0) recent_recording$classification[pp] = "Laag Actief"
colnames(recent_recording) = c("ID", "Datum start meeting", "Kans op laag actief (%)",
"Gemiddelde beweging t.o.v. referentie groep (z-score)",
"Uren slaap per nacht", "Gem. duur slaapperioden (uren)",
"Slaap efficiency binnen slaapperioden (%)", "Gemiddelde slaap start t.o.v. middernacht",
"Gemiddelde slaap einde t.o.v. middernacht", "Klassificatie") #,"feedback")
varnames = c("Datum start meeting", "Klassificatie", "Kans op laag actief (%)",
"Gemiddelde beweging t.o.v. referentie groep (z-score)")
recent_recording = recent_recording[,c(1:4,10,5:9)]
cat("\n")
part5_daysummary_file = grep(dir(paste0(outputdir,"/results"), full.names = TRUE),
pattern = "part5_daysummary", value = T)
part5_daysummary = read.csv(file = part5_daysummary_file, sep = sep, dec = dec)
part2_daysummary_file = grep(dir(paste0(outputdir,"/results"), full.names = TRUE),
pattern = "part2_daysummary", value = T)
part2_daysummary = read.csv(file = part2_daysummary_file, sep = sep, dec = dec)
# Load sleep/wake times from both GGIR and diary
part4_nightsummaryfull_file = grep(dir(paste0(outputdir,"/results/QC"), full.names = TRUE),
pattern = "part4_nightsummary_sleep_full.csv", value = T)
part4_nightsummary = read.csv(part4_nightsummaryfull_file, sep = sep, dec = dec)
days = c("Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday", "Sunday")
dagen = c("maandag", "dinsdag", "woensdag", "donderdag", "vrijdag", "zaterdag", "zondag")
for (ki in 1:length(days)) {
part5_daysummary$weekday = gsub(x = part5_daysummary$weekday,
pattern = days[ki], replacement = dagen[ki])
part2_daysummary$weekday = gsub(x = part2_daysummary$weekday,
pattern = days[ki], replacement = dagen[ki])
part4_nightsummary$weekday = gsub(x = part4_nightsummary$weekday,
pattern = days[ki], replacement = dagen[ki])
}
for (ID in unique(recent_recording$ID)) {
# Get rows for this ID
P5D = part5_daysummary[which(part5_daysummary$ID == ID),]
P2D = part2_daysummary[which(part2_daysummary$ID == ID),]
P4N = part4_nightsummary[which(part4_nightsummary$ID == ID),]
# Put date in data format
P5D$calendar_date = as.Date(P5D$calendar_date, "%Y-%m-%d")
P2D$calendar_date = as.Date(P2D$calendar_date, "%Y-%m-%d")
P4N$calendar_date = as.Date(P4N$calendar_date, "%d/%m/%Y")
# Merge time spent in two acceleration range to sleep object, to ease plotting later on
# Combine P4 and P2 output to ease plotting
P4N = base::merge(P4N, P2D[,c("calendar_date", MVPAdefinition, "step_count_sum_0.24hr")],
by = c("calendar_date"), all.x = TRUE)
colnames(P4N)[which(colnames(P4N) %in% MVPAdefinition == TRUE)] = c("MVPA") #"SB", "LIPA",
# Remove problematic nights from the sleep results P4
missingSleeplog = which(P4N$guider != "sleeplog")
if (length(missingSleeplog) > 0) {
is.na(P4N[missingSleeplog, grep(pattern = "guider", x = colnames(P4N), invert = FALSE)]) = TRUE
}
sleepProblem = which(P4N$cleaningcode > 1 & P4N$guider == "sleeplog")
if (length(sleepProblem) > 0) {
is.na(P4N[sleepProblem, grep(pattern = "calendar|ID|night|guider|weekday|page|filename|cleaningcode",
x = colnames(P4N), invert = TRUE)]) = TRUE
}
insufficientValidData4 = which(P4N$fraction_night_invalid > 0.33)
if (length(insufficientValidData4) > 0) {
is.na(P4N[insufficientValidData4, grep(pattern = "calendar|ID|night|guider|weekday|page|filename|cleaningcode",
x = colnames(P4N), invert = TRUE)]) = TRUE
}
#---------------------------------------------
# add empty rows for missing days
P5D = P5D[order(P5D$calendar_date),]
dates = P5D$calendar_date
dates_theoretical = seq(min(dates), max(dates), by = 1)
missing_dates = dates_theoretical[which(dates_theoretical %in% dates == FALSE)]
if (length(missing_dates) > 0) {
newvalues = (nrow(P5D) + 1):(nrow(P5D) + length(missing_dates))
P5D[newvalues,] = NA
P5D$calendar_date[newvalues] = format(as.Date(missing_dates), "%Y-%m-%d")
P5D$weekday[newvalues] = weekdays(missing_dates, abbreviate = FALSE)
if (length(which(P5D$weekday %in% days == TRUE)) > 0) {
for (ki in 1:length(days)) {
P5D$weekday = gsub(x = P5D$weekday,
pattern = days[ki], replacement = dagen[ki])
}
}
}
P5D = P5D[order(P5D$calendar_date),] # order again now missing day has been added
P5D$calendar_date = format(as.Date(P5D$calendar_date), "%Y-%m-%d")
P5D$calendar_date = format(as.Date(P5D$calendar_date), "%d")
# Now shorten weekday
dagen_kort = c("ma", "di", "wo", "do", "vr", "za", "zo")
# dagen = c("maandag", "dinsdag", "woensdag", "donderdag", "vrijdag", "zaterdag", "zondag")
dagen = c("maandag", "dinsdag", "woensdag", "donderdag", "vrijdag", "zaterdag", "zondag")
for (ki in 1:length(days)) {
P5D$weekday = gsub(x = P5D$weekday,
pattern = dagen[ki], replacement = dagen_kort[ki])
P4N$weekday = gsub(x = P4N$weekday,
pattern = dagen[ki], replacement = dagen_kort[ki])
}
# Extract person summary for part 5
P5P = part5_summary[which(part5_summary$ID == ID),]
pdffile = paste0(outputdir, "/results/Beweeg_en_slaap_rapport_",ID,".pdf")
pdf(file = pdffile, paper = "a4r", width = 11.69, height = 8.27)
par(mfrow = c(4,1), mar = c(4, 4, 3, 0.5), oma = c(0,0,0,0)) # las = 3
img <- png::readPNG(system.file("extdata/Amsterdam_UMC_Logo2.png", package = "ActChronicFatigue"))
# Key facts
keystats = t(recent_recording[which(recent_recording$ID == ID), varnames])
colnames(keystats)[1] = paste0("ID: ", ID)
CXmain = 1.0
CXdays = 1.0
relprop = as.numeric(keystats[3])
if (keystats[2] != "Laag Actief") relprop = 100 - relprop
IDtitle = paste0("ID: ", ID, " | Start: ", keystats[1])
title = paste0("Klassificatie: ", keystats[2], " (zekerheid ", relprop,"%) | ",
"Beweging (z-score) = ", keystats[4])
CL = 1.0
CXdots = 1.5
#========================================
# Plot activity time series (based on GGIR part 5 output)
NBars = 3
TS = P5D$ACC_day_mg
TS = c(TS, rep(NA, NBars))
plotlogo = function() {
grid::grid.raster(img, x = .91, y = 0.98, width = .18, interpolate = FALSE) # print homer in ll conrner
grid::grid.text(x = 0.91, y = 0.93,
gp = grid::gpar(fontsize = 7, col = "black", font = 2),
label = "Afdeling Medische Psychologie, NKCV") #, xpd = NA, col ="red", cex= 1, pos = 4)
}
plot(TS, type = "b", pch = 20, ylim = range(c(model_threshold * 2, P5D$ACC_day_mg), na.rm = T),
main = "", xlab = "", ylab = "Beweging per dag", bty = "l", axes = FALSE, cex.main = 1.5,
cex.lab = CL, cex = CXdots)
mtext(title, side = 3, adj = 0, line = 1.2, cex = 0.9, font = 2);
axis(side = 1, at = 1:(length(TS) - NBars), labels = c(paste0(P5D$weekday, " ", P5D$calendar_date)), #, rep("", 3)
cex.axis = CXdays)
lines(x = c(1,length(TS)), y = rep(model_threshold, 2), col = "black", lty = 2, lwd = 1.3)
text(x = length(TS) - 3, y = model_threshold + 10, labels = "Grenswaarde laag actief", pos = 4, cex = 1)
plotlogo()
#========================================
# Plot sleep duration (based on GGIR part 4 output)
if (nrow(P4N) > nrow(P5D) & nrow(P4N) > 1) {
P4N = P4N[2:nrow(P4N),]
}
if (sleepwindowType == "TimeInBed") {
A = P4N[,c("SleepDurationInSpt" , "SptDuration", "sleeplatency", "guider_inbedDuration", "sleepefficiency")]
title = "Slapen en wakker liggen (met slaap efficientie (%) aangeduid in wit)" #Slaap (donker), Wakker na start slaap (licht), Slaap latency (lichter), en
maxvalue = max(c(A[,2] + A[,3]), na.rm = T)
} else {
A = P4N[,c("SleepDurationInSpt" , "SptDuration")] #/ 60
title = "Slaap" #Slaap (donker), Wakker na start slaap (licht)
maxvalue = max(A[,2], na.rm = T)
}
A[,2] = A[,2] - A[,1]
if (sleepwindowType == "TimeInBed") {
B = A[,c(1,2, 3)]
} else {
B = A[,c(1,2)]
}
B[(nrow(B) + 1):(nrow(B) + NBars),] <- NA
brpos = barplot(t(as.matrix(B)), space = rep(0, nrow(B)), ylab = "Tijd in uren", beside = FALSE, #space = c(0, 0.2),
ylim = c(0, maxvalue + 6), cex.axis = 0.9,
names.arg = c(P4N$weekday, rep("",nrow(B) - nrow(P4N))),
cex.names = CXdays, cex.lab = CL, cex.main = CXmain,
legend.text = c("Slaap", "Wakker na in slaap vallen", "Wakker voor in slaap vallen"),
main = "",
args.legend = list(ncol = 1, x = nrow(B), cex = 1, bty = "n"))
mtext(title, side = 3, adj = 0, line = 1.2, cex = 0.9, font = 2);
if (sleepwindowType == "TimeInBed") {
text(brpos, rep(1, length(A$guider_inbedDuration)),
labels = round(100 * (A$sleepefficiency), digits = 0),
cex = 0.8,
col = "white")
}
#========================================
# Niet gedragen tijd
A = P5D[,c("nonwear_perc_day", "nonwear_perc_spt")]
A[(nrow(A) + 1):(nrow(A) + NBars),] <- NA
barplot(t(as.matrix(A)), ylab = "Percentage (%)", beside = TRUE, space = c(0, 0.2),
names.arg = c(P5D$weekday, rep("",NBars)), cex.names = CXdays, cex.lab = CL, cex.main = CXmain,
main = "", legend.text = c("overdag", "nacht"),
args.legend = list(x = (nrow(A)*2) + 1, ncol = 1, cex = 1, bty = "n"), ylim = c(0,100))
mtext("Beweegmeter niet gedragen?", side = 3, adj = 0, line = 1.2, cex = 0.9, font = 2);
#=======================================
# MVPA
maxvalue = max(c(P4N$MVPA, 30), na.rm = T) + 10
boutdur = unlist(strsplit(unlist(strsplit(MVPAdefinition, "B"))[2], "M"))[1]
title = paste0("Matig to Zwaar intensief gedrag in blokken van minimaal ", boutdur, " minuten")
PASB = P4N[,c("MVPA")]
PASB = c(PASB, matrix(NA,3, 1))
brpos = barplot(t(as.matrix(PASB)),
space = rep(0, length(PASB)),
ylab = "Tijd in minuten",
ylim = c(0, maxvalue), cex.axis = 0.9,
names.arg = c(P4N$weekday, rep("",NBars)),
cex.names = CXdays, cex.lab = CL, cex.main = CXmain,
main = "", col = "grey")
lines(x = c(0, length(P4N$weekday) + 3), y = c(30, 30), type = "l", lty = 2, lwd = 1.3)
mtext(title, side = 3, adj = 0, line = 1.2, cex = 0.9, font = 2);
text(x = length(P4N$weekday), y = 35, labels = "Beweegrichtlijn (30 min. per dag)", cex = 1, font = 1, pos = 4);
mtext(text = paste0("Software versie: ", packageVersion("ActChronicFatigue")),
side = 1, col = "black", cex = 0.7, line = 2.5, font = 2, adj = 0)
mtext(text = IDtitle,
side = 1, col = "black", cex = 0.7, line = 2.5, font = 2, adj = 1)
#========================================
# Slaap waak tijden
sleeponset_ts = P4N$sleeponset_ts # c("22:00", "23:20", "01:10")
wakeup_ts = P4N$wakeup_ts # c("6:55", "8:30", "10:40")
if (sleepwindowType == "SPT") {
sleeponset_log_ts = P4N$guider_onset_ts
wakeup_log_ts = P4N$guider_wakeup_ts
} else if (sleepwindowType == "TimeInBed") {
sleeponset_log_ts = P4N$guider_inbedStart_ts
wakeup_log_ts = P4N$guider_inbedEnd_ts
}
WV = which((is.na(wakeup_ts) == FALSE & is.na(sleeponset_ts) == FALSE) |
(is.na(wakeup_log_ts) == FALSE & is.na(sleeponset_log_ts) == FALSE))
sleeponset_ts = as.POSIXlt(sleeponset_ts, format = "%H:%M")
wakeup_ts = as.POSIXlt(wakeup_ts, format = "%H:%M")
sleeponset_log_ts = as.POSIXlt(sleeponset_log_ts, format = "%H:%M")
wakeup_log_ts = as.POSIXlt(wakeup_log_ts, format = "%H:%M")
for (j in 1:length(WV)) {
if (!is.na(sleeponset_ts[j])) {
if (sleeponset_ts[j] < as.POSIXlt("12:00", format = "%H:%M")) {
sleeponset_ts[j] = sleeponset_ts[j] + (24 * 3600)
}
}
if (!is.na(sleeponset_log_ts[j])) {
if (sleeponset_log_ts[j] < as.POSIXlt("12:00", format = "%H:%M")) {
sleeponset_log_ts[j] = sleeponset_log_ts[j] + (24 * 3600)
}
}
if (!is.na(wakeup_ts[j]) & !is.na(sleeponset_ts[j])) {
if (wakeup_ts[j] < sleeponset_ts[j] |
(wakeup_ts[j] > as.POSIXlt("00:00", format = "%H:%M")) &
(wakeup_ts[j] < as.POSIXlt("12:00", format = "%H:%M"))) {
wakeup_ts[j] = wakeup_ts[j] + (24 * 3600)
}
}
if (!is.na(wakeup_log_ts[j]) & !is.na(sleeponset_log_ts[j])) {
if (wakeup_log_ts[j] < sleeponset_log_ts[j] |
(wakeup_log_ts[j] > as.POSIXlt("00:00", format = "%H:%M")) &
(wakeup_log_ts[j] < as.POSIXlt("12:00", format = "%H:%M"))) {
wakeup_log_ts[j] = wakeup_log_ts[j] + (24 * 3600)
}
}
if (!is.na(wakeup_ts[j])) wakeup_ts[j] = wakeup_ts[j] - (24 * 3600)
if (!is.na(wakeup_log_ts[j])) wakeup_log_ts[j] = wakeup_log_ts[j] - (24 * 3600)
}
YLIM = c(as.POSIXlt("12:00", format = "%H:%M") - 12 * 3600, (as.POSIXlt("12:00", format = "%H:%M") + 24 * 3600))
timeaxis = c(as.numeric(YLIM[1]) + (c(0:36) * (3600)))
timeaxislabels = format(as.POSIXlt(timeaxis, origin = "1970-1-1"), "%H:%M")
waki = 1:(max(WV) - 1)
onsi = 2:max(WV)
# Check whether onset occurs before waking up, if yes then add 24 hour to onset
correctionNeeded = which(sleeponset_ts[onsi] < wakeup_ts[waki])
if (length(correctionNeeded) > 0) {
sleeponset_ts[onsi][correctionNeeded] = sleeponset_ts[onsi][correctionNeeded] + (24 * 3600)
}
COL = c("blue", "red", "black")
CL = 1
par(mfrow = c(1,1), oma = c(0,0,0,0))
plot(WV[waki], wakeup_ts[waki], type = "b", pch = 16, xlab = "", ylab = "",
ylim = as.numeric(YLIM), axes = FALSE, main = "Opsta- en bed-tijden", cex.lab = CL, cex = CXdots, col = COL[1])
for (ii in -24:24) {
abline(h = as.numeric(as.POSIXlt("24:00", format = "%H:%M")) + (ii * 3600), col = "grey", lty = 3)
}
axis(side = 1, at = 1:length(onsi), labels = P4N$weekday[onsi], cex.axis = CXdays)
abline(v = waki, col = "grey", lty = 3 ) # vertical lines
par(las = 1)
axis(side = 2, at = timeaxis, labels = as.character(timeaxislabels), cex.axis = 1.0)
lines(WV[waki], sleeponset_ts[onsi], type = "b", lty = 1, pch = 17, cex = CXdots, col = COL[1])
lines(WV[waki], wakeup_log_ts[waki], type = "b", lty = 1, pch = 16, cex = CXdots, col = COL[2])
lines(WV[waki], sleeponset_log_ts[onsi], type = "b", lty = 1, pch = 17, cex = CXdots, col = COL[2])
CXdots = 1.3
legend("top", legend = c("beweegmeter opstatijden", "beweegmeter slaaptijden",
"dagboek opstatijden", "dagboek bedtijden"), #, "meest inactieve 5 uur"),
col = COL[c(1, 1, 2, 2)], lty = rep(1,4), pch = c(16, 17, 16, 17), ncol = 3, cex = 0.8, bg = "white")
par(mfrow = c(2,1), mar = c(4, 4, 3, 0.5), oma = c(0,0,0,0)) # las = 3
YMAX = max(P4N$step_count_sum_0.24hr, na.rm = TRUE)
YMAX = ceiling(YMAX / 5000) * 5000 # ifelse(test = max(P4N$step_count_sum_0.24hr) > 10000, yes = 20000, no = 15000)
brpos = barplot(t(as.matrix(P4N$step_count_sum_0.24hr)),
space = rep(0, length(P4N$step_count_sum_0.24hr)),
ylab = "",
ylim = c(0, YMAX), cex.axis = 0.9,
names.arg = c(P4N$weekday), #, rep("", NBars)),
cex.names = CXdays, cex.lab = CL, cex.main = CXmain,
main = "", col = adjustcolor("lightblue", alpha.f = 0.2))
grid(nx = NA, ny = NULL,
lty = 2, # Grid line type
col = "gray", # Grid line color
lwd = 1) # Grid line width
mtext("Stappen per dag", side = 3, adj = 0, line = 1.2, cex = 0.9, font = 2);
dev.off()
cat(paste0("\nDe PDF rapporten zijn nu opgeslagen in ",pdffile,"."))
}
}
}
vincentvanhees/ActChronicFatigue documentation built on Aug. 31, 2024, 11:55 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions summarise
Documented in summarise
#' summarise
#'
#' @param outputdir ath to output directory generated by GGIR
#' @param part5_summary Content of GGIR part 5 personsummary file
#' @param model_threshold Acceleration threshold corresponding to the model, which can be derived as negative of the first coefficient divided by the second coefficient.
#' @param referentiewaarden Mean and standard deviation ENMO_fullRecordingMean in a reference population.
#' @param sleepwindowType See, GGIR documentation.
#' @param MVPAdefinition MVPA variable name from the GGIR part2 daysummary output
#' @param threshold_wrist Threshold to classify active or passive
#' @param sep separator used by csv files stored by GGIR
#' @param dec dec separator used by csv files stored by GGIR
#' @return no object is returned, only a summary is printed to the screen
#' @export
#' @importFrom grDevices dev.off pdf adjustcolor
#' @importFrom graphics abline axis barplot par grid
#' @importFrom graphics legend lines text mtext rasterImage
#' @importFrom utils packageVersion
summarise = function(outputdir,
part5_summary,
model_threshold = 75,
referentiewaarden,
sleepwindowType="SPT",
MVPAdefinition = "MVPA_E5S_B10M80%_T100_BFEN_0-24hr",
threshold_wrist = 0.5,
sep = sep, dec = dec) {
part2_summary_file = grep(dir(paste0(outputdir,"/results"), full.names = TRUE),
pattern = "part2_summary", value = T)
part2_summary = read.csv(file = part2_summary_file, sep = sep, dec = dec)
extractid = function(x) {
tmp = unlist(strsplit(x," "))
if (length(tmp) > 1) {
tmp = tmp[1]
} else {
tmp = x
}
return(tmp)
}
if (is.character(part2_summary$ID) == TRUE) {
part2_summary$ID = sapply(part2_summary$ID, FUN = extractid)
}
part2_summary = part2_summary[,c("ID", "ENMO_fullRecordingMean")] # "BFEN_fullRecordingMean" "ENMO_fullRecordingMean"
colnames(part2_summary) = c("ID","Activity_zscore")
part2_summary$Activity_zscore = round((part2_summary$Activity - referentiewaarden[1]) / referentiewaarden[2], digits = 1) # - 23.64) / 15.5 for ENMO
part5_summary = merge(part5_summary, part2_summary, by.x = "ID2", by.y = "ID")
# check for which IDs reports have already been generated:
filesinresults = dir(paste0(outputdir, "/results"))
report_names = grep(pattern = "eweeg_en_slaap_rapport_", x = filesinresults, value = TRUE)
IDdone = c()
if (length(report_names) > 0) {
extractID = function(x) {
return(unlist(strsplit(x, "_rapport_|[.]pd"))[2])
}
IDdone = unlist(lapply(X = report_names, FUN = extractID))
}
if (length(IDdone) > 0) {
most_recent_recordings = which(part5_summary$ID2 %in% IDdone == FALSE)
} else {
most_recent_recordings = 1:nrow(part5_summary)
}
if (length(most_recent_recordings) == 0) {
cat("\nDe bestanden waren al eerder genalyseerd en de rapporten staan in de resultaten. Mocht je de rapporten opnieuw willen creeeren, verwijder dan eerst de oude pdf bestanden.")
} else {
cat(paste0("Klassificaties voor meest recente ", length(most_recent_recordings), " metingen: "))
part5_summary = part5_summary[order(as.Date(part5_summary$calendar_date)),]
if ("sleep_efficiency_after_onset_pla" %in% colnames(part5_summary)) {
colnames(part5_summary)[which(colnames(part5_summary) == "sleep_efficiency_after_onset_pla")] = "sleep_efficiency_pla"
}
recent_recording = part5_summary[most_recent_recordings,
c("ID2", "calendar_date", "prop_perv_passive", "Activity_zscore",
"dur_spt_sleep_min_pla", "dur_spt_min_pla", "sleep_efficiency_pla",
"sleeponset_pla", "wakeup_pla")]
recent_recording$dur_spt_sleep_min_pla = round(recent_recording$dur_spt_sleep_min_pla / 60, digits = 2)
recent_recording$dur_spt_min_pla = round(recent_recording$dur_spt_min_pla / 60, digits = 2)
recent_recording$sleep_efficiency_pla = round(recent_recording$sleep_efficiency_pla * 100, digits = 0)
recent_recording$sleeponset_pla = round(recent_recording$sleeponset_pla - 24, digits = 2)
recent_recording$wakeup_pla = round(recent_recording$wakeup_pla - 24, digits = 2)
recent_recording$prop_perv_passive = round(recent_recording$prop_perv_passive * 100)
recent_recording$classification = "Actief"
pp = which(recent_recording$prop_perv_passive >= threshold_wrist * 100)
if (length(pp) > 0) recent_recording$classification[pp] = "Laag Actief"
colnames(recent_recording) = c("ID", "Datum start meeting", "Kans op laag actief (%)",
"Gemiddelde beweging t.o.v. referentie groep (z-score)",
"Uren slaap per nacht", "Gem. duur slaapperioden (uren)",
"Slaap efficiency binnen slaapperioden (%)", "Gemiddelde slaap start t.o.v. middernacht",
"Gemiddelde slaap einde t.o.v. middernacht", "Klassificatie") #,"feedback")
varnames = c("Datum start meeting", "Klassificatie", "Kans op laag actief (%)",
"Gemiddelde beweging t.o.v. referentie groep (z-score)")
recent_recording = recent_recording[,c(1:4,10,5:9)]
cat("\n")
part5_daysummary_file = grep(dir(paste0(outputdir,"/results"), full.names = TRUE),
pattern = "part5_daysummary", value = T)
part5_daysummary = read.csv(file = part5_daysummary_file, sep = sep, dec = dec)
part2_daysummary_file = grep(dir(paste0(outputdir,"/results"), full.names = TRUE),
pattern = "part2_daysummary", value = T)
part2_daysummary = read.csv(file = part2_daysummary_file, sep = sep, dec = dec)
# Load sleep/wake times from both GGIR and diary
part4_nightsummaryfull_file = grep(dir(paste0(outputdir,"/results/QC"), full.names = TRUE),
pattern = "part4_nightsummary_sleep_full.csv", value = T)
part4_nightsummary = read.csv(part4_nightsummaryfull_file, sep = sep, dec = dec)
days = c("Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday", "Sunday")
dagen = c("maandag", "dinsdag", "woensdag", "donderdag", "vrijdag", "zaterdag", "zondag")
for (ki in 1:length(days)) {
part5_daysummary$weekday = gsub(x = part5_daysummary$weekday,
pattern = days[ki], replacement = dagen[ki])
part2_daysummary$weekday = gsub(x = part2_daysummary$weekday,
pattern = days[ki], replacement = dagen[ki])
part4_nightsummary$weekday = gsub(x = part4_nightsummary$weekday,
pattern = days[ki], replacement = dagen[ki])
}
for (ID in unique(recent_recording$ID)) {
# Get rows for this ID
P5D = part5_daysummary[which(part5_daysummary$ID == ID),]
P2D = part2_daysummary[which(part2_daysummary$ID == ID),]
P4N = part4_nightsummary[which(part4_nightsummary$ID == ID),]
# Put date in data format
P5D$calendar_date = as.Date(P5D$calendar_date, "%Y-%m-%d")
P2D$calendar_date = as.Date(P2D$calendar_date, "%Y-%m-%d")
P4N$calendar_date = as.Date(P4N$calendar_date, "%d/%m/%Y")
# Merge time spent in two acceleration range to sleep object, to ease plotting later on
# Combine P4 and P2 output to ease plotting
P4N = base::merge(P4N, P2D[,c("calendar_date", MVPAdefinition, "step_count_sum_0.24hr")],
by = c("calendar_date"), all.x = TRUE)
colnames(P4N)[which(colnames(P4N) %in% MVPAdefinition == TRUE)] = c("MVPA") #"SB", "LIPA",
# Remove problematic nights from the sleep results P4
missingSleeplog = which(P4N$guider != "sleeplog")
if (length(missingSleeplog) > 0) {
is.na(P4N[missingSleeplog, grep(pattern = "guider", x = colnames(P4N), invert = FALSE)]) = TRUE
}
sleepProblem = which(P4N$cleaningcode > 1 & P4N$guider == "sleeplog")
if (length(sleepProblem) > 0) {
is.na(P4N[sleepProblem, grep(pattern = "calendar|ID|night|guider|weekday|page|filename|cleaningcode",
x = colnames(P4N), invert = TRUE)]) = TRUE
}
insufficientValidData4 = which(P4N$fraction_night_invalid > 0.33)
if (length(insufficientValidData4) > 0) {
is.na(P4N[insufficientValidData4, grep(pattern = "calendar|ID|night|guider|weekday|page|filename|cleaningcode",
x = colnames(P4N), invert = TRUE)]) = TRUE
}
#---------------------------------------------
# add empty rows for missing days
P5D = P5D[order(P5D$calendar_date),]
dates = P5D$calendar_date
dates_theoretical = seq(min(dates), max(dates), by = 1)
missing_dates = dates_theoretical[which(dates_theoretical %in% dates == FALSE)]
if (length(missing_dates) > 0) {
newvalues = (nrow(P5D) + 1):(nrow(P5D) + length(missing_dates))
P5D[newvalues,] = NA
P5D$calendar_date[newvalues] = format(as.Date(missing_dates), "%Y-%m-%d")
P5D$weekday[newvalues] = weekdays(missing_dates, abbreviate = FALSE)
if (length(which(P5D$weekday %in% days == TRUE)) > 0) {
for (ki in 1:length(days)) {
P5D$weekday = gsub(x = P5D$weekday,
pattern = days[ki], replacement = dagen[ki])
}
}
}
P5D = P5D[order(P5D$calendar_date),] # order again now missing day has been added
P5D$calendar_date = format(as.Date(P5D$calendar_date), "%Y-%m-%d")
P5D$calendar_date = format(as.Date(P5D$calendar_date), "%d")
# Now shorten weekday
dagen_kort = c("ma", "di", "wo", "do", "vr", "za", "zo")
# dagen = c("maandag", "dinsdag", "woensdag", "donderdag", "vrijdag", "zaterdag", "zondag")
dagen = c("maandag", "dinsdag", "woensdag", "donderdag", "vrijdag", "zaterdag", "zondag")
for (ki in 1:length(days)) {
P5D$weekday = gsub(x = P5D$weekday,
pattern = dagen[ki], replacement = dagen_kort[ki])
P4N$weekday = gsub(x = P4N$weekday,
pattern = dagen[ki], replacement = dagen_kort[ki])
}
# Extract person summary for part 5
P5P = part5_summary[which(part5_summary$ID == ID),]
pdffile = paste0(outputdir, "/results/Beweeg_en_slaap_rapport_",ID,".pdf")
pdf(file = pdffile, paper = "a4r", width = 11.69, height = 8.27)
par(mfrow = c(4,1), mar = c(4, 4, 3, 0.5), oma = c(0,0,0,0)) # las = 3
img <- png::readPNG(system.file("extdata/Amsterdam_UMC_Logo2.png", package = "ActChronicFatigue"))
# Key facts
keystats = t(recent_recording[which(recent_recording$ID == ID), varnames])
colnames(keystats)[1] = paste0("ID: ", ID)
CXmain = 1.0
CXdays = 1.0
relprop = as.numeric(keystats[3])
if (keystats[2] != "Laag Actief") relprop = 100 - relprop
IDtitle = paste0("ID: ", ID, " | Start: ", keystats[1])
title = paste0("Klassificatie: ", keystats[2], " (zekerheid ", relprop,"%) | ",
"Beweging (z-score) = ", keystats[4])
CL = 1.0
CXdots = 1.5
#========================================
# Plot activity time series (based on GGIR part 5 output)
NBars = 3
TS = P5D$ACC_day_mg
TS = c(TS, rep(NA, NBars))
plotlogo = function() {
grid::grid.raster(img, x = .91, y = 0.98, width = .18, interpolate = FALSE) # print homer in ll conrner
grid::grid.text(x = 0.91, y = 0.93,
gp = grid::gpar(fontsize = 7, col = "black", font = 2),
label = "Afdeling Medische Psychologie, NKCV") #, xpd = NA, col ="red", cex= 1, pos = 4)
}
plot(TS, type = "b", pch = 20, ylim = range(c(model_threshold * 2, P5D$ACC_day_mg), na.rm = T),
main = "", xlab = "", ylab = "Beweging per dag", bty = "l", axes = FALSE, cex.main = 1.5,
cex.lab = CL, cex = CXdots)
mtext(title, side = 3, adj = 0, line = 1.2, cex = 0.9, font = 2);
axis(side = 1, at = 1:(length(TS) - NBars), labels = c(paste0(P5D$weekday, " ", P5D$calendar_date)), #, rep("", 3)
cex.axis = CXdays)
lines(x = c(1,length(TS)), y = rep(model_threshold, 2), col = "black", lty = 2, lwd = 1.3)
text(x = length(TS) - 3, y = model_threshold + 10, labels = "Grenswaarde laag actief", pos = 4, cex = 1)
plotlogo()
#========================================
# Plot sleep duration (based on GGIR part 4 output)
if (nrow(P4N) > nrow(P5D) & nrow(P4N) > 1) {
P4N = P4N[2:nrow(P4N),]
}
if (sleepwindowType == "TimeInBed") {
A = P4N[,c("SleepDurationInSpt" , "SptDuration", "sleeplatency", "guider_inbedDuration", "sleepefficiency")]
title = "Slapen en wakker liggen (met slaap efficientie (%) aangeduid in wit)" #Slaap (donker), Wakker na start slaap (licht), Slaap latency (lichter), en
maxvalue = max(c(A[,2] + A[,3]), na.rm = T)
} else {
A = P4N[,c("SleepDurationInSpt" , "SptDuration")] #/ 60
title = "Slaap" #Slaap (donker), Wakker na start slaap (licht)
maxvalue = max(A[,2], na.rm = T)
}
A[,2] = A[,2] - A[,1]
if (sleepwindowType == "TimeInBed") {
B = A[,c(1,2, 3)]
} else {
B = A[,c(1,2)]
}
B[(nrow(B) + 1):(nrow(B) + NBars),] <- NA
brpos = barplot(t(as.matrix(B)), space = rep(0, nrow(B)), ylab = "Tijd in uren", beside = FALSE, #space = c(0, 0.2),
ylim = c(0, maxvalue + 6), cex.axis = 0.9,
names.arg = c(P4N$weekday, rep("",nrow(B) - nrow(P4N))),
cex.names = CXdays, cex.lab = CL, cex.main = CXmain,
legend.text = c("Slaap", "Wakker na in slaap vallen", "Wakker voor in slaap vallen"),
main = "",
args.legend = list(ncol = 1, x = nrow(B), cex = 1, bty = "n"))
mtext(title, side = 3, adj = 0, line = 1.2, cex = 0.9, font = 2);
if (sleepwindowType == "TimeInBed") {
text(brpos, rep(1, length(A$guider_inbedDuration)),
labels = round(100 * (A$sleepefficiency), digits = 0),
cex = 0.8,
col = "white")
}
#========================================
# Niet gedragen tijd
A = P5D[,c("nonwear_perc_day", "nonwear_perc_spt")]
A[(nrow(A) + 1):(nrow(A) + NBars),] <- NA
barplot(t(as.matrix(A)), ylab = "Percentage (%)", beside = TRUE, space = c(0, 0.2),
names.arg = c(P5D$weekday, rep("",NBars)), cex.names = CXdays, cex.lab = CL, cex.main = CXmain,
main = "", legend.text = c("overdag", "nacht"),
args.legend = list(x = (nrow(A)*2) + 1, ncol = 1, cex = 1, bty = "n"), ylim = c(0,100))
mtext("Beweegmeter niet gedragen?", side = 3, adj = 0, line = 1.2, cex = 0.9, font = 2);
#=======================================
# MVPA
maxvalue = max(c(P4N$MVPA, 30), na.rm = T) + 10
boutdur = unlist(strsplit(unlist(strsplit(MVPAdefinition, "B"))[2], "M"))[1]
title = paste0("Matig to Zwaar intensief gedrag in blokken van minimaal ", boutdur, " minuten")
PASB = P4N[,c("MVPA")]
PASB = c(PASB, matrix(NA,3, 1))
brpos = barplot(t(as.matrix(PASB)),
space = rep(0, length(PASB)),
ylab = "Tijd in minuten",
ylim = c(0, maxvalue), cex.axis = 0.9,
names.arg = c(P4N$weekday, rep("",NBars)),
cex.names = CXdays, cex.lab = CL, cex.main = CXmain,
main = "", col = "grey")
lines(x = c(0, length(P4N$weekday) + 3), y = c(30, 30), type = "l", lty = 2, lwd = 1.3)
mtext(title, side = 3, adj = 0, line = 1.2, cex = 0.9, font = 2);
text(x = length(P4N$weekday), y = 35, labels = "Beweegrichtlijn (30 min. per dag)", cex = 1, font = 1, pos = 4);
mtext(text = paste0("Software versie: ", packageVersion("ActChronicFatigue")),
side = 1, col = "black", cex = 0.7, line = 2.5, font = 2, adj = 0)
mtext(text = IDtitle,
side = 1, col = "black", cex = 0.7, line = 2.5, font = 2, adj = 1)
#========================================
# Slaap waak tijden
sleeponset_ts = P4N$sleeponset_ts # c("22:00", "23:20", "01:10")
wakeup_ts = P4N$wakeup_ts # c("6:55", "8:30", "10:40")
if (sleepwindowType == "SPT") {
sleeponset_log_ts = P4N$guider_onset_ts
wakeup_log_ts = P4N$guider_wakeup_ts
} else if (sleepwindowType == "TimeInBed") {
sleeponset_log_ts = P4N$guider_inbedStart_ts
wakeup_log_ts = P4N$guider_inbedEnd_ts
}
WV = which((is.na(wakeup_ts) == FALSE & is.na(sleeponset_ts) == FALSE) |
(is.na(wakeup_log_ts) == FALSE & is.na(sleeponset_log_ts) == FALSE))
sleeponset_ts = as.POSIXlt(sleeponset_ts, format = "%H:%M")
wakeup_ts = as.POSIXlt(wakeup_ts, format = "%H:%M")
sleeponset_log_ts = as.POSIXlt(sleeponset_log_ts, format = "%H:%M")
wakeup_log_ts = as.POSIXlt(wakeup_log_ts, format = "%H:%M")
for (j in 1:length(WV)) {
if (!is.na(sleeponset_ts[j])) {
if (sleeponset_ts[j] < as.POSIXlt("12:00", format = "%H:%M")) {
sleeponset_ts[j] = sleeponset_ts[j] + (24 * 3600)
}
}
if (!is.na(sleeponset_log_ts[j])) {
if (sleeponset_log_ts[j] < as.POSIXlt("12:00", format = "%H:%M")) {
sleeponset_log_ts[j] = sleeponset_log_ts[j] + (24 * 3600)
}
}
if (!is.na(wakeup_ts[j]) & !is.na(sleeponset_ts[j])) {
if (wakeup_ts[j] < sleeponset_ts[j] |
(wakeup_ts[j] > as.POSIXlt("00:00", format = "%H:%M")) &
(wakeup_ts[j] < as.POSIXlt("12:00", format = "%H:%M"))) {
wakeup_ts[j] = wakeup_ts[j] + (24 * 3600)
}
}
if (!is.na(wakeup_log_ts[j]) & !is.na(sleeponset_log_ts[j])) {
if (wakeup_log_ts[j] < sleeponset_log_ts[j] |
(wakeup_log_ts[j] > as.POSIXlt("00:00", format = "%H:%M")) &
(wakeup_log_ts[j] < as.POSIXlt("12:00", format = "%H:%M"))) {
wakeup_log_ts[j] = wakeup_log_ts[j] + (24 * 3600)
}
}
if (!is.na(wakeup_ts[j])) wakeup_ts[j] = wakeup_ts[j] - (24 * 3600)
if (!is.na(wakeup_log_ts[j])) wakeup_log_ts[j] = wakeup_log_ts[j] - (24 * 3600)
}
YLIM = c(as.POSIXlt("12:00", format = "%H:%M") - 12 * 3600, (as.POSIXlt("12:00", format = "%H:%M") + 24 * 3600))
timeaxis = c(as.numeric(YLIM[1]) + (c(0:36) * (3600)))
timeaxislabels = format(as.POSIXlt(timeaxis, origin = "1970-1-1"), "%H:%M")
waki = 1:(max(WV) - 1)
onsi = 2:max(WV)
# Check whether onset occurs before waking up, if yes then add 24 hour to onset
correctionNeeded = which(sleeponset_ts[onsi] < wakeup_ts[waki])
if (length(correctionNeeded) > 0) {
sleeponset_ts[onsi][correctionNeeded] = sleeponset_ts[onsi][correctionNeeded] + (24 * 3600)
}
COL = c("blue", "red", "black")
CL = 1
par(mfrow = c(1,1), oma = c(0,0,0,0))
plot(WV[waki], wakeup_ts[waki], type = "b", pch = 16, xlab = "", ylab = "",
ylim = as.numeric(YLIM), axes = FALSE, main = "Opsta- en bed-tijden", cex.lab = CL, cex = CXdots, col = COL[1])
for (ii in -24:24) {
abline(h = as.numeric(as.POSIXlt("24:00", format = "%H:%M")) + (ii * 3600), col = "grey", lty = 3)
}
axis(side = 1, at = 1:length(onsi), labels = P4N$weekday[onsi], cex.axis = CXdays)
abline(v = waki, col = "grey", lty = 3 ) # vertical lines
par(las = 1)
axis(side = 2, at = timeaxis, labels = as.character(timeaxislabels), cex.axis = 1.0)
lines(WV[waki], sleeponset_ts[onsi], type = "b", lty = 1, pch = 17, cex = CXdots, col = COL[1])
lines(WV[waki], wakeup_log_ts[waki], type = "b", lty = 1, pch = 16, cex = CXdots, col = COL[2])
lines(WV[waki], sleeponset_log_ts[onsi], type = "b", lty = 1, pch = 17, cex = CXdots, col = COL[2])
CXdots = 1.3
legend("top", legend = c("beweegmeter opstatijden", "beweegmeter slaaptijden",
"dagboek opstatijden", "dagboek bedtijden"), #, "meest inactieve 5 uur"),
col = COL[c(1, 1, 2, 2)], lty = rep(1,4), pch = c(16, 17, 16, 17), ncol = 3, cex = 0.8, bg = "white")
par(mfrow = c(2,1), mar = c(4, 4, 3, 0.5), oma = c(0,0,0,0)) # las = 3
YMAX = max(P4N$step_count_sum_0.24hr, na.rm = TRUE)
YMAX = ceiling(YMAX / 5000) * 5000 # ifelse(test = max(P4N$step_count_sum_0.24hr) > 10000, yes = 20000, no = 15000)
brpos = barplot(t(as.matrix(P4N$step_count_sum_0.24hr)),
space = rep(0, length(P4N$step_count_sum_0.24hr)),
ylab = "",
ylim = c(0, YMAX), cex.axis = 0.9,
names.arg = c(P4N$weekday), #, rep("", NBars)),
cex.names = CXdays, cex.lab = CL, cex.main = CXmain,
main = "", col = adjustcolor("lightblue", alpha.f = 0.2))
grid(nx = NA, ny = NULL,
lty = 2, # Grid line type
col = "gray", # Grid line color
lwd = 1) # Grid line width
mtext("Stappen per dag", side = 3, adj = 0, line = 1.2, cex = 0.9, font = 2);
dev.off()
cat(paste0("\nDe PDF rapporten zijn nu opgeslagen in ",pdffile,"."))
}
}
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.