fitbit-wrangling.Rmd

---
title: "fitbit analysis for CSEP project"
author: "Shelby Sturrock"
date: "14/01/2022"
output: 
  html_document:
    keep_md: true
---

```{r setup, include=FALSE, results=FALSE, echo=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(dplyr)
library(tidyr)
library(tidyverse)
library(foreign)
library(readxl)
library(gridExtra)
library(zoo)
knitr::opts_knit$set(root.dir = '~/Desktop/CSEP fitbit')
getwd()
```

This report includes (a) cleaning and exploration of intra-day and day-level Fitbit activity data, (b) creation of a day-level "table of power", including day-level data from Fitbit (vs. minute-level data aggregated at the level of the day), and (c) preliminary analyses of day-level MVPA vs. day-level engagement/clicks. 

# 1. Read in Data

## A. FitBit Data

### Read in and consolidate data

There are 18 csv files from each of four time periods. Read in all data from each of the four time periods (folder). 

```{r, results="hide"}

# sample data
  sample<-list.files("sample_data", pattern="*.csv", full.names=TRUE)
  sample_csv<-lapply(sample, read.csv)
  names(sample_csv)<-paste(gsub("sample_data/|.csv","",sample),".1",sep="")
  list2env(sample_csv,envir=.GlobalEnv)
  rm(sample,sample_csv)

# nov 08 - nov 21
  nov08<-list.files("fitbit_data/nov08-nov21", pattern="*.csv", full.names=TRUE)
  nov08_csv<-lapply(nov08, read.csv)
  names(nov08_csv)<-paste(gsub("fitbit_data/nov08-nov21/|.csv","",nov08),".2",sep="")
  list2env(nov08_csv,envir=.GlobalEnv)
  rm(nov08,nov08_csv)
  
# nov 22 - dec 05
  nov22<-list.files("fitbit_data/nov22-dec05", pattern="*.csv", full.names=TRUE)
  nov22_csv<-lapply(nov22, read.csv)
  names(nov22_csv)<-paste(gsub("fitbit_data/nov22-dec05/|.csv","",nov22),".3",sep="")
  list2env(nov22_csv,envir=.GlobalEnv)
  rm(nov22,nov22_csv)

# dec 06 - dec 12
  dec06<-list.files("fitbit_data/dec06-dec12", pattern="*.csv", full.names=TRUE)
  dec06_csv<-lapply(dec06, read.csv)
  names(dec06_csv)<-paste(gsub("fitbit_data/dec06-dec12/|.csv","",dec06),".4",sep="")
  list2env(dec06_csv,envir=.GlobalEnv)
  rm(dec06,dec06_csv)
    
```

Now consolidate data from four time periods. This will leave us with 18 unique dataframes. 

```{r}

# generate list of unique datasets (there are 18)
  files<-unique(gsub("\\.1|\\.2|\\.3|\\.4","",ls()))
    
# for each unique type of data, rbind data from the three periods together
  for(i in 1:length(files)) {
      
      pattern<-files[i]
      X<-mget(ls(pattern=pattern)) %>% bind_rows()
      distinct(X)
      assign(gsub("_merged","",pattern),X)
      rm(list = ls()[grepl(pattern, ls())])
      rm(X)
      i<-i+1
      
    }
    
# clean environment
  rm(i,pattern,files)
    
```

The 18 data types are as follows: 

```{r}
  ls()
```

## B. Read in demographics and ID

```{r}
  key <- foreign::read.spss("demo_key/Codes & Demographics.sav", to.data.frame=TRUE)
  key$PAC_Email<-tolower(key$PAC_Email)
  key<-key %>% rename(email=PAC_Email,
                      Id=Fitbit_ID)
  key$email<-gsub(" ","",key$email)
```

## C. Read in engagement data

```{r}
  engage<-read_xlsx("engagement_data/17.01.22.24hmg_daily_ps_appengagements.xlsx")
```

# 2. Explore FitBit data

## A. Intra-day data

Intraday datasets including minute-level data: 
  1. heart rate (heartrate_1min)
  2. calories (minuteCaloriesNarrow)
  3. intensity (minuteIntensitiesNarrow)
  4. METs (minuteMETsNarrow)
  5. steps (minuteStepsNarrow)
  6. sleep stage (minuteSleep)

Sleep stage information is also available measured at 30-second increments (30secondSleepStages). 

### Explore intra-day data: 

#### i. Minute-level steps:  

Number of steps per minute
```{r}
  str(minuteStepsNarrow) # steps per minute
```

0 minutes have missing data for steps -- must be set to 0 during non-wear hours
```{r}
  length(minuteStepsNarrow[is.na(minuteStepsNarrow$Steps),]$Id)
```

Plot of step frequency
```{r,echo=FALSE}
  summary(minuteStepsNarrow$Steps)
  ggplot(minuteStepsNarrow,aes(Steps)) + 
         geom_histogram() + 
         theme_classic() + 
         labs(y="Steps per minute")
```


#### ii. Minute-level calories

Number of calories per minute. 
```{r}
  str(minuteCaloriesNarrow) # calories burned each minute
```

Very few records have 0 calories. Fitbit appears to apply a baseline number of calories per minute based on estimated basal metabolic rate. 
```{r}
  length(minuteCaloriesNarrow[minuteCaloriesNarrow$Calories==0,]$Id)
```

Plot of calorie/min frequency
```{r,echo=FALSE}
  summary(minuteCaloriesNarrow$Calories)
  ggplot(minuteCaloriesNarrow,aes(Calories)) + geom_histogram() +   theme_classic() + labs(x="Calories per minute")
```


#### iii. Minute-level activity intensity

Minute-level intensity
```{r}
  str(minuteIntensitiesNarrow) # intensity of activity during each minute -- 0, 1, 2 and 3
```

Plot of intensity frequency
```{r,echo=FALSE}
  summary(minuteIntensitiesNarrow$Intensity)
  ggplot(minuteIntensitiesNarrow,aes(Intensity)) + geom_histogram() + theme_classic() + labs(x="Minute-level intensity")
```

