R/Wprime.R

In trackeR: Infrastructure for Running, Cycling and Swimming Data from GPS-Enabled Tracking Devices

#' W' expended.
#'
#' Calculate W' expended, i.e., the work capacity above critical
#' power/speed which has been depleted and not yet been replenished.
#'
#' @param object Univariate \code{\link[zoo]{zoo}} object containing
#'     the time stamped power output or speed values. (Power should be
#'     in Watts, speed in meters per second.)
#' @param w0 Inital capacity of W', as calculated based on the
#'     critical power model by Monod and Scherrer (1965).
#' @param cp Critical power/speed, i.e., the power/speed which can be
#'     maintained for longer period of time.
#' @param version How should W' be replenished? Options include
#'     \code{'2015'} and \code{'2012'} for the versions presented in
#'     Skiba et al. (2015) and Skiba et al. (2012), respectively. See
#'     Details.
#' @param meanRecoveryPower Should the mean of all power outputs below
#'     critical power be used as recovery power? See Details.
#'
#' @details
#'
#' Skiba et al. (2015) and Skiba et al. (2012) both describe an
#' exponential decay of \eqn{W'} expended over an interval
#' \eqn{[t_{i-1}, t_i)} if the power output during this interval is
#' below critical power:
#'
#' \deqn{W_exp (t_i) = W_exp(t_{i-1}) * exp(nu * (t_i - t_{i-1}))}
#'
#' However, the factor nu differs: Skiba et al. (2012) describe it as
#' \eqn{1/\tau} with \eqn{\tau} estimated as
#'
#' \deqn{tau = 546 * exp(-0.01 * (CP - P_i)) + 316}
#'
#' Skiba et al. (2015) use \eqn{(P_i - CP) / W'_0}.  Skiba et
#' al. (2012) and Skiba et al. (2015) employ a constant recovery power
#' (calculated as the mean over all power outputs below critical
#' power). This rationale can be applied by setting the argument
#' \code{meanRecoveryPower} to \code{TRUE}. Note that this uses
#' information from all observations with a power output below
#' critical power, not just those prior to the current time point.
#'
#' @references
#'
#' Monod H, Scherrer J (1965). 'The Work Capacity of a Synergic
#' Muscular Group.' Ergonomics, 8(3), 329--338.
#'
#' Skiba PF, Chidnok W, Vanhatalo A, Jones AM (2012). 'Modeling the
#' Expenditure and Reconstitution of Work Capacity above Critical
#' Power.' Medicine & Science in Sports & Exercise, 44(8), 1526--1532.
#'
#' Skiba PF, Fulford J, Clarke DC, Vanhatalo A, Jones AM
#' (2015). 'Intramuscular Determinants of the Abilility to Recover
#' Work Capacity above Critical Power.' European Journal of Applied
#' Physiology, 115(4), 703--713.
#'
#' @export
Wexp <- function(object, w0, cp, version = c("2015", "2012"), meanRecoveryPower = FALSE) {
    version <- match.arg(version, c("2015", "2012"))

    time <- index(object)
    power <- coredata(object)

    dp <- power - cp
    delta <- as.numeric(difftime(c(time[-1], time[length(time)]), time, units = "secs"))
    ind <- dp > 0

    wna <- is.na(power)
    ## if all power values are NA, return NA
    if (all(wna)) {
        ret <- cbind(coredata(object), power)
        colnames(ret) <- c("movement", "wprime")
        ret <- zoo(ret, order.by = index(object))
        return(ret)
    }
    ## otherwise carry on
    time <- time[!wna]
    power <- power[!wna]
    dp <- dp[!wna]
    delta <- delta[!wna]
    ind <- ind[!wna]

    ## Done in papers (Skiba et al (2012) and Skiba et al (2015)) but not necessary with
    ## assumption of constant power over intervals between observations: use mean power over
    ## all instances with power below CP as recovery power (here weighted by length of time
    ## interval)
    if (meanRecoveryPower) {
        recovery <- power <= cp
        rPower <- power[recovery]
        rDelta <- delta[recovery]
        rDeltaSum <- sum(rDelta)
        dp[recovery] <- sum(rPower * rDelta/sum(rDeltaSum)) - cp
    }

    ## recovery factor
    if (version == "2015") {
        recoveryFactor <- dp/w0
    }
    else {
        ## version 2012
        tau <- 546 * exp(0.01 * dp) + 316
        recoveryFactor <- -1/tau
    }

    ## W'expended is W' previously depleted and not yet recovered
    Wexp <- rep(NA, length(ind))
    Wexp[1] <- if (ind[1])
        dp[1] * delta[1] else 0
    for (i in 2:length(ind)) {
        if (ind[i]) {
            Wexp[i] <- Wexp[i - 1] + dp[i] * delta[i]
        } else {
            Wexp[i] <- Wexp[i - 1] * exp(recoveryFactor[i] * delta[i])
        }
    }

    ret <- rep(NA, length(object))
    ret[!wna] <- Wexp
    ret <- cbind(coredata(object), ret)
    colnames(ret) <- c("movement", "wprime")
    ret <- zoo(ret, order.by = index(object))

    return(ret)
}

