ACTman: Actigraphy manager

Documented in ACTman

#############################################################################
### ACTman package                                                        ###
### Script authors: Yoram Kunkels, Stefan Knapen, & Ando Emerencia        ###
### Most recent Update: 16-04-2018                                        ###
### Supported devices: Actiwatch 2 Respironics & MW8                      ###
###=======================================================================###
### Revision History:                                                     ###
### 16-04-2018: Added Actogram Functionality.                             ###
###~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~###

#' ACTman - Actigraphy Manager
#'
#' ACTman manages actigraphy data whilst offering pre-processing and analyses options.
#' This initial version supports the 'Actiwatch 2' and 'MotionWatch 8' actigraphy devices,
#' whilst allowing for both sleep and circadian rhythm analyses.
#'
#' @param workdir The working directory of the script.
#' @param sleepdatadir An optional vector specifying the directory for actogram and sleep analysis data.
#' @param myACTdevice Name of the input device used. Should be either 'Actiwatch2' or 'MW8'.
#' @param iwantsleepanalysis Boolean value indicating whether sleep analysis should be performed.
#' @param plotactogram Value indicating if and what kind of actogram has to be plotted. Can be either '48h', '24h', or FALSE.
#' @param selectperiod Boolean value indicating whether a specific period has to be selected.
#' @param startperiod An optional vector specifying single or multiple period starts. Should be in the format "2016-10-03 00:00:00".
#' @param daysperiod An optional vector specifying the length in days of the period.
#' @param endperiod An optional argument that is a date string (format: "2016-10-03 00:00:00"), denoting the end of the data subset to be analyzed. Only used if daysperiod is not specified.
#' @param movingwindow Boolean value indicating whether a moving window should be utilised.
#' @param movingwindow.size An optional vector specifying the length in days of the moving window. Default is 14 days.
#' @param movingwindow.jump An optional vector specifying the length of the jumps with which the moving window is shifted each iteration. Default is 1 day.
#' @param circadian_analysis Boolean value indicating whether non-parametric circadian rhythm analysis should be performed.
#' @param nparACT_compare Boolean value indicating that comparison with another actigraphy R package should be performed. If TRUE, the values for IS, IV, RA, L5, L5_starttime, M10, and M10_starttime of the nparACT_base_loop function are recorded in the returned overview variable.
#' @param na_omit Boolean value indicating whether NA's should be omitted.
#' @param na_impute Boolean value indicating whether NA's should be imputed.
#' @param missings_report Boolean value indicating whether missings promt should appear.
#' @param lengthcheck Boolean value. If TRUE, the dataset is shortened to the start date plus 14 days, and observations more than 14 days after the start date are removed.
#'
#' @return if iwantsleepanalysis, this returns the sleepdata overview, else if movingwindow, it returns the moving window results, and otherwise it returns the actdata overview.
#' @examples
#' \dontrun{
#' View(ACTman::ACTman(workdir = "C:/Bibliotheek/Studie/PhD/Publishing/ACTman/R-part/mydata2",
#'                     myACTdevice = "Actiwatch2",
#'                     iwantsleepanalysis = FALSE,
#'                     plotactogram = FALSE,
#'                     selectperiod = FALSE,
#'                     startperiod = "2016-10-03 00:00:00",
#'                     daysperiod = 14, movingwindow = FALSE, circadian_analysis = TRUE))
#' }
#' @importFrom stats na.omit
#' @importFrom utils View
#' @importFrom utils read.csv
#' @importFrom utils read.table
#' @importFrom utils tail
#' @importFrom utils write.table
#'
#' @export
ACTman <- function(workdir = "C:/Bibliotheek/Studie/PhD/Publishing/ACTman/R-part/mydata",
                   sleepdatadir = paste("C:/Bibliotheek/Studie/PhD/Publishing/",
                                        "ACTman/R-part/Actogram & Sleep analysis", sep = ""),
                   myACTdevice = "Actiwatch2", iwantsleepanalysis = FALSE, plotactogram = FALSE,
                   selectperiod = FALSE, startperiod = NULL, daysperiod = FALSE, endperiod = NULL,
                   movingwindow = FALSE, movingwindow.size = 14, movingwindow.jump = 1,
                   circadian_analysis = TRUE, nparACT_compare = FALSE, na_omit = FALSE, na_impute = FALSE,
                   missings_report = TRUE, lengthcheck = TRUE, i_want_EWS = FALSE) {

  ## Step 1: Basic Operations-----------------------------------------------------------------

  ## Set working directory:
  setwd(workdir)

  ## List actigraphy files and seperate out sleeplogs:
  pattern_file <- ".csv"
  ACTdata.files <- sort(list.files(getwd(), pattern = pattern_file))
  if (any((grep(pattern = "sleeplog", x = ACTdata.files)))) {
    ACTdata.files <- ACTdata.files[-(grep(pattern = "sleeplog", x = ACTdata.files))] # Remove any sleeplogs from data listing
  }
  if (any((grep(pattern = "markers", x = ACTdata.files)))) {
    ACTdata.files <- ACTdata.files[-(grep(pattern = "markers", x = ACTdata.files))] # Remove any markers files from data listing
  }

  ## Initialise overview:
  ACTdata.overview <- data.frame("filename" = ACTdata.files, "start" = NA, "end" = NA, "end2" = NA,
                                 "numberofobs" = NA, "numberofobs2" = NA, "recordingtime" = NA,
                                 "recordingtime2" = NA, "summertime.start" = NA, "summertime.end" = NA, "missings" = NA,
                                 "missings_perc" = NA, "IS" = NA, "IV" = NA, "RA" = NA, "L5" = NA, "L5_starttime" = NA, "M10" = NA,
                                 "M10_starttime" = NA, "r2.IS" = NA, "r2.IV" = NA, "r2.RA" = NA, "r2.L5" = NA,
                                 "r2.L5_starttime" = NA, "r2.M10" = NA, "r2.M10_starttime" = NA, "last5act.active" = NA,
                                 "lengthcheck" = NA)

  ## Initiate parameters:
  i <- 1 # set i
  secshour <- 60 * 60 # Seconds per hour.
  secsday <- 24 * secshour # Seconds per day.
  secs14day <- secsday * 14 # Seconds in 14 days.
  minsaday <- (secsday / 60) # Minutes per day.

  ## Semantic checks:
  if (selectperiod && length(startperiod) != length(ACTdata.files)) {
    stop(paste("The number of start periods does not match the number of data files found:", startperiod, length(ACTdata.files)))
  }

  if (".mtn" %in% substr(list.files(getwd()), nchar(list.files(getwd())) - 4 + 1, nchar(list.files(getwd())))) {
    message(paste("There is at least 1 unsupported Actigraphy file format present in the working directory!
                  Please convert, stash, or remove these unsupported files before rerunning:"))
    print(list.files(getwd())[grep(".mtn", list.files(getwd()))])
    stop()
  }


  ## Step 2: Main ACTman Loop------------------------------------------------------------------------
  ## Description: ...

  for (i in 1:length(ACTdata.files)) {

    ## Test for user-defined arguments:
    if (myACTdevice != "MW8" && myACTdevice != "Actiwatch2") {
      stop(paste("Unknown value for myACTdevice (should be MW8, Actiwatch2):", myACTdevice))
    }

    print(paste("*** Start of Dataset", i, "***"))
    print("")
    print(paste("Dataset Name:", ACTdata.overview[i, "filename"]))
    print("")


    ## Step 2.1.: Reading Data-----------------------------------------------------------------------

    ## Reading in .CSV data, plus some device-specific Data Management:
    ## Device-specific Data Management is required as raw data format differs between
    ## devices, in e.g., headers and raw data locations.
    if (myACTdevice == "Actiwatch2") {

      ACTdata.1 <- read.csv(paste(ACTdata.files[i]), header = FALSE)
      ACTdata.1.sub <- ACTdata.1[, c(4, 5, 6)]
      colnames(ACTdata.1.sub) <- c("Date", "Time", "Activity")

    } else if (myACTdevice == "MW8") {

      ACTdata.1 <- read.csv(paste(ACTdata.files[i]), header = FALSE, fill = TRUE, stringsAsFactors = FALSE, col.names = c("A", "B", "C"))

      if (all(is.na(ACTdata.1$B)) && all(is.na(ACTdata.1$C))) { ## Note! If TRUE ACTdata.1 is (suddenly) tab-seperated!

        ACTdata.1 <- read.csv(paste(ACTdata.files[i]), header = FALSE, fill = TRUE, stringsAsFactors = FALSE, col.names = c("A", "B", "C"), sep = "\t")

      }




      if (any(ACTdata.1[, 1] == "Raw data:")) {
        ACTdata.1 <- as.data.frame(ACTdata.1[((which(ACTdata.1[, 1] == "Raw data:")) + 2):nrow(ACTdata.1), ])
      } else {
        ACTdata.1 <- read.csv(paste(ACTdata.files[i]), header = TRUE, fill = TRUE, stringsAsFactors = FALSE, col.names = c("A", "B", "C"))
      }

      ## Make a copy of the original data to work on:
      ACTdata.1.sub <- ACTdata.1
      colnames(ACTdata.1.sub) <- c("Date", "Time", "Activity")
      ACTdata.1.sub$Activity <- as.numeric(ACTdata.1.sub$Activity)

      ## Test for 30 sec. bins:
      if (any(grepl(pattern = ":30", x = ACTdata.1$B[1:2]))) {
        print("Detecting Epoch Length.......")
        print("Warning: 30 sec. Epoch's Detected!")
        print("Action: Binning 30 sec. Epochs in 60 sec. Epochs")
        print("")

        ## If epochs are 30 sec. instead of 60 sec., bin them together to form 60 sec. epochs.
        ACTdata.TEMP <- ACTdata.1[(grepl(pattern = ".*(?<!:30)$", x = ACTdata.1$B, perl = TRUE)), ]

        halfminute_data <- as.numeric(ACTdata.1[(grepl(pattern = ".*(?<=:30)$", x = ACTdata.1$B, perl = TRUE)), ]$C)
        if (length(grep(pattern = ":30$", x = ACTdata.1$B[1]))) {
          # if it starts wtih :30, throw first val away
          halfminute_data <- tail(halfminute_data, n = -1)
        }
        if (length(grep(pattern = ":00$", x = ACTdata.1$B[length(ACTdata.1$B)]))) {
          # if it ends on a full minute, add a zero value
          halfminute_data <- c(halfminute_data, 0)
        }
        ACTdata.TEMP$C <- as.numeric(ACTdata.TEMP$C) + halfminute_data
        # #! Workaround for aforementioned issue
        # ACTdata.TEMP <- ACTdata.TEMP[ - which(grepl("00:30", ACTdata.TEMP$B)), ]

        ## Write binned data in ACTdata.TEMP to workable data object ACTdata.1.sub:
        ACTdata.1.sub <- ACTdata.TEMP
        colnames(ACTdata.1.sub) <- c("Date", "Time", "Activity")
        ACTdata.1.sub$Activity <- as.numeric(ACTdata.1.sub$Activity)
        rm(ACTdata.TEMP) # Remove temporary data object
        rm(halfminute_data)

      } else {
        ## Make no changes if 60 sec. bins
        print("Detecting Epoch Length.......")
        print("Normal 60 sec. Epochs detected")
        print("No changes made")
        print("")
      }
    }

    ## Step 2.2: Managing the Data---------------------------------------------------------------

    ## Reformat dates and times into required format for further processing:
    ACTdata.1.sub$Date <- gsub(pattern = "/20", replacement = "/", x = ACTdata.1.sub$Date) # Take only last two year digits
    ACTdata.1.sub$Date <- paste(ACTdata.1.sub$Date, ACTdata.1.sub$Time) # Merge Date Time

    if (grepl("-", ACTdata.1.sub$Date[1])) { # Reformat Date
      ACTdata.1.sub$Date <- strptime(ACTdata.1.sub$Date, "%Y-%m-%d %H:%M:%S")
    } else {
      ACTdata.1.sub$Date <- strptime(ACTdata.1.sub$Date, "%d/%m/%y %H:%M:%S")
    }

    ACTdata.1.sub$Time <- NULL # Remove empty Time variable


    ## Check for empty first row and if so, remove it:
    if (all(is.na(ACTdata.1.sub[1, ]))) { ACTdata.1.sub <- ACTdata.1.sub[-1, ] }


    ## Period Selection (user-defined option):
    ## Start by obtaining the row location of the user-defined start of period (startperiod)
    ## (format: "2016-10-03 00:00:00").
    if (selectperiod) {
      startperiod.loc <- which(ACTdata.1.sub$Date == startperiod[i])
      ## Then, if number of days from startperiod (daysperiod) is given, take only data from
      ## start (startperiod.loc) untill start + specified number of days (startperiod.loc + (daysperiod*minsaday)).
      if (daysperiod) {
        ACTdata.1.sub <- ACTdata.1.sub[(startperiod.loc:(startperiod.loc + (daysperiod*minsaday))), ]
        ## Else, if no daysperiod if given, see if user-defined end of period (endperiod) is given.
        ## If so, take only data from user-defined start of period to end of period.
      } else if (endperiod %in% ACTdata.1.sub$Date) {
        endperiod.loc <- which(ACTdata.1.sub$Date == endperiod)
        ACTdata.1.sub <- ACTdata.1.sub[(startperiod.loc:endperiod.loc), ]
        ## Else, take only data from start of period to end of dataset.
      } else {
        ACTdata.1.sub <- ACTdata.1.sub[(startperiod.loc:(nrow(ACTdata.1.sub))), ]
      }
    }


    ## Write values to overview file:
    ## Start- and end dates and times of actigraph data:
    start_date <- ACTdata.1.sub$Date[1] # Get start date
    end_date <- ACTdata.1.sub$Date[nrow(ACTdata.1.sub)] # Get end date.
    nr_obs <- nrow(ACTdata.1.sub) # Get number of observations.
    ACTdata.overview[i, "start"] <- as.character(start_date) # Write start date to overview.
    ACTdata.overview[i, "end"] <- as.character(end_date) # Write end date to overview.
    ## Recordingtime (time difference between start and end date):
    ACTdata.overview[i, "recordingtime"] <- round((as.POSIXct(start_date) - as.POSIXct(end_date)), 2) # write recordingtime to overview
    ## Number of observations:
    ACTdata.overview[i, "numberofobs"] <- nr_obs


    ## Identify Last Whole 24 Hour component and its Position:
    ACTdata.1.sub.lastwhole24h <- ACTdata.1.sub[tail(grep("00:00:00", ACTdata.1.sub$Date), 2), "Date"]
    ACTdata.1.sub.lastwhole24h <- ACTdata.1.sub.lastwhole24h[1]


    ## Add 14 days in a way that respects daylight savings time changes:
    #! Number of days needs to be made dynamic!! Needs to correspond to number of days in dataset!!
    ACTdata.1.sub.14day <- increase_by_days(ACTdata.1.sub$Date[1], 14)


    ## Lengthcheck:

    if (lengthcheck) {
      ## If Dataset is longer than Start Date plus 14 days, Remove data recorded thereafter:
      print("Task 2: Detecting if Dataset is longer than Start Date plus 14 full days")
      if (ACTdata.1.sub[nrow(ACTdata.1.sub), "Date"] > ACTdata.1.sub.14day) {
        print("Warning: Dataset is longer than Start Date plus 14 full days!")

        ACTdata.1.sub <- ACTdata.1.sub[1:((secs14day/60) + 1), ]

        print("Action: Observations recorded after Start Date plus 14 full days were removed.")
        print("Task 2 DONE: Dataset shortened to Start Date plus 14 days. Tasks 3, 4, and 5 also OK")
        print("")
        ACTdata.overview$lengthcheck[i] <- TRUE
      } else {
        print("Task 2 OK: Dataset not longer than Start Date plus 14 days.")
        print("")
      }

      # Update overview after removal
      ACTdata.overview[i, "numberofobs2"] <- nrow(ACTdata.1.sub)
      ACTdata.overview[i, "recordingtime2"] <- round(as.POSIXct(ACTdata.1.sub$Date[1]) - as.POSIXct(ACTdata.1.sub$Date[nrow(ACTdata.1.sub)]), 2)
      ACTdata.overview[i, "end2"] <- as.character(ACTdata.1.sub$Date[nrow(ACTdata.1.sub)]) # write end date to overview
    }



    ## Handling of Missing Data (NA's):
    ## Write Missings values (total number of missings and percentage of total data missing)
    ## to overview file:
    ACTdata.overview[i, "missings"] <- sum(is.na(ACTdata.1.sub$Activity))  # write missings to overview
    ACTdata.overview[i, "missings_perc"] <- round(ACTdata.overview[i, "missings"] / ACTdata.overview[i, "numberofobs"], 3) # write missings percentage to overview
    ## Report missings in console:
    print("Task: Reporting NA's")
    print(paste("Number of NA's in this Dataset:", ACTdata.overview[i, "missings"]))
    print(paste("This is:", round(ACTdata.overview[i, "missings"] / ACTdata.overview[i, "numberofobs"], 3), "% of the total number of observations!"))
    print("")
    ## If user-defined argument "na_omit" is TRUE, then use na.omit{stats} to row-wise delete NA's:
    if (na_omit) {
      print("Row-wise removal of NA's as user defined na.omit = TRUE")
      ACTdata.1.sub <- na.omit(ACTdata.1.sub) #! Creates error because many NA's!!!!
      print("All NA's removed!")
    }
    ## If user-defined argument "na_impute" is TRUE, then use mice{mice} to impute missings through
    ## Multivariate Imputation by Chained Equations (MICE). This installs the 'mice' package and dependencies.
    if (na_impute) {
      ## Impute Missings
      # if (!require('mice')) { install.packages('mice', dep = TRUE)}; library('mice')
      # if (!require('pscl')) { install.packages('pscl', dep = TRUE)}; library('pscl')
      # if (!require('accelmissing')) { install.packages('accelmissing', dep = TRUE)}; library('accelmissing')
      tempData <- mice::mice(matrix(data = c(ACTdata.1.sub$Activity, rep.int(x = 0, times = (ACTdata.overview[i, "numberofobs"]))), ncol = 2), m = 5, maxit = 50, meth = "pmm", seed = 500)
      tempData2 <- mice::complete(tempData, 1)
      ACTdata.1.sub$Activity <- tempData2$V1
    }

    ## User-control over Analysis if too much Missings:
    ## Too much is specified in this case as > 0.01% missing of total dataset.
    #! This 0.01% is arbitrarily chosen, as of now no suitable validated criterion is yet found.
    ## Initialise exception handling for when there are no missings:
    number_of_missings <- ifelse(is.na(ACTdata.overview[i, "missings"]), 0, ACTdata.overview[i, "missings"])
    ## If a missings-prompt is required (default is missings_report = TRUE), ..
    if (missings_report && !na_impute) {
      ## .. see if number of missings exceeds 0.01% of total number of data.
      if ((number_of_missings / ACTdata.overview[i, "numberofobs"]) > 0.01) {
        ## If so, explain situation to user via text prompt, and give them the choice to
        ## either continue or stop with the analyses.
        message("\nMore than 0.01% of data is missing!\nAnalysis results might deviate from true values!")
        message("Do you want to continue?")
        missings_prompt_answer <- readline(prompt = "Enter 'y' for Yes or 'n' for No:")
        ## If user decides to stop the analyses, stop the analyses and report stop message.
        if (missings_prompt_answer == "n") {
          message("Stopped by user!")
          ({break()})
        }
        ## If user decides to continue the analyses, continue and report continue message.
        if (missings_prompt_answer == "y") {
          print("Continue analysis with > 0.01% missings")
          print("")
        } else {
          message("Unknown input!")
          ({break()})
        }
      }
    }


    # ## Set NA to 0
    # #! Debug and/or possible future functionality if required
    # ACTdata.1.sub[is.na(ACTdata.1.sub)] <- 0


    ## Check if there is activity in the tail of the dataset. As sometimes at the end of the study
    ## the actigraph is handed over by the participant, but not immediately stopped.
    ## If Activity in Last 5 observations is on average zero, Skip to Last Activity:
    ACTdata.1.sub.last5act <- ACTdata.1.sub$Activity[(nrow(ACTdata.1.sub) - 4):nrow(ACTdata.1.sub)] # Last 5 activity counts in dataset
    ACTdata.1.sub.last5act.active <- sum(ACTdata.1.sub.last5act, na.rm = T) >= (5 * length(ACTdata.1.sub.last5act)) # Is there on average more than 5 counts per obs?
    print("Task: Checking for Activity in Last 5 observations")
    if (ACTdata.1.sub.last5act.active == FALSE) {
      print("Warning: No Activity in Last 5 observations!")
      print("Last 5 Activity Counts, before Correction:")
      print(ACTdata.1.sub.last5act)
      ACTdata.1.sub <- ACTdata.1.sub[1:max(which(ACTdata.1.sub$Activity >= (5 * length(ACTdata.1.sub.last5act)))), ] # Shortens data untill reached last activity
      ACTdata.1.sub.last5act <- ACTdata.1.sub$Activity[(nrow(ACTdata.1.sub) - 4):nrow(ACTdata.1.sub)] # Last 5 activity counts in dataset
      ACTdata.overview$last5act.active[i] <- FALSE
      print("Last 5 Activity Counts, after Correction:")
      print(ACTdata.1.sub.last5act)
      print("Task DONE: Dataset Skipped to last Activity.")
      print("")
    } else {
      print("Task OK: Dataset contained Activity in Last 5 observations.")
      print("")
    }


    ## Update overview file:
    ## Write number of observations, total recording time, and end date/time, as these values
    ## might have been altered after previous data-management steps.
    ACTdata.overview[i, "numberofobs2"] <- nrow(ACTdata.1.sub)
    ACTdata.overview[i, "recordingtime2"] <- round(as.POSIXct(ACTdata.1.sub$Date[1]) - as.POSIXct(ACTdata.1.sub$Date[nrow(ACTdata.1.sub)]), 2)
    ACTdata.overview[i, "end2"] <- as.character(ACTdata.1.sub$Date[nrow(ACTdata.1.sub)]) # write end date to overview


    ## Step 2.3: Write managed data to file for analyses-------------------------------------------
    ## Create a new directory in working directory for writing managed data files.
    dir.create(file.path(workdir, "Managed Datasets"), showWarnings = FALSE)
    setwd(file.path(workdir, "Managed Datasets"))
    ## Set parameters and go to new directory
    wd <- getwd()
    name <- paste(gsub(pattern = ".csv", replacement = "", x = ACTdata.files[i]))
    newdir <- paste(wd, name, sep = "/")
    dir.create(newdir, showWarnings = FALSE)
    setwd(newdir)
    ## Write managed data:
    ACTdata.1.sub$Date = format(ACTdata.1.sub$Date, "%Y-%m-%d %H:%M:%S")
    write.table(ACTdata.1.sub, quote = FALSE, row.names = FALSE, 
      col.names = FALSE, file = paste(gsub(pattern = ".csv", 
        replacement = "", x = ACTdata.files[i]), "MANAGED.txt"))

    ## Step 2.4: Initialising analyses and funtionalities--------------------------------------
    ## Description: .....

    ## Read managed dataset for analyses and functionalities:
    CRV.data <- read.table(file = file.path(newdir, paste(gsub(pattern = ".csv", replacement = "", x = ACTdata.files[i]), "MANAGED.txt")),
                           stringsAsFactors = FALSE)
    colnames(CRV.data) <- c("Date", "Time", "Activity")


    sleepdata.overview <- NULL
    rollingwindow.results <- NA
    ## Moving/Rolling Window
    ## Check first if Moving Window is required, as this requires it's own analysis calls.
    #! Add Sleep-analysis for Rolling Window!
    #! Add possibility to change 'jump-length' of rolling window (now 1 day) to multiple days!
    if (movingwindow) {
      rollingwindow <- function(x, window, jump) {
        ## Initialise parameters:
        out <- data.frame()
        n <- nrow(x)
        rollingwindow.results <- as.data.frame(matrix(nrow = (floor(((n - window) / jump))), ncol = 21))
        ## Set number of iterations at number of rows of (data - windowsize) / (minutes per day (1440) * jump)
        for (i in 1:((floor(((n - window) / jump))) + 1)) {
          ## Take data as 1 till windowsize for first iteration, for further iterations take
          ## data as starting at ((iteration - 1) * (minutes per day * jump)), and ending at
          ## ((iteration - 1) * (minutes per day * jump)) plus windowsize.
          if (i == 1) {
            out <- x[i:window, ]
          } else {
            out <- x[((i - 1) * jump):(((i - 1) * jump) + window), ]
          }
          ## Write selected period to dataset (CRV.data) and add relevant column names:
          CRV.data <- out
          if (ncol(CRV.data) > 2) {
            colnames(CRV.data) <- c("Date", "Time", "Activity")
          } else if (ncol(CRV.data) == 2) {
            colnames(CRV.data) <- c("Date", "Activity")
          }

          ## Use the nparcalc{ACTman} function to calculate circadian rhythm variables over
          ## the selected period within the window. Write results to rollingwindow.results.
          r2 <- nparcalc(myACTdevice = myACTdevice, movingwindow = movingwindow, CRV.data = CRV.data, ACTdata.1.sub = ACTdata.1.sub, out = out)
          rollingwindow.results[i, 1] <- as.character(strftime(CRV.data[1, "Date"], format = "%Y-%m-%d %H:%M:%S"))
          rollingwindow.results[i, 2] <- as.character(strftime(CRV.data[nrow(CRV.data), "Date"], format = "%Y-%m-%d %H:%M:%S"))
          rollingwindow.results[i, 3] <- r2$IS
          rollingwindow.results[i, 4] <- r2$IV
          rollingwindow.results[i, 5] <- round(r2$RA, 2)
          rollingwindow.results[i, 6] <- round(r2$L5, 2)
          rollingwindow.results[i, 7] <- r2$L5_starttime
          rollingwindow.results[i, 8] <- round(r2$M10, 2)
          rollingwindow.results[i, 9] <- r2$M10_starttime
          rollingwindow.results[i, 10] <- r2$Mean
          rollingwindow.results[i, 11] <- r2$Variance
          rollingwindow.results[i, 12] <- r2$SD
          rollingwindow.results[i, 13] <- r2$CoV
          rollingwindow.results[i, 14] <- r2$Skewness
          rollingwindow.results[i, 15] <- r2$Kurtosis
          rollingwindow.results[i, 16] <- r2$Autocorr
          rollingwindow.results[i, 17] <- r2$Autocorr_lag2
          rollingwindow.results[i, 18] <- r2$Autocorr_lag3
          rollingwindow.results[i, 19] <- r2$Autocorr_lag60
          rollingwindow.results[i, 20] <- r2$Autocorr_lag120
          rollingwindow.results[i, 21] <- r2$Time_to_Recovery
          colnames(rollingwindow.results) <- c("starttime", "endtime", "IS", "IV", "RA", "L5", "L5_starttime",
                                               "M10", "M10_starttime", "Mean", "Variance", "SD",
                                               "Coeff_of_Var", "Skewness", "Kurtosis", "Autocorr_lag1",
                                               "Autocorr_lag2", "Autocorr_lag3", "Autocorr_lag60",
                                               "Autocorr_lag120", "Time_to_Recovery")

          ## Report the calculated circadian results in console:
          print("---------------------------------------------------------------------------------")
          print(paste("Rolling window CRV analysis output - Window step:", (i - 1)))
          print(paste("Begin time:", CRV.data[1, "Date"]))
          print(paste("End time:", CRV.data[nrow(CRV.data), "Date"]))
          print(paste("nOBS:", nrow(CRV.data)))
          print("")
          print("Circadian Rhythm Variables")
          print(paste("IS: ", r2$IS))
          print(paste("IV: ", r2$IV))
          print(paste("RA: ", round(r2$RA, 2)))
          print(paste("L5: ", round(r2$L5, 2)))
          print(paste("L5_starttime: ", r2$L5_starttime))
          print(paste("M10: ", round(r2$M10, 2)))
          print(paste("M10_starttime: ", r2$M10_starttime))
          print("")
          print("Early-Warning Signals")
          print(paste("Mean: ", r2$Mean))
          print(paste("Variance: ", r2$Variance))
          print(paste("SD: ", r2$SD))
          print(paste("Coefficient of Variation: ", r2$CoV))
          print(paste("Skewness: ", r2$Skewness))
          print(paste("Kurtosis: ", r2$Kurtosis))
          print(paste("Autocorr at-lag-1: ", r2$Autocorr))
          print(paste("Autocorr at-lag-2: ", r2$Autocorr_lag2))
          print(paste("Autocorr at-lag-3: ", r2$Autocorr_lag3))
          print(paste("Autocorr at-lag-60: ", r2$Autocorr_lag60))
          print(paste("Autocorr at-lag-120: ", r2$Autocorr_lag120))
          print(paste("Time_to_Recovery: ", r2$Time_to_Recovery))
          print("---------------------------------------------------------------------------------")

        }
        rollingwindow.results # Needed for output in .CSV

        rollingwindow.results <- rollingwindow.results
      }
      ## Assign results from rolling window:
      rollingwindow.results <- rollingwindow(x = CRV.data, window = (1440 * (movingwindow.size)), jump = (1440 * (movingwindow.jump)))

      ## Initialise normal circadian rhythm analysis without moving window:
    } else {
      if (circadian_analysis) {
        ## Use the nparcalc{ACTman} function to calculate circadian rhythm variables over
        ## the whole period:

        r2 <- nparcalc(myACTdevice = myACTdevice, movingwindow = movingwindow, CRV.data = CRV.data, ACTdata.1.sub = ACTdata.1.sub)

        ## Attach r2 output to overview
        ACTdata.overview[i, "r2.IS"] <- r2$IS
        ACTdata.overview[i, "r2.IV"] <- r2$IV
        ACTdata.overview[i, "r2.RA"] <- round(r2$RA, 2)
        ACTdata.overview[i, "r2.L5"] <- r2$L5
        ACTdata.overview[i, "r2.L5_starttime"] <- r2$L5_starttime
        ACTdata.overview[i, "r2.M10"] <- round(r2$M10, 2)
        ACTdata.overview[i, "r2.M10_starttime"] <- r2$M10_starttime
      }
    }
    ## If a comparison with another actigraphy R package is required, run nparACT_base_loop{nparACT}:
    if (nparACT_compare) {

      ## Use nparACT Package to calculate circadian rhythm variables:
      r <- nparACT::nparACT_base_loop(path = newdir, SR = 1/60, fulldays = T, plot = T)

      ## Attach nparACT output to overview
      ACTdata.overview[i, "IS"] <- r$IS
      ACTdata.overview[i, "IV"] <- r$IV
      ACTdata.overview[i, "RA"] <- r$RA
      ACTdata.overview[i, "L5"] <- r$L5
      ACTdata.overview[i, "L5_starttime"] <- r$L5_starttime
      ACTdata.overview[i, "M10"] <- r$M10
      ACTdata.overview[i, "M10_starttime"] <- r$M10_starttime
    }

    ## Set wd back to main workdir
    setwd(workdir)


    ## Create "Results" directory for rolling window results:
    workdir_temp <- getwd()
    dir.create(file.path(workdir_temp, "Results"), showWarnings = FALSE)
    setwd(file.path(workdir_temp, "Results"))

    ## Write rollingwindow.results to .CSV
    # if (movingwindow) {
    #   write.table(rollingwindow.results,
    #               file = paste("rollingwindow-results", i, ".csv", sep = ""),
    #               row.names = F, sep = ",")
    # }

    if (movingwindow) {
      write.table(rollingwindow.results,
                  file = paste(substr(ACTdata.files[i], 1, (nchar(ACTdata.files[i]) - 4)),
                               "-rollingwindow-results.csv", sep = ""),
                  row.names = F, sep = ",")
    }


    paste(substr(ACTdata.files[i], 1, (nchar(ACTdata.files[i]) - 4)),
          "-sleepresults.csv", sep = "")

    ## Set working directory back to main working directory:
    setwd(workdir_temp)
    rm(workdir_temp)


    ## Sleep Analysis:
    ## Use the sleepdata_overview{ACTman} function to calculate sleep variables over
    ## the whole period.
    if (iwantsleepanalysis) {
      # sleepdata.overview <- sleepdata_overview(workdir = sleepdatadir, actdata = ACTdata.1.sub, i = i, lengthcheck = lengthcheck)
      sleepdata.overview <- sleepdata_overview(workdir = workdir, actdata = ACTdata.1.sub, i = i, lengthcheck = lengthcheck, ACTdata.files = ACTdata.files)
    }

    ## Actogram:
    ## Use the plot_actogram{ACTman} function to plot an Actogram of the whole period.
    if (plotactogram != FALSE) {
      plot_actogram(workdir = workdir, ACTdata.1.sub = ACTdata.1.sub, i = i, plotactogram = plotactogram,
                    rollingwindow.results = rollingwindow.results, i_want_EWS = i_want_EWS)
    }


    ## Report progress in console:
    print(paste("--------------------------------------", "END OF DATASET", i, "---", "@",
                round(i * (100 / length(ACTdata.files))), "% DONE",  "--------------------------------------"))
  }


  ## Step 3: After loop processing-----------------------------------------------------------------

  ## Transform negative recordingtimes to positive:
  ACTdata.overview$recordingtime <- ((ACTdata.overview$recordingtime) ^ 2) ^ (1 / 2)
  ACTdata.overview$recordingtime2 <- ((ACTdata.overview$recordingtime2) ^ 2) ^ (1 / 2)

  ## Assign zero to missings column in overview when there are no missings:
  ACTdata.overview[is.na(ACTdata.overview[, "missings"]), "missings"] <- 0

  ## Update overview if comparison to nparACt is not required:
  if (!nparACT_compare) {
    ACTdata.overview["IS"] <- NULL
    ACTdata.overview["IV"] <- NULL
    ACTdata.overview["RA"] <- NULL
    ACTdata.overview["L5"] <- NULL
    ACTdata.overview["L5_starttime"] <- NULL
    ACTdata.overview["M10"] <- NULL
    ACTdata.overview["M10_starttime"] <- NULL
    colnames(ACTdata.overview) <- gsub(pattern = "r2.", x = colnames(ACTdata.overview), replacement = "")
  }

  ## Subset experimental variables
  ACTdata.1.sub.expvars <- ACTdata.overview[c("IS", "IV", "RA", "L5", "L5_starttime", "M10", "M10_starttime")]
  # colnames(ACTdata.1.sub.expvars) <- c("IS", "IV", "RA", "L5", "L5 Start time", "M10", "M10 Start time")


  ## Create "Results" directory:
  workdir_temp <- getwd()
  dir.create(file.path(workdir_temp, "Results"), showWarnings = FALSE)
  setwd(file.path(workdir_temp, "Results"))

  ## Write results of circadian analysis to .CSV
  if (circadian_analysis) {
    write.table(ACTdata.1.sub.expvars, file = "ACTdata_circadian_res.csv", sep = ",", row.names = F)
  }

  ## Write ACTdata.overview to .CSV
  write.table(ACTdata.overview, file = "ACTdata_overview.csv", sep = ",", row.names = F)

  ## Set working directory back to main working directory:
  setwd(workdir_temp)
  rm(workdir_temp)


  ## Returned result.
  if (iwantsleepanalysis) {
    sleepdata.overview
  } else if (movingwindow) {
    rollingwindow.results
  } else {
    ACTdata.overview
  }
}

#' Data included in the package
#'
#' @name ACTdata.1
#' @docType data
#' @author Yoram K. Kunkels \email{y.k.kunkels@umcg.nl}
#' @references \url{exampledata.com}
#' @keywords data
NULL