Extracting the data

# remotes::install_github("mathesong/rwithings")

library(tidyverse)
library(rwithings)

Now that we have everything loaded, we can begin. First, we must authenticate ourselves and get a token.

client_id <- MY_CLIENT_ID            # Secret!
client_secret <- MY_CLIENT_SECRET    # Secret!

token <- withings_auth(client_id, client_secret)

And now we can use that to download some data using the rwithings package.

sleepdat <- getsleepsummary(token, startdate = "2018-10-01", 
                            enddate = "2018-12-10")

Cleaning and Analysis

And now we can pull out the bits that we want: how long I spent awake and in light and deep sleep. I also found a nifty function in the chron package to say whether a day was a weekend day or not in case that was a factor.

sleepsummary <- sleepdat$body$series %>%
  select(date, ends_with("duration")) %>%
  rename(Date=date,
         Awake = data.wakeupduration,
         LightSleep = data.lightsleepduration,
         DeepSleep = data.deepsleepduration) %>%
  mutate(Date=as.Date(Date),
         Weekend = chron::is.weekend(Date))

I noticed some dates duplicated in this data, and this corresponded with naps.

sleepsummary %>% 
  filter(Date == "2018-11-11")

##         Date Awake LightSleep DeepSleep Weekend
## 1 2018-11-11  1680      12960     14040    TRUE
## 2 2018-11-11   300       5460      6420    TRUE

I solved this by removing the second duplicate of each day’s sleep.

sleepsummary <- sleepsummary[!duplicated(sleepsummary$Date, fromLast = F),]

sleepsummary %>% 
  filter(Date == "2018-11-11")

##         Date Awake LightSleep DeepSleep Weekend
## 1 2018-11-11  1680      12960     14040    TRUE

I wanted to look at proportional sleep in each of the different states, so I looked at proportions of the total sleep duration for each night.

sleepprops <- sleepsummary %>%
  mutate(TotalSleep = Awake + LightSleep + DeepSleep) %>%
  gather(SleepType, Minutes, -Date, -TotalSleep, -Weekend) %>%
  mutate(SleepProportion = Minutes/TotalSleep,
         SleepType = factor(SleepType, levels = c("Awake",
                                                  "LightSleep",
                                                  "DeepSleep")))

Plotting

I decided to use this oppotunity to play with the hrbrthemes themes that I’ve seen so much talk of, which just look gorgeous! I also picked some colours using the Color Calculator.

# remotes::install_github("hrbrmstr/hrbrthemes")
library(hrbrthemes)

sleep_proportions <- ggplot(sleepprops, aes(x=Date, y=SleepProportion, 
                                            colour=SleepType)) +
  geom_point(aes(shape=Weekend), size=3) +
  geom_line() +
  geom_vline(xintercept = as.Date("2018-11-30"),
             linetype="dotted", size=1) +
  geom_smooth(data=filter(sleepprops, 
                          Date<=as.Date("2018-11-30"))) +
  facet_grid(SleepType ~ ., scales = "free") +
  scale_shape_manual(values=c(19, 1)) +
  theme_ipsum_rc() +
  labs(x="Date", y="Sleep Stage Proportions",
       title="My Sleeping Patterns Before my PhD Defence",
       subtitle="Nightly proportional sleep stages in the two months prior",
       caption="Recordings were made using a Withings Steel HR") +
  scale_color_manual(values=c("#ffd631", "#31d6ff", "#d12869"))

sleep_proportions

Conclusions

So it seems that the proportion of hours spent awake or in light sleep were increasing, while my deep sleep was decreasing over the course of the months preceding my defence. This is exactly what one might have predicted! I must admit that I didn’t expect this to actually come out, as my two-year-old daughter decided to start waking up in the middle of the night from time to time in the weeks preceding the defence, but this didn’t get any better or worse in the last days preceding the defence itself, so I don’t think it explains the change, but might contribute to the trend being a little noisier.

Extracting the data

I first extracted the data as a CSV, and read it into R using the rwithings package function for munging this particular data.

hr <- read_csv_startdurval("../../static/data/20181213_ThesisDefence/raw_tracker_hr.csv") %>% 
  rename(Start = start, Duration = duration, Freq = value,
         Date = date, End = end, Mid = mid)

## Date in ISO8601 format; converting timezone from UTC to "Europe/Berlin".

During the defence