#' W': work capacity above critical power/speed.
#'
#' @aliases trackeRWprime
#'
#' Based on the critical power model for cycling (Monod and Scherrer,
#' 1965), W' (read W prime) describes the finite work capacity above
#' critical power (Skiba et al., 2012).  While W' is depleted during
#' exercise above critical power, it is replenished during exercise
#' below critical power. Thus, it is of interest how much of this work
#' capacity has been depleted and not yet been replinished again,
#' named W' expended, or how much of this work capacity is still
#' available, named W' balance.  This principal is applied to runners
#' by subsituting power and critical power with speed and critical
#' speed, respectively (Skiba et al., 2012).
#'
#' @param object A \code{\link{trackeRdata}} object.
#' @param session A numeric vector of the sessions to be used,
#'     defaults to all sessions.
#' @param quantity Should W' \code{'expended'} or W' \code{'balance'}
#'     be returned?
#' @inheritParams Wexp
#' @param parallel Logical. Should computation be carried out in
#'     parallel? If \code{TRUE} computation is performed in parallel
#'     using the backend provided to \code{\link{foreach}}. Default is
#'     \code{FALSE}.
#' @param ... Currently not used.
#'
#' @return An object of class \code{trackeRWprime}.
#'
#' @details
#'
#' #' Skiba et al. (2015) and Skiba et al. (2012) both describe an
#' exponential decay of \eqn{W'} expended over an interval
#' \eqn{[t_{i-1}, t_i)} if the power output during this interval is
#' below critical power:
#'
#' \deqn{W_exp (t_i) = W_exp(t_{i-1}) * exp(nu * (t_i - t_{i-1}))}
#'
#' However, the factor nu differs: Skiba et al. (2012) describe it as
#' \eqn{1/\tau} with \eqn{\tau} estimated as
#'
#' \deqn{tau = 546 * exp(-0.01 * (CP - P_i)) + 316}
#'
#' Skiba et al. (2015) use \eqn{(P_i - CP) / W'_0}.  Skiba et
#' al. (2012) and Skiba et al. (2015) employ a constant recovery power
#' (calculated as the mean over all power outputs below critical
#' power). This rationale can be applied by setting the argument
#' \code{meanRecoveryPower} to \code{TRUE}. Note that this uses
#' information from all observations with a power output below
#' critical power, not just those prior to the current time point.
#'
#' @references
#'
#' Monod H, Scherrer J (1965). 'The Work Capacity of a Synergic
#' Muscular Group.' Ergonomics, 8(3), 329--338.
#'
#' Skiba PF, Chidnok W, Vanhatalo A, Jones AM (2012). 'Modeling the
#' Expenditure and Reconstitution of Work Capacity above Critical
#' Power.' Medicine & Science in Sports & Exercise, 44(8), 1526--1532.
#'
#' Skiba PF, Fulford J, Clarke DC, Vanhatalo A, Jones AM
#' (2015). 'Intramuscular Determinants of the Abilility to Recover
#' Work Capacity above Critical Power.' European Journal of Applied
#' Physiology, 115(4), 703--713.
#'
#' @export
#' @examples
#' \dontrun{
#' data('runs', package = 'trackeR')
#' wexp <- Wprime(runs, session = c(11,13), cp = 4, version = '2012')
#' plot(wexp)
#' }
Wprime <- function(object, session = NULL, quantity = c("expended", "balance"), w0, cp,
    version = c("2015", "2012"), meanRecoveryPower = FALSE, parallel = FALSE,
    ...) {
    ## prep args
    quantity <- match.arg(quantity, c("expended", "balance"))
    if (quantity == "balance" & missing(w0))
        stop("Please provide w0 or choose quantity = 'expended'.")
    version <- match.arg(version, c("2015", "2012"))
    if (version == "2015" & missing(w0))
        stop("Please provide w0 or choose version = '2012'.")
    if (missing(w0))
        w0 <- NA

    ## units
    units <- getUnits(object)

    ## select sessions
    if (is.null(session))
        session <- 1:length(object)
    object <- object[session]
    sports <- get_sport(object)

    ## FIXME: change here to accommodate multisport environment for
    ## now, do not allow computation of Wprime if user does not
    ## understand what they are trying to do
    if (sum(c("running", "cycling") %in% na.omit(unique(sports))) != 1) {
        stop("Wprime applies only for running or only for cycling sessions")
    }

    cycling <- all(sports == "cycling")
    ps <- ifelse(cycling, "power", "speed")
    if (cycling) {
        if (units$unit[units$variable == "power" & units$sport == "cycling"] != "W") {
            object <- change_units(object, variable = "power", unit = "W", sport = unique(sports))
            units <- getUnits(object)
            conversion <- match.fun(paste(units$unit[units$variable == "power" & units$sport == "cycling"], "W", sep = "2"))
            cp <- conversion(cp)
        }
    }
    else {
        if (units$unit[units$variable == "speed" & units$sport == "running"] != "m_per_s") {
            object <- change_units(object, variable = "speed", unit = "m_per_s", sport = unique(sports))
            units <- getUnits(object)
            conversion <- match.fun(paste(units$unit[units$variable == "speed" & units$sport == "running"], "m_per_s",
                sep = "2"))
            cp <- conversion(cp)
        }
    }


    ## get W'
    getW <- function(j) {
        sess <- object[[j]]
        W <- Wexp(sess[, ps], w0 = w0, cp = cp, version = version, meanRecoveryPower = meanRecoveryPower)
        if (quantity == "balance")
            W$wprime <- w0 - W$wprime
        return(W)
    }

    foreach_object <- eval(as.call(c(list(quote(foreach::foreach),
                                          j = seq.int(nsessions(object))))))

    if (parallel) {
        setup_parallel()
        ret <- foreach::`%dopar%`(foreach_object, getW(j))
    }
    else {
        ret <- foreach::`%do%`(foreach_object, getW(j))
    }

    ## class and return
    attr(ret, "quantity") <- quantity
    attr(ret, "w0") <- if (missing(w0))
        NA else w0
    attr(ret, "cp") <- cp

    ## FIXME: fixme here to accommodate for multi-sport environment
    attr(ret, "cycling") <- cycling
    attr(ret, "sport") <- sports
    attr(ret, "unit") <- units[units$variable == ps & units$sport == unique(sports), ]
    class(ret) <- "trackeRWprime"
    return(ret)
}

