This is an exploratory data analysis of the mHealth data from:

Digital Health: Tracking Physiomes and Activity Using Wearable Biosensors Reveals Useful Health-Related Information Xiao Li , Jessilyn Dunn , Denis Salins , Gao Zhou, Wenyu Zhou, Sophia Miryam Schüssler-Fiorenza Rose, Dalia Perelman, Elizabeth Colbert, Ryan Runge, Shannon Rego, Ria Sonecha, Somalee Datta, Tracey McLaughlin, Michael P. Snyder Published: January 12, 2017 https://doi.org/10.1371/journal.pbio.2001402

The argument that I’m trying to make is that the SeeClickFix or crime data are similar to those from the Snyder data. They have seasonal variability and lots of things changing, so it can be difficult to isolate the effect of an intervention. Let’s take the example of an exercise intervention, where individuals are asked to increase their activity and monitor the effect on some other outcome, like heart rate. You can use the variability to which different individuals followed the intervention as a means to identify the signal from that intervention.

The two that I expect to relate are MOVES and heart rate.

Activity is likely to have an effect on heart rate, but are there some activities that have a larger effect than others? Can we figure out which is most effective at reducing heart rate by comparing how much the different activities are performed by different people?

This is similar to SeeClickFix and Crime, in that there are many different categories of posts, and they may effect crime to varying degrees.

I suppose I should just check a standard regression after an intervention and see if the amount of seeclickfix use was significant.

These are the devices for which data are available on the paper’s website:

HR and SpO2 data collected using four devices (Scanadu Scout, iHealth-finger, Masimo, and Basis) were very close to that of the WA instrument over a wide range of values using the Bland–Altman method of comparison [16,17]

So I’ll compare the MOVES data, as an analog of SeeClickFix, to heart rate, collected by any of four devices. First I’ll check out the MOVES data.

files <- list.files("wearables/Subject1_rawdata/MOVES/", full.names = TRUE)
head(files)
## [1] "wearables/Subject1_rawdata/MOVES//activities.csv"
moves <- read.csv(files)

There are a series of different activites with dates and duration time stamps.

summary(moves)
##        Date             Activity           Group      
##  3/18/15 :   48   walking   :14011            :    7  
##  7/24/14 :   43   transport : 1874   cycling  :  880  
##  8/23/14 :   43   cycling   :  880   running  :  184  
##  3/24/15 :   40   running   :  184   transport: 1923  
##  1/5/16  :   39   airplane  :   49   walking  :14011  
##  12/27/15:   39   basketball:    3                    
##  (Other) :16753   (Other)   :    4                    
##                        Start                              End       
##  2014-09-18T12:32:55-07:00:    2   2014-03-21T15:26:27-07:00:    1  
##  2014-03-21T15:26:22-07:00:    1   2014-03-21T16:04:14-07:00:    1  
##  2014-03-21T16:02:33-07:00:    1   2014-03-21T16:26:22-07:00:    1  
##  2014-03-21T16:05:20-07:00:    1   2014-03-21T16:30:20-07:00:    1  
##  2014-03-21T16:26:43-07:00:    1   2014-03-21T17:15:32-07:00:    1  
##  2014-03-21T17:15:04-07:00:    1   2014-03-21T18:19:28-07:00:    1  
##  (Other)                  :16998   (Other)                  :16999  
##     Duration          Distance            Steps           Calories    
##  Min.   :    0.0   Min.   :   0.000   Min.   :   0.0   Min.   :  0.0  
##  1st Qu.:   90.0   1st Qu.:   0.040   1st Qu.:  36.0   1st Qu.:  2.0  
##  Median :  203.0   Median :   0.138   Median : 173.0   Median :  8.0  
##  Mean   :  513.9   Mean   :  19.219   Mean   : 334.2   Mean   : 17.6  
##  3rd Qu.:  434.0   3rd Qu.:   0.367   3rd Qu.: 403.0   3rd Qu.: 18.0  
##  Max.   :50207.0   Max.   :7572.527   Max.   :8170.0   Max.   :567.0  
##

