In humanfactors/FIPS: Fatigue Impairment Prediction Suite (FIPS)

library(knitr)
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Introduction

This vignette demonstrates how to generate a sleep times dataframe and then proceed to process this for modeling with the Three Process Model of Fatigue. The code here demonstrates the usage from a higher-level user perspective.

Note that in future releases, we may seek to create functions that explicitly generate these sleep times for you. Furthermore, it is possible to simply create the sleeptimes in Excel or another program and transform to the format shown here. Therefore, the sleep generation features are presented here for practical purposes.

library(FIPS)
library(dplyr)
library(tibble)
library(ggplot2)
library(lubridate)
library(colorspace)
library(tidyr)

Data Preparation

Sleep Times Format

The dataframe below shows a prototypical FIPS 'sleep times' dataframe. This form of dataframe is intended for individuals who are manually inputting sleep history (e.g., from pencil forms to a spreadsheet). Below, we generate this data structure ourselves for convenience. Specifically, we want to simulate a scenario where an individual obtains exactly 7.5 hours of sleep per night for 6 nights. There are likely multiple ways to achieve this, but the generation function below allows sufficient flexibility to have continuously offset sleeptimes (e.g., by changing by = '24 hours' to another value).

The sleeptimes should correspond to only one participant (multiple people are not currently supported in FIPS without map functions). The default column names for this dataframe are sleep.id, sleep.start, and sleep.end. It is recommended that you explicitly specify timezones in all functions to avoid any silent errors relating to time zone specifications.

example.sleeptimes <- tibble::tibble(
  sleep.start = seq(
    from = lubridate::ymd_hms('2018-05-01 23:00:00', tz = "Australia/Perth"), 
    to = lubridate::ymd_hms('2018-05-07 17:00:00', tz = "Australia/Perth"),
    by = '24 hours'),
  sleep.end = sleep.start + lubridate::dhours(7.5),
  sleep.id = rank(sleep.start))

print(example.sleeptimes)

Prior to actually conducting the simulation, you must convert the sleep times to the continuous "FIPS_df" format. This format is a continuous time series style dataframe that contains calculated variables to be interpreted by the FIPS model functions. This dataframe contains the following headings: datetime, sleep.id, wake_status, wake_status_int, change_point, switch_direction, status_duration, total_prev, time, day, sim_hours. Information regarding these will be presented in the section below, but first let's quickly run through how to generate the FIPS dataframe from sleep times.

The parse_sleeptimes function from FIPS will takes in 'sleep times' generated previously, and note the arguments below.

sleeptimes — The sleep times dataframe generated previously
series.start — This is the start datetime of the entire simulation series.
series.end — This will be the datetime of the entire simulation series
sleep.start.col — This is the name of your sleep.start column in sleep times if changed from default above
sleep.end.col — This is the name of your sleep.end column in sleep times if changed from default above
sleep.id.col — This is the name of your sleep.id column in sleep times if changed from default above
roundvalue — The epoch of the series (i.e., rounding value). Please read further information below!

The setting of roundvalue is critically important. It determines the epoch or spacing between observations in your dataset (in minutes). At a value of 1, the simulation is updated every 1 minute. Consequently, this will increase the size (in rows) of your dataset by a factor of 5 (relative to 5 minutes). In most cases, 5, 10 or even 15 minutes should be sufficient.

Moreover, note that all sleep observations will have their datetime rounded to this value. For example, with a roundvalue = 5 the datetime 2018-05-07 21:02:02 would be rounded to 2018-05-07 21:00:00. For this reason, it is ideal if your series.start and series.end are rounded to the same epoch value as you request. Seconds should never be included in your datetime (biomathematical models are not sensitive to this time resolution anyway).

# Simulation start date time (i.e., when you want first predictions to begin)
simulation.start = lubridate::ymd_hms('2018-05-01 07:00:00', tz = "Australia/Perth")

# Simulation end date time (i.e., when you want predictions to end)
# In this case it ends 
simulation.end = lubridate::ymd_hms('2018-05-07 21:00:00', tz = "Australia/Perth")

# The Continuous FIPS_df dataframe format
# This creates the format ready for simulation
simulated.dataframe = parse_sleeptimes(
  sleeptimes = example.sleeptimes,
  series.start = simulation.start,
  series.end = simulation.end,
  sleep.start.col = "sleep.start",
  sleep.end.col = "sleep.end",
  sleep.id.col = "sleep.id",
  roundvalue = 5)

print(simulated.dataframe)

Bitvector Sequence Format

It is common to represent sleep history information as a bitvector (i.e., a sequence of 1's and 0's). In bitvectors, sleep and wake times are represented by 1's or 0's, with each bit representing an equal time duration (e.g., 1 minute in that status). The bitvector sequence also must be relative to a start datetime. This format is commonly outputted from actigraphy devices (a form of wearable sleep tracker) and corresponding sleep detection algorithms (e.g., Cole-Kripke). Other proprietary software packages (e.g., SAFTE-FAST) require imported data to be in bit vector form (though the exact required formats do vary).

FIPS expects bitvectors to repesent sleep as 0 and wake as 1. The parse_sleepwake_sequence function can transform a bit vector sequence to a compliant FIPS_df. Below, an example of this data and steps to transform are provided, however, note that we do not use this dataframe again within this vignette.

# Simulation start date time (i.e., when you want first predictions to begin)
simulation.start = lubridate::ymd_hms('2018-05-01 10:00:00', tz = "Australia/Perth")