#### iv. Minute-level METs

Minute-level METs
```{r}
  str(minuteMETsNarrow) # MET value measued each minute
```

Plot of strange MET values -- values are much higher than expected
```{r,echo=FALSE}
  summary(minuteMETsNarrow$METs)
  ggplot(minuteMETsNarrow,aes(METs)) + geom_histogram() + theme_classic() + labs(x="minute-level METs")
```


#### v. minute-level heart rate (fewer observations than above)

Minute-level heart rate
```{r}
  str(heartrate_1min) # heart rate measured each minute -- fewer measurements than above
```

Plot of minute-level heart rate
```{r,echo=FALSE}
  summary(heartrate_1min$Value)
  ggplot(heartrate_1min,aes(Value)) + geom_histogram() + theme_classic() + labs(x="minute-level heart rate")
```

### Merge minute-level activity data into one dataframe

First, merge minute-level calories, intensities, METs and steps
```{r}
  minute <- minuteCaloriesNarrow %>% full_join(minuteIntensitiesNarrow, by = c("ActivityMinute","Id")) %>%
                                     full_join(minuteMETsNarrow, by = c("ActivityMinute","Id")) %>%
                                     full_join(minuteStepsNarrow, by = c("ActivityMinute","Id")) 
  rm(minuteCaloriesNarrow,minuteIntensitiesNarrow,minuteMETsNarrow,minuteStepsNarrow)
```

Next, merge minute-level heart rate data to other minute-level activity data
```{r}
  heartrate_1min$ActivityMinute<-heartrate_1min$Time
  heartrate_1min<-heartrate_1min %>% select(-Time)
  minute<-merge(minute,heartrate_1min,by=c("ActivityMinute","Id"),all.x=TRUE)
  rm(heartrate_1min)
```

Generate dateTime and date variables with common naming convention 
```{r}
  minute$dateTime<-lubridate::mdy_hms(minute$ActivityMinute)
  minute$date<-as.Date(minute$dateTime)
  minute<-minute %>% select(-ActivityMinute)
```

Assess time period of data
```{r}
  min(minute$date)
  max(minute$date)
```

### Explore and clean intraday sleep data

Sleep observations can be measured on the half or full minute (one or the other within a given bout of sleep)
```{r}
  minuteSleep$test1<-ifelse(grepl(":30 ",minuteSleep$date),":30",":00")
  table(minuteSleep$test1)
```

However, all other minute-level data are measured on the full minute. 
```{r}
  minute$test1<-ifelse(grepl(":30 ",minute$dateTime),":30",":00")
  table(minute$test1)
  minute<-minute %>% select(-test1)
```

In order to merge with the other minute-level data, I will update the date of all observations to the full minute:
```{r}
# since sleep is either measured at :00 or :30 within a bout of sleep,
# we can update all records to ":00 " without introducing any duplicates
# this is necessary to merge sleep records with the other minute-level data (always measured on full minute)
  minuteSleep$date<-gsub(":30 ",":00 ",minuteSleep$date)
  minuteSleep$dateTime<-lubridate::mdy_hms(minuteSleep$date)
  minuteSleep$date<-as.Date(minuteSleep$dateTime)
  minuteSleep$test2<-ifelse(grepl(":30 ",minuteSleep$date),":30",":00")
  table(minuteSleep$test2)
```

Check the date range of the sleep data
```{r}
  summary(minuteSleep$date)
```

Filter to the same date range as the other minute-level data. 
```{r}
  minuteSleep<-minuteSleep[minuteSleep$date<"2021-12-13",]
```

Another issue: There are many duplicate sleep records (multiple rows corresponding to the same participant-minute):
```{r}
  test<-minuteSleep %>% select(Id,dateTime) %>% group_by(Id,dateTime) %>% mutate(length=length(Id))
  table(test$length) # this is unrelated to updating the date/time from :30 to :00 above
  # length(unique(test[test$length==2,]$Id)) 
    # almost every participant has duplicate records for sleep
```

Importantly, all but 6 observations (3 minutes for 1 participants) are exact duplicates (same Id, date AND value)
```{r,warnings=FALSE}
  test3<-minuteSleep %>% select(Id,dateTime,value) %>% group_by(Id,dateTime) %>% 
                                                       mutate(length=length(Id)) %>%
                                                       mutate(unique=length(unique(value)))
# however, only 6 records (3 minutes for 1 participant) have more than one unique "value" within the same minute
# the rest of the duplicate rows are exact duplicates (same "value" (sleep stage) across both observations)
  table(test3[test3$length==2,]$unique) 
  rm(test,test3)
```

I will remove exact duplicates (across all columns in the data frame):
```{r}
# create de-duplicated version of the minuteSleep dataframe
# need to decide how to handle the 3 minutes with two unique "values"
  minuteSleep<-minuteSleep %>% select(-test1,-test2)
  minuteSleep<-distinct(minuteSleep) 
```

There are still over 800 duplicate records (same participant-minute):
```{r}
  test4<-minuteSleep %>% select(Id,dateTime,value) %>% group_by(Id,dateTime) %>% 
                                                       mutate(length=length(Id)) %>%
                                                       mutate(unique=length(unique(value)))
  table(test4$length)
```

It appears as though all 800+ duplicate rows are within the same participant -- the time and value are the same, but the logId differs between observations: 
```{r}
  minuteSleep<-distinct(minuteSleep, across(-"logId"),.keep_all = TRUE) 
  test4<-minuteSleep %>% select(Id,dateTime,value) %>% group_by(Id,dateTime) %>% 
                                                       mutate(length=length(Id)) %>%
                                                       mutate(unique=length(unique(value)))
  length(unique(test4[test4$length==2,]$Id))
  table(test4$length)
  duplicates<-test4[test4$length==2,]
```

The remaining 6 duplicates represent 3 minutes from one participant where the value also differs between the two observations. I will recode value to NA and then de-duplicate again: 
```{r}
  minuteSleep$value<-ifelse(minuteSleep$Id %in% duplicates$Id & 
                            minuteSleep$dateTime %in% duplicates$dateTime,NA,minuteSleep$value)
  minuteSleep<-distinct(minuteSleep,across(-"logId"),.keep_all = TRUE) 
  rm(duplicates)
```


