#############################################################################
### nparcalc{ACTman} ###
### Script authors: Yoram Kunkels, Stefan Knapen, & Ando Emerencia ###
### Most recent Update: 27-02-2018 ###
###~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~###
#' nparcalc
#'
#' Calculate non-parametrical circadian rhythm variables:
#' IS, IV, RA, L5, L5_starttime, M10, and M10_starttime.
#'
#' @param myACTdevice Name of the input device used. Should be either Actiwatch2 or MW8.
#' @param movingwindow A boolean indicating whether moving window is used.
#' @param CRV.data CRV data
#' @param ACTdata.1.sub Managed data set
#' @param out Optional. When movingwindow is TRUE, this is the current window of data.
#'
#' @return A list with the result values IS, IV, RA, L5, L5_starttime, M10, and M10_starttime.
#'
#' @importFrom stats na.pass
#' @importFrom utils tail
#' @importFrom stats aggregate
#' @importFrom moments skewness
#' @importFrom moments kurtosis
nparcalc <- function(myACTdevice, movingwindow, CRV.data, ACTdata.1.sub, out = NULL) {
## Step 1: Basic Operations-------------------------------------------------------------------
## Initialise results list:
result <- list()
## Read Data, and combine seperate date and time variables when needed:
## If there are more than 2 columns, combine date and time variables, and remove extra column.
if (ncol(CRV.data) > 2) {
CRV.data$Date <- paste(CRV.data$Date, " ", CRV.data$Time)
CRV.data <- CRV.data[, -2]
## Else, just assign correct column names.
} else {
colnames(CRV.data) <- c("Date", "Activity")
}
## Prune data untill only full 24h days are obtained
CRV.data.wholehours <- CRV.data[grep("00:00", CRV.data[, "Date"]), ]
CRV.data.start <- which(CRV.data$Date == CRV.data.wholehours[1, "Date"])
## Device- and functionality specific identification of dataset end:
#! Note: There is repetition in this loop > can be written more efficiently!
## If device is MW8, then take dataset end at last "00:00:00" of dataset.
if (myACTdevice == "MW8") {
CRV.data.end <- tail(grep("00:00:00", ACTdata.1.sub$Date), 2)[1]
## Else, if moving window is required take dataset end at last "00:00:00" of "out" dataset.
} else {
if (movingwindow) {
CRV.data.end <- which(out == "00:00:00")[length(which(out == "00:00:00"))]
## Else, assume device to be Actiwatch2 and take dataset end at last "00:00:00" of dataset.
} else {CRV.data.end <- tail(grep("00:00:00", ACTdata.1.sub$Date), 2)[1]}
}
## Assign dataset for circadian rhythm analysis to be dataset "CRV.data" with a period of the
## aforementioned start and end date/times:
CRV.data <- CRV.data[CRV.data.start:CRV.data.end, ]
## Step 2: Calculate Circadian Rhythm Variables-------------------------------------------------
## IS: Interdaily stability
## Source: Van Someren, E., Swaab, D., Colenda, C., Cohen, W., McCall, W. and Rosenquist, P. (1999).
## Bright Light Therapy: Improved Sensitivity to Its Effects on Rest-Activity Rhythms in Alzheimer Patients
## by Application of Nonparametric Methods. Chronobiology International, 16(4), pp.505-518.
## Calculate Xi (which represent the individual datapoints / consecutive hourly means)
## by aggregation with hourly breaks:
xi <- aggregate(CRV.data[, "Activity"],
list(hour = cut(as.POSIXct(CRV.data[, "Date"]), breaks = "hour")),
mean, na.action = na.pass, na.rm = TRUE)
## Exception for moving window to retain last observation for correct number of observtions:
#! Wat mij niet helemaal duidelijk was is waarom je het laatste meetpunt in xi weggooit als er geen moving window is?
if (!movingwindow) {
xi <- xi[1:(nrow(xi) - 1), ]
} else {
xi <- xi[1:(nrow(xi)), ]
}
## Assign only activity data to xi:
xi <- xi$x
## Calculate X (or X_bar; the mean of all the data):
X <- mean(xi, na.rm = T)
## Calculate the Denominator part of the IS formula (van Someren et al., 1999):
xi_X <- xi - X # Differences between consecutive hourly means and overall mean.
sq.xi_X <- xi_X ^ 2 # Square of differences.
sum.sq.xi_X <- sum(sq.xi_X, na.rm = T) # Sum of squares.
n <- sum(!is.na(xi)) # Get number of hours (should be 168 for 7 day intervals (7*24)).
sum.sq.xi_X.perhour <- sum.sq.xi_X / n # Sum of squares per hour
## Calculate the Nominator part of the IS formula (van Someren et al., 1999):
## Workaround to assure xi length is multiple of the number of rows [24]
xi_sub <- xi[1:(24 * floor(length(xi) / 24))]
Xh <- rowMeans(matrix(xi_sub, nrow = 24), na.rm = T) # Calculate hourly means.
Xh_X <- Xh - X # difference 24 hour means and overall mean.
sum.sq.Xh_X <- sum(Xh_X ^ 2, na.rm = T) # sum of squares.
sum.sq.Xh_X.perhour <- sum.sq.Xh_X / 24 # sum of squares per hour.
IS <- sum.sq.Xh_X.perhour / sum.sq.xi_X.perhour # hourly variance divided by overall variance.
IS <- round(x = IS, digits = 2) # rounding down to 2 decimals.
result$IS <- IS # Write IS to results
##---------------------------------------------------------------------------------------------------------
## IV: Interdaily Variability
## Source: Van Someren, E., Swaab, D., Colenda, C., Cohen, W., McCall, W. and Rosenquist, P. (1999).
## Bright Light Therapy: Improved Sensitivity to Its Effects on Rest-Activity Rhythms in Alzheimer Patients
## by Application of Nonparametric Methods. Chronobiology International, 16(4), pp.505-518.
Xi_diffXi <- diff(xi) # difference Xi and previous Xi (the difference between all successive hours) (!!! Warning because no previous hour for 1st obs.)
sum.sq.Xi_diffXi <- sum(Xi_diffXi ^ 2, na.rm = T) # sum of squares
sum.sq.Xi_diffXi.perhour <- sum.sq.Xi_diffXi / (n - 1) # sum of squares per hour
Xi_X <- xi - X # difference Xi and overall mean
sum.sq.Xi_X <- sum(Xi_X ^ 2, na.rm = T) # sum of squares
# sum.sq.Xi_X.perhour <- sum.sq.Xi_X / (n - 1) #! Delen door n niet (n - 1)!
sum.sq.Xi_X.perhour <- sum.sq.Xi_X / n # sum of squares per hour minus 1
IV <- sum.sq.Xi_diffXi.perhour / sum.sq.Xi_X.perhour
IV <- round(x = IV, digits = 2) # rounding down to 2 decimals.
result$IV <- IV # Write IS to results
##---------------------------------------------------------------------------------------------------------
## L5: Average of the 5 Lowest Hourly Means (within each day(!))
## Source: Witting, W., Kwa, I., Eikelenboom, P., Mirmiran, M., & Swaab, D. (1990).
## Alterations in the circadian rest-activity rhythm in aging and Alzheimer's disease.
## Biological Psychiatry, 27(6), 563-572. http://dx.doi.org/10.1016/0006-3223(90)90523-5
## Calculate mean activity per minute over the whole period.
## The idea is to calculate an average day to calculate L5 and M10 over.
averageday <- matrix(c(substr(CRV.data[1:1440, "Date"], 14, 22), rep(NA, 1440)), nrow = 1440, ncol = 2)
selection_mat <- matrix(FALSE, 1440)
for (aa in 1:1440) {
selection_mat[aa, ] <- TRUE
mean(CRV.data[selection_mat, "Activity"])
averageday[aa, 2] <- mean(CRV.data[selection_mat, "Activity"], na.rm = TRUE) # added "na.rm = TRUE" for moving window
selection_mat <- matrix(FALSE, 1440)
}
averageday[, 2] <- as.numeric(averageday[, 2])
#! DEBUG to allow for window to cross midnight point
averageday <- rbind(averageday, averageday)
averageday_loc_L5 <- matrix(NA, nrow(averageday))
## Loop for calculating mean activity in 300 minute intervals ~ 5 hours
for (gg in 1:(nrow(averageday) - 300)) {
averageday_loc_L5[gg] <- mean(as.numeric(averageday[c(gg:(299 + gg)), 2]))
}
L5_starttime <- averageday[which.min(averageday_loc_L5), 1]
L5 <- round(averageday_loc_L5[which.min(averageday_loc_L5)], 2)
result$L5 <- L5
result$L5_starttime <- L5_starttime
##---------------------------------------------------------------------------------------------------------
## M10: Average of the 10 Highest Hourly Means (within each day(!))
## Source: Witting, W., Kwa, I., Eikelenboom, P., Mirmiran, M., & Swaab, D. (1990).
## Alterations in the circadian rest-activity rhythm in aging and Alzheimer's disease.
## Biological Psychiatry, 27(6), 563-572. http://dx.doi.org/10.1016/0006-3223(90)90523-5
## Calculate mean activity per minute over the whole period.
## The idea is to calculate an average day to calculate L5 and M10 over.
averageday_loc_M10 <- matrix(NA, nrow(averageday))
## Loop for calculating mean activity in 600 minute intervals ~ 10 hours
for (hh in 1:(nrow(averageday) - 600)) {
averageday_loc_M10[hh] <- mean(as.numeric(averageday[c(hh:(599 + hh)), 2]), na.rm = TRUE)
}
M10_starttime <- averageday[which.max(averageday_loc_M10), 1]
M10 <- round(averageday_loc_M10[which.max(averageday_loc_M10)], 2)
result$M10 <- M10
result$M10_starttime <- M10_starttime
##---------------------------------------------------------------------------------------------------------
## RA: Relative Amplitude
## Source: Van Someren, E., Swaab, D., Colenda, C., Cohen, W., McCall, W. and Rosenquist, P. (1999).
## Bright Light Therapy: Improved Sensitivity to Its Effects on Rest-Activity Rhythms in Alzheimer Patients
## by Application of Nonparametric Methods. Chronobiology International, 16(4), pp.505-518.
## Calculate Amplitude by subtracting M10 from L5
Amp <- M10 - L5
## Calculate Relative Amplitude by dividing Amplitude by the sum of L5 and M10, and assign to result
RA <- Amp / (L5 + M10)
result$RA <- RA
##---------------------------------------------------------------------------------------------------------
## EWS Part
## Mean
Mean <- round(mean(CRV.data[, "Activity"], na.rm = TRUE), 2)
result$Mean <- Mean
## Variance
Variance <- round(var(CRV.data[, "Activity"], na.rm = TRUE), 2)
result$Variance <- Variance
## Standard Deviation
SD <- round(sd(CRV.data[, "Activity"], na.rm = TRUE), 2)
result$SD <- SD
## Coefficient of Variation
CoV <- round(((sd(CRV.data[, "Activity"], na.rm = TRUE) / mean(CRV.data[, "Activity"], na.rm = TRUE)) * 100), 2)
result$CoV <- CoV
## Skewness
Skewness <- moments::skewness(CRV.data[, "Activity"], na.rm = TRUE)
result$Skewness <- round(Skewness, 2)
## Kurtosis
Kurtosis <- moments::kurtosis(CRV.data[, "Activity"], na.rm = TRUE)
result$Kurtosis <- round(Kurtosis, 2)
## Autocorrelation at lag-1, lag-2, lag-3, etc.
Autocorr <- acf(x = CRV.data[, "Activity"], lag.max = 120, na.action = na.pass, plot = FALSE)
result$Autocorr <- round(Autocorr$acf[2], 2)
result$Autocorr_lag2 <- round(Autocorr$acf[3], 2)
result$Autocorr_lag3 <- round(Autocorr$acf[4], 2)
result$Autocorr_lag60 <- round(Autocorr$acf[61], 2)
result$Autocorr_lag120 <- round(Autocorr$acf[121], 2)
## Time to recovery (if not recovered within timeframe assign max_lag))
if (length(which(Autocorr$acf < 0.2)) != 0) {
result$Time_to_Recovery <- which(Autocorr$acf < 0.2)[1]
} else {
result$Time_to_Recovery <- 120
}
## Assign data to results
result$CRV_data <- CRV.data
## Return the list with all result values.
result
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.