Each date has ~40 entries, comprised of different activities like cycling and running. I’ll collapse all of the different activities into a single entry per day.

moves1 <- moves %>%
  ddply(c("Date", "Group"), summarize, 
          duration = sum(Duration),
          distance = sum(Distance),
          calories = sum(Calories))

# I'm going to exclude this for now, since "transport" is also in the list, and 
# the total amount of moving + transport time doesn't really make sense. I should do the total active time.
movesActive <- moves1[-grep("transport", moves1$Group),]

movesActive <- movesActive %>%
  ddply("Date", summarize,
          duration = sum(duration),
          distance = sum(distance),
          calories = sum(calories))
  
movesActive$Group <- "active"

moves2 <- rbind(moves1, movesActive)

moves2$Date <- as.Date(as.character(moves2$Date), format = "%m/%d/%y")

moves2 <- moves2[order(moves2$Date),]

moves2 %>%
  ggplot(aes(x = Date, y = duration)) +
    geom_line(aes(group = Group, color = Group))

Why is that so spikey? Maybe I have to smooth by collapsing to weeks? Also, what’s that red line without a label?

table(moves2$Group)
## 
##             cycling   running transport   walking    active 
##         6       321       167       529       740       740
moves2[grep("^$", moves2$Group),]
##            Date Group duration distance calories
## 1647 2014-09-15           5100        0      351
## 1662 2014-09-18           2100        0      199
## 1666 2014-09-19           3900        0      323
## 167  2014-10-12           1200        0      140
## 553  2014-12-31           1800        0       93
## 1660 2015-09-17           2100        0      199

Those have a huge duration! Since I don’t know what they are I’ll remove them.

moves2 <- moves2[-grep("^$", moves2$Group),]

moves2 %>%
  ggplot(aes(x = Date, y = duration)) +
    geom_line(aes(group = Group, color = Group))

Still too spikey. I’ll aggregate to week.

moves2$week <- format(moves2$Date, format = "%Y-%U")

moves3 <- moves2 %>%
  ddply(c("week", "Group"), summarize,
          duration = sum(duration),
          distance = sum(distance),
          calories = sum(calories))

moves3 %>%
  ggplot(aes(x = week, y = duration)) +
    geom_line(aes(group = Group, color = Group)) +
    theme(
      axis.text.x = element_blank()
    )

What if I exclude transport?

moves_notrans <- moves3[-grep("transport", moves3$Group),]

pmoves <- moves_notrans %>%
  ggplot(aes(x = week, y = duration)) +
    geom_line(aes(group = Group, color = Group)) +
    theme(
      axis.text.x = element_blank()
    ) +
  labs(title = "MOVES data",
         subtitle = "Activity over time") 
pmoves

How about by month?

moves2$month <- format(moves2$Date, format = "%Y-%m")

moves4 <- moves2 %>%
  ddply(c("month", "Group"), summarize,
          duration = sum(duration),
          distance = sum(distance),
          calories = sum(calories))

moves4 %>%
  ggplot(aes(x = month, y = duration)) +
    geom_line(aes(group = Group, color = Group)) +
    theme(
      axis.text.x = element_blank()
    ) +
  labs(title = "MOVES data",
         subtitle = "Activity over time")

I like the weeks best.

Ok, now I need the heart rate information. First I’ll check masimo.

files <- list.files("wearables/Subject1_rawdata/Masimo/", full.names = TRUE)
head(files)
## [1] "wearables/Subject1_rawdata/Masimo//Subject1_MPSnyder_masimohealth_010516_to_011916.csv"
## [2] "wearables/Subject1_rawdata/Masimo//Subject1_MPSnyder_masimohealth_010516_to_012816.csv"
## [3] "wearables/Subject1_rawdata/Masimo//Subject1_MPSnyder_masimohealth_012816_to_021116.csv"
## [4] "wearables/Subject1_rawdata/Masimo//Subject1_MPSnyder_masimohealth_021216_to_021216.csv"
## [5] "wearables/Subject1_rawdata/Masimo//Subject1_MPSnyder_masimohealth_021716_to_021716.csv"
## [6] "wearables/Subject1_rawdata/Masimo//Subject1_MPSnyder_masimohealth_021916_to_030316.csv"
masimo <- lapply(files, function(x) read.csv(x, skip = 1, header = FALSE, as.is = TRUE))