During the defence, I set my watch to “workout mode”, which means that it records heart rate data at a higher frequency. For this analysis, I thought a good strategy would be to calculate the mean heart rate for each minute, weighted by the duration of the measurement interval. What this means is that I first want to limit the data set to the time of the defence, and then calculate minute means.

library(lubridate)
tz <- Sys.timezone()

dday <- hr %>%
  filter(Date=="2018-11-30") %>%
  filter(Start > as.POSIXct("2018-11-30 09:15:00", tz = tz)) %>%
  filter(Start < as.POSIXct("2018-11-30 14:15:00", tz = tz)) %>%
  mutate(Minute = floor_date(Mid, unit="minute"))

minsummary <- dday %>%
  group_by(Minute) %>%
  summarise(MinuteFreq = weighted.mean(Freq, w=Duration))

dday <- left_join(dday, minsummary)

labels <- tibble(
          Time = c("2018-11-30 09:33:00", "2018-11-30 09:37:00",
                   "2018-11-30 09:58:00", "2018-11-30 10:43:00",
                   "2018-11-30 12:03:00", "2018-11-30 12:05:00", 
                   "2018-11-30 12:14:00", "2018-11-30 12:48:00", 
                   "2018-11-30 13:17:00", "2018-11-30 13:41:00",
                   "2018-11-30 13:56:00", "2018-11-30 13:58:00"),
         Event = c("Chairperson introduces\nthe defence",
                   "Opponent's presentation\nbegins", 
                   "My presentation\nbegins",
                   "Opponent's questions\nbegin", 
                   "Opponent's questions\nend",
                   "Short break", 
                   "Committee questions\nbegin",
                   "Committee questions\nend", 
                   "Walk to\nreception\n(15 min)", 
                   "Awaiting\ndecision",
                   "Committee\nreturns", 
                   "I passed!")) %>%
  mutate(Freq = 120,  # This is just for the annotation position
         MinuteFreq = Freq,
         Time = as.POSIXct(Time),
         Minute = Time) 

library(ggrepel)

dday_minuteplot <- ggplot(dday, aes(x=Minute, y=MinuteFreq)) +
  geom_segment(data=labels, aes(x=Minute, xend=Minute, 
                                y=120, yend=25), size=0.4) +
  geom_segment(aes(x=Minute, xend=Minute, y=MinuteFreq, yend=25,
                   colour=MinuteFreq), size=0.8) +
  scale_colour_viridis_c("Frequency\n(bpm)", option = "A") +
  scale_x_datetime(breaks = "1 hour", date_labels = "%I:%M %p") +
  theme_ipsum_rc() +
  labs(x="Time", y="Heart Rate (bpm)",
       title="Heart rate during my PhD Defence",
       subtitle="Weighted mean frequencies recorded for each minute",
       caption="Recordings were made using a Withings Steel HR") +
  geom_point(data=labels) +
  geom_text_repel(
    angle        = 30,
    segment.size = 0.2,
    ylim = c(125,NA),
    data=labels, aes(label=Event),
    size=3
  ) +
  ylim(25, 140)

dday_minuteplot

Heart Rate Distributions

Next, I wanted to look at the daily distributions of my heart rate in the days before and after. There are several complications to this:

• Waking and sleeping heart rates are quite different
• Heart rates are typically right-skewed distributions, so a mean doesn’t cut it, and the median depends on the amount of time spent exercising etc
• The watch does not record heart rate at a constant sampling rate, and this is particularly different between “workout mode” (during the defence, as well as during exercise sessions) and normal tracking, so I can’t just use a summary statistic on the raw data, since highly-sampled periods are also likely to be those with higher heart rate

What I aimed to do was to separate waking and sleeping times from the sleep data, show the whole distribution of heart rate (and cut it in the figure at 120 during waking hours because above there is pretty much only exercise) as a ridge plot. I decided to use medians as a summary statistic: while mode would probably be more appropriate physiologically, it’s just not easy to calculate in a robust fashion. Getting a physiologically meaningful summary was not the goal, but rather having a number to make comparison easier, for which the median should do just fine. In order to deal with the different sampling rates, I opted to linearly interpolate heart rate over the day into increments of 5 seconds so that highly-sampled and sparsely-sampled times get as much as importance/coverage as one another. I also decided to only look at a month before in the plot in order not to make it too busy.

