R/plotLZISD.R

In ppernot/ErrViewLib: Library for ErrView

#' Estimate statistics of the inverse of SD of a sample
#'
#' @param Z (vector) a data sample
#' @param method (string) one of 'bootstrap' (default) and 'cho'
#' @param CImethod (string) one of 'bca' (default), 'perc' and 'basic'
#'   for the CI estimation algorithm from bootstrap sample
#' @param nBoot (integer) number of bootstrap repeats
#' @param level (numeric) a probability level for CI (default: 0.95)
#' @param parallel (string) one of 'no' (default) and 'multicore'
#'
#' @return A list containing the mean, sd and ci for the SD^-1
#'     of the Z sample
#' @export
#'
rsdZCI = function (
  Z,
  method = c('bootstrap','cho'),
  CImethod = c('bca','perc','basic'),
  nBoot = 1500,
  level = 0.95,
  parallel = c('no','multicore')
) {

  method   = match.arg(method)
  CImethod = match.arg(CImethod)
  parallel = match.arg(parallel)

  wtd.sd = function(x,weights=NULL,normwt=TRUE)
    1 / sqrt(Hmisc::wtd.var(x,weights=weights,normwt=normwt))

  switch (
    method,
    bootstrap = {
      bst = boot::boot(Z, wtd.sd, R = nBoot, stype = 'f',
                       parallel = parallel, normwt = TRUE)
      bci = boot::boot.ci(bst, conf = level, type = CImethod )
      ci = switch(
        CImethod,
        bca   = bci$bca[1, 4:5],
        perc  = bci$perc[1, 4:5],
        basic = bci$basic[1, 4:5]
      )
      list(
        mean   = bci$t0,
        sd     = sd(bst$t),
        ci     = ci,
        method = method,
        CImethod = CImethod,
        level  = level
      )
    },
    cho = {
      # LUP formula for y = 1/sqrt(Var(X))
      V  = var(Z)
      S  = 1 / sqrt(V)
      SD = sqrt(ErrViewLib::varvar(Z)) / (2*V^(3/2))
      list(
        mean   = S,
        sd     = SD,
        ci     = S + qnorm((1 + level) / 2) * c(-1, 1) * SD,
        method = method,
        level  = level
      )
    }
  )
}
#' Plot local inverse of z-score standard deviation 1/SD(Z)
#'
#' @param X (vector) abscissae of the Z values
#' @param Z (vector) set of z-score values to be tested
#' @param aux (vector) auxilliary vector to resolve ties in X sorting
#' @param sdZ (numeric) target value for 1/SD(Z) (default `1`)
#' @param logX (logical) log-transform X
#' @param method (string) method used to estimate 95 percent CI on 1/SD(Z)
#' @param BSmethod (string) bootstrap variant
#' @param nBoot (integer) number of bootstrap replicas
#' @param intrv (object) intervals generated by `genIntervals` (default: `NULL`)
#' @param nBin (integer) number of intervals for local coverage stats
#' @param slide (logical) use sliding window for subsetting (X,Z)
#' @param equiPop (logical) generate intervals  with equal bin counts
#'   (default: `TRUE`)
#' @param popMin (integer) minimal bin count in an interval
#' @param logBin (logical) if `equiPop = FALSE`, one can choose between
#'   equal range intervals, or equal log-range intervals
#'   (default `TRUE`)
#' @param ylim (vector) limits of the y axis
#' @param title (string) a title to display above the plot
#' @param label (integer) index of letter for subplot tag
#' @param gPars (list) graphical parameters
#' @param plot (logical) plot the results
#' @param score (logical) estimate calibration stats (default: `FALSE`)
#' @param xlab (string) X axis label
#' @param xlim (vector) min and max values of X axis
#' @param add (logical) add to previous graph ?
#' @param col (integer) color index of curve to add
#'
#' @return Invisibly returns a list of LZISD results. Mainly used
#'   for its plotting side effect.
#' @export
#'
#' @examples
#' \donttest{
#'   uE  = sqrt(rchisq(1000, df = 4))  # Re-scale uncertainty
#'   E   = rnorm(uE, mean=0, sd=uE)    # Generate errors
#'   plotLZISD(uE, E/uE, method = 'cho', ylim = c(0,2))
#' }
plotLZISD = function(
  X, Z,
  aux       = NULL,
  sdZ       = 1,
  logX      = FALSE,
  nBin      = NULL,
  intrv     = NULL,
  equiPop   = TRUE,
  popMin    = 30,
  logBin    = TRUE,
  plot      = TRUE,
  slide     = FALSE,
  method    = c('cho','bootstrap'),
  BSmethod  = c('bca','perc','basic'),
  nBoot     = 500,
  xlab      = 'Calculated value',
  xlim      = NULL,
  ylim      = NULL,
  title     = '',
  score     = FALSE,
  add       = FALSE,
  col       = 5,
  label     = 0,
  gPars     = ErrViewLib::setgPars()
) {

  method   = match.arg(method)
  BSmethod = match.arg(BSmethod)

  N = length(Z)

  if(is.null(aux))
    ord  = order(X)
  else
    ord  = order(X,aux)
  xOrd = X[ord]
  zOrd = Z[ord]

  # Design local areas
  if(is.null(intrv)) {
    if(is.null(nBin))
      nBin  = max(min(floor(N/150),15),2)
    if(nBin <= 0)
      stop('>>> nBin should be > 0')
    Xin = N
    if(!equiPop)
      Xin = xOrd
    intrv = ErrViewLib::genIntervals(Xin, nBin, slide, equiPop, popMin, logBin)
  }
  nBin0 = nBin
  nBin = intrv$nbr

  # if(is.null(slide))
  #   slide = intrv$nbr <= 4

  # LZV values
  mV = loV = upV = mint = c()
  for (i in 1:nBin) {
    sel = intrv$lwindx[i]:intrv$upindx[i]
    M      = length(sel)
    zLoc   = zOrd[sel]
    zs     = rsdZCI(zLoc,
                    nBoot = max(nBoot, M),
                    method = method,
                    CImethod = BSmethod)
    mV[i]  = zs$mean
    loV[i] = zs$ci[1]
    upV[i] = zs$ci[2]
    mint[i] = mean(range(xOrd[sel])) # Center of interval
  }
  zs = rsdZCI(
    zOrd,
    nBoot = max(nBoot, N),
    method = method,
    CImethod = BSmethod
  )
  mV0 = zs$mean
  loV0 = zs$ci[1]
  upV0 = zs$ci[2]

  if(plot) {
    if(length(gPars) == 0)
      gPars = ErrViewLib::setgPars()

    for (n in names(gPars))
      assign(n, rlist::list.extract(gPars, n))

    if(is.null(xlim))
      xlim = range(xOrd)

    if(!add) {
      par(
        mfrow = c(1, 1),
        mar = c(mar[1:3],3), # Reserve of right margin space
        mgp = mgp,
        pty = 's',
        tcl = tcl,
        cex = cex,
        lwd = lwd,
        cex.main = 1
      )

      if(is.null(ylim))
        ylim = range(c(loV, upV, sdZ))

      matplot(
        mint,
        mV,
        xlab = xlab,
        ylab = 'Local Z-score Inverse SD',
        xlim = xlim,
        xaxs = 'i',
        ylim = ylim,
        type = 'b',
        lty  = 3,
        pch  = 16,
        lwd  = lwd,
        cex  = ifelse(slide,0.5,1),
        col  = cols[col],
        main = title,
        log  = ifelse(logX,'xy','y')
      )
      grid(equilogs = FALSE)
      abline(h   = sdZ,
             lty = 2,
             col = cols[1],
             lwd = lwd)
      mtext(text = paste0(signif(sdZ,2),' -'),
            side = 2,
            at = sdZ,
            col = cols[1],
            cex = 0.75*cex,
            las = 1,
            adj = 1,
            font = 2)

      # Mean variance header
      ypos = 10^par("usr")[4]
      mtext(text = ' Average',
            side = 4,
            at =  ypos,
            col = cols[1],
            cex = 0.75*cex,
            las = 1,
            font = 2)

    } else {

      lines(
        mint,
        mV,
        type = 'b',
        lty  = 3,
        pch  = 16,
        lwd  = lwd,
        cex  = ifelse(slide,0.5,1),
        col  = cols[col]
      )
    }

    # Mean + CI
    pm = round(mV0,2)
    mtext(text = paste0('- ',pm),
          side = 4,
          at   = mV0,
          col  = cols[col],
          cex  = 0.75*cex,
          las  = 1,
          font = 2)
    segments(
      xlim[2], loV0,
      xlim[2], upV0,
      col  = cols[col],
      lwd  = 6 * lwd,
      lend = 1
    )

    if(slide) {
      ipl = seq(1,length(mint),length.out=nBin0)
      polygon(
        c(mint,rev(mint)),
        c(loV, rev(upV)),
        col = cols_tr[col],
        border = NA)
      segments(
        mint[ipl], loV[ipl],
        mint[ipl], upV[ipl],
        col  = cols[col],
        lwd  = 1.5 * lwd,
        lend = 1)

    } else {
      segments(
        mint, loV,
        mint, upV,
        col  = cols[col],
        lwd  = 1.5 * lwd,
        lend = 1)
    }

    box()

    if(label > 0)
      mtext(
        text = paste0('(', letters[label], ')'),
        side = 3,
        adj = 1,
        cex = cex,
        line = 0.3)

  }

  ZISDE = ZISDEUp = ZISDM = ZISDMs = NA
  if(score) {
    scores = abs(log(mV))
    # Max deviation
    im = which.max(scores)
    ZISDM = exp( sign(log(mV[im])) * scores[im] )
    # Significant ?
    ZISDMs = FALSE
    if(sdZ < loV[im] | sdZ > upV[im])
      ZISDMs = TRUE
    # Mean deviation
    ZISDE = exp(mean(scores))
    scores = pmax(log(upV/mV),log(mV/loV))
    ZISDEUp = exp(mean(scores))
  }

  invisible(
    list(
      mint   = mint,
      lwindx = intrv$lwindx,
      upindx = intrv$upindx,
      pc     = mV,
      pcl    = loV,
      pcu    = upV,
      meanP  = mV0,
      meanPl = loV0,
      meanPu = upV0,
      ZISDE  = ZISDE,
      ZISDEUp= ZISDEUp,
      ZISDM  = ZISDM,
      ZISDMs = ZISDMs
    )
  )
}