R/pava.R

In targeted: Targeted Inference

#' @description Pooled Adjacent Violators Algorithm
#' @title Pooled Adjacent Violators Algorithm
#' @param y response variable
#' @param x (optional) predictor vector (otherwise y is assumed
#' to be a priori sorted according to relevant predictor)
#' @param weights weights (optional) weights
#' @return List with index (idx) of jump points and values (value)
#' at each jump point.
#' @export
#' @author Klaus K. Holst
#' @aliases pava isoreg isoregw
#' @examples
#' x <- runif(5e3, -5, 5)
#' pr <- lava::expit(-1 + x)
#' y <- rbinom(length(pr), 1, pr)
#' pv <- pava(y, x)
#' plot(pr ~ x, cex=0.3)
#' with(pv, lines(sort(x)[index], value, col="red", type="s"))
pava <- function(y, x=numeric(0), weights=numeric(0)) {
    if (length(x)) {
        ord <- order(x)
        y <- y[ord]
        if (length(weights)==length(y)) weights <- weights[ord]
    }
    if (length(weights)!=length(y)) weights <- numeric(0)
    val <- .pava(y, x, weights)
    return(val)
}

#' @export
isoregw <- function(x, y, weights=NULL, ...) {
  if (NCOL(x) > 1) stop("Expect only one predictor variable.")

  x <- as.vector(x) # in case x is a n x 1 matrix
  ord <- order(x)
  if (is.null(weights)) {
      pv <- pava(y[ord])
  } else {
      ## replication weights
      pv <- try(pava(y[ord], x=numeric(0), weights[ord]), silent=TRUE)
      if (inherits(pv, "try-error")) return(base::identity)
  }
  suppressWarnings(f <- tryCatch(approxfun(x[ord[pv$index]], pv$value,
                                            rule=2, method="constant"),
                                  error=function(...) base::identity))
  return(f)
}