Testing correlation between data gathered with Firstbeat device against data from Oura cloud (from the ring) and Google Fit cloud (multiple devices)

firstbeat_analyses.R

# It's fake
# Data from cloud but already on computer
# Testing self measurement devices
# CC Lauri Ahonen 2018

require(scales)
library(tidyverse)
library(lubridate)
library(jsonlite)

Sys.setenv(TZ="Europe/Helsinki")


# --- funcs for firstbeat ------------------------------------------------------

read_firstbeat_file <- function(file) {
    ibis <- read_csv(file, skip=5, col_names = FALSE,
                col_types = cols(col_integer()))[[1]]
    sig_start <- read_delim(file,delim = "=", n_max = 1, skip = 2,
                        col_names = FALSE, col_types = cols('c', 'c'))[[2]]
    sig_start <- dmy_hms(sig_start, tz="Europe/Helsinki")
    ts = c(0, cumsum(ibis[1:(length(ibis)-1)]) / 1000) + sig_start
    sig <- as_tibble(list(t = ts, ibi = ibis))
    sig
}

rm_artefacts <- function(data) {
    filt <- data$ibi < 2000
    filtdata <- data[filt,]
    filtdata
}

avg_hrv <- function(data, seg = 300) {
    seg.len = seg
    start <- first(data$t); minute(start) <- 0; second(start) <- 0
    end <- last(data$t); minute(end) <- 0; second(end) <- 0
    hour(end) <- hour(end) + 1
    seg <- seq(from = start, to = end - seg.len, by = seg.len)
    segmat <- matrix(data = seg, nrow = length(seg), ncol = 2, byrow = FALSE)
    segmat[,2] <- segmat[,2] + seg.len
    tnum <- data$t
    ibis <- data$ibi
    sig <- list()
    sig$t <- seg
    sig$ibi <- apply(segmat, 1, function(x)
                    ibi_mean(ibis[tnum > x[1] & tnum < x[2]]))
    sig$hr <- apply(segmat, 1, function(x)
                    ibi_mean(ibis[tnum > x[1] & tnum < x[2]], type = "hr"))
    sig$rmssd <- apply(segmat, 1, function(x)
                    ibi_rmssd(ibis[tnum > x[1] & tnum < x[2]]))

    as_tibble(sig)
}

ibi_rmssd <- function(ibi) {
    sqrt(sum(diff(ibi)^2) / (length(ibi) - 1))
}

ibi_mean <- function(ibi, type = "ibi") {
    if (type == "ibi")
        res <- mean(ibi)
    if (type == "hr")
        res <- mean(6e4 / ibi)
    res
}

# --- funcs for google fit data ------------------------------------------------

read_daily_gdata <- function(file) { # vaati aika paljon purkkaa
    aggtable <- read_csv(file, col_types = cols(.default = "n",
                            'Start time' = "c", 'End time' = "c"))
    date <- ymd(tools::file_path_sans_ext(basename(file)))
    aggtable[['Start time']] <- as.POSIXct(hms(aggtable[['Start time']])+date)
    aggtable[['End time']] <- as.POSIXct(hms(aggtable[['End time']])+date)
    aggtable
}


# --- funcs for ouraring -------------------------------------------------------

trans_odata <- function(ouradata) {
    ouradata %>% by_row(get_timestamps, .to = "t") %>% unnest %>%
    filter(hr_5min > 0)
}

get_timestamps <- function(data) {
    start <- ymd_hms(data$bedtime_start)
    datalen <- length(data$hr_5min[[1]])
    round_date(seq(start, by = 300, length = datalen), "5 minutes")
}

### READING THE DATA ###

# --- read the first beat data -------------------------------------------------

files <- list.files(pattern = '.sdf', full.names = TRUE,
        path = '/Users/laurivaltteri/gdrive/quantified_me/20180923_firstbeat/')

fb_raw_data <- files %>% map_df(read_firstbeat_file) %>% rm_artefacts

fbdata_15min <- avg_hrv(fb_raw_data, 900)
fbdata_5min <- avg_hrv(fb_raw_data, 300)