First, I interpolate the data into 5 second intervals. I do this by calculating the number of seconds, dividing it by 5, and then interpolating into that many points. I’ll first set everything as “Waking”, and then change the sleeping times to “Sleeping” later.

hr_ridgedat <- filter(hr, Date >= as.Date("2018-11-01"))

n_periods <- abs(as.numeric(
  difftime(
    hr_ridgedat$Mid[1],   # First
    hr_ridgedat$Mid[nrow(hr_ridgedat)], # Last
    units = "s")))/5  # no. seconds divided by 5 

hr_interp <- tibble(Time = seq(hr_ridgedat$Mid[1],  # First
                               hr_ridgedat$Mid[nrow(hr_ridgedat)],  # Last
                               length.out = n_periods),
                    Waking = "Waking",
                    Date = as.Date(Time))

hr_interp$Freq <- approx(x = hr_ridgedat$Mid,  # interpolate
                         y = hr_ridgedat$Freq,
                         xout = hr_interp$Time, method = "linear")$y

Next, I’ll put together a set of sleep times. Frustratingly, if I went to bed too late (after ~02h00), the watch tended to assume I didn’t sleep at all. I did in fact sleep every night before the defence, and never after 03h00, so I just set my sleeping times for those days in which it didn’t record sleep as being from 03h00. For waking up, however, I can be pretty sure it was around 06h00 each morning, because that’s when my 2 year old daughter wakes up like clockwork!

sleeptimes <- sleepdat$body$series %>% 
  select(Start=startdate, End=enddate, Date=date) %>%
  mutate(Date=as.Date(Date)) %>% 
  filter(Date >= as.Date("2018-11-01"))

dateseq <- tibble::tibble(Date=seq(sleeptimes$Date[1], Sys.Date(),by = 1) )

sleeptimes <- full_join(dateseq, sleeptimes)

## Joining, by = "Date"

measuredsleep <- filter(sleeptimes, !is.na(Start))

missingsleep <- filter(sleeptimes, is.na(Start)) %>%
  mutate(Start = as.POSIXct(strptime(paste(Date, "3:00"),
                   format = "%Y-%m-%d %H:%M"), tz = Sys.timezone()),
         End = as.POSIXct(strptime(paste(Date, "6:00"),
                   format = "%Y-%m-%d %H:%M"), tz = Sys.timezone()))

sleeptimes <- bind_rows(measuredsleep, missingsleep)

Now I can use this data to define the sleeping times. Also, I need to assign sleeping heart rate to the correct days, which does not always correspond to the time recorded: if I was sleeping at 23h00 on day 1, the sleep should be assigned to day 2 along with the rest of my sleep until 06h00 on day 2. This I can also do using the sleeptimes data.

I couldn’t readily think of a tidy way of doing this operation, so I used a for loop. Sadface.

for(i in 1:nrow(sleeptimes)) {
  
  hr_interp$Waking <- ifelse(hr_interp$Time > sleeptimes$Start[i] & 
                               hr_interp$Time < sleeptimes$End[i],
                             yes = "Sleeping", no = hr_interp$Waking)
  
  hr_interp$Date <- ifelse(hr_interp$Time > sleeptimes$Start[i] & 
                             hr_interp$Time < sleeptimes$End[i],
                             yes = sleeptimes$Date[i], no = hr_interp$Date)
  
}

saveRDS(hr_interp, file = "../../static/data/20181213_ThesisDefence/hr_interp.rds")

Now, we want to divide the data into waking and sleeping. For some reason, the Date column has now become a double, so we’ll also convert it to date again.

hr_interp <- readRDS("../../static/data/20181213_ThesisDefence/hr_interp.rds")

hr_interp <- hr_interp %>%
  mutate(Date = as.Date(Date, origin="1970-01-01", tz = Sys.timezone()))

hr_interp_waking <- filter(hr_interp, Waking=="Waking")

hr_interp_sleep <- filter(hr_interp, Waking=="Sleeping")

Let’s quickly make some median summaries for both.

hr_interp_waking_summary <- hr_interp_waking %>%
  group_by(Date) %>%
  summarise(MedFreq = median(Freq),
            MedFreq_round = formatC(MedFreq,1,format="f"))

hr_interp_sleep_summary <- hr_interp_sleep %>%
  group_by(Date) %>%
  summarise(MedFreq = median(Freq),
            MedFreq_round = formatC(MedFreq,1,format="f"))

