Kubios Time-Varying Results

---

R Programs I used

library(readr) # read CSV files
library(tidyr) #  data manipulation
library(dplyr)# data manipulation
library(lubridate)# timestamp data
library(ggplot2) # plots
library(plotly) # interactive plots
library(viridis) # color palettes
library(viridisLite)

Load Time-Varying Results

# Time-Varying Results
Time.Varing.Results <- read_csv("20190307-7007_hrv.csv", 
                               skip = 148)
# Columns removed
Time.Varing.Results <- Time.Varing.Results[-c(1, 41)]
# Rows kept
Time.Varing.Results  <- Time.Varing.Results[c(2:142890),]
# Renamed variables
Time.Varing.Results<-Time.Varing.Results %>% 
  rename(Artifact=`Artifact percentage`, PNS_index= `PNS index`, SNS_index=`SNS index`,Stress_index=`Stress index`,
         Mean_RR=`Mean RR`, STD_RR=  `STD RR`, Mean_HR=`Mean HR`, STD_HR= `STD HR`,Min_HR= `Min HR`, Max_HR=`Max HR`, 
         VLF_peak=`VLF peak`, LF_peak= `LF peak`,HF_peak= `HF peak`,VLF_power.ms2= `VLF power`,LF_power.ms2=`LF power`,
         HF_powerms2=`HF power`, VLF_power.log= `VLF power_1`, LF_power.log=`LF power_1`,HF_power.log=  `HF power_1`,VLF_power.per=`VLF power_2`,       
         LF_power.per= `LF power_2`,HF_power.per=`HF power_2`, LF_power.nu=`LF power_3`,HF_power.nu= `HF power_3`,LF_HF= `LF/HF ratio`,SD2_SD= `SD2/SD1`,
         DFA_a1=    `DFA a1`,DFA_a2=`DFA a2`)               

# Check if data are coded as characters
sapply(Time.Varing.Results , class) 

         Time      Artifact     PNS_index     SNS_index  Stress_index       Mean_RR        STD_RR       Mean_HR        STD_HR        Min_HR        Max_HR 
  "character"   "character"     "numeric"     "numeric"     "numeric"   "character"   "character"   "character"   "character"   "character"   "character" 
        RMSSD          NNxx         pNNxx         HRVti          TINN      VLF_peak       LF_peak       HF_peak VLF_power.ms2  LF_power.ms2   HF_powerms2 
  "character"   "character"   "character"     "numeric"   "character"   "character"   "character"   "character"   "character"   "character"   "character" 
VLF_power.log  LF_power.log  HF_power.log VLF_power.per  LF_power.per  HF_power.per   LF_power.nu   HF_power.nu         LF_HF           EDR           SD1 
  "character"   "character"   "character"   "character"   "character"   "character"   "character"   "character"     "numeric"   "character"   "character" 
          SD2        SD2_SD          ApEn        SampEn        DFA_a1        DFA_a2 
  "character"     "numeric"     "numeric"     "numeric"     "numeric"     "numeric" 

Time & HRV data

Fundamentally, there are two types of data when using time-varying results: (1) HRV and (2) Time. The following code breaks down the file into these two types and then combines them into a single data frame. My usual program for timestamp data (lubridate) did not work. The time data did not parse time data into hours, minutes, and seconds for every row. Consequently, I made my time calculations of seconds and hours.

## HRV Data
HRV<-Time.Varing.Results %>%
  select(Artifact:DFA_a2) %>% 
  mutate_if(is.character,as.numeric)

# Time Data
Time<-Time.Varing.Results %>%
  select(Time)
# First,I created a vector called "Seconds." Because there are 142,889 observations and the sampling rate is in seconds,there are 142,889 seconds in this sample.
Seconds<-1:142889
Time<-cbind(Time,Seconds)

# Hours
Time<-Time %>% 
  mutate(Hours=Seconds/3600)

# Combined HRV & Time Data
Time.Varying.Data<-cbind(Time,HRV)

Heart Rate Data

At the root of all HRV variability statistics are the RR intervals. You can convert these RR intervals into beat-per-minute (bpm). I strongly recommend using this metric because it’s intuitive to most people and often used in many scientific disciplines.

