{"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"4.0.5"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"library(tidyverse) # metapackage of all tidyverse packages\nlibrary(arrow) # for reading parquet files\nlibrary(ggthemes) # for ggplot themes\nlibrary(tictoc) # for measuring command processing times\nlibrary(pryr) # for memory usage data\nlist.files(path = \"../input/child-mind-institute-detect-sleep-states\")\noptions(warn=-1)\nread_train_parquet_imp <- FALSE # flag set to read modified train_series.parquet instead of original.  This saves ~ 45 minutes of manipulation to the original file.\n","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","execution":{"iopub.status.busy":"2023-09-26T23:02:24.656932Z","iopub.execute_input":"2023-09-26T23:02:24.658714Z","iopub.status.idle":"2023-09-26T23:02:24.709685Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Read the train events\ntrain_events <- read.csv(\"../input/child-mind-institute-detect-sleep-states/train_events.csv\")\n\n# Read the the sample submission\nsample_submission <- read.csv(\"../input/child-mind-institute-detect-sleep-states/sample_submission.csv\")","metadata":{"execution":{"iopub.status.busy":"2023-09-26T23:02:24.712798Z","iopub.execute_input":"2023-09-26T23:02:24.714427Z","iopub.status.idle":"2023-09-26T23:02:24.778123Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## From the competition material, we learn the train_events.csv contains \"the sleep logs for series in the training set recording onset and wake events.\"\n\n* series_id - Unique identifier for each series of accelerometer data in train_series.parquet.\n* night - An enumeration of potential onset / wakeup event pairs. At most one pair of events can occur for each night.\n* event - The type of event, whether onset or wakeup.\n* step and timestamp - The recorded time of occurence of the event in the accelerometer series.","metadata":{}},{"cell_type":"code","source":"train_events <- train_events %>% mutate(timestamp = ymd_hms(timestamp))\n\ntrain_events %>% head()","metadata":{"execution":{"iopub.status.busy":"2023-09-26T23:02:24.781109Z","iopub.execute_input":"2023-09-26T23:02:24.782819Z","iopub.status.idle":"2023-09-26T23:02:24.855932Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## We also learn from the competition material that the train_series.parquet file \"is a continuous recording of accelerometer data for a single subject spanning many days.\"\n\n* series_id - Unique identifier for each accelerometer series.\n* step - An integer timestep for each observation within a series.\n* timestamp - A corresponding datetime with ISO 8601 format %Y-%m-%dT%H:%M:%S%z.\n* anglez - As calculated and described by the GGIR package, z-angle is a metric derived from individual accelerometer components that is commonly used in sleep detection, and refers to the angle of the arm relative to the vertical axis of the body\n* enmo - As calculated and described by the GGIR package, ENMO is the Euclidean Norm Minus One of all accelerometer signals, with negative values rounded to zero. While no standard measure of acceleration exists in this space, this is one of the several commonly computed features","metadata":{}},{"cell_type":"code","source":"if(!read_train_parquet_imp){\n   # Read train & test parquet files\n    options(arrow.skip_nul = TRUE) # required to deal with null values, otherwise errors out during parquet read\n    train_series.parquet <- read_parquet(\"../input/child-mind-institute-detect-sleep-states/train_series.parquet\", as_data_frame = TRUE)\n    test_series.parquet <- read_parquet(\"../input/child-mind-institute-detect-sleep-states/test_series.parquet\", as_data_frame = TRUE)\n    \n    # Convert timestamp column (initially a character) to datetime format\n    train_series.parquet <- train_series.parquet %>% mutate(timestamp = ymd_hms(timestamp))\n\n    # If we bring back the first 15 rows, we see that the timestamp entries are in 5 second increments.\n    train_series.parquet %>% head(15)\n\n    train_series.parquet %>% nrow()\n    # 127946340 rows!!\n}else{\n    train_series.parquet <- read_parquet('/kaggle/working/train_series.parquet.imp')\n    \n}","metadata":{"execution":{"iopub.status.busy":"2023-09-26T23:02:24.859062Z","iopub.execute_input":"2023-09-26T23:02:24.860670Z","iopub.status.idle":"2023-09-26T23:05:12.745981Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# So what do these sleep patterns look like?\n### Let's take a look at 5 consecutive days worth of measurements for one series_id.  Note the absence of a onset / awake for 2018-08-19, which we will explore later.","metadata":{}},{"cell_type":"code","source":"# Grab the first series_id \nfirst_series_id <- train_events %>% pull(series_id) %>% unique() %>% head(1)\n\n# Separate out onset & wakeup events in the train_events dataset\nonset_events <- train_events %>% filter(series_id %in% first_series_id, event == 'onset')\nwakeup_events <- train_events %>% filter(series_id %in% first_series_id, event == 'wakeup')\n\n# Adjust plot size\noptions(repr.plot.width = 15, repr.plot.height =6)\n\n# Plot enmo data vs. time, with vertical ref. lines for onset events (pink) and wakeup_events(green)\ntrain_series.parquet %>% filter(series_id %in% first_series_id) %>% ggplot(aes(x=timestamp)) + geom_line(aes(y=enmo), color = 'steelblue') + \nscale_x_datetime(limit=c(as_datetime(\"2018-08-18\"),as_datetime(\"2018-08-23\")), breaks = '1 day') + \ngeom_vline(xintercept = onset_events$timestamp, color = 'pink') + geom_vline(xintercept = wakeup_events$timestamp, color = 'green') +theme_clean()","metadata":{"execution":{"iopub.status.busy":"2023-09-26T23:05:12.749345Z","iopub.execute_input":"2023-09-26T23:05:12.751083Z","iopub.status.idle":"2023-09-26T23:05:25.180616Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Let's do the same for anglez\n\n# Adjust plot size\noptions(repr.plot.width = 15, repr.plot.height =6)\n\n# Plot anglez data vs. time, with vertical ref. lines for onset events (pink) and wakeup_events(green)\ntrain_series.parquet %>% filter(series_id %in% first_series_id) %>% ggplot(aes(x=timestamp)) + geom_line(aes(y=anglez), color = 'steelblue') + \nscale_x_datetime(limit=c(as_datetime(\"2018-08-18\"),as_datetime(\"2018-08-23\")), breaks = '1 day') + \ngeom_vline(xintercept = onset_events$timestamp, color = 'pink') + geom_vline(xintercept = wakeup_events$timestamp, color = 'green') +theme_clean()","metadata":{"execution":{"iopub.status.busy":"2023-09-26T23:05:25.183968Z","iopub.execute_input":"2023-09-26T23:05:25.185679Z","iopub.status.idle":"2023-09-26T23:05:36.939476Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# The \"concrete instructions\" from the competition hosts gives us some ideas to help focus EDA.  I've annotated some things we'll explore for the first few:\n\n1. A single sleep period must be at least 30 minutes in length\n * ***1.a What is the average sleep period length? What is the distribution of sleep period lengths?***\n * ***1.b What are the average onset / awake times? What are their distributions?***   \n2. If no valid sleep window is identifiable, neither an onset nor a wakeup event is recorded for that night.\n * ***2.a How often is a sleep window not identifiable?***\n3. A single sleep period can be interrupted by bouts of activity that do not exceed 30 consecutive minutes\n * **3.aWhat do the sleep periods look like in the .parquet data?**","metadata":{}},{"cell_type":"markdown","source":"# (1a) What is the average sleep period? What is the distribution of sleep period lengths?","metadata":{}},{"cell_type":"code","source":"# To answer these, we can use the train_events.csv.  Since the sleep period doesn't end until the wakeup event, we'll add a sleep_duration value for those rows (in minutes) that have wakeup events.\n\n# Adjust plot size\noptions(repr.plot.width = 15, repr.plot.height =6)\n\n# The lag function works perfectly for these kinds of operations, with an ifelse to determine which row to add the value to.  We use NA as the value for onset rows.\ntrain_events <- train_events %>% mutate(sleep_duration = ifelse(event == 'wakeup', difftime(timestamp, lag(timestamp), units = 'mins'), NA)) \ntrain_events %>% head()\n\n# With that established, the mean sleep_duration in minutes is straightforward:\nmean_sleep <- train_events %>% pull(sleep_duration) %>% mean(na.rm = TRUE)\n\n# As well as the distribution of sleep_duration, with the mean plotted as a vertical red line\ntrain_events %>% ggplot(aes(x=sleep_duration,)) + geom_bar(color = 'steelblue') + geom_vline(xintercept = mean_sleep, color= 'red') + theme_clean()","metadata":{"execution":{"iopub.status.busy":"2023-09-26T23:05:36.943050Z","iopub.execute_input":"2023-09-26T23:05:36.944680Z","iopub.status.idle":"2023-09-26T23:05:37.493041Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# One of the first things that jump out are the outliers in sleep_duration\n# Let's take a look at sleep durations 1000+ to start, which shows two series_ids for night = 20\ntrain_events %>% filter(sleep_duration > 1000)","metadata":{"execution":{"iopub.status.busy":"2023-09-26T23:05:37.496180Z","iopub.execute_input":"2023-09-26T23:05:37.497790Z","iopub.status.idle":"2023-09-26T23:05:37.535082Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# If we filter on series_id = '148471991ffb', night == 19, 20 21, we see both the onset & wakeup entries that produce this large value in sleep_duration\ntrain_events %>% filter(series_id %in% c('148471991ffb'), night %in% c(19, 20, 21))\n\n# Is 1133 a valid sleep_duration, considering this individual only slept a little over 4 hours the night before?","metadata":{"execution":{"iopub.status.busy":"2023-09-26T23:05:37.539368Z","iopub.execute_input":"2023-09-26T23:05:37.541229Z","iopub.status.idle":"2023-09-26T23:05:37.580448Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# We do the same for series_id = 'c3072a759efb', night == 19, 20 21.\ntrain_events %>% filter(series_id %in% c('c3072a759efb'), night %in% c(19,20,21))\n\n# These results are more concerning as it appears we are missing an onset entry for the night = 20.  This has caused the lag function to fail, incorrectly using the previous wakeup row instead of an onset row.\n# We will need to scrub the data a little more or decide to discard orphan onset/wakeup entries.\ndifftime(as_datetime('2018-03-17 13:52:00'), as_datetime('2018-03-18 13:44:00'), units = 'mins') # this shows the incorrect difftime calculation for the two wakeup events.","metadata":{"execution":{"iopub.status.busy":"2023-09-26T23:05:37.584146Z","iopub.execute_input":"2023-09-26T23:05:37.586244Z","iopub.status.idle":"2023-09-26T23:05:37.643465Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# (1b) What are the average onset / awake times? What are their distributions?","metadata":{}},{"cell_type":"code","source":"# When gathering the onset times, we see there are NAs in the dataset.  We'll remove them for now and make a note to decide what to do with them later.\n\n# Next, we want to calculate the minutes_since_midnight, get their average.\navg_onset_mins_since_midnight <- train_events %>% filter(event == 'onset' & !is.na(timestamp)) %>% \nmutate(onset_minutes_since_midnight = as.integer(format(as.POSIXct(timestamp), format = \"%H\")) * 60 + as.integer(format(as.POSIXct(timestamp), format = \"%M\"))) %>% \npull(onset_minutes_since_midnight) %>% mean() \n\n# Now we convert back to a time object, and get our average onset time.\navg_onset <- format(as.POSIXct('1900-1-1') + avg_onset_mins_since_midnight, '%H:%M:%S')\n\n# Then we repeat the process for wakeup\navg_awake_mins_since_midnight <- train_events %>% filter(event == 'wakeup' & !is.na(timestamp)) %>% \nmutate(wakeup_minutes_since_midnight = as.integer(format(as.POSIXct(timestamp), format = \"%H\")) * 60 + as.integer(format(as.POSIXct(timestamp), format = \"%M\"))) %>% \npull(wakeup_minutes_since_midnight) %>% mean() \n\n# Now we convert back to a time object, and get our average wakeup time.\navg_onset_time <- hm(paste(avg_onset_mins_since_midnight %/% 60, floor(avg_onset_mins_since_midnight %% 60), sep=\":\"))\navg_awake_time <- hm(paste(avg_awake_mins_since_midnight %/% 60, floor(avg_awake_mins_since_midnight %% 60), sep=\":\"))\n\npaste('Average onset time is: ',avg_onset_time)\npaste('Average awake time is: ', avg_awake_time)","metadata":{"execution":{"iopub.status.busy":"2023-09-26T23:05:37.647405Z","iopub.execute_input":"2023-09-26T23:05:37.649665Z","iopub.status.idle":"2023-09-26T23:05:37.723294Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# (2a) How often is a sleep window not identifiable?","metadata":{}},{"cell_type":"code","source":"# From our plot of series_id == '038441c925bb' earlier, we noticed the absence of onset/awake events for 2018-08-19.  This is what it looks like in tabular form:\ntrain_events %>% filter(series_id == '038441c925bb') %>% select(-sleep_duration) %>% head(12)\n\n# If we simply sum NAs, we get an odd number and we don't quite get what we are looking for because we identified earlier that there are orphan onset/awake pairs in the dataset:\nsum(is.na(train_events$timestamp))","metadata":{"execution":{"iopub.status.busy":"2023-09-26T23:05:37.727127Z","iopub.execute_input":"2023-09-26T23:05:37.729094Z","iopub.status.idle":"2023-09-26T23:05:37.790170Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# The orphan onset/awake pair we identified earlier looks like this in data -- it simply is missing:\ntrain_events %>% filter(series_id == 'c3072a759efb', night %in% c(19,20,21)) %>% select(-sleep_duration) \n\n# If we want to identify the orphan onset/awake pairs we can use the following code:\ntrain_events %>% group_by(series_id, night) %>% summarize(totals = n()) %>% filter(totals != 2)\n\n# As you can see from the output, not only do we have one series_id that is missing an event entry, but we have three that have an additional one (more than two)\n# Choosing one, it looks like this in the data:\ntrain_events %>% filter(series_id == '655f19eabf1e', night %in% c(2,3,4)) %>% select(-sleep_duration)","metadata":{"execution":{"iopub.status.busy":"2023-09-26T23:05:37.794527Z","iopub.execute_input":"2023-09-26T23:05:37.796712Z","iopub.status.idle":"2023-09-26T23:05:39.813538Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# (3a) What do the sleep periods look like in the .parquet data?  This is a little tricky, we first have to indentify the 'onset' and 'wakeup' timestamps from train_events, then we need to iterate through the remaining rows and determine if they are within a sleep interval or an awake interval.","metadata":{}},{"cell_type":"code","source":"if(!read_train_parquet_imp){\n    \n    # In order to identify sleep periods in the training data, we want to merge the train_events on series_id, timestamp, and step to capture the onset & awake fields\n    train_series.parquet <- left_join(train_series.parquet, train_events %>% filter(!is.na(timestamp)), by=c('series_id', 'timestamp', 'step')) # left_join on series_id, timestamp, & step\n}else{\n    \n    # This will create a dataframe with an event column that is mostly NA, but will have onset & awake populated where appropriate from the train_events.  We can see that by filtering on the event column:\n    train_series.parquet %>% filter(!is.na(event)) %>% head()\n}","metadata":{"execution":{"iopub.status.busy":"2023-09-26T23:05:39.816666Z","iopub.execute_input":"2023-09-26T23:05:39.818312Z","iopub.status.idle":"2023-09-26T23:06:25.462019Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if(!read_train_parquet_imp){\n    tic('Start flagging awake / sleeping')\n    unique_series_ids <- train_events %>% filter(!is.na(timestamp)) %>% pull(series_id) %>% unique()\n\n    for(id in unique_series_ids){\n\n        # Get the onset & awake events separately\n        onset_events <- train_events %>% filter(series_id == id & event == 'onset' & !is.na(timestamp))\n        wakeup_events <- train_events %>% filter(series_id == id & event == 'wakeup'& !is.na(timestamp))\n\n        # Join the data on itself to create intervals of sleep time\n        sleep_interval <- left_join(onset_events, wakeup_events, by=c('series_id', 'night')) \n        sleep_interval <- sleep_interval %>% select(series_id, timestamp.x, timestamp.y) %>% mutate(sleep_int = interval(timestamp.x, timestamp.y)) %>% pull(sleep_int) %>% as.list()\n\n        train_series.parquet[which(train_series.parquet$series_id == id),'event'] <- train_series.parquet[which(train_series.parquet$series_id == id),] %>% mutate(event = ifelse(ymd_hms(timestamp) %within% sleep_interval & is.na(event), 'sleeping', event)) %>% select(event)   \n\n        # Repeat the process to create intervals of awake time. We need to create a new column called 'new_night' that will align the nights correctly, accounting for holes in the dataset.\n        if(nrow(onset_events) > 1){ # we remove one row from the onset_events df, so if there is only one, we won't be able to make awake determinations\n\n                wakeup_events$new_night <- 1:nrow(wakeup_events)\n                onset_events <- onset_events[-1,]\n                onset_events$new_night <- 1:nrow(onset_events)\n\n                awake_interval <- left_join(wakeup_events,onset_events, by=c('series_id', 'new_night'))\n                awake_interval <- awake_interval %>% select(series_id, timestamp.x, timestamp.y) %>% mutate(awake_int = interval(timestamp.x, timestamp.y)) %>% pull(awake_int) %>% as.list() \n\n                train_series.parquet[which(train_series.parquet$series_id == id),'event'] <- train_series.parquet[which(train_series.parquet$series_id == id),] %>% mutate(event = ifelse(ymd_hms(timestamp) %within% awake_interval & is.na(event), 'awake', event)) %>% select(event)  \n        }\n    }\n    toc()\n}","metadata":{"execution":{"iopub.status.busy":"2023-09-26T23:06:25.464873Z","iopub.execute_input":"2023-09-26T23:06:25.466366Z","iopub.status.idle":"2023-09-26T23:54:55.264446Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_series.parquet %>% head()","metadata":{"execution":{"iopub.status.busy":"2023-09-26T23:55:14.078258Z","iopub.execute_input":"2023-09-26T23:55:14.080074Z","iopub.status.idle":"2023-09-26T23:55:14.118220Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if(!read_train_parquet_imp){\n    # Write the parquet file for later use\n    write_parquet(train_series.parquet, 'train_series.parquet.imp')\n}","metadata":{"execution":{"iopub.status.busy":"2023-09-26T23:55:32.210315Z","iopub.execute_input":"2023-09-26T23:55:32.213011Z","iopub.status.idle":"2023-09-26T23:56:27.684115Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Grab the first series_id \nfirst_series_id <- train_events %>% pull(series_id) %>% unique() %>% head(1)\n\nonset_events <- train_events %>% filter(series_id %in% first_series_id, event == 'onset')\nwakeup_events <- train_events %>% filter(series_id %in% first_series_id, event == 'wakeup')\n\n# Adjust plot size\noptions(repr.plot.width = 15, repr.plot.height =6)\n\n# Plot enmo data vs. time, with vertical ref. lines for onset events (pink) and wakeup_events(green)\ntrain_series.parquet %>% filter(series_id %in% first_series_id) %>% ggplot(aes(x=timestamp)) + geom_line(aes(y=enmo, color = event)) + \nscale_x_datetime(limit=c(as_datetime(\"2018-08-18\"),as_datetime(\"2018-08-23\")), breaks = '1 day') + \ngeom_vline(xintercept = onset_events$timestamp, color = 'pink') + geom_vline(xintercept = wakeup_events$timestamp, color = 'green') + theme_clean()","metadata":{"execution":{"iopub.status.busy":"2023-09-26T23:57:14.274901Z","iopub.execute_input":"2023-09-26T23:57:14.277155Z","iopub.status.idle":"2023-09-26T23:57:26.661831Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# The remaining \"concrete instructions\" notes from the competition hosts will be covered in future versions -- thanks!\n\n* No sleep windows can be detected unless the watch is deemed to be worn for the duration (elaborated upon, below)\n* The longest sleep window during the night is the only one which is recorded\n* Sleep events do not need to straddle the day-line, and therefore there is no hard rule defining how many may occur within a given period. However, no more than one window should be assigned per night. For example, it is valid for an individual to have a sleep window from 01h00–06h00 and 19h00–23h30 in the same calendar day, though assigned to consecutive nights\n* There are roughly as many nights recorded for a series as there are 24-hour periods in that series.\n","metadata":{}},{"cell_type":"code","source":"# This is the format for our submission\nsample_submission %>% head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# I'll continue to add as time permits -- hope it was helpful!","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}