##############################################################################
### ACTman{plot_actogram} ###
### Script authors: Yoram Kunkels, Stefan Knapen, & Ando Emerencia ###
### Most recent Update: 21-12-2017 ###
### Goal: plot activity counts over time, with 24 hour lines (or 48 hours) ###
### ###
### To add: Grey part before start, Title, ###
### Export to PDF function, Days more clearly written ###
###========================================================================###
### Revision History: ###
### 16-04-2018: Added x-axis hours & 24h Plotting ###
###~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@###
#' plot_actogram
#'
#' Function to plot 48 hour Actograms.
#'
#' @param workdir the working directory as supplied to ACTman.
#' @param ACTdata.1.sub The managed data set
#' @param i The index of the current file in ACTdata.files
#' @param plotactogram Value indicating if and what kind of actogram has to be plotted. Can be either '48h', '24h', or FALSE
#'
#'
#' @importFrom grDevices dev.off
#' @importFrom grDevices pdf
#' @importFrom graphics axis
#' @importFrom graphics barplot
#' @importFrom graphics par
#' @importFrom stats na.omit
plot_actogram <- function(workdir, ACTdata.1.sub, i, plotactogram, rollingwindow.results, i_want_EWS) {
### Part 1: Basic Operations ----------------------------------------------------------------------------
# Set and save working directory
workdir.save <- getwd()
setwd(workdir)
act_data <- ACTdata.1.sub # Copy data for editing required for plotting
act_data <- within(act_data, Date <- as.character(act_data$Date)) # Date as character for plotting
ndays.plot <- round(abs(as.numeric(round(as.POSIXct(act_data$Date[1]) - as.POSIXct(act_data$Date[nrow(act_data)]), 2))))
### Part 2: 1st day Selection & Midnight Detection 2nd day ---------------------------------------------------
## Midnight Detection 2nd day
day2start <- ((which(substr(act_data$Date, start = 12, stop = 19) == "00:00:00")[1]) - 1) # Start of day2 at midnight.
## 1st Day selection
day1 <- act_data[1:day2start, ]
day1.rest <- 1440 - day2start
## Filling in period before 1st day selection to get full 24h in 1st day
day1.rest.mat <- matrix(data = 0, nrow = day1.rest, ncol = 2) # Create empty matrix
colnames(day1.rest.mat) <- colnames(day1) # Equalise colnames
# day1 <- within(day1, Date <- as.POSIXct(as.character(Date), format = "%Y-%m-%d %H:%M:%S")) # Do POSIXct magic (http://grokbase.com/t/r/r-help/112jk8eyqp/r-problem-with-rbind-when-data-frame-contains-an-date-time-variable-posixt-posixlt)
day1 <- rbind(day1.rest.mat, day1) # Combine day1 data with empty pre-day1 matrix
### Part 3: Loop for assigning subsequent days (14 days max for plot) & setting ylim's ----------------------
## Assign other days
for (i.plot in 2:(ndays.plot)) {
if (i.plot == 2) {
assign(paste("day", i.plot, sep = ""), act_data[((day2start + (1440 ^ (i.plot - 2))):(day2start + (1440 * (i.plot - 1)))), ])
} else {
assign(paste("day", i.plot, sep = ""), act_data[((day2start + ((1440 * (i.plot - 2)) + 1)):(day2start + (1440 * (i.plot - 1)))), ])
}
}
## Define the range so all plots will be equal height (ylim)
ylimit_TEMP <- ls()[setdiff(grep("day", ls()), c(grep("rest", ls()), grep("start", ls()), grep("plot", ls())))]
ylimit <- range(na.omit(eval(parse(text = ylimit_TEMP))[, "Activity"]))
# ylimit <- range(na.omit(c(day1$Activity, day2$Activity, day3$Activity, day4$Activity, day5$Activity,
# day6$Activity, day7$Activity, day8$Activity, day9$Activity, day10$Activity,
# day11$Activity, day12$Activity, day13$Activity, day14$Activity))) #! Obsolete
### Part 4: Combining Days & 48 hour Doubleplot -----------------------------------------------------------
## Combine days
for (j.plot in 1:(ndays.plot - 1)) {
assign(paste("day.", j.plot, ".", (j.plot + 1), sep = ""),
rbind(eval(parse(text = paste("day", j.plot, sep = ""))), eval(parse(text = paste("day", (j.plot + 1), sep = "")))))
}
### Work in Progress ------------------------------------------------------------------------------
#! TEST for moving window plot
## combine all days in 1 plot
## Runner
# ACTman::ACTman(workdir = "C:/Bibliotheek/Studie/PhD/Publishing/ACTman/R-part/TEST-RESULTS",
# myACTdevice = "MW8", lengthcheck = F, nparACT_compare = F, iwantsleepanalysis = F,
# circadian_analysis = T, plotactogram = "24h", movingwindow = TRUE, movingwindow.size = 7)
if (i_want_EWS == TRUE && is.na(rollingwindow.results)) {
message("Cannot create EWS plot without rolling window results")
message("Please make sure that both 'i_want-EWS' AND 'movingwindow' arguments are TRUE")
message("Quitting...")
stop()
}
if (i_want_EWS == TRUE) { # ## Initialise empty matrix for timestamps and activity counts
LOLkat <- matrix(NA, nrow = (1440 * ndays.plot), ncol = 2)
## Assign timestamps and activity counts to matrix
for (k.plot in 1:ndays.plot) {
if (k.plot == 1) {
LOLkat[k.plot:1440, 1] <- eval(parse(text = paste("day", k.plot, "[ , 1]", sep = "")))
LOLkat[k.plot:1440, 2] <- eval(parse(text = paste("day", k.plot, "[ , 2]", sep = "")))
} else {
LOLkat[((((k.plot - 1) * 1440) + 1):(k.plot * 1440)), 1] <- eval(parse(text = paste("day", k.plot, "[ , 1]", sep = "")))
LOLkat[((((k.plot - 1) * 1440) + 1):(k.plot * 1440)), 2] <- eval(parse(text = paste("day", k.plot, "[ , 2]", sep = "")))
}
}
# plotme <- "Time_to_Recovery"
plotme <- colnames(rollingwindow.results)[10:21]
for (EWS_count in 1:length(plotme)) {
## Initialise .PDF plot in A4 size (11.7 x 8.3 inches)
png(paste("Actigraphy EWS Plot - ", plotme[EWS_count], ".png"), width = 842, height = 595, units = "px")
## Create barplot
bp2 <- barplot(as.numeric(LOLkat[, 2]), plot = FALSE)
## Obtain barplot range
# bp2_ylim <- range(as.numeric(LOLkat[, 2]))
bp2_ylim <- range(na.omit(as.numeric(LOLkat[, 2])))
roundup_power_10 <- function(x) 10 ^ ceiling(log10(x))
bp2_ylim_upper <- roundup_power_10(max(bp2_ylim))
## Plot barplot
barplot(as.numeric(LOLkat[, 2]), ylim = c(0, bp2_ylim_upper))
## Create labels
x_labels2 <- substr(LOLkat[, 1], 1, 10)
l.plot_n <- length(unique(substr(LOLkat[, 1], 1, 10)[!substr(LOLkat[, 1], 1, 10) == "0"]))
x_labels_pos2_start <- matrix(NA, nrow = l.plot_n, ncol = 1)
l.plot <- 0
## Obtain label positions
for (l.plot in 1:l.plot_n) {
## Assure l.plot is not NA, otherwise assign last non-NA value
if (is.na(unique(substr(LOLkat[, 1], 1, 10)[!substr(LOLkat[, 1], 1, 10) == "0"])[l.plot])) {
l.plot <- max(which(is.na(unique(substr(LOLkat[, 1], 1, 10)[!substr(LOLkat[, 1], 1, 10) == "0"])) == FALSE))
}
x_labels_pos2_start[l.plot, ] <- which.max(x_labels2 == unique(substr(LOLkat[, 1], 1, 10)[!substr(LOLkat[, 1], 1, 10) == "0"])[l.plot])
}
## Assign start and end dates
LOLkat_startdate <- x_labels2[x_labels_pos2_start][1] # Activity data start date
LOLkat_enddate <- x_labels2[x_labels_pos2_start][length(x_labels2[x_labels_pos2_start])] # Activity data end date
## Create axis and vertical day lines
axis(side = 1, at = bp2[1 + x_labels_pos2_start], labels = x_labels2[x_labels_pos2_start], las = 2, cex.axis = 0.8)
abline(v = bp2[1 + x_labels_pos2_start], col = "blue")
#! debug
rollingwindow.results <<- rollingwindow.results
LOLkat <<- LOLkat
## Assign moving window results (remove first obs to account for non-24h day)
rollingwindow.results <- rollingwindow.results[2:nrow(rollingwindow.results), ]
## Asssign occurences when days from plotted actigraphy data co-occur with days from the moving window results
matched_dates <- paste(x_labels2[x_labels_pos2_start], "00:00:00") %in% rollingwindow.results$endtime
## Assign x-coordinates for moving window results
plot_points_x <- bp2[1 + x_labels_pos2_start][matched_dates]
# # plotme <- "Time_to_Recovery"
# plotme <- colnames(rollingwindow.results)[10:21]
#
# for(EWS_count in 1:length(plotme)) {
## Plot results for moving window on existing barplot at designated x-coordinates
#! N.b. Scaling should be made dynamic!
# plotme <- "Time_to_Recovery"
scaling_var <- (bp2_ylim_upper / max(rollingwindow.results[(1:(nrow(rollingwindow.results) - 1)), plotme[EWS_count]]))
## Plotting "plotme" EWS over actogram
## Create right axis for EWS measure
## Exception for when length(at) != length(labels)
#! Added na.omit() functions for when "to is not finite" error.
if (length(seq(0, round(max(na.omit(rollingwindow.results[(1:(nrow(rollingwindow.results) - 1)), plotme[EWS_count]]))),
by = (round(max(na.omit(rollingwindow.results[(1:(nrow(rollingwindow.results) - 1)), plotme[EWS_count]]))) / 2)))
== length(c(0, (bp2_ylim_upper / 2), bp2_ylim_upper))) {
axis(side = 4,
labels = seq(0, round(max(na.omit(rollingwindow.results[(1:(nrow(rollingwindow.results) - 1)), plotme[EWS_count]]))),
by = (round(max(na.omit(rollingwindow.results[(1:(nrow(rollingwindow.results) - 1)), plotme[EWS_count]]))) / 2)),
at = c(0, (bp2_ylim_upper / 2), bp2_ylim_upper))
} else {
axis(side = 4,
labels = c(0, (bp2_ylim_upper / 2), bp2_ylim_upper),
at = c(0, (bp2_ylim_upper / 2), bp2_ylim_upper))
}
## Plot EWS over actogram (exceptions for when length of x- and y-values differ)
if (length(plot_points_x) < length((rollingwindow.results[(1:(nrow(rollingwindow.results))), plotme[EWS_count]]
* scaling_var))) {
# print("DEBUG--------------------------------------------------------------------------")
# print(length(plot_points_x))
# print(length((rollingwindow.results[(1:(nrow(rollingwindow.results))), plotme] * scaling_var)))
points(x = plot_points_x, type = "l", col = "red",
y = (rollingwindow.results[(1:(nrow(rollingwindow.results))), plotme[EWS_count]]
* scaling_var)[1:length(plot_points_x)])
} else if (length(plot_points_x) > length((rollingwindow.results[(1:(nrow(rollingwindow.results))), plotme] * scaling_var))) {
points(x = plot_points_x[1:length(rollingwindow.results[(1:(nrow(rollingwindow.results))), plotme[EWS_count]])],
type = "l", col = "red", y = (rollingwindow.results[(1:(nrow(rollingwindow.results))), plotme[EWS_count]]
* scaling_var))
} else if (length(plot_points_x) == length((rollingwindow.results[(1:(nrow(rollingwindow.results))),
plotme[EWS_count]] * scaling_var))) {
points(x = plot_points_x, type = "l", col = "red",
y = (rollingwindow.results[(1:(nrow(rollingwindow.results))), plotme[EWS_count]] * scaling_var))
}
## Create Title
mtext(paste("Total Activity data with", plotme[EWS_count], "moving window"), line = 1, cex = 1.0)
mtext(paste("start date:", LOLkat_startdate, "end date:", LOLkat_enddate), line = 0, cex = 0.8)
dev.off()
}
# dev.off()
}
### /Work in Progress ------------------------------------------------------------------------------
if (plotactogram == "48h") {
## Initialise .PDF plot in A4 size (11.7 x 8.3 inches)
pdf(paste("Actigraphy Data - 48h Plot", i, ".pdf"), width = 11.7, height = 8.3)
## Plot initialisation & parameters
par(mfrow = c(14, 1)) # Set plot parameters
par(mar = c(0.5, 4, 0.5, 4)) # Set margins
# par(mai = c(1, 1, 1, 1))
bp <- barplot(day.1.2$Activity, ylim = ylimit, ylab = "Day 1", plot = FALSE)
barplot(day.1.2$Activity, ylim = ylimit, ylab = "Day 1")
# abline(v = day2start + 75, lty = 2) # Set start moment
x_labels <- substr(day.1.2$Date, nchar(day.1.2$Date) - 8 + 1, nchar(day.1.2$Date))
x_labels_pos <- grep("00:00", x_labels)
x_labels <- x_labels[x_labels_pos]
x_labels <- substr(x_labels, start = 1, stop = 5)
axis(side = 1, at = bp[1 + x_labels_pos[ c(TRUE, FALSE)]], labels = x_labels[ c(TRUE, FALSE)])
axis(side = 1, at = bp[1 + x_labels_pos[ c(FALSE, TRUE)]], labels = FALSE, col.ticks = "red")
## Filling in plot with loop for other days
for (k.plot in 2:(ndays.plot - 1)) {
barplot(eval(parse(text = paste("day.", k.plot, ".", (k.plot + 1), "$Activity", sep = ""))),
ylim = ylimit, ylab = paste("Day", k.plot))
x_labels <- substr(day.2.3$Date, nchar(day.2.3$Date) - 8 + 1, nchar(day.2.3$Date))
x_labels_pos <- grep("00:00", x_labels)
x_labels <- x_labels[x_labels_pos]
x_labels <- substr(x_labels, start = 1, stop = 5)
axis(side = 1, at = bp[1 + x_labels_pos[ c(TRUE, FALSE)]], labels = x_labels[ c(TRUE, FALSE)])
axis(side = 1, at = bp[1 + x_labels_pos[ c(FALSE, TRUE)]], labels = FALSE, col.ticks = "red")
}
dev.off()
}
### Part 5: 24 hour plot ------------------------------------------------------------------------
## 24 hour plot
if (plotactogram == "24h") {
## Initialise .PDF plot in A4 size (11.7 x 8.3 inches)
pdf(paste("Actigraphy Data - 24h Plot", i, ".pdf"), width = 11.7, height = 8.3)
par(mfrow = c(14, 1)) # Set parameters for plots
par(mar = c(0.5, 4, 0.5, 4)) # Set margins
bp <- barplot(day1$Activity, ylim = ylimit,
ylab = substr(day1$Date[(which(day1$Date == "0")[length(which(day1$Date == "0"))] + 1)], 6, 10),
plot = F)
barplot(day1$Activity, ylim = ylimit,
ylab = substr(day1$Date[(which(day1$Date == "0")[length(which(day1$Date == "0"))] + 1)], 6, 10),
plot = T)
# abline(v = day2start+75, lty = 2) # Add start line (+75 = marge)
x_labels <- substr(day2$Date, nchar(day2$Date) - 8 + 1, nchar(day2$Date))
x_labels_pos <- grep("00:00", x_labels)
x_labels <- x_labels[x_labels_pos]
axis(side = 1, at = bp[1 + x_labels_pos], labels = x_labels)
## Filling in plot with loop for other days
for (k.plot in 2:(ndays.plot - 1)) {
barplot(eval(parse(text = paste("day", k.plot, "$Activity", sep = ""))),
ylim = ylimit,
ylab = substr((eval(parse(text = paste("day", k.plot, "$Date", sep = ""))))[1], 6, 10)
)
# axis(side = 1, at = (x_labels_pos), labels = x_labels)
axis(side = 1, at = bp[1 + x_labels_pos], labels = x_labels)
}
dev.off()
}
### Part 6: Finishing Operations -------------------------------------------------------------------
setwd(workdir.save)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.