### Merge minute-level sleep and activity data

Merge minute-level sleep and activity data
```{r}
  minute <- minute %>% full_join(minuteSleep, by = c("dateTime","Id","date"))
  rm(minuteSleep)
```

Rename columns in minute-level data (so it is clear which columns are minute- vs. day-level)
```{r}
  minute<-minute %>% rename(minuteCalories=Calories,
                            minuteHeartRate=Value,
                            minuteIntensity=Intensity,
                            minuteMETs=METs,
                            minuteSteps=Steps,
                            minuteSleepStage=value,
                            minuteSleepLogId=logId) %>%
                     relocate(dateTime,.after=Id) %>%
                     relocate(date,.after=dateTime) %>%
                     arrange(Id,dateTime)
```

### Identify non-wear time

Will set time as non-wear if every value (except calories) is 0 (and intensity is 1) for 10 minutes in a row. 
```{r}
  minute$nonWearFlag<-ifelse(minute$minuteSteps==0 & minute$minuteIntensity==0 & is.na(minute$minuteHeartRate),1,0)
  minute<-minute %>% arrange(Id,dateTime) %>% group_by(Id) %>% mutate(nonWearMean = zoo::rollapply(nonWearFlag, width=10, 
                                                                                                   align="right", mean, fill=NA))
  minute$nonWearMean<-ifelse(is.na(minute$nonWearMean),minute$nonWearFlag,minute$nonWearMean) 
    # first 10 observations will be non-wear time for all participants
```

### Reformat minute-level intensity level, considering non-wear time

Reformat minute-level intensity 
```{r}
  minute$minuteIntensity1<-ifelse(minute$minuteIntensity==0,"minuteSedentary",
                               ifelse(minute$minuteIntensity==1,"minuteLightlyActive",
                                      ifelse(minute$minuteIntensity==2,"minuteFairlyActive",
                                             ifelse(minute$minuteIntensity==3,"minuteVeryActive",NA))))
  minute$minuteIntensity1<-ifelse(minute$nonWearMean==1,"non-wear",minute$minuteIntensity1)
# remove rows where nonWearMean is NA (the first 9 rows per participant)
  minute<-minute[!is.na(minute$nonWearMean),]
  summary(as.factor(minute$minuteIntensity1))
  minute$minuteIntensityValue<-1
  minute<-pivot_wider(minute,names_from="minuteIntensity1",values_from="minuteIntensityValue")
```


## B. Day-level data

### Explore day-level data

#### i. Day-level calories

```{r}
  str(dailyCalories)
```

Distribution: Daily calories
```{r,echo=FALSE}
  summary(dailyCalories$Calories)
  ggplot(dailyCalories,aes(Calories)) + geom_histogram() + theme_classic() + labs(x="Daily Calories")
```

#### ii. Day-level intensities

```{r}
  str(dailyIntensities)
```

#### iii. Day-level steps

```{r}
  str(dailySteps)
```

Distribution: daily steps
```{r,echo=FALSE}
  summary(dailySteps$StepTotal)
  ggplot(dailySteps,aes(StepTotal)) + geom_histogram() + theme_classic() + labs(x="Daily Steps")
```

#### iv. Day-level sleep 

```{r}
  str(sleepDay)
```

Distribution: daily total time in bed
```{r,echo=FALSE}
  summary(sleepDay$TotalTimeInBed)
  ggplot(sleepDay,aes(TotalTimeInBed)) + geom_histogram() + theme_classic() + labs(x="Daily Time in Bed")
```

#### v. Day-level sleep stages

```{r}
  str(sleepStagesDay) # same variables as above (Id, SleepDay, TotalSleepRecords, TotalMinutesAsleep, TotalTimeInBed) plus total minutes awake and spent in each sleep stage (light/deep/REM)
  # will use sleepStagesDay instead of sleepDay because it contains all of the same and more information
  rm(sleepDay)
```

### Merge day-level activity data

```{r}
  daily<-dailyCalories %>% full_join(dailyIntensities, by = c("ActivityDay","Id")) %>%
                           full_join(dailySteps, by = c("ActivityDay","Id")) 
  rm(dailyCalories,dailyIntensities,dailySteps)
```

Generate date variable (named to align with date column in minute-level datasets)
```{r}
  str(daily$ActivityDay)
  daily$date<-as.Date(daily$ActivityDay,"%m/%d/%Y")
```

### Merge day-level activity data and day-level sleep data 

Generate date in the day-level sleep data 
```{r}
  sleepStagesDay$date<-as.Date(lubridate::mdy_hms(sleepStagesDay$SleepDay))
```

Merge day-level activity data and day-level sleep data
```{r}
  daily<-daily %>% full_join(sleepStagesDay, by = c("date","Id"))
  daily<-daily %>% select(-ActivityDay,-SleepDay)
  rm(sleepStagesDay)
```

Rename columns in day-level data (so it is clear which columns are minute- vs. day-level)
```{r}
  daily<-daily %>% rename(TotalCalories=Calories) %>%
                   relocate(date,.after=Id) 
```

Create intervention flag in day-level data
```{r}
  daily$intervention<-ifelse(daily$date < "2021-10-25","pre",
                      ifelse(daily$date >= "2021-12-05","post","intervention"))
  daily$intervention<-factor(daily$intervention,levels=c("pre","intervention","post"))
  table(daily$intervention)
```


### Merge minute- and day-level dataframes

```{r}
  data<-merge(daily,minute,by=c("Id","date"))
```

### C. Explore remaining dataframes

#### Battery (and syncEvents) data

Clean columns in battery data
```{r}
  battery$dateTime<-as.character(lubridate::mdy_hms(battery$DateTime))
  substr(battery$dateTime,18,19)<-"00"
  battery$dateTime<-as.POSIXct(battery$dateTime,tz="UTC")
  battery$lastSync<-lubridate::mdy_hms(battery$LastSync)
  battery<-battery %>% select(Id,dateTime,lastSync,DeviceName,BatteryLevel)
  unique(battery$BatteryLevel)
  length(battery[battery$BatteryLevel=="Empty",]$BatteryLevel)
```

