Quantified Self Experiments / Data Visualization - Oura ring sleep patterns

Clustering techniques may provide interesting insights into patterns of our biomarkers. Here is simple example on clustering oura ring 2.0 sleep data.

Some of oura sleep data seems to be not accurate, so i'm focusing only on a few biomarkers: breathing rate, nightly HR and HRV, temperature deviation, time in bed. Analyzed dataset contains 667 nights.

Here i begin with 2 simple algorithmsfor clustering: K-means (KNN) clustering and Hierarchical clustering (HC). Both algorithms have different settings, but number of clusters and point distance method made a powerful influence on results. Knowledge about clustering i've got from this beatiful data analysis textbook (p. 516)

I've started with 3 clusters and pearson correlation as method for distance and got interesting results. Radar plot is intuitive way to visualize clusters if there is 3-8 variables of interest. The first run gave me interesting results:

Left: Hierarchical clustering. Right: K-means clustering. Numbers are max (1st line) and min (2nd line)

Red and green clusters looks same for both HC and KNN:

Looking at KNN results we can see a clear patters of sick night on one side and restorative night on other side and usual night beetween them. Running analysis with more clusters or euclidean distance method did not produce additional meaningful clusters, so i've finished with 3 clusters and correlation as distance.

As we identified 3 night patterns, lets describe them in numbers.

Restorative (173) Usual (466) Sick (28 nights)
HR, bpm 46.4 [44.1,47.4] 51.2 [49.9,52.8] 57 [55.4,60]
HRV, ms 69 [62,75] 49 [45,53] 40 [35,45]
temperature delta, C -0.09 [-0.25,0.05] 0 [-0.1,0.13] 1.02 [0.63,1.41]
breathing rate, rpm 12.6 [12.4,12.9] 13.1 [12.9,13.4] 14.6 [13.6,15.3]
time in bed, h 8.5 [8.1,9] 8.8 [8.2,9.4] 9.5 [8.6,10.1]

Data represented as median and IQR [q25,q75].

On the sick night my temperature raise by more than +0.63C, i sleep ~9.5 hours with high HR>55 and low HRV <45. Restorative night not differ from usual night on sleep duration (~8.5h) and temperature (<0.15C) and breathing rate (~12.9). The main difference in restorative night compared to usual night is high HRV (+25ms) and lower HR (-5 bmp).

These results may seen as obvious, but i think it's always a good idea to move from subjective feelings of a good night to objective measurements and classification. Also these results provide objective thresholds which can help categorise 5 different measures into 3 simple categories: restorative, usual and sick.

These results should not be generalized and might be person-specific. Soon i'm going to analyse my sleep stages patterns by using Dreem 2 EEG data.

Data availability & Information

Welcome for questions, suggestions and critics in comments below.

Original unmodified (exported) oura raw data for oura is here.

R Code (raw data download make take a while, file is ~10MB):

library(mise); mise()
RNGkind(sample.kind = "Rejection")
options(scipen=999)
describe.data <- function(df, digits = 2) {
  for(i in 1:length(colnames(df))) {
    x <- df[,colnames(df)[i]]
    if(is.numeric(x)) {
      print(paste0(colnames(df)[i], "(", length(x) , "): ",round(mean(x),digits),"±",round(sd(x),digits)," ", round(median(x),digits),"[",round(quantile(x,.25),digits),",",round(quantile(x,.75), digits),"] ", round(min(x),digits),"-",round(max(x),digits))) 
    }
  }
}
radar.graph <- function(df, n.cols, title)
{
  areas <- NULL
  set.seed(7)
  for(i in 1:nrow(df)) areas <- cbind(areas, rgb(rand(1), rand(1), rand(1), 0.33))
  library(fmsb)
  radarchart(df, maxmin = F, title = title, cglty = 2, cglcol = "gray", pcol = 1:n.cols, plwd = 2, pfcol = areas)
}
library(jsonlite)
url <- 'https://blog.kto.to/uploads/'
cache.filename <- 'oura_data_2022-08-19T03-36-46.json'
rawoura <- fromJSON(paste0(url, cache.filename))
sleep<-flatten(data.frame(rawoura$sleep))
cols <- c("temperature_delta", "rmssd", "breath_average", "hr_average", "duration")
library(dplyr)
ouratemp <- sleep %>% select(all_of(cols))
df <- na.omit(ouratemp)
df$rmssd <- as.numeric(df$rmssd)
df$duration <- df$duration / 3600
df.ori <- df
set.seed(1)
#1 decide if features should be standardized: usually when units are differs
scale.data <- T; scale.text <- "Unscaled features"
if(scale.data) {
  scale.text <- "Scaled features"
  for(i in 1:length(cols)) {
    df[, cols[i]] <- scale(df[, cols[i]])
  }
}
# 2 decide dissimilarity measure: 
# 2.1 euclidean distance: focuses on magnitude
# 2.2 correlation: correlation between observations for 1-n variables, focuses on shapes of observation profile rather than magnitude 
dist.method <- "pearson"
#dist.method <- "euclidean"
library(factoextra)
df.dist <- get_dist(x = df, method = dist.method)
# 3 decide for linkage of dendrogram: complete & average often prefered over single.
hc.methods <- c("average", "complete") #, "single")
#4 decide where to cut dendrogram and how many clusters for KNN
n.clusters <- 3
#make a matrix for KNN
mx <- as.matrix(df)
km.out <- kmeans (mx, n.clusters, nstart = 50)
km.out$tot.withinss
for(i in 1:length(hc.methods)) {
  hc.final <- hclust(df.dist, method = hc.methods[i])
  hc.clusters <- cutree(hc.final, n.clusters); table(hc.clusters)
}
df.ori$knn.cluster <- km.out$cluster
df.ori$hc.cluster <- hc.clusters
for(i in 1:n.clusters) {
  print(paste0("Cluster #", i, ", KNN"))
  print(describe.data(df.ori[df.ori$knn.cluster == i, cols]))
}
for(i in 1:n.clusters) {
  print(paste0("Cluster #", i, ", HC"))
  print(describe.data(df.ori[df.ori$hc.cluster == i, cols]))
}
df.r.hc <- aggregate(cbind(df.ori[, cols[1:length(cols)]]), by = list(cluster = df.ori$hc.cluster), FUN = median)
df.r.hc <- df.r.hc[order(df.r.hc$temperature_delta, decreasing = T), ]
df.r.hc$cluster <- NULL
par(mfrow = c(1, 2))
radar.graph(df.r.hc, length(cols), "HC Clusters")
df.r.knn <- aggregate(cbind(df.ori[, cols[1:length(cols)]]), by = list(cluster = df.ori$knn.cluster), FUN = median)
df.r.knn <- df.r.knn[order(df.r.knn$temperature_delta, decreasing = T), ]
df.r.knn$cluster <- NULL
radar.graph(df.r.knn, length(cols), "KNN Clusters")

Statistical analysis

RStudio version 2022.07.1 and R version 4.2.1.