I also created a new variable called “Threshold.” This variable helps determine when there was a deviation in heart rate. I defined this deviation empirically: 2 standard deviations above the mean. Feel free to use a different method if you find information in the literature that is more informative. Importantly, we can determine how long the event occurred by knowing when heart rate was above a defined threshold.

# Mean HR Histogram
Time.Varying.Data %>%  
  ggplot(aes(x=Mean_HR))+
  geom_histogram(bins = 100,color = "black")+
  theme_bw()+ 
  labs(
    title    = "Mean HR", 
    y        = "Count",
    x        = "bpm")

# Descriptive Statistics
Descriptives<-Time.Varying.Data  %>% 
  summarise(mean=mean(Mean_HR),median=median(Mean_HR), 
            min=min(Mean_HR), max=max(Mean_HR),SD= sd(Mean_HR), n =n()) %>% 
  mutate(SD2= mean+(2*SD))
  
# Threshold 
Time.Varying.Data<-Time.Varying.Data  %>% 
  mutate(Threshold= ifelse(Mean_HR >=Descriptives$SD2,"Above","Normal"))

Ploting time-varying data using ggplot2

I strongly recommend plotting your heart rate (x = time, y = heart rate). You will be able to determine if a significant event occurred.

I prefer using ggplot2 program in R for data visualization. There are many online tutorials out there that are easy to follow. Plus, it is part of the tiddyverse (https://www.tidyverse.org) set of programs.

# Scatter plot: Heart Rate
Time.Varying.Data%>%
  ggplot(aes(x=Hours, y=Mean_HR))+
  geom_hline(yintercept=Descriptives$SD2, color ="red", linetype="dashed", alpha =.75)+
  geom_hline(yintercept = Descriptives$mean, color ="green")+ # Mean
  geom_line()+
  theme_bw()+
  annotate("text", label = "2 SD Above Mean", x = 3, y = 130, size = 4, color = "red")+
  annotate("text", label = "Mean", x = 3, y = 60, size = 4, color = "darkgreen")+
  ylim(40,200)+
  labs(                         
    title    = "Time & Mean HR", 
    y        = "bpm",
    x        = "Hours")

# Scatter Plot SNS Index
Time.Varying.Data%>%
  ggplot(aes(x=Hours, y=SNS_index))+
  geom_line()+
  theme_bw()+
  labs(                         
    title    = "Time & SNS Index", 
    y        = "SNS Index",
    x        = "Hours")

  
# Heart Rate and SNS Index
Time.Varying.Data %>% 
  ggplot(aes(x=Hours, y=Mean_HR, color =SNS_index))+
  geom_hline(yintercept=Descriptives$SD2, color ="red", linetype="dashed", alpha =.75)+
  geom_hline(yintercept = Descriptives$mean, color ="green")+ # Mean
  scale_color_viridis(option = "B")+ # Inferno Color Palette
  geom_line(size =1)+
  theme_dark()+
  annotate("text", label = "2 SD Above Mean", x = 3, y = 130, size = 4, color = "red")+
  annotate("text", label = "Mean", x = 3, y = 60, size = 4, color = "green")+
  ylim(40,200)+
  labs(
    color = "SNS Index",
    title    = "Time & Mean HR", 
    y        = "bpm",
    x        = "Hours")

Ploting time-varying data using Plotly

When I explore my data, I sometimes like to use this program called “Plotly.” The advantage of this program is it can be interactive. That is, if you hover your cursor over the values, you can see the x- and y- values. This feature is helpful if you are quickly trying to see when an event occurred. Another interesting feature with Plotly: you can create 3D figures.

# 2D 
Time.Varying.Data%>% 
plot_ly(x = ~Hours, y = ~Mean_HR,type = 'scatter', mode = 'lines', name ='bpm') %>% 
  add_segments(x = 0, xend = 40, y = 118.9955, yend = 118.9955, name='2 SD Above Mean') %>% 
  add_segments(x = 0, xend = 40, y = 70.00502, yend = 70.00502, name='Mean') %>% 
  layout(
    title = 'Time & Mean Heart Rate',
    xaxis = list(title = 'Hours'),
    yaxis = list(title = 'bpm'))

# 3D  
Time.Varying.Data%>% 
  plot_ly(x = ~Hours, y = ~Mean_HR,z = ~SNS_index,type = 'scatter3d', mode = 'lines')