Sometimes there are multiple sync's per minute, but they always have the same battery level
```{r}
# sometimes there are multiple syncs within the same minute
  test<-battery %>% group_by(Id,dateTime) %>% mutate(length=length(Id),
                                                     unique=length(unique(Id)))
  table(test$unique)
  rm(test)
```

Take list of distinct battery measurements (one per minute maximum)
```{r}
  battery<-distinct(battery,across(-"lastSync"))
```


syncEvents dataframe does not contain any additional info above and beyond battery data
```{r}
  names(syncEvents)
# will use battery data instead (has the same and more information as syncEvents)
  rm(syncEvents)
```

Merge battery data with activity and sleep data
```{r}
  data<-merge(data,battery,by=c("Id","dateTime"),all=TRUE)
  rm(battery)
```

#### sleepLogInfo

sleepLogInfo contains other information re: sleep, including the time they went to sleep, how efficient the sleep was, minutes to fall asleep, time in bed, awakeCount, AwakeDuration, RestlessCount, Restless Duration
```{r}
# may not be necessary for this study? will omit for now and revisit
  str(sleepLogInfo)
  rm(sleepLogInfo)
```

#### sleepStageLogInfo

This seems to contain more information relating to each individual sleep bout of sleep, including the time spent in each intensity
```{r}
# however, primary dataframe (data) already contains much of the high-level summary statistics re: sleep
  str(sleepStageLogInfo)
  rm(sleepStageLogInfo)
```

#### 30 second sleep stages

We already have minute-level sleep stages (incorporated into the dataframe); however, this has data regarding the sleep level, as well as any interruptions (short wakes) and what the sleep stage was during that short wake
```{r}
  str(`30secondSleepStages`)
  rm(`30secondSleepStages`)
```

#### heartRateZones

Long dataset with day-level heart rate zone information with four rows per person, one per each zone (out of range, fat burn, cardio, peak), and additional columns defining the minimum and maximum HR values for each zone. 
```{r}
# structured this way to have a min and max column defining the heart rate range corresponding to the label (range will differ between participants)
  str(heartRateZones)
```

Re-structure heartRateZones before merging
```{r}

# create a "range column" for each (can always parse into min/max columns if necessary?)
  heartRateZones$range<-paste(heartRateZones$Min,"-",heartRateZones$Max,sep="")

# pivot wider
  HRrange<-heartRateZones %>% select(Id,ActivityDay,ZoneName,range)
  HRrange<-pivot_wider(data=HRrange,id_cols=c("Id","ActivityDay"),names_from="ZoneName",values_from="range")
  HRrange<-HRrange %>% rename(range_outOfRange=`Out of Range`,
                              range_fatBurn=`Fat Burn`,
                              range_cardio=`Cardio`,
                              range_peak=`Peak`)
  
# generate minutes in each HR zone
  HRzoneMinutes<-heartRateZones %>% select(Id,ActivityDay,ZoneName,Minutes)
  HRzoneMinutes<-pivot_wider(data=HRzoneMinutes,id_cols=c("Id","ActivityDay"),names_from="ZoneName",values_from="Minutes")
  HRzoneMinutes<-HRzoneMinutes %>% rename(total_outOfRange=`Out of Range`,
                                          total_fatBurn=`Fat Burn`,
                                          total_cardio=`Cardio`,
                                          total_peak=`Peak`)
  
# merge back together
  HRzones<-HRrange %>% full_join(HRzoneMinutes,by=c("Id","ActivityDay"))
  
# remove dataframes I don't need
  rm(HRzoneMinutes,heartRateZones,HRrange)
```

Generate date column in day-level heartrate zone data
```{r}
  HRzones$date<-as.Date(HRzones$ActivityDay,"%m/%d/%Y")
  HRzones<-HRzones %>% select(-ActivityDay)
```

Merge with data
```{r}
  data<-merge(data,HRzones,by=c("Id","date"))
  rm(HRzones)
```

#### activitylogs

```{r}
# similar to sleepLogInfo and sleepStageLogInfo, but for workouts
# one row per workout
  str(activitylogs)
```

Can merge workout data at the start time?
```{r}
  activitylogs$dateTime<-paste(activitylogs$Date,activitylogs$StartTime,sep=" ")
  substr(activitylogs$dateTime,18,19)<-"00"
  activitylogs$dateTime<-lubridate::mdy_hms(activitylogs$dateTime)
```

Select and rename columns 
```{r}
  activitylogs$date<-as.Date(activitylogs$dateTime)
  activitylogs<-activitylogs %>% select(-Date,-StartTime,-ActivityType) %>%
                                 rename(activity_duration=Duration,
                                        activity=Activity,
                                        activity_logType=LogType,
                                        activity_steps=Steps,
                                        activity_distance=Distance,
                                        activity_elevationGain=ElevationGain,
                                        activity_calories=Calories,
                                        activity_sedentaryMins=SedentaryMinutes,
                                        activity_lightlyActiveMins=LightlyActiveMinutes,
                                        activity_fairlyActiveMins=FairlyActiveMinutes,
                                        activity_veryActiveMins=VeryActiveMinutes,
                                        activity_avgHeartRate=AverageHeartRate,
                                        activity_outOfRangeMins=OutOfRangeHeartRateMinutes,
                                        activity_fatBurnMins=FatBurnHeartRateMinutes,
                                        activity_cardioMins=CardioHeartRateMinutes,
                                        activity_peakMins=PeakHeartRateMinutes) %>%
                                  relocate(date,.before=activity_duration) %>%
                                  relocate(dateTime,.after=date)
```


