R/time_since_event.R

In doBy: Groupwise Statistics, LSmeans, Linear Estimates, Utilities

.timeSinceEvent <- function(yvar, tvar=seq_along(yvar)){

  if (!(is.numeric(yvar) | is.logical(yvar))){
    stop("yvar must be either numeric or logical")
  }

  yvar[is.na(yvar)] <- 0

  event.idx <- which(yvar==1)
  if (length(event.idx)==0){
    return(NULL)
  }

  res <- vector("list", length(event.idx))
  for (kk in seq_along(event.idx)){
    jj <- event.idx[kk]
    res[[kk]] <- tvar-tvar[jj]
  }

  res      <-do.call(rbind, res)
  #print(res)
  abs.tse  <-apply(abs(res),2,min)
  #print(abs.tse)

  sgn 	<- rep.int(0, length(yvar))
  sss 	<- subSeq(yvar,item=0)
  for (ii in 1:nrow(sss)){
    idx  <- sss[ii,1]:sss[ii,2]
    idx1 <- seq_along(idx)
    uuu  <- abs.tse[idx]
    fff <- .find.inc.dec(uuu)
    sgn[idx] <- fff
  }



  #print(sgn)
  sign.tse <- abs.tse * sgn
  eventWindow <- rep(NA, length(abs.tse))

  #print(sign.tse)

  curr.state <- 1
  eventWindow[1] <- curr.state
  for (jj in 2:(length(eventWindow)-1)){
    if (sign.tse[jj] <= 0 & sign.tse[jj-1] >0 ){
      curr.state <- curr.state + 1
    }
    eventWindow[jj] <- curr.state
  }
  eventWindow[length(eventWindow)] <- curr.state
  #print(eventWindow)

  run <- cumsum(yvar)
  un <- unique(run)
  tlist <- list()
  for (ii in 1:length(un)){
    vv <- un[ii]
    yy <- yvar[run==vv]
    tt <- tvar[run==vv]
    tt <- tt - tt[1]
    tlist[[ii]] <- tt
  }
  timeAfterEvent <- unlist(tlist)
  timeAfterEvent[run==0] <- NA

  yvar2 <- rev(yvar)
  tvar2 <- rev(tvar)

  run2 <- cumsum(yvar2)
  un2 <- unique(run2)
  tlist2 <- list()
  for (ii in 1:length(un2)){
    vv <- un2[ii]
    yy <- yvar2[run2==vv]
    tt <- tvar2[run2==vv]
    tt <- tt - tt[1]
    tlist2[[ii]] <- tt
  }
  timeBeforeEvent <- unlist(tlist2)
  timeBeforeEvent[run2==0] <- NA

  timeBeforeEvent <- rev(timeBeforeEvent)
  run[run==0]<-NA

  #aux <- cbind(yvar,tvar, run, timeAfterEvent, timeBeforeEvent)
  aux <- cbind(run, tae=timeAfterEvent, tbe=timeBeforeEvent)
  ans <- cbind(data.frame(yvar=yvar, tvar=tvar), abs.tse, sign.tse, ewin=eventWindow, aux)
  return(ans)
}






.find.inc.dec <- function(ttt){
  ans <- rep.int(0,length(ttt))
  max.ttt <- max(ttt)
  wm <- which(ttt==max.ttt)
  if (max(wm)==length(ttt)){ ## The last element is maximal
    mmm <- max(wm)
    ans[] <- 1
  } else {
    if (min(wm)==1){ ## The first element is maximal
      mmm <- 1
      ans[] <- -1
    } else {
      mmm <- ceiling(mean(wm))
      ans[1:(mmm-1)] <- 1
      ans[mmm:length(ans)] <- -1
    }
  }
  ans
}
# ttt <- c(1,1,1,2,2,3,4,4,4,5,5,5,3,2,2,2,1,1)
# ttt <- c(1,1,1,2,2,3,4,4,4,5,5,4)
# ttt <- rev(ttt)



