R/heatmaps.R

In wadpac/sleepsight-analytics-pipeline: Wearable Sensor and Survey Processing Pipeline

#' heatmaps
#'
#' @param Dshort data.frame generated by agg.sleepsight
#' @param Dlong data.frame generated by agg.sleepsight
#' @param heatmapsfile png.file to which the plot should be written to
#' @param heatmapsfile_steps png.file to which the plot should be written to
#' @param simplify.behavioralclasses TRUE or FALSE
#' @param Dsurvey data.frame generated by agg.sleepsight
#' @param startDate Character,  optional argument to specify the startDate of the recording, e.g. 2019-01-01
#' @param endDate Character, optional argument to specify the startDate of the recording, e.g. 2019-01-01
#' @param desiredtz timezone (character) in Europe/London format
#' @return no output, just a file is stored
#' @export
#' @importFrom graphics plot
#' @importFrom utils timestamp
#' @importFrom cowplot get_legend plot_grid
#' @import ggplot2
# #' @importFrom ggplot2  aes element_text geom_raster geom_tile ggplot labs scale_colour_manual scale_fill_continuous scale_fill_manual scale_x_discrete theme theme_bw xlab ylab
heatmaps = function(Dshort, Dlong, heatmapsfile, heatmapsfile_steps,
                    simplify.behavioralclasses, Dsurvey = c(),
                    startDate = c(), endDate=c(), desiredtz) {
  cat("\n* plot heatmaps")
  # Declare functions:
  addDoublePlot = function(x) {
    doubleplot = x
    doubleplot$hour_in_day = doubleplot$hour_in_day + 24
    doubleplot$date = doubleplot$date - 1
    x = rbind(x, doubleplot)
    x = x[order(x$hour_in_day,x$date),]
    # x = x[order(x$date, x$hour_in_day),]
    x$date = as.character(x$date)
    return(x)
  }
  #-------------------------------------------
  if (length(startDate) > 0 & length(endDate) > 0) {
    startDateNum = as.numeric(as.Date(startDate)) * 3600*24
    endDateNum = as.numeric(as.Date(endDate)) * 3600*24
    timeNumLong = as.numeric(as.POSIXlt(Dlong$time,"1970-01-01",tz=desiredtz))
    timeNumShort = as.numeric(as.POSIXlt(Dshort$time,"1970-01-01",tz=desiredtz))
    validdatesLong = which(timeNumLong > startDateNum & timeNumLong < endDateNum)
    validdatesShort = which(timeNumShort > startDateNum & timeNumShort < endDateNum)
    if (length(validdatesLong) != 0) Dlong = Dlong[validdatesLong,]
    if (length(validdatesShort) != 0) Dshort = Dshort[validdatesShort,]
  }
  # color coding:
  coldic <- c("active" = "#D55E00",
              "inactive" = "#F0E442",
              "sleep" = "#56B4E9",
              "sustained inactive" = "#009E73",
              "inconsluvie" = "#CC79A7",
              "no data" = "#999999")
  coldic2 = c("bed" = "black","rise" = "purple")
  
  if (simplify.behavioralclasses == TRUE) {
    # simplify status classes to only active and inactive:
    Dshort$status[which(Dshort$status %in% c("sustained inactive","sleep") == TRUE)] = "inactive"
  }
  dates_on_x_axis = unique(Dshort$date)
  dates_on_x_axis = dates_on_x_axis[seq(1,length(dates_on_x_axis),round(length(dates_on_x_axis) / 25))]
  Dshort$date = as.Date(Dshort$date)

  # Add clock times to Dsurvey for risetime
  Dsurvey = Dsurvey[!duplicated(Dsurvey),]
  Dsurvey$date = as.character(as.Date(Dsurvey$risetime))
  Dsurvey$hour_rise = Dsurvey$risetime$hour
  Dsurvey$min_rise = Dsurvey$risetime$min
  Dsurvey$hour_in_day_rise = Dsurvey$hour_rise + (Dsurvey$min_rise/60)
  Dsurvey$month_rise = paste0(format(Dsurvey$risetime,"%m"),"-",format(Dsurvey$risetime,"%Y"))
  Dsurvey$day_rise = format(Dsurvey$risetime,"%d")
  Dsurvey$weekday_rise = weekdays(Dsurvey$risetime)
  # Add clock times to Dsurvey for bedtime
  Dsurvey$hour_bed = Dsurvey$bedtime$hour
  Dsurvey$min_bed = Dsurvey$bedtime$min
  Dsurvey$hour_in_day_bed = Dsurvey$hour_bed + (Dsurvey$min_bed/60)
  Dsurvey$month_bed = paste0(format(Dsurvey$bedtime,"%m"),"-",format(Dsurvey$bedtime,"%Y"))
  Dsurvey$day_bed = format(Dsurvey$bedtime,"%d")
  Dsurvey$weekday_bed = weekdays(Dsurvey$bedtime)
  
  Dsurvey2 = data.frame(date= Dsurvey$date, hr = Dsurvey$hour_in_day_rise, action = rep("rise",nrow(Dsurvey)),
                        col = rep(as.character(coldic2["bed"]),nrow(Dsurvey)), shape = rep(17,nrow(Dsurvey))) #
  Dsurvey3 = data.frame(date= Dsurvey$date, hr = Dsurvey$hour_in_day_bed, action = rep("bed",nrow(Dsurvey)),
                        col = rep(as.character(coldic2["rise"]),nrow(Dsurvey)), shape = rep(17,nrow(Dsurvey))) #
  DS = rbind(Dsurvey2, Dsurvey3)
  
  Dshort = addDoublePlot(Dshort)
  
  textsize.y = 16
  legend.size = 14
  y.breaks = c(0,12,24,36,48)
  # We want to plot a heatmap with on top of that the bedtimes and rising times according to the survey
  # and a separate legend for each.
  # This plot was difficult to create in ggplot directly, but with the help of 
  # package cowplot we are able to put it together.
  
  # First focus on getting the plot itself, getting both legends is not important now:
  p0 = ggplot() +
    geom_tile(data=Dshort,
              aes(date, hour_in_day, fill=status)) + # Add tiles (heatmap)
    geom_point(data=Dsurvey,
               aes(x=date, y=hour_in_day_rise),
               color = as.character(coldic2["rise"])) +  # Add the line plot with risetime
    geom_point(data=Dsurvey,
               aes(x=date, y=hour_in_day_bed),
               color = as.character(coldic2["bed"])) +  # Add the line plot with bedtime
    xlab("") +
    # ylab("Hour in day") +
    theme_bw() +
    scale_colour_manual(name="Passive data", values = coldic) +
    scale_fill_manual(name="Passive data", values = coldic) +
    scale_x_discrete(breaks = dates_on_x_axis) +
    scale_y_continuous(name="Hour in day", breaks=y.breaks) +
    theme(axis.text.x = element_text(angle = 45),
          axis.text.y = element_text(size = textsize.y),
          axis.title.y = element_text(size = textsize.y),
          legend.text = element_text(size = legend.size),
          legend.title = element_text(size = legend.size)) #
    
  
  # Now create new plot for the points only and focus on getting a good legend for it:
  p2 = ggplot() +
    geom_point(data=DS,
               aes(x=date, y= hr, color=action),
               shape = DS$shape) +  # Add the line plot with risetime
      xlab("") + ylab("") +
    # ylab("Hour in day") +
    theme_bw() +
    scale_colour_manual(name="Sleep Diary", values = coldic2) +
    scale_x_discrete(breaks = dates_on_x_axis) +
    scale_y_continuous(name="Hour in day", breaks=y.breaks) +
    theme(axis.text.x = element_text(angle = 45),
          axis.text.y = element_text(size = textsize.y),
          axis.title.y = element_text(size = textsize.y),
          legend.text = element_text(size = legend.size),
          legend.title = element_text(size = legend.size)) #

  # Extract main legends that we need
  leg0 <- get_legend(p0)
  leg2 <- get_legend(p2)
  # Combine the legends:
  leg12 <- plot_grid(leg0, leg2,  ncol = 1)
  # remove legend from p0:
  p0 = p0 +  theme(legend.position="none") 
  # Combine the legends with the main panel
  final_p <- plot_grid(p0,
                       leg12,
                       nrow = 1,
                       align = "l",
                       axis = "t",
                       rel_widths = c(1, 0.2))
        # theme(legend.text = element_text(size = legend.size))
  # write to png file
  png(filename = heatmapsfile ,width = 15, height = 7,units = "in",res = 400)                       
  print(final_p)
  dev.off()
  
  # Heatmap of steps
  dates_on_x_axis = unique(Dlong$date)
  dates_on_x_axis = dates_on_x_axis[seq(1,length(dates_on_x_axis),round(length(dates_on_x_axis) / 25))]

  
  Dlong$date = as.Date(Dlong$date)
  Dlong = addDoublePlot(Dlong)
  Dlong$steps = as.numeric(Dlong$steps)
  png(filename = heatmapsfile_steps ,width = 15, height = 7,units = "in",res = 400)
  myplot = ggplot(Dlong, aes(date, hour_in_day)) +
    geom_raster(aes(fill=steps)) +
    xlab("") +
    ylab("Hour in day") +
    theme_bw() +
    labs(fill = "steps per 30 minutes") +
    scale_fill_continuous(low="black", high="green",na.value = "#000000") + 
    scale_x_discrete(breaks = dates_on_x_axis) +
    theme(axis.text.x = element_text(angle = 45))
  print(myplot)
  dev.off()
}