Generate count of workouts per day, unique types of workouts per day, total duration of workouts
```{r}
  activitylogs$activity_duration<-activitylogs$activity_sedentaryMins+
                                  activitylogs$activity_lightlyActiveMins+
                                  activitylogs$activity_fairlyActiveMins+
                                  activitylogs$activity_veryActiveMins
  activitylogs<-activitylogs %>% group_by(Id,date) %>% mutate(day_activity_count=length(Id),
                                                              day_activity_duration=sum(activity_duration),
                                                              day_activity_unique=length(unique(activity))) %>%
                                                       relocate(day_activity_count,.after=dateTime) %>%
                                                       relocate(day_activity_duration,.after=day_activity_count) %>%
                                                       relocate(day_activity_unique,.after=day_activity_duration) 
```


Merge with data
```{r}
  data<-merge(data,activitylogs[,!(colnames(activitylogs) %in% c("date"))],by=c("Id","dateTime"),all.x=TRUE)
  rm(activitylogs)
```

## D. Compare minute-level and day-level data

Aggregate at the level of the day
```{r}
  assess<-data %>% group_by(Id,date) %>% summarise(Id=identity(Id),
                                                   date=identity(date),
                                                   TotalCalories=identity(TotalCalories),
                                                   minuteCalories=as.numeric(format(trunc(sum(minuteCalories),0),nsmall=0)),
                                                   nonWearMinutes=sum(`non-wear`,na.rm=TRUE),
                                                   SedentaryMinutes=identity(SedentaryMinutes),
                                                   minuteSedentary=sum(minuteSedentary,na.rm=TRUE),
                                                   LightlyActiveMinutes=identity(LightlyActiveMinutes),
                                                   minuteLightlyActive=sum(minuteLightlyActive,na.rm=TRUE),
                                                   FairlyActiveMinutes=identity(FairlyActiveMinutes),
                                                   minuteFairlyActive=sum(minuteFairlyActive,na.rm=TRUE),
                                                   VeryActiveMinutes=identity(VeryActiveMinutes),
                                                   minuteVeryActive=sum(minuteVeryActive,na.rm=TRUE),
                                                   nonWearMinutes=sum(nonWearMinutes,na.rm=TRUE),
                                                   TotalSleps=identity(StepTotal),
                                                   minuteSteps=sum(minuteSteps),
                                                   TotalTimeInBed=identity(TotalTimeInBed),
                                                   TotalTimeAsleep=identity(TotalMinutesAsleep),
                                                   TotalTimeAwake=identity(TotalTimeAwake)) %>% distinct()
# merge Id, date and nonWearMinutes (per Id and date) with daily data
  daily<-merge(daily,assess[,colnames(assess) %in% c("Id","date","nonWearMinutes")],by=c("Id","date"))
```

Explore distribution of nonWearMinutes
```{r}
  nonWearDist<-ggplot(daily,aes(x=nonWearMinutes)) + 
                          geom_bar(fill="#619CFF") + 
                          theme_classic() +
                          labs(x="Daily minutes of non-wear time",
                               y="Frequency",
                               title="Non-wear time") +
                          theme(axis.title.x = element_text(face="bold"),
                                axis.title.y = element_text(face="bold"),
                                plot.title = element_text(face="bold",
                                                          hjust=0.5))
  nonWearDist
  
  nonWearOverTime<-ggplot(daily,aes(x=date,y=nonWearMinutes)) + 
                          geom_jitter(color="#619CFF",alpha=0.1) + 
                          theme_classic() +
                          labs(x="Date",
                               y="Daily non-wear minutes",
                               title="Non-wear time") +
                          theme(axis.title.x = element_text(face="bold"),
                                axis.title.y = element_text(face="bold"),
                                plot.title = element_text(face="bold",
                                                          hjust=0.5))
  nonWearOverTime
```


How many minutes of data are accounted for accross the various day-level intensity categories?
Calculate as the sum of SedentaryMinutes, LightlyActiveMinutes, FairlyActiveMinutes, VeryActivitMinutes, TotalTimeInBed, nonWearMinutes
```{r}
  assess$sum_minutes<-assess$SedentaryMinutes+
                      assess$LightlyActiveMinutes+
                      assess$FairlyActiveMinutes+
                      assess$VeryActiveMinutes+
                      assess$TotalTimeInBed
  assess<-distinct(assess)
  summary(assess$sum_minutes)
  length(assess[assess$sum_minutes!=1440,]$Id) # 362 people-days with either more or less than 1440 minutes/day of day-level data (under 10% of days)
```

Create a dataframe that calculates the difference between totals from the day-level data vs. totals aggregated up from the minute-level data
```{r}
# compare
  compare<-assess %>% mutate(diffCalories=TotalCalories-minuteCalories,
                             diffSed=SedentaryMinutes-minuteSedentary, # day-level minute minute-level 
                             diffLight=LightlyActiveMinutes-minuteLightlyActive,
                             diffFairly=FairlyActiveMinutes-minuteFairlyActive,
                             diffVery=VeryActiveMinutes-minuteVeryActive,
                             diffSteps=TotalSleps-minuteSteps,
                             nonWearMinutes=identity(nonWearMinutes))
```

Difference in calories: 
```{r}
  compare<-distinct(compare)
  summary(compare$diffCalories) 
      # agree for ~2300 participant-days; 
      # disagrees for just over 100 participant-days, most by a relatively small amount, 12 by more than 100 calories 
  ggplot(compare,aes(diffCalories)) + geom_histogram() + theme_classic() + labs(x="Difference in daily calories")
```
Difference in minutes of sedentary time: 
```{r}
  # diffSed: issues with day-level measure of sedentary time -- appears to take another category into account?
  ggplot(compare,aes(diffSed)) + geom_histogram() + theme_classic() + labs(x="Difference in daily sedentary time")
  ggplot(compare,aes(diffSed+TotalTimeInBed)) + geom_histogram() + theme_classic() + labs(x="Difference in daily sedentary time")
```

Difference in minutes of light activity: 
```{r}
  # diffSed: issues with day-level measure of sedentary time -- appears to take another category into account?
  summary(compare$diffLight) # disagrees for 4 people-days, by no more than 12 minutes
  ggplot(compare,aes(diffLight)) + geom_histogram() + theme_classic() + labs(x="Difference in daily light activity")
```