# Create a single dataframe
df <- do.call(rbind, masimo)

# Change column names
names(df) <- c("Session", "Index", "Timestamp", "Day", "Date", "Year", "Time",
                "SpO2", "Pulse Rate")

# Convert to long format and visualize
df[, c("Timestamp", "SpO2", "Pulse Rate")] %>%
  gather(measurement, value, SpO2, "Pulse Rate", -Timestamp) %>%
  ggplot(aes(x = Timestamp, y = value)) +
    geom_line(aes(group = measurement, color = measurement)) +
    labs(title = "Masimo data",
         subtitle = "Heart rate and Blood Oxygen Saturation over time")

Wow there are a lot of dropouts, I should remove outliers and view the mean per hour, or day or something like that. First I’ll have to wrangle the time into an easier format. I wonder if this would be easier if I used one of the other device’s data? I’ll check a few others for easier date measurements.

files <- list.files("wearables/Subject1_rawdata/ihealth-finger/", full.names = TRUE)
head(files)
## [1] "wearables/Subject1_rawdata/ihealth-finger//ihealth.txt"
x <- read.table(files, sep = "\t", header = TRUE, as.is = TRUE)
head(x)
##      Date  Time     TotalTime SpO2   Pulse X
## 1 6/19/15 14:44 6/19/15 14:44  97% 101 bpm  
## 2 6/19/15 14:45 6/19/15 14:45  96%  90 bpm  
## 3 6/19/15 15:56 6/19/15 15:56  96%  98 bpm  
## 4 6/19/15 15:58 6/19/15 15:58  97%  93 bpm  
## 5 6/19/15 16:08 6/19/15 16:08  97%  93 bpm  
## 6 6/19/15 16:20 6/19/15 16:20  98%  95 bpm
dim(x)
## [1] 5151    6
# Reformat columns
x$time <- as.POSIXct(x$TotalTime, tz = "America/Los_Angeles", format = "%m/%d/%y %H:%M")
x$SpO2 <- as.numeric(gsub("%", "", x$SpO2))
x$Pulse <- as.numeric(gsub(" bpm", "", x$Pulse))

x[, c("time", "SpO2", "Pulse")] %>%
  gather(measurement, value, SpO2, Pulse, -time) %>%
  ggplot(aes(x = time, y = value)) +
    geom_line(aes(group = measurement, color = measurement)) +
    labs(title = "ihealth-finger Data",
         subtitle = "Heart rate and Blood Oxygen Saturation over time")

There are about 100 fold fewer observations for ihealth finger than for Masimo. I’ll go back to Masimo and clean it up a bit.

# Format the date
df$datetime <- paste(df$Date, df$Year, df$Time, sep = ", ") %>%
            gsub(" Pacific Standard Time", "", .) %>%
            as.POSIXct(datetime, tz = "America/Los_Angeles", format = " %B %d, %Y, %I:%M:%S %p")  

dfp <- df[, c("datetime", "SpO2", "Pulse Rate")]

# Drop zeros
# sd3p <- which(dfp$Pulse < mean(dfp$Pulse) - sd(dfp$Pulse)*4)
sd3p <- which(dfp$Pulse < 20)
dfp$Pulse[sd3p] <- NA

# sd3s <- which(dfp$SpO2 < mean(dfp$SpO2) - sd(dfp$SpO2)*4)
sd3s <- which(dfp$SpO2 < 20)
dfp$SpO2[sd3s] <- NA

dfp %>%
  gather(measurement, value, SpO2, "Pulse Rate", -datetime, na.rm = TRUE) %>%
  ggplot(aes(x = datetime, y = value)) +
    geom_line(aes(group = measurement, color = measurement)) +
    labs(title = "Masimo data",
         subtitle = "Heart rate and Blood Oxygen Saturation over time")