# Example bitvector sequence. This typically would be imported directly via a textfile.
# Here we generate, though typically this would be returned by a ReadLines/read.delim/read.table
bv.sleep.sequence = rep(rep(c(1,0), 6), sample(20:40, 12))

bv.sim.dataframe = parse_sleepwake_sequence(
  seq = bv.sleep.sequence,
  series.start = simulation.start,
  epoch = 15)

print(bv.sim.dataframe)

Modelling and the Simulation Dataframe

Now that you have generated the FIPS format, you can now apply the FIPS simulation functions to this to actually run BMM predictions on the series. In the example below, we will run a Three Process Model simulation over the FIPS_df series we just created. To do this, we will use the FIPS_simulation function, which takes in three arguments: a FIPS_df object, a specification of a modeltype (see help for model types currently implemented), and a pvec which is a vector of parameters for the model.

The parameter vectors provided by default in the package are those reported by Ingre et al. (2014). Please see the help files for citations and further information. These defaults are customisable, and if you are using Rstudio or Emacs ESS, you will get full autocompletion with descriptions of each parameter (See help("TPM_make_pvec", "FIPS")).

We will use the default parameter vector for the Three process model created by FIPS::TPM_make_pvec(), which takes on the following values:

pdf <- data.frame("Parameter" = names(TPM_make_pvec()),
                  "Value" = unname(TPM_make_pvec()))

knitr::kable(list(pdf[1:7,], pdf[8:15,]), row.names=FALSE, digits=3)

Calling FIPS_simulate() will the produces a FIPS_df with model predictions and predicted time-varying process variables (e.g., s, c). The returned FIPS_df will also now inherit the FIPS_simulation class. A FIPS_simulation object has attributes containing the parameter vector, the modeltype string, and the pvec used, and several other values used internally. A custom print function of the object will reveal all this information. Note that this print function does mask some of the columns for ease of reading.

# Run a simulation with the three process model
TPM.simulation.results = FIPS_simulate(
  FIPS_df = simulated.dataframe, # The FIPS_df
  modeltype = "TPM",             # three process model
  pvec = TPM_make_pvec()       # parameter vector
  )

TPM.simulation.results

datetime = vector of datetime stamps separated at equidistance intervals.
sleep.id = a supplementary variable indicating sleep episode identifier.
wake_status = Awake (T) or asleep (F) at that epoch interval/epoch
wake_status_int = Awake (1) or asleep (0) at that epoch interval/epoch
change_point = Whether the individual changed wake status at that interval/epoch.
switch_direction = Whether switch was to sleep or to wake
status_duration = How long individual has been in status at that current time point
total_prev = If a switch has occured, how long were that in the previous status.
time = time of day in decimal hour
day = days into simulation
sim_hours = hours simulation has run for total
s — Estimate of S (homeostatic) process at that time
l — Estimate of l (lower asmytope) process at that time
c — Estimate of C (24-hour circadian) process at that time
w — Estimate of W (sleep intertia) process at that time
u — Estimate of U (12-hour circadian) process at that time
alertness — Currently just s + c + u
KSS — This is equal to 10.6 + -0.6 * (alertness)

Plotting

Plot 1: Parameter Plots

FIPS provides a default plot method to visualize model predictions and time-varying process estimates over time, discussed in detail in the plotting vignette. These plots aid with debugging and understanding how the different model parameters contribute to predictions.

# Plot the whole time series
plot(TPM.simulation.results, plot_stat=c("alertness", "s", "c"))

#Narrow in on the first 90 hours
plot(TPM.simulation.results, plot_stat=c("alertness", "s", "c"),
      from= '2018-05-01 23:00:00',
      to= as_datetime('2018-05-01 23:00:00') + hours(90))

Plot 2: Heatmap Prediction Plots

There are many other ways that visualizations could help inform theory and practice. The below heat plot is an example of how dangerous times in a mission, according to the TPM predictions, could be visualized. In the example below, the regions in deep orange indicate a greatly increased risk of fatigue. However, the fact that the increased fatigue occurs at the start of the mission indicates the initalisation parameters (e.g., S0) are to blame, so this isn't a cause for significant concern.

if (!requireNamespace("colorspace", quietly = TRUE)) {
    stop("Package \"colorspace\" needed for this example to work. Please install it.",
      call. = FALSE)
  }

plot3_w = TPM.simulation.results %>%
  # Change epoch to 30 minutes for this plot (optional)
  mutate(time = strftime(datetime, format = "%H:%M")) %>% 
  separate(time, into = c("plot_hour", "plot_minute"), ":") %>%
  mutate(halfhour = if_else(plot_minute < 30, 0, 0.5) + as.numeric(plot_hour)) %>% 
  group_by(day, plot_hour, halfhour) %>% 
  summarise(KSS = mean(KSS), .groups = "keep") %>% 
  filter(day > 1) %>% 
  # Start Plot
  ggplot(aes(day, halfhour, fill = KSS)) +
  geom_tile(aes(fill = KSS), color = "white") +
  scale_x_continuous(breaks = seq(0,7,1)) +
  scale_y_continuous(breaks = seq(0,23,1)) +
  labs(y = "Hour of Day", x = "Day of Mission") +
  colorspace::scale_fill_continuous_divergingx(name = "KSS", palette = "Zissou 1", limits = c(1,9), mid = 5) +
  ggtitle(label = "Three Process Model (KSS)", subtitle = "Heatmap of TPM predicted performance by 0.5hr") +
  theme(axis.title = element_text(size = 12)) 

plot(plot3_w)