#' Plot W'.
#'
#' @param x An object of class \code{trackeRWprime} as returned by \code{\link{Wprime}}.
#' @param session A numeric vector of the sessions to be plotted, defaults to all sessions.
#' @param dates Logical. Should the date of the session be used in the panel header?
#' @param scaled Logical. Should the W' be scaled to the movement variable (power or speed)
#'     which is then plotted in the background?
#' @param ... Currently not used.
#' @export
#' @examples
#' \dontrun{
#' data('runs', package = 'trackeR')
#' wexp <- Wprime(runs, session = 1:3, cp = 4, version = '2012')
#' plot(wexp, session = 1:2)
#' }
plot.trackeRWprime <- function(x, session = NULL, dates = TRUE, scaled = TRUE, ...) {
    quantity <- attr(x, "quantity")
    sports <- attr(x, "sport")
    cp <- attr(x, "cp")
    cycling <- attr(x, "cycling")
    Wunit <- if (cycling) "[W]" else "[m]"
    mylabels <- c(paste0(ifelse(cycling, "Power", "Speed"),
                         " [",
                         prettifyUnits(attr(x, "unit")$unit), "]"),
                  paste("W'", quantity, "[scaled]"))

    ## select sessions
    if (is.null(session))
        session <- seq_along(x)

    x <- x[session]

    ## transform W' to match power/speed scale
    if (scaled) {
        sdMov <- stats::sd(unlist(lapply(x, function(z) z$movement)), na.rm = TRUE)
        mMov <- mean(unlist(lapply(x, function(z) z$movement)), na.rm = TRUE)

        x <- lapply(x, function(z) {
            w <- (coredata(z$wprime) - mean(coredata(z$wprime), na.rm = TRUE))/stats::sd(coredata(z$wprime),
                na.rm = TRUE)
            w <- w * sdMov  #sd(coredata(z$movement), na.rm = TRUE)
            z$wprime <- w + mMov
            return(z)
        })
    }

    ## class and get data
    attr(x, "sport") <- sports[session]
    class(x) <- "trackeRWprime"
    df <- fortify(x, melt = TRUE)



    ## prepare session id for panel header
    if (dates) {
        df$SessionID <- format(session[df$SessionID])
        df$SessionID <- gsub(" ", "0", df$SessionID)
        df$SessionID <- paste0(paste(df$SessionID, df$Sport, sep = ": "), "\n", format(df$Index, "%Y-%m-%d"))
    }
    else {
        df$SessionID <- paste0(paste(df$SessionID, df$Sport, sep = ": "))
        ## factor(df$SessionID, levels = seq_along(session), labels = session)
    }
    df$Series <- factor(df$Series)

    ## check that there is data to plot
    for (l in levels(df$Series)) {
        if (all(is.na(subset(df, Series == l, select = "Value"))))
            df <- df[!(df$Series == l), ]
    }

    ## make facets
    singleSession <- length(session) == 1L
    facets <- if (singleSession)
        NULL else ". ~ SessionID"

    if (scaled) {
        ## basic plot
        p <- ggplot(data = df, mapping = aes_(x = quote(Index), y = quote(Value))) +
            ylab("") + xlab("Time")
        ## lines for power/speed and W'
        p <- p + geom_line(aes_(group = quote(Series), col = quote(Series)),
            na.rm = TRUE) + scale_colour_manual(name = "", labels = mylabels,
            values = c("gray", "black"))
        ## add line for cp
        p <- p + geom_hline(data = data.frame(cp = cp), aes(yintercept = cp),
            col = gray(0.25))
    }
    else {
        ## basic plot
        p <- ggplot(data = subset(df, Series == "wprime"), mapping = aes_(x = quote(Index),
            y = quote(Value))) + ylab(paste("W'", quantity, Wunit)) + xlab("Time")
        ## lines for W'
        p <- p + geom_line(na.rm = TRUE)
    }
    ## add facet if necessary
    if (!is.null(facets)) {
        p <- p + facet_grid(facets, scales = "free")
    }
    ## add bw theme
    p <- p + theme_bw() + theme(legend.position = "top", axis.text.x = element_text(angle = 50,
        hjust = 1))

    return(p)
}