Maybe the Withings data would be better? Weight vs activity would be easier for people to connect with?

files <- list.files("wearables/Subject1_rawdata/Withings/", full.names = TRUE)
head(files)
## [1] "wearables/Subject1_rawdata/Withings//Withings_activity_and_calories.csv"
## [2] "wearables/Subject1_rawdata/Withings//Withings_weight.xlsx"
withings <- readxl::read_xlsx(files[2])

Oh man that looks so clean. The date is already formatted correctly.

withings %>%
  ggplot(aes(x = Date, y = Weight)) +
    geom_line() +
    theme_bw() +
    labs(title = "Withings data",
         subtitle = "Weight over time")

There are a few extremely low values – I’ll remove those

withings2 <- withings[!(withings$Weight < 120 | withings$Weight > 153), ]

pwithings <- withings2 %>%
  ggplot(aes(x = Date, y = Weight)) +
    geom_line() +
    theme_bw() +
    labs(title = "Withings data",
         subtitle = "Weight over time")
pwithings  

gridExtra::grid.arrange(pmoves, pwithings, ncol = 1)

files <- list.files("wearables/Subject1_rawdata/Withings", full.names = TRUE)
head(files)
## [1] "wearables/Subject1_rawdata/Withings/Withings_activity_and_calories.csv"
## [2] "wearables/Subject1_rawdata/Withings/Withings_weight.xlsx"
withingso <- read_csv(files[1])
## Parsed with column specification:
## cols(
##   Date = col_character(),
##   Steps = col_integer(),
##   `Active calories` = col_integer()
## )
names(withingso)
## [1] "Date"            "Steps"           "Active calories"
withingso$Date <- as.Date(withingso$Date, format = "%m/%d/%y")

withingso %>%
  gather(measurement, value, -Date) %>%
  ggplot(aes(x = Date, y = value)) +
    geom_line(aes(group = measurement, color = measurement)) +
    labs(title = "Withings data",
         subtitle = "Steps and Active Calories over time")
## Warning: Removed 116 rows containing missing values (geom_path).

files <- list.files("wearables/Subject1_rawdata/Scanadu/", full.names = TRUE)
head(files)
## [1] "wearables/Subject1_rawdata/Scanadu//ScanaduData.tsv"
scanadu <- read.table(files, sep = "\t", header = TRUE)
names(scanadu)
## [1] "Data" "Time" "HR"   "SpO2"
# Plot the number of readings per day for the devices reporting HR and SpO2
rbind(data.frame(table(scanadu$Data), device = "Scanadu"),
              data.frame(table(df$Date), device = "Masimo"),
              data.frame(table(x$Date), device = "iHealth-finger")) %>%
  ggplot(aes(Freq)) +
    geom_density(stat = "bin", aes(fill = device), alpha = 0.5) +
    labs(x = "Mesurements per day",
         fill = "Device",
         title = "Number of recorded measurements per day",
         subtitle = "The three devices for HR and SpO2") +
    theme_bw()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

files <- list.files("wearables/Subject1_rawdata/Radiation/", full.names = TRUE)
head(files)
## [1] "wearables/Subject1_rawdata/Radiation//Stanford_Wearables_Radiation_06072016.xlsx"
rad <- readxl::read_xlsx(files)
radlab <- rad[1,2]
rad <- rad[-1,]
names(rad) <- c("datetime", "rad")

rad$rad <- as.numeric(rad$rad)
rad$datetime <- as.POSIXct(rad$datetime, tz = "America/Los_Angeles", format = "%Y-%m-%d %H")

ggplot(rad, aes(x = datetime, y = rad)) +
  geom_line() +
  theme_bw() +
  labs(title = "RadTarge radiation measurements",
       y = radlab,
       x = "Date")