#p <- ggplot() + geom_line(data=fbdata_5min, aes(t, hr, colour = "5 min")) +
#                geom_line(data=fbdata_15min, aes(t, hr, colour = "15 min")) +
#                scale_x_datetime('time', date_breaks = "6 hours",
#                labels = date_format("%H:%M", tz= 'Europe/Helsinki')) + theme_bw()

p <- ggplot()

## --- read the google fit data ------------------------------------------------

#start <- date(first(firstb_15min$t))
#end <- date(last(firstb_15min$t))
#days <- seq(start, end, by = "day")
f1 <- '/Users/laurivaltteri/gdrive/Takeout/Takeout/Fit/Daily Aggregations/2018-09-23.csv'
f2 <- '/Users/laurivaltteri/gdrive/Takeout/Takeout/Fit/Daily Aggregations/2018-09-24.csv'
f3 <- '/Users/laurivaltteri/gdrive/Takeout/Takeout/Fit/Daily Aggregations/2018-09-25.csv'
f4 <- '/Users/laurivaltteri/gdrive/Takeout/Takeout/Fit/Daily Aggregations/2018-09-26.csv'
files <- c(f1, f2, f3, f4)

gdata <- files %>% map_df(read_daily_gdata)
gdata['t'] <- with_tz(gdata['Start time'], tz = "Europe/Helsinki") # add matching column for time

gdata_with_fb <- left_join(fbdata_15min,gdata)

p <- p + geom_line(data=gdata_with_fb, aes(t, hr, colour = "Firstbeat 15 min average")) +
geom_point(data=gdata_with_fb, aes(t, `Average heart rate (bpm)`,colour = "Google Fit"))

# --- read the oura data -------------------------------------------------------

file <- "/Users/laurivaltteri/gdrive/Takeout/Oura/20180923.json"
ouradata <- fromJSON(file)$sleep
odata <- trans_odata(ouradata)
odata$t <- with_tz(odata$t, tz = "Europe/Helsinki") # some stylizing
odata <- odata %>% rename(rmssd_mean = rmssd)

odata_with_fb <- left_join(fbdata_5min,odata, by = 't')

p <- p + geom_point(data=odata_with_fb, aes(t, hr_5min,colour = "Oura")) +
geom_line(data=odata_with_fb, aes(t, hr, colour = "Firstbeat 5 min average"))


# --- play with the data -------------------------------------------------------

# hr figure to illustrate
p + scale_x_datetime('time', date_breaks = "6 hours",
labels = date_format("%H:%M", tz= 'Europe/Helsinki')) + theme_bw()

# rmssd figure to take back to reality
g <- ggplot() + geom_point(data=odata_with_fb, aes(t, rmssd_5min,colour = "Oura")) +
geom_line(data=odata_with_fb, aes(t, rmssd, colour = "Firstbeat 5 min average")) +
scale_x_datetime('time', date_breaks = "6 hours", labels = date_format("%H:%M",
tz= 'Europe/Helsinki')) + theme_bw()
g

# some correlation testing
filtdata <- filter(odata_with_fb,rmssd < 125 & rmssd > 0 & rmssd_5min > 0)
ggplot(filtdata, aes(hr_5min,hr)) + geom_point() + geom_smooth(method=lm)
cor.test(filtdata$hr_5min, filtdata$hr, method = "pearson", conf.level = 0.95)

filtdata <- filter(odata_with_fb,rmssd < 125 & rmssd > 0 & rmssd_5min > 0)
ggplot(filtdata, aes(rmssd_5min,rmssd)) + geom_point() + geom_smooth(method=lm)
cor.test(filtdata$rmssd_5min, filtdata$rmssd, method = "pearson", conf.level = 0.95)

# for gdata as well
filtdata <- filter(gdata_with_fb,rmssd < 125 & rmssd > 0)
ggplot(filtdata, aes(`Average heart rate (bpm)`,hr)) + geom_point() + geom_smooth(method=lm)
cor.test(filtdata[["Average heart rate (bpm)"]], filtdata$hr, method = "pearson", conf.level = 0.95)