Difference in minutes of fairly active: 
```{r}
  # diffSed: issues with day-level measure of sedentary time -- appears to take another category into account?
  summary(compare$diffFairly) # disagrees for 3 people-days, by no more than 7 minutes
  str(compare$diffFairly)
  ggplot(compare,aes(diffFairly)) + geom_histogram() + theme_classic() + labs(x="Difference in fairly active minutes")
```

Difference in minutes of very active: 
```{r}
  summary(compare$diffVery) # disagrees for 2 people-days, by no more than 5 minutes
  ggplot(compare,aes(diffVery)) + geom_histogram() + theme_classic() + labs(x="Difference in very active minutes")
```

Difference in steps: 
```{r}
  summary(compare$diffSteps) # disagrees for ~150 people-days, 12 people/days differ by more than 100 steps (9 by >1000 steps)
  ggplot(compare,aes(diffSteps)) + geom_histogram() + theme_classic() + labs(x="Difference in daily steps")
```

```{r}
  rm(assess,compare,data)
```


# 3. Add Demographic and Daily Engagement Data

## Merge Demographic and Engagement Data using email

Check for IDs that are in engagement and not key
```{r}
# check for IDs that are in engagement and not key
  unique(engage[!(engage$email %in% key$email),]$email)
```

Check for IDs that are in key and not engagement
```{r}
# check for IDs that are in the key but not the engagement data
  unique(key[!(key$email %in% engage$email),]$email)
```

Note: There are three participants missing the letter location identifier in their engagement data email address (missing an A specifically). The email address in the key dataframe contains the A, so the records did not merge. All three participants look to be from Queens, and should have an A in their email. 

I will manually correct these email addresses and then merge the data again
```{r}
  engage$email<-ifelse(engage$email=="pacappstudy+08@gmail.com","pacappstudy+a08@gmail.com",
                       ifelse(engage$email=="pacappstudy+03@gmail.com","pacappstudy+a03@gmail.com",
                              ifelse(engage$email=="pacappstudy+28@gmail.com","pacappstudy+a28@gmail.com",engage$email)))
```

Try merging again with corrected email addresses
```{r}
  engage<-merge(engage,key,by=c("email"))
  rm(key)
```

Note that all records from the engagement dataset have been retained after this correction. 

Reformat date so it aligns with the daily fitbit data (as.Date vs. as.POSIXct)
```{r}
  engage$date<-as.Date(engage$date)
```

Trim whitespace from Id columns in engage data
```{r}
  engage$Id<-gsub(" ","",engage$Id)
```

There are some participants who have a row indicating "no_engagement" (during the study) but one or more records indicating that they had engaged with a specific article on a specific date: 
```{r}
  any_engage<-engage[,colnames(engage) %in% c("Id","any_engagement")]
  any_engage<-distinct(any_engage)
  any_engage<-any_engage %>% group_by(Id) %>% mutate(n_rows=length(Id))
  length(unique(engage[engage$Id %in% any_engage[any_engage$n_rows==2,]$Id,]$Id))
```

Need to check with Geri to see how this variable was calculated. For now, I will generate a new version (ever_engage) that = YES if there is at least one row where date and title are not missing, NO otherwise
```{r}
  engage<-engage %>% group_by(Id) %>% mutate(ever_engage=max(click))
  table(engage$any_engagement,engage$ever_engage)
  any_engage<-engage[,colnames(engage) %in% c("Id","ever_engage")]
  any_engage<-distinct(any_engage)
```

Remove rows where ever_engage = 1 but any_engagement = no_engagement
```{r}
  engage<-engage[!(engage$any_engagement=="no_engagement" & engage$ever_engage==1),]
```

## Merge daily fitbit data with engagement data

Split time-varying and non-time-varying factors
```{r}
  names(engage)
  engage_participant<-engage %>% select(email,Id,ever_engage,Age,Year_HS_Graduation,Gender,Ethnicity,University,Undergrad_Program,
                                        Current_Living_Situation,PACapp_Use,Tracking_Movement_Behaviour,
                                        Questionnaire_ID_1,Questionnaire_ID_2)
  engage_time_varying<-engage %>%
                       select(Id,date,title,intervention_content,featured_content,pubdate,weekofpublish,type,click,engage)
```

Merge daily data and variables that are not time/date-dependent
```{r}
  daily<-distinct(merge(daily[daily$Id!="Study Coordinator QU",],
                        engage_participant,by=c("Id")))
```

Merge time-varying engagement variables
```{r}
  daily<-distinct(merge(daily,
                        engage_time_varying,by=c("Id","date"),all.x=TRUE))
```

<!-- ## Exploratory: Compare FitBit intensity categorization with alternate approaches.  -->

<!-- ### Merge age from enagement data with minute-level heart rate data to calcualte various intensity metrics -->

<!-- Merge minute-level HR with ages -->
```{r}
  # minSimple<-minute[minute$Id!="Study Coordinator QU",
  #                              colnames(minute) %in%
  #                              c("Id","dateTime","date","minuteHeartRate","minuteIntensity","minuteSteps")]
  # minHR<-merge(minSimple,engage_participant[,colnames(engage_participant) %in% c("Id","Age")],by=c("Id"))
  # minHR<-distinct(minHR)
  # rm(minSimple)
```

<!-- ### Generate alternative intensity categories using demographics -->

<!-- #### minuteIntensity1: fitbit categories -->

<!-- This is a four-level variable. I am assuming that 0 = sedentary, 1 = lightly active, 2 = fairly active, 3 = very active. I am also assuming that fairly active is approximately moderate intensity, very active is approximately vigorous intensity.  -->
```{r}
  # minHR$minuteIntensity1<-ifelse(minHR$minuteIntensity %in% c("2","3"),"MVPA","non-MVPA")
  # table(minHR$minuteIntensity1)
```

<!-- #### minuteIntensity2: 100 steps/minute -->