files <- list.files("wearables/Subject1_rawdata/Basis/Basis_Peak/", full.names = TRUE)
head(files)
## [1] "wearables/Subject1_rawdata/Basis/Basis_Peak//bodymetrics_to_03-31-16.csv"
## [2] "wearables/Subject1_rawdata/Basis/Basis_Peak//sleep_to_03-31-16.csv"
body <- read.csv(files[1])
body$date <- gsub("Z$", "", body$date)
head(body$date)
## [1] "2015-05-20 18:50" "2015-05-20 18:51" "2015-05-20 18:52"
## [4] "2015-05-20 18:53" "2015-05-20 18:54" "2015-05-20 18:55"
body$date <- as.POSIXct(body$date, tz = "America/Los_Angeles", format = "%Y-%m-%d %H:%M")

body[,c("date", "heart.rate", "skin.temp")] %>%
  gather(measurement, value, -date) %>%
  ggplot(aes(date, value)) +
    geom_line(aes(group = measurement, color = measurement)) +
    theme_bw() +
    labs(title = "Basis Data",
         subtitle = "Heart Rate and Skin Temperature")
## Warning: Removed 128 rows containing missing values (geom_path).

body %>%
  ggplot(aes(date, gsr)) +
    geom_line() +
    theme_bw() +
    labs(title = "Basis Data",
         subtitle = "Galvanic Skin Response")
## Warning: Removed 64 rows containing missing values (geom_path).

sleep <- read.csv(files[2], as.is = TRUE)
summary(sleep)
##  local_start_time   local_end_time     heart_rate_avg      calories       
##  Length:711         Length:711         Min.   : 51.15   Min.   :   1.964  
##  Class :character   Class :character   1st Qu.: 62.21   1st Qu.:  49.100  
##  Mode  :character   Mode  :character   Median : 67.35   Median : 147.900  
##                                        Mean   : 67.65   Mean   : 162.758  
##                                        3rd Qu.: 72.32   3rd Qu.: 244.700  
##                                        Max.   :106.18   Max.   :1331.600  
##  actual_minutes  light_minutes     deep_minutes     rem_minutes    
##  Min.   : 20.0   Min.   : 20.00   Min.   :  0.00   Min.   :  0.00  
##  1st Qu.: 50.5   1st Qu.: 45.50   1st Qu.:  0.00   1st Qu.:  0.00  
##  Median :150.0   Median : 91.00   Median : 21.00   Median : 33.00  
##  Mean   :163.6   Mean   : 94.93   Mean   : 30.53   Mean   : 37.78  
##  3rd Qu.:248.0   3rd Qu.:131.50   3rd Qu.: 52.00   3rd Qu.: 65.00  
##  Max.   :632.0   Max.   :398.00   Max.   :168.00   Max.   :145.00  
##  interruption_minutes unknown_minutes  interruptions    toss_and_turn  
##  Min.   : 0.000       Min.   : 0.000   Min.   :0.0000   Min.   : 0.00  
##  1st Qu.: 0.000       1st Qu.: 0.000   1st Qu.:0.0000   1st Qu.: 4.00  
##  Median : 0.000       Median : 0.000   Median :0.0000   Median :10.00  
##  Mean   : 1.778       Mean   : 0.121   Mean   :0.3446   Mean   :11.74  
##  3rd Qu.: 0.000       3rd Qu.: 0.000   3rd Qu.:1.0000   3rd Qu.:17.00  
##  Max.   :28.000       Max.   :23.000   Max.   :4.0000   Max.   :63.00  
##  start_timestamp     start_time_iso     start_time_timezone
##  Min.   :1.433e+09   Length:711         Length:711         
##  1st Qu.:1.440e+09   Class :character   Class :character   
##  Median :1.445e+09   Mode  :character   Mode  :character   
##  Mean   :1.446e+09                                         
##  3rd Qu.:1.451e+09                                         
##  Max.   :1.459e+09                                         
##  start_time_offset end_timestamp       end_time_iso      
##  Min.   :-480.0    Min.   :1.433e+09   Length:711        
##  1st Qu.:-480.0    1st Qu.:1.440e+09   Class :character  
##  Median :-420.0    Median :1.445e+09   Mode  :character  
##  Mean   :-344.6    Mean   :1.446e+09                     
##  3rd Qu.:-300.0    3rd Qu.:1.451e+09                     
##  Max.   : 120.0    Max.   :1.459e+09                     
##  end_time_timezone  end_time_offset     X          
##  Length:711         Min.   :-480.0   Mode:logical  
##  Class :character   1st Qu.:-480.0   NA's:711      
##  Mode  :character   Median :-420.0                 
##                     Mean   :-343.8                 
##                     3rd Qu.:-300.0                 
##                     Max.   : 120.0
time <- strsplit(sleep$local_start_time, split = " ")
time <- data.frame(do.call(rbind, time))
time[,2] <- toupper(time[,2])
y <- apply(time, 1, function(x) paste(x, collapse = " "))
sleep$time <- as.POSIXct(y, tz = "America/Los_Angeles", format = "%d %b %Y %H:%M:%S")
names(sleep)
##  [1] "local_start_time"     "local_end_time"       "heart_rate_avg"      
##  [4] "calories"             "actual_minutes"       "light_minutes"       
##  [7] "deep_minutes"         "rem_minutes"          "interruption_minutes"
## [10] "unknown_minutes"      "interruptions"        "toss_and_turn"       
## [13] "start_timestamp"      "start_time_iso"       "start_time_timezone" 
## [16] "start_time_offset"    "end_timestamp"        "end_time_iso"        
## [19] "end_time_timezone"    "end_time_offset"      "X"                   
## [22] "time"
sleep[,c("time", "actual_minutes", "light_minutes", "deep_minutes", "rem_minutes", "interruption_minutes", "unknown_minutes")] %>%
  gather(type, minutes, -time) %>%
  ggplot(aes(x = time, y = minutes)) +
    geom_line(aes(group = type, color = type))