#' Fortify a trackeRWprime object for plotting with ggplot2.
#'
#' @param model The \code{trackeRWprime} object as returned by \code{\link{Wprime}}.
#' @param data Ignored.
#' @param melt Logical. Should the data be melted into long format
#'     instead of the default wide format?
#' @param ... Ignored.
#' @export
fortify.trackeRWprime <- function(model, data, melt = FALSE, ...) {
    ret <- list()
    sports <- get_sport(model)
    for (i in seq_along(model)) {
        ret[[i]] <- zoo::fortify.zoo(model[[i]], melt = melt)
        ret[[i]]$SessionID <- i
        ret[[i]]$Sport <- sports[i]
    }
    ret <- do.call("rbind", ret)
    return(ret)
}

#' @rdname nsessions
#' @export
nsessions.trackeRWprime <- function(object, ...) {
    length(object)
}


#' @rdname get_sport
#' @export
get_sport.trackeRWprime <- function(object, ...) {
    attr(object, "sport")
}


#' Change the units of the variables in an \code{\link{trackeRWprime}} object
#'
#' @param object An object of class \code{\link{trackeRWprime}}.
#' @param variable A vector of variables to be changed.
#' @param unit A vector with the units, corresponding to variable.
#' @param ... Currently not used.
#' @export
change_units.trackeRWprime <- function(object,
                                       variable,
                                       unit,
                                       ...) {
    ## get current unit
    current <- getUnits(object)
    if (missing(variable))
        variable <- ifelse(attr(object, "cycling"), "power", "speed")
    if (missing(unit) & !missing(variable)) {
        unit <- variable
        variable <- ifelse(attr(object, "cycling"), "power", "speed")
    }
    if (attr(object, "cycling")) {
        if (variable != "power")
            stop("can only change measurement units for power.")
    } else {
        if (variable != "speed")
            stop("can only change measurement units for speed.")
    }

    ## change units
    for (i in variable) {
        currentUnit <- current$unit[current$variable == i]
        newUnit <- unit[which(variable == i)]
        if (currentUnit != newUnit) {
            conversion <- match.fun(paste(currentUnit, newUnit, sep = "2"))
            ## change data
            for (session in seq_along(object)) {
                object[[session]][, "movement"] <- conversion(object[[session]][, "movement"])
            }
            ## change units attribute
            current$unit[current$variable == i] <- newUnit
        }
    }

    ## update attributes and return
    attr(object, "units") <- current
    return(object)
}