<!-- based on research from Catrine Tudor-Locke, 100 steps/minute and 130 steps/minute can be used as thresholds for moderate and vigorous physical activity (https://ijbnpa.biomedcentral.com/articles/10.1186/s12966-019-0769-6)  -->
```{r}
  # minHR$minuteIntensity2<-ifelse(minHR$minuteSteps>100,"MVPA","non-MVPA")
  # table(minHR$minuteIntensity2) # fewer MVPA minutes using this vs. the fitbit categories
```

<!-- #### minuteIntensity3: Karvonen formula -->

<!-- Karvonen formula estimates: calculate heart rate resevere by taking maximum - resting heart rate. We then calculate the threshold HR (approximately corresponding to MVPA) by taking resting heart rate + heart rate reserve (maximum minus resting), multiplied by 0.4 (40%). Max is calculated by taking 220 minus age, and resting is estimated as the median HR among records under 80 bpm for that particular participant (https://bcmj.org/articles/science-exercise-prescription-martti-karvonen-and-his-contributions) -->
```{r}
# heart rate reserve = max - resting (how to determine resting heart rate?) (mean/median of values less than 80 as resting heart rate)
# resting + (max-resting)*0.4
# 60 + (140)*0.4 = 116
  # estRHR<-minHR[!is.na(minHR$minuteHeartRate) & minHR$minuteHeartRate<80,]
  # estRHR<-estRHR %>% select(Id,minuteHeartRate)
  # estRHR<-estRHR %>% group_by(Id) %>% mutate(meanRHR=format(round(mean(minuteHeartRate),0),nsmall=0),
  #                                            medianRHR=median(minuteHeartRate)) %>% 
  #                                     select(Id,meanRHR,medianRHR)
  # estRHR<-distinct(estRHR)
  # minHR<-merge(minHR,estRHR,by="Id")
  # minHR$thresholdHR<-(minHR$medianRHR+(((220-minHR$Age)-minHR$medianRHR)*0.4))
  # minHR$minuteIntensity3<-ifelse(minHR$minuteHeartRate >= minHR$thresholdHR,"MVPA","non-MVPA")
  # table(minHR$minuteIntensity3) # fewer using this method
```

<!-- #### comparison of various intensity categorizations -->

```{r}
# fitbit vs. 100 steps/min
  # test<-table(minHR$minuteIntensity1,minHR$minuteIntensity2)
  # chisq.test(test)
```

<!-- There is not a lot of agreement here... more minutes categorized as 2/3 intensity were flagged as non-MVPA than MVPA.  -->
```{r}
# fitbit vs. hr one steps/min
  # test2<-table(minHR$minuteIntensity1,minHR$minuteIntensity3)
  # chisq.test(test2)
```

<!-- Next steps: Replace chi-square with non-parametric ranked statistic (to revisit). For now, I will proceed using the day-level intensity data provided by Fitbit (minutes spent in each intensity category)  -->

<!-- Remove engagement data frames -->
```{r}
  # rm(engage,engage_participant,engage_time_varying)
```


# 4. Create click and engagement metrics (including count of clicks and engagements per person-day, person-period and person, and mean daily clicks and engagements per person-period and person-day)

Recode ever_engage for figures
```{r}
  daily$ever_engage<-ifelse(daily$ever_engage==1,"Clicked or engaged at least once","Never engaged or clicked")
```

Add week start date and day of week for plotting and analysis
```{r}
  daily$weekNum<-format(daily$date,"%W")
  daily<-daily %>% group_by(weekNum) %>% mutate(weekStart=min(date))
  daily$weekday<-factor(weekdays(daily$date),
                        levels=c("Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"))
  daily$mon_wed<-ifelse(daily$weekday %in% c("Monday","Wednesday"),"yes","no")
```

Calculate daily minutes of MVPA by adding fairly active and very active minutes
```{r}
  daily$totalMVPA<-daily$FairlyActiveMinutes+daily$VeryActiveMinutes
```

Remove participant-days where totalMVPA is na
```{r}
  daily<-daily[!is.na(daily$totalMVPA),]
```

Generate new clickInt and engageInt variables that are binary flags for clicks and engagement with intervention content specifically
```{r}
# set click and engage to 0 if click is NA
  daily$click<-ifelse(is.na(daily$click),0,daily$click)
  daily$engage<-ifelse(is.na(daily$engage),0,daily$engage)
  
# create new dummy variables for clicks/engagement with intervention content
  daily$clickInt<-ifelse(daily$click==1 & daily$intervention_content=="yes" & !is.na(daily$intervention_content),1,0)
  summary(as.factor(daily$clickInt)) # there are two clicks where the title and intervention_content are unknown but click == 1
  daily$engageInt<-ifelse(daily$engage==1 & daily$intervention_content=="yes",1,0)
  summary(as.factor(daily$engageInt))
```

#### Summarize clicks and engagement: Counts of clicks, engagement over different time periods (overall, by intervention period, by day)
(1) per participant-day (dayCountAny, dayEngageAny, dayEngageInt, dayCountInt): group by date and Id, sum the binary/dummy variable for clicks, engagements, clicks on intervention content, and engagement on intervention content to get counts per day
(2) per participant-period (periodCountAny, periodEngageAny, periodEngageInt, periodCountInt): group by Id and intervention period, sum the binary/dummy variable for clicks, engagements, clicks on intervention content, and engagement on intervention content to get counts per intervention period
(3) per participant (personCountAny, personEngageAny, personEngageInt, personCountInt): group by Id, sum the binary/dummy variable for clicks, engagements, clicks on intervention content, and engagement on intervention content to get counts per intervention period
```{r}
# set click and engage to 0 if na
  daily<-daily %>% group_by(Id,date) %>% mutate(dayCountEngageAny=sum(engage),
                                                dayCountClickAny=sum(click),
                                                dayCountEngageInt=sum(engageInt),
                                                dayCountClickInt=sum(clickInt)) #%>%
                 # ungroup() %>% group_by(Id,intervention) %>% mutate(periodCountClickAny=sum(click),
                 #                                                    periodCountEngageAny=sum(engage),
                 #                                                    periodCountClickInt=sum(clickInt),
                 #                                                    periodCountEngageInt=sum(engageInt)) %>% 
                 # ungroup() %>% group_by(Id) %>% mutate(personCountClickAny=sum(click),
                 #                                       personCountEngageAny=sum(engage),
                 #                                       personCountClickInt=sum(clickInt),
                 #                                       personCountEngageInt=sum(engageInt))
```

#### Summarize clicks and engagement: Mean daily enagement over different periods (overall, by intervention period)
<!-- Calculate average daily clicks and engagements per participant, and per participant and period (pre, intervention, post) -->
```{r}
  # daily<-daily %>% group_by(Id) %>% mutate(personMeanDailyClickAny=mean(dayCountClickAny),
  #                                          personMeanDailyEngageAny=mean(dayCountEngageAny),
  #                                          personMeanDailyClickInt=mean(dayCountClickInt),
  #                                          personMeanDailyEngageInt=mean(dayCountEngageInt),
  #                                          personMeanDailyMVPA=mean(totalMVPA),
  #                                          personMedianDailyMVPA=median(totalMVPA)) %>%
  #                  group_by(Id,intervention) %>% mutate(periodMeanDailyClickAny=mean(dayCountClickAny),
  #                                                       periodMeanDailyEngageAny=mean(dayCountEngageAny),
  #                                                       periodMeanDailyClickInt=mean(dayCountClickInt),
  #                                                       periodMeanDailyEngageInt=mean(dayCountEngageInt),
  #                                                       periodMeanDailyMVPA=mean(totalMVPA),
  #                                                       periodMedianDailyMVPA=median(totalMVPA))
```


# 5. Clean and export TOP

## Remove participant-day with < 10 hours of data per day (across the 5 movement categories) or > 24 hours of data
```{r}
  daily2<-daily
  daily$sumMinutes<-daily$SedentaryMinutes+
                    daily$LightlyActiveMinutes+
                    daily$totalMVPA+
                    daily$TotalTimeInBed
  hist(daily$sumMinutes)
  daily<-daily[daily$sumMinutes<=1440 & daily$sumMinutes >=600,]
```

## Remove participant-days where TotalTimeAsleep == 0 and SedentaryMinutes == 1440 (or greater?)

First, explore participant-days with 0 time in bed 
```{r}
  length(daily[daily$TotalTimeInBed==0,]$Id) # 486 participant-days with 0 time in bed
    summary(daily[daily$TotalTimeInBed==0,]$SedentaryMinutes) # almost half had 0 1440 sedentary minutes, most 3/4 had at least 1134... could the balance be time when the device was dead? so 1440 = device was off, <1440 = off all day but dead part of the day 
    summary(daily[daily$TotalTimeInBed==0,]$StepTotal)
    summary(daily[daily$TotalTimeInBed==0,]$TotalMinutesAsleep) # all participant-days also have 0 minutes asleep
```

Next, explore participant-days with 0 minutes asleep 
```{r}
  length(daily[daily$TotalMinutesAsleep==0,]$Id) # 565 particiapnt-days with 0 time asleep -- same as the number with 0 time in bed
    summary(daily[daily$TotalMinutesAsleep==0,]$TotalTimeInBed)
```

Also explore participant-days with > or = 1440 minutes of sedentary minutes
```{r}
  length(daily[daily$SedentaryMinutes>=1440,]$Id) # 242 participant-days with 1440 or more sedentary minutes
  test<-daily[daily$SedentaryMinutes>=1440,]
    unique(daily[daily$SedentaryMinutes>=1440,]$TotalTimeInBed) # all 242 participant-days with 1440 or more sedentary minutes also had 0 minutes in bed
```

...and participant-days with 0 steps
```{r}
  length(daily[daily$StepTotal==0,]$Id)
    unique(daily[daily$StepTotal==0,]$TotalTimeInBed) # all 233 participant-days with 0 steps also had 0 minutes in bed
```

...and participant-days with 0 minutes of light intensity activity
```{r}
  length(daily[daily$LightlyActiveMinutes==0,]$Id)
      unique(daily[daily$LightlyActiveMinutes==0,]$TotalTimeInBed) # all 143 participant-days with 0 minutes of LPA also had 0 minutes in bed
```


Conclusion: Will delete all participant-days with 0 minutes in bed/asleep -- resolves missingness in other variables (e.g., sedentary time == 1440 (no other movement categories captured) and steps (0 steps within a day))
```{r}
  daily<-daily[daily$TotalTimeInBed!=0,]
```

## Remove participant-days with 0 minutes of sedentary time
```{r}
  length(daily[daily$SedentaryMinutes == 0,]$Id)
  daily<-daily[daily$SedentaryMinutes != 0,]
```

## Remove participant-days with < 10 hours of wear time (nonWearMinutes >= 840 minutes per day)
Remove days with < 10 hours of wear time (14+ hours of non-wear time)
```{r}
  length(daily[daily$nonWearMinutes >= 840,]$Id)
  daily<-daily[daily$nonWearMinutes < 840,]
```


## Remove the first day of data for all participants
Explore the mean (and max) minutes of non-wear time by day. Numbers on the first day may be influenced by low wear time (e.g., very high sedentary time, low light, low sleep, etc.). Mean non-wear time is an order of magnitude higher on the first day vs. the second highest day (towards the end of the study/during the post-intervention period). As a resul, will also deleted the first day of data
```{r}
  nonWear<-daily %>% select(date,nonWearMinutes) %>% group_by(date) %>% arrange(date) %>%
                             mutate(meanNonWear=mean(nonWearMinutes),
                                    maxNonWear=max(nonWearMinutes)) %>% 
                             select(-nonWearMinutes,-Id) %>% distinct() %>% arrange(desc(meanNonWear))
  nonWear
  daily<-daily[daily$date>min(daily$date),]
```

## Write CSV file for TOP

```{r}
  names(daily)
  daily<-daily %>% select(-weekNum)
  daily$pubdate<-as.character(daily$pubdate)
  daily[is.na(daily)]<-""
  write.csv(daily,"daily-TOP-CSEP.csv",row.names = FALSE)
```