## Warning: Removed 6 rows containing missing values (geom_path).

sleep$day <- format(sleep$time, format = "%Y-%m-%d")
sleep$day <- as.Date(sleep$time, format = "%Y-%m-%d")

sleep[,c("day", "actual_minutes", "light_minutes", "deep_minutes", "rem_minutes", "interruption_minutes", "unknown_minutes")] %>%
  ddply("day", summarize,
        actual = sum(actual_minutes),
        light = sum(light_minutes),
        deep = sum(deep_minutes),
        rem = sum(rem_minutes),
        interruption = sum(interruption_minutes),
        unknown = sum(unknown_minutes)) %>%
  gather(type, minutes, -day) %>%
  ggplot(aes(x = day, y = minutes)) +
    geom_line(aes(group = type, color = type)) +
    theme_bw() +
    labs(title = "BASIS sleep data")
## Warning: Removed 6 rows containing missing values (geom_path).

files <- list.files("wearables/Subject1_rawdata/Basis/MPS_B1/", full.names = TRUE)
head(files)
## [1] "wearables/Subject1_rawdata/Basis/MPS_B1//mpsnyder_20140523.csv"
## [2] "wearables/Subject1_rawdata/Basis/MPS_B1//mpsnyder_20140524.csv"
## [3] "wearables/Subject1_rawdata/Basis/MPS_B1//mpsnyder_20140525.csv"
## [4] "wearables/Subject1_rawdata/Basis/MPS_B1//mpsnyder_20140526.csv"
## [5] "wearables/Subject1_rawdata/Basis/MPS_B1//mpsnyder_20140527.csv"
## [6] "wearables/Subject1_rawdata/Basis/MPS_B1//mpsnyder_20140528.csv"
basis <- lapply(files[1], read.csv)
basis <- do.call(rbind, basis)
# Remove a big pile of NAs
basis <- basis[!is.na(basis$hr),]
time <- gsub("..7:00", "", basis$time)
basis$time <- as.POSIXct(time, tz = "America/Los_Angeles", format = "%Y-%m-%d %H:%M:%S")


ggplot(basis, aes(x = time, y = accel_magnitude)) +
  geom_line() +
  labs(title = "Basis data",
       subtitle = "Acceleration") +
  theme_bw()
## Warning: Removed 99 rows containing missing values (geom_path).