Waking Heart Rate Distributions

So, first, let’s look at my heart rate each day, including the defence day. I’ll make a couple of annotations for the plot. And I’ll also cut everything off at 120 bpm. And we can show the medians on the left.

library(ggridges)

## 
## Attaching package: 'ggridges'

## The following object is masked from 'package:ggplot2':
## 
##     scale_discrete_manual

annotations_waking <- tribble(
  ~Date,   ~Label,
  "2018-11-30", "Defence day!") %>%
  mutate(Date = as.Date(Date))


hr_ridgeplot <- ggplot(hr_interp_waking, 
                       aes(x=Freq, y=Date, group=Date, fill=..x..)) +
  geom_density_ridges_gradient(scale = 1.7, rel_min_height = 0.01, 
                               gradient_lwd = 1.) +
  xlim(40, 120) +
  scale_fill_viridis_c(name = "Heart\nRate\n(bpm)", option = "A") +
  labs(title = 'Waking Heart Rate',
       subtitle = "Distributions and medians of interpolated heart rate during waking hours",
       caption = "Recordings were made using a Withings Steel HR",
       x = "Heart Rate (bpm)") +
  theme_ridges(font_size = 13, grid = TRUE) + 
  theme(axis.title.y = element_blank()) +
  geom_point(data = hr_interp_waking_summary, aes(x = MedFreq, y=Date+0.25), 
             colour="white", size=1.5) +
  theme_ipsum_rc() +
  geom_text(data=hr_interp_waking_summary, 
            aes(x = 45, y = Date+0.5, label=MedFreq_round), 
            size=3, color = "red") +
  geom_text(data=annotations_waking, aes(x = 58, y = Date+0.5, label=Label), 
            size=3, color = "red")

hr_ridgeplot

## Picking joint bandwidth of 1.72

Conclusions

I must admit, I was quite surprised by how clear the increase for the defence day is compared to the rest! The median is clearly higher than the mode in almost all cases, but as a summary statistic for comparing daily heart rates, it does the job.

Sleeping Heart Rate Distributions

Next, let’s look at my heart rate during sleep. For each day, the sleep is that which ends on that particular day.

annotations_sleep <- tribble(
  ~Date,   ~Label,
  "2018-11-02", "Thesis sent to printers",
  "2018-11-30", "Night before",
  "2018-12-01", "It's over!",
  "2018-12-02", "Recovery"
) %>%
  mutate(Date = as.Date(Date))


hr_ridgeplot_sleep <- ggplot(hr_interp_sleep, 
                             aes(x=Freq, y=Date, group=Date, fill=..x..)) +
  geom_density_ridges_gradient(scale = 1.7, rel_min_height = 0.01, 
                               gradient_lwd = 1.) +
  xlim(40, 90) +
  scale_fill_viridis_c(name = "Heart\nRate\n(bpm)", option = "A") +
  labs(title = 'Sleeping Heart Rate',
       subtitle = "Distributions and medians of interpolated heart rate during sleeping hours",
       caption = "Recordings were made using a Withings Steel HR") +
  theme_ridges(font_size = 13, grid = TRUE) + 
  theme(axis.title.y = element_blank()) +
  geom_point(data = hr_interp_sleep_summary, aes(x = MedFreq, y=Date+0.25), 
             colour="white", size=1.5) +
  theme_ipsum_rc() +
  geom_text(data=hr_interp_sleep_summary, 
            aes(x = 42, y = Date+0.5, label=MedFreq_round), 
            size=3, color = "red") +
  geom_text(data=annotations_sleep, aes(x = 85, y = Date+0.5, label=Label), 
            size=3, color = "red")

hr_ridgeplot_sleep

## Picking joint bandwidth of 0.816

Conclusions

I couldn’t believe this when I saw it! Even looking at sleeping heart rate was a bit of an afterthought. The night before can obviously be attributed to nerves, the night of the defence was presumably the result of masses of adrenaline coursing through my body throughout the day having a lasting impact into the night, and I love that the next night is one of the lowest of the nights recorded, in recovery. I was a little surprised about the night on November 2nd, and I had to look it up: it turned out it was the night after I had sent the final version of my thesis to the printers to be printed: I had nightmares all night that I’d forgotten something critical, misspelled a name of a committee member, or that I’d left something out of the main description pages. I must say, I was really surprised by how well the watch was able to pick this all up.