#' @title Calculate "time since event" in a vector.
#' 
#' @description Calculate "time since event" in a vector.
#'
#' @details Events are coded as 1 in numeric vector (and non-events are
#'     coded with values different from 1). \code{timeSinceEvent} will give the
#'     time since event (with and without sign). In a logical vector, events are
#'     coded as TRUE and all non-events as FALSE.
#' 
#' @param yvar A numerical or logical vector specifying the events
#' @param tvar An optional vector specifying time
#' @return A dataframe with columns 'yvar', 'tvar', 'abs.tse' (absolute time
#'     since nearest event), 'sign.tse' (signed time since nearest event) and
#'     'run' (indicator of the time window around each event).
#' @note NA's in yvar are converted to zeros.
#' @author Søren Højsgaard, \email{sorenh@@math.aau.dk}
#' @seealso \code{\link{subSeq}}, \code{\link{rle}}
#' @keywords utilities
#' @examples
#' 
#' ## Events:
#' yvar <- c(0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0,
#'           0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0)
#' 
#' ## Plot results:
#' tse<- timeSinceEvent(yvar)
#' plot(sign.tse~tvar, data=tse, type="b")
#' grid()
#' rug(tse$tvar[tse$yvar==1], col=4, lwd=4)
#' points(scale(tse$run), col=tse$run, lwd=2)
#' lines(abs.tse + .2 ~ tvar, data=tse, type="b", col=3)
#' 
#' ## Find times for which time since an event is at most 1:
#' tse$tvar[tse$abs<=1]
#' 
#' yvar <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 
#' 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
#' )
#' tvar <- c(207, 208, 208, 208, 209, 209, 209, 209, 210, 210, 211, 211, 
#' 211, 212, 213, 213, 214, 214, 215, 216, 216, 216, 216, 217, 217, 
#' 217, 218, 218, 219, 219, 219, 219, 220, 220, 221, 221, 221, 221, 
#' 222, 222, 222)
#' 
#' timeSinceEvent(yvar, tvar)
#' 
#' 
#' 
#' @export timeSinceEvent
timeSinceEvent <- function(yvar, tvar=seq_along(yvar)){

  if (!(is.numeric(yvar) | is.logical(yvar))){
    stop("yvar must be either numeric or logical")
  }

  yvar[is.na(yvar)] <- 0

  event.idx <- which(yvar==1)
  if (length(event.idx)==0){
    return(NULL)
  }
  n.event <- length(event.idx)

  ## find event times
  event.time <- tvar[event.idx]

  ## get time difference to each event
  rrr <-  do.call(rbind, lapply(event.idx, function(ii) tvar-tvar[ii]))
  abs.tse <- apply(abs(rrr),2,min)


  ## get the event windows (~ symmetrical around event time)
  ewin<-rep.int(NA, length(yvar))
  if (n.event>1){
    ff<-event.time[1:(n.event-1)]+diff(event.time)/2
    ewin[tvar<=ff[1]] <- 1
    for (ii in 2:(length(ff)-0)){
      ewin[tvar>ff[ii-1] & tvar<=ff[ii] ] <- ii
    }
    ewin[tvar>ff[length(ff)]] <- n.event
  } else {
    ewin[] <- n.event
  }

  ## get the signs
  ggg <- list()
  for (ii in 1:(length(event.idx))){
    ggg[[ii]] <- rrr[ii,ewin==ii]
  }
  ggg <- unlist(ggg)
  sign.tse <- sign(ggg)*abs.tse

  run <- cumsum(yvar)

  un <- unique(run)
  tlist <- list()
  for (ii in 1:length(un)){
    vv <- un[ii]
    yy <- yvar[run==vv]
    tt <- tvar[run==vv]
    tt <- tt - tt[1]
    tlist[[ii]] <- tt
  }
  tae <- unlist(tlist)
  tae[run==0] <- NA

  yvar2 <- rev(yvar)
  tvar2 <- rev(tvar)

  run2 <- cumsum(yvar2)
  un2 <- unique(run2)
  tlist2 <- list()
  for (ii in 1:length(un2)){
    vv <- un2[ii]
    yy <- yvar2[run2==vv]
    tt <- tvar2[run2==vv]
    tt <- tt - tt[1]
    tlist2[[ii]] <- tt
  }
  tbe <- unlist(tlist2)
  tbe[run2==0] <- NA

  tbe <- rev(tbe)
  run[run==0]<-NA

  ans <- cbind(data.frame(yvar=yvar, tvar=tvar), abs.tse, sign.tse, ewin=ewin,
               run, tae=tae, tbe=tbe)
  ans
}