.packageName <- "NADA"
###
# NADA for R by Lopaka Lee.
#
# Version 1.2-1
# Copyright (2004, 2005) Lopaka Lee
#
# A S-language software module based on 
# methodologies described by Dennis R. Helsel in his book 
# Nondetects and Data Analysis: Statistics for Censored Environmental Data.
#
# NADA is free software; you can redistribute it and/or modify it under
# the terms of the GNU General Public License as published by the Free
# Software Foundation; either version 2 of the License, or (at your option)
# any later version.
# 
# NADA is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
# for more details.  You should have received a copy of the GNU General
# Public License along with NADA; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
###

### Definitions common to All sections in this package

## Globals

.onLoad = function(lib, pkg) require(methods)

NADAprobs = c(0.05,0.10,0.25,0.50,0.75,0.90,0.95)

## Generics

setGeneric("plot", function(x, y, ...) standardGeneric("plot"))

setGeneric("print", function(x, ...) standardGeneric("print"))

setGeneric("summary", function(object, ...) standardGeneric("summary"))

setGeneric("mean", function(x, ...) standardGeneric("mean"))

setGeneric("sd", function(x, na.rm=FALSE) standardGeneric("sd"))

setGeneric("median", function(x, na.rm=FALSE) standardGeneric("median"))

setGeneric("quantile", function(x, ...) standardGeneric("quantile"))

setGeneric("predict", function(object, ...) standardGeneric("predict"))

setGeneric("lines", function(x, ...) standardGeneric("lines"))

## Classes

setClass("NADAlist", "list")

## Methods

setMethod("print", signature("NADAlist"), function(x, ...)
{
    tag = names(x)
    for (i in 1:length(x))
      {
        cat(tag[i], "\n")
        print(x[[i]])
        cat("\n")
      }
})

#-->> BEGIN general utility functions

##
# splitQual extracts qualifed and unqualifed vectors from a vector
# containing concatenated qualifiying characters and numeric values
# like "<0.5".  Only handles one kind of censoring character/symbol.
splitQual =
function(v, qual.symbol = "<")
{
    qual.index = grep(qual.symbol, x=as.character(v))  

    qual.chars = as.character(v[qual.index])
    qual = as.numeric(sub(qual.symbol, "", qual.chars))

    unqual.index = -1 * qual.index
    unqual = as.numeric(as.character(v[unqual.index]))
    
    return(list(qual         = qual, 
                unqual       = unqual, 
                qual.index   = qual.index,
                unqual.index = unqual.index))
}

## pct.censored  -- Simple function to save some typing
pct.censored =
function(obs, censored)
{
    if (!is.logical(censored)) 
      {
        stop("censored indicator must be logical vector!\n")
      }

    return(100*(length(obs[censored])/length(obs)))
}

#-->> END general utility functions

#-->> BEGIN survival analysis based functions -- "cencen" section 

## Generics

setGeneric("cenfit", 
           function(obs, censored, groups, ...) standardGeneric("cenfit"))

setGeneric("cendiff", 
           function(obs, censored, groups, ...) standardGeneric("cendiff"))

setGeneric("flip", function(x) standardGeneric("flip"))

## Classes 

setOldClass("Surv")
setOldClass("survfit")

setClass("Cen", representation(Surv="Surv", flipFactor="numeric"))
setClass("cenfit", representation(survfit="survfit"))

## Methods

# Begin cencen utility functions

# LCL() and UCL() return string representations of 
# the lower and upper conf limits of cenfit object ("0.95LCL")
LCL =
function(x) { paste(x@survfit$conf.int, "LCL", sep='') }

UCL =
function(x) { paste(x@survfit$conf.int, "UCL", sep='') }

# stepfind() -- More applicable version of stats::stepfun(). 
# Used in predict.cenfit() and quantile.cenfit().
# Given decreasingly ordered x and y vectors of a step function, 
# find the y value associated with any given x value.  
# Optionally right or left looking on the number line. 
stepfind =
function(x, y, val, right=TRUE) 
{
    findStep =
    function(x, y, val, right)
    {
        i = length(x)
        if (val >= max(x)) return(max(y))
        if (val <= min(x)) return(min(y))

        while (as.logical(i)) 
          {
            if (x[i] > val || identical(all.equal(x[i], val), TRUE)) break
            i = i - 1
          }
        return(ifelse(right, y[i], y[i+1]))
    }

    sapply(val, findStep, x=x, y=y, right=right)
}

## End cencen utility functions

# Cen() is analgous to Surv() in survival package.
# However, Cen() can be used in the flip() routines below to
# rescale left-censored data into right-censored data 
# for use in the survival package routines.
Cen =
function(obs, censored, type="left")
{
    new("Cen", Surv=Surv(obs, !censored, type=type), flipFactor = max(obs))
}

setMethod("print", signature(x="Cen"), function(x, ...) print(x@Surv))

## flip() methods for rescaling input for use in survival routines.

# Note that flip()ing a Cen object results in a Surv object.
setMethod("flip", signature(x="Cen"), function(x)
{
    surv = x@Surv
    # Subtract the flipFactor from all observations
    surv[,1] = x@flipFactor - surv[,1]

    # Update the censoring type
    if (attr(surv, "type") == "right") attr(surv, "type") = "left"
    else if (attr(surv, "type") == "left") attr(surv, "type") = "right"
    else stop("Can only flip \"left\" or \"right\" censored Cen objects")

    return(surv)
})

# flip()ing a formula just symbolically updates the 
# response (which should be a Cen object). 
# Result is like: flip(Cen(obs, cen))~groups
setMethod("flip", signature(x="formula"), function(x) update(x, flip(.)~.))

## Begin cencen.* functions
# These routines are allow cenfit and cendiff methods
# to be used with Cen objects, vectors, or formulas as input.
# Note they all convert input to a formula and call the formula method.

cencen.Cen =
function(obs, censored, groups, ...)
{
    cl = match.call()
    f = as.formula(substitute(a~1, list(a=cl[[2]])))
    environment(f) = parent.frame()
    callGeneric(f, ...)
}

cencen.vectors =
function(obs, censored, groups, ...)
{
    cl = match.call()
    f = as.formula(substitute(Cen(a, b)~1, list(a=cl[[2]], b=cl[[3]])))
    environment(f) = parent.frame()
    callGeneric(f, ...)
}

cencen.vectors.groups =
function(obs, censored, groups, ...)
{
    cl = match.call()
    f = substitute(Cen(a, b)~g, list(a=cl[[2]], b=cl[[3]], g=cl[[4]]))
    f = as.formula(f)
    environment(f) = parent.frame()
    callGeneric(f, ...)
}
## End cencen.* routines

# cenfit for formulas
setMethod("cenfit", 
          signature(obs="formula", censored="missing", groups="missing"), 
          function(obs, censored, groups, ...)
{
    sf = survfit(flip(obs), ...)

    cenObj = eval(obs[[2]], environment(obs))
    sf$time = cenObj@flipFactor - sf$time 

    new("cenfit", survfit=sf)
})

setMethod("cenfit", 
          signature(obs="Cen", censored="missing", groups="missing"), 
          cencen.Cen)

setMethod("cenfit", 
          signature(obs="numeric", censored="logical", groups="missing"), 
          cencen.vectors)

setMethod("cenfit", 
          signature(obs="numeric", censored="logical", groups="factor"), 
          cencen.vectors.groups)

# cendiff for formulas
setMethod("cendiff", 
          signature(obs="formula", censored="missing", groups="missing"), 
          function(obs, censored, groups, rho=1, ...)
{
    x = survival::survdiff(flip(obs), rho=rho, ...)
    # To do: modify the call object to reflect cenfit call
    # for now, just nullify it -- users typically remember the call.
    x$call = NULL
    return(x)
})

#setMethod("cendiff", 
#          signature(obs="Cen", censored="missing", groups="missing"), 
#          cencen.Cen)

setMethod("cendiff", 
          signature(obs="numeric", censored="logical", groups="factor"), 
          cencen.vectors.groups)

# Indexing cenfit objects.
# When this happens, we throw away the strata information.
setMethod("[", signature(x="cenfit", i="numeric", j="missing"), 
          function(x, i, drop=FALSE) 
{
    s = x@survfit

    if (is.null(s$strata)) stop("can't index; object contains only one ECDF")

    s$n = s$strata.all[i]

    s = s[i]

    s$strata        = NULL
    s$strata.all    = NULL
    s$ntimes.strata = NULL

    new("cenfit", survfit=s)
})

setMethod("plot", signature(x="cenfit"),
          function(x, y, log  = 'x', axLimitFactor = 0.8, 
                   ylab = "Probability", xlab = "Value", 
                   lty  = seq(1,6), ...)
{
    s = x@survfit

    firstx = (min(s$time)*axLimitFactor)

    plot(s, log=log, firstx=firstx, ylab=ylab, xlab=xlab, lty=lty, ...)

    # Draw a vertical line at the minimum observation -- this is
    # the lower extent of the step curve or ECDF.
    abline(v=min(s$time))

    abline(h=1.0)
})

setMethod("lines", "cenfit", function(x, ...) lines(x@survfit, ...))

# Private predict method for cenfit objects -- does not handle strata!
.predict.cenfit =
function(x, newdata, conf.int=FALSE) 
{
    ret  = NULL
    s = x@survfit
    pred = stepfind(s$time, s$surv, newdata)

    if (!conf.int) ret = pred
    else
      {
        pred.l = stepfind(s$time, s$lower, newdata)
        pred.u = stepfind(s$time, s$upper, newdata)

        ret = data.frame(newdata, pred, pred.l, pred.u)

        names(ret) = c("obs", "prob", LCL(x), UCL(x))
      }

    return(ret)
}

# Public predict method for cenfit objects
setMethod("predict", signature(object="cenfit"),
          function(object, newdata, conf.int=FALSE) 
{
    ret  = NULL
    s = object@survfit

    if (is.null(s$strata)) ret = .predict.cenfit(object, newdata, conf.int)
    else
      {
        for (i in 1:length(s$strata))
          {
            ret[[i]] = .predict.cenfit(object[i], newdata, conf.int)
          }
        names(ret) = names(s$strata)
        class(ret) = "NADAlist"
      }
      
    return(ret)
})

# Private quantile method for cenfit objects -- does not handle strata!
.quantile.cenfit =
function(x, newdata, conf.int=FALSE) 
{
    ret  = NULL
    s = x@survfit
    quan = stepfind(s$surv, s$time, newdata, FALSE)

    if (!conf.int) 
      {
        ret = quan
        names(ret) = 
          paste(formatC(100 * newdata, format = "fg", wid = 1), "%", sep='')
      }
    else
      {
        quan.l = stepfind(s$surv, s$lower, newdata, FALSE)
        quan.u = stepfind(s$surv, s$upper, newdata, FALSE)

        ret = data.frame(newdata, quan, quan.l, quan.u)

        names(ret) = c("quantile", "obs", LCL(x), UCL(x))
      }

    return(ret)
}

# Public quantile method for cenfit objects
setMethod("quantile", signature(x="cenfit"),
          function(x, probs = NADAprobs, conf.int = FALSE, ...)
{
    ret = NULL
    s = x@survfit

    if (is.null(s$strata)) ret = .quantile.cenfit(x, probs, conf.int)
    else
      {
        for (i in 1:length(s$strata))
          {
            ret[[i]] = .quantile.cenfit(x[i], probs, conf.int)
          }
        names(ret) = names(s$strata)
        class(ret) = "NADAlist"
      }

    return(ret)
})

# Private mean method for cenfit objects -- does not handle strata!
# Given a cenfit or survfit object returns the calculated mean and se(mean)
.mean.cenfit =
function(x)
{
    x = x@survfit

    stime   = x$time
    surv    = x$surv
    n.risk  = x$n.risk
    n.event = x$n.event

    mean = NULL
    varmean = NULL

    min.stime = min(stime)
    min.time = min(0, min.stime)
    n = length(stime)

    hh = c(ifelse((n.risk[-n] - n.event[-n]) == 0, 
                   0, 
                   n.event[-n]/(n.risk[-n] * (n.risk[-n] - n.event[-n]))), 0)

    ndead = sum(n.event)
    dif.time = c(diff(c(min.time, stime)), 0)

    if (!is.matrix(surv)) 
      {
        mean = dif.time * c(1, surv)
        varmean = sum(rev(cumsum(rev(mean))^2)[-1] * hh)
      }
    else 
      {
        n = nrow(surv)
        mean = dif.time * rbind(1, surv)
        if (n == 1) 
            temp = mean[2, , drop = FALSE]
        else temp = (apply(mean[(n + 1):2, , drop = FALSE], 2, 
            cumsum))[n:1, , drop = FALSE]
        varmean = c(hh %*% temp^2)
      }

    ret = c(sum(mean), sqrt(varmean))
    names(ret) = c("rmean", "se(rmean)")

    return(ret)
}

# Public mean method for cenfit objects
setMethod("mean", signature(x="cenfit"), function(x, ...)
{
    ret = NULL
    s = x@survfit

    if (is.null(s$strata)) ret = .mean.cenfit(x)
    else
      {
        for (i in 1:length(s$strata)) ret[[i]] = .mean.cenfit(x[i])
        names(ret) = names(s$strata)
        class(ret) = "NADAlist"
      }

    return(ret)
})

setMethod("median", signature(x="cenfit"), function(x, na.rm=FALSE)
{
    # To do: remove NAs?
    ret = as.vector(quantile(x, 0.5))
    names(ret) = names(x@survfit$strata)
    return(ret)
})

setMethod("sd", signature(x="cenfit"), function(x, na.rm=FALSE)
{
    # To do: remove NAs?
    ret = NULL
    s = x@survfit

    if (is.null(s$strata)) ret = (sqrt(s$n) * as.numeric(mean(x)[2]))
    else
      {
        for (i in 1:length(s$strata))
          {
            n = as.numeric(s$strata.all[i])
            ret[[i]] = (sqrt(n) * as.numeric(mean(x[i])[2]))
          }
        names(ret) = names(s$strata)
      }

    return(ret)
})

setMethod("print", signature(x="cenfit"), function(x, ...)
{
    s = x@survfit

    summaryVec =
    function(x)
    {
      s = x@survfit

      n      = s$n
      cen    = s$n - sum(s$n.event)
      median = median(x)
      mean   = mean(x)

      return(c(n, cen, median, mean))
    }

    ret = NULL
    tag = c("n", "n.cen", "median", "mean", "se(mean)")

    if (is.null(s$strata))
      {
        ret = summaryVec(x)
        names(ret) = tag
      }
    else
      {
        ret = summaryVec(x[1])
        for (i in 2:length(s$strata))
          {
            ret = rbind(ret, summaryVec(x[i]))
          }
        colnames(ret) = tag
        rownames(ret) = names(s$strata)
      }

    print(ret, ...)
    invisible(ret)
})

setMethod("summary", signature(object="cenfit"),
          function(object, times, censored=FALSE, scale=1, ...)
{
    strataSummary =
    function(x)
    {
      s = x@survfit

      # Reverse vectors -- we have leftist views ;)
      s = lapply(s, rev)
      ret = data.frame(s$time, s$n.risk, s$n.event, 
                       s$surv, s$std.err, s$lower, s$upper)
      names(ret) = c("obs", "n.risk", "n.event", 
                     "prob", "std.err", LCL(x), UCL(x))
      return(ret)
    }

    s   = object@survfit
    ret = NULL

    if (is.null(s$strata)) ret = strataSummary(object)
    else
      {
        for (i in 1:length(s$strata)) ret[[i]] = strataSummary(object[i])
        names(ret) = names(s$strata)
        class(ret) = "NADAlist"
      }

    return(ret)
})

#-->> END survival analysis based methods -- "cencen" section
#-->> BEGIN Regression on Order Statistics (ROS) section

## Generics

setGeneric("ros", 
           function(obs, censored, forwardT = "log", reverseT = "exp") 
             standardGeneric("ros"))

## Classes

setClass("ros", "list")

## Methods

# ros -- Regression on Order Statistics (ROS).
# An implementation of ROS for left-censored data (less-thans) 
# containing one to multiple censoring thresholds. See man page.
setMethod("ros", signature(obs="numeric", censored="logical"),
          function(obs, censored, forwardT = "log", reverseT = "exp")
{
    if (is.null(forwardT) || is.null(reverseT)) {
        forwardT = reverseT = ".trueT"
    }
    else if (!exists(forwardT)) {
         stop("Can not find Forward Transformation function: ", forwardT, "\n")
    }
    else if (!exists(reverseT)) {
        stop("Can not find Reverse Transformation function: ", reverseT, "\n")
    }

    if ( (length(obs[censored])/length(obs)) > 0.8 ) {
        warning("Input > 80% censored -- Results are tenuous.\n")
    }

    ix = order(obs)
    obs = obs[ix]
    censored = censored[ix]

    pp = hc.ppoints(obs, censored)

    pp.nq = qnorm(pp[!censored])
    obs.transformed = get(forwardT)(obs[!censored])
    hc = lm(obs.transformed ~ pp.nq)

    oldClass(hc) = c("ros", "lm")
    hc$obs      = obs
    hc$modeled  = obs 
    hc$pp       = pp
    hc$censored = censored 
    hc$reverseT = reverseT
    hc$forwardT = forwardT

    hc$modeled[censored] = predict(hc, qnorm(pp[censored]))

    return(hc)
})

# ROS is technically "Linear Regression on Order Statistics" 
# or LROS.  Version 1.0-1 had funtions and objects as "lros".  
# As of Version 1.0-2, I've renamed things to "ros" to keep consistent
# with the literature on this method which calls it just ROS.
# This is provided for backward compatibility.
lros = ros

#  .trueT is provided so that ros() can be used with no transforms.
#  (It is a quick hack -- rather than recoding ros and family)
.trueT =
function(x)
{
    return(x)
}

setMethod("print", signature(x="ros"), function(x, ...)
{
    cat("\nMultiply-Censored ROS Model\n\n")

    uncen.n = length(x$modeled[!x$censored])
    cen.n   = length(x$modeled[x$censored])

    cat("           N:", (uncen.n + cen.n), "\n")
    cat("    Censored:", cen.n, "\n")
    cat("  % Censored:", format((cen.n/(cen.n+uncen.n))*100, digits=4), "\n")
    cat("\n");
    
    cat("        Mean:", format(mean(x), digits=4), "\n")
    cat("      StdDev:", format(sd(x), digits=4), "\n")
    cat("\n");
    cat("   Quantiles:\n")
    print(quantile(x))
    cat("\n");
    cat("Use summary() to view the linear regression model\n\n")
})

setMethod("summary", "ros", function(object, plot=FALSE, ...)
{
    ret = summary.lm(object, ...)
    if (plot)
      {
        oldClass(object) = "lm"
        plot.lm(object, ...)
      }
    return(ret)
})

# S3 Method -- ros conversion to data.frame discards all linear model info
as.data.frame.ros = 
function (x, row.names = NULL, optional = FALSE)
{
    x = list(obs=x$obs, censored=x$censored, pp=x$pp, modeled=x$modeled)
    as.data.frame(x, row.names, optional)
}

setMethod("mean", signature(x="ros"), function(x, ...) mean(x$modeled, ...))

setMethod("sd", signature(x="ros"), function(x, na.rm=FALSE)
{
    sqrt(var(x=x$modeled))
})

setMethod("median", "ros", function(x, na.rm=FALSE) median(x$modeled))

## Query and prediction functions for ROS objects

setMethod("quantile", signature(x="ros"), function(x, probs=NADAprobs,...)
{
    quantile(x$modeled, probs, ...)
})

# predict method -- given new normal quantiles of plotting positions,
# returns the corresponding modeled values.
setMethod("predict", "ros", function(object, newdata, ...)
{
    predicted = 
      as.vector(predict.lm(object, newdata=data.frame(pp.nq=newdata), ...))

    return(get(object$reverseT)(predicted))
})

## Generic plotting functions

setMethod("lines", signature(x="ros"), function (x, ...)
{
    minmax.nq = qnorm(c(min(x$pp), max(x$pp)))
    lines(x=minmax.nq, y=predict(x, minmax.nq), ...)
})


# Constructs a prob-plot representation of a ROS model
setMethod("plot", signature(x="ros", y="missing"),
          function(x, 
                   plot.censored = FALSE, 
                   lm.line = TRUE, 
                   grid    = TRUE, 
                   ylab    = "Value", 
                   pch     = 16,
                   ... )
{
    ##
    # To do:
    #   Refactor this routine -- it is long and ugly.
    #   Constrain y ticks to factors of 10?
    #   Draw in grid lines at y ticks
    #
    uncen = x$modeled[!x$censored]
    cen   = x$modeled[x$censored]

    pp.uncen.nq = qnorm(x$pp[!x$censored])
    pp.cen.nq   = qnorm(x$pp[x$censored])

    ymin = min(c(uncen, cen))
    ymax = max(c(uncen, cen))
    xmin = min(c(pp.uncen.nq, pp.cen.nq))
    xmax = max(c(pp.uncen.nq, pp.cen.nq))

    if (x$forwardT == "log") {
        plot(y    = uncen, 
             x    = pp.uncen.nq, 
             ylim = c(ymin, ymax), 
             xlim = c(xmin, xmax), 
             ylab = ylab,
             xlab = "Normal Quantiles",
             log  = "y",
             yaxt = "n",
             pch  = 16,
             ...
        )
     }
     else {
        plot(y    = uncen, 
             x    = pp.uncen.nq, 
             ylim = c(ymin, ymax), 
             xlim = c(xmin, xmax), 
             ylab = ylab,
             xlab = "Normal Quantiles",
             yaxt = "n",
             pch  = 16,
             ...
        )
     }

    axis(2, axTicks(2))

    # Draw a line through the predicted xmin and xmax points 
    if (lm.line) lines(x)

    # Draw top "Chance of Exceedance" axis
    labels = c("5","10","25","50","75","90","95")
    atv = qnorm(c(.05, .1, .25, .50, .75, .90, .95))
    axis(3, at=atv, labels=rev(labels), las=2)
    mtext("Percent Chance of Exceedance", side=3, line=3)

    # Plot the synthetic censored points -- if requested
    if (plot.censored) points(y=cen, x=pp.cen.nq)

    # Draw in grid at major divisions  -- if requested
    if (grid)
      {
        abline(v=atv, lty="dotted")
        abline(h=axTicks(2), lty="dotted")
      }
})


## Routines for calculating Helsel-Cohn style plotting positions

# hc.ppoints calculates computes Wiebull-type plotting postions of data
# containing mixed uncensored and left-censored data. See man page.
# If there are no censored values, the plotting postitions are calculated
# using the standard ppoints() function.
hc.ppoints = 
function(obs, censored)
{    
    if (!is.logical(censored)) 
      {
        stop("censored indicator must be logical vector!\n")
      }

    pp = numeric(length(obs))

    if (!any(censored)) pp = ppoints(obs)
    else
      {
        uncen = obs[!censored]
        cen   = obs[censored]

        cn = .cohnN(uncen, cen)
        pp[!censored] = hc.ppoints.uncen(cn=cn)
        pp[censored]  = hc.ppoints.cen(cn=cn)
      }

    return(pp)
}

# .cohnN Calculates 'Cohn' numbers -- quantities described by 
# Helsel and Cohn's (1988) reformulation of a prob-plotting formula
# described by Hirsch and Stedinger (1987).
.cohnN =
function(uncen, cen)
{
    alldata = c(uncen, cen)
    A = B = C = P = numeric()

    limit = sort(unique(cen))

    a = length(uncen[uncen < limit[1]])
    if (a > 0) { limit = c(0, limit) }

    i = length(limit)

    A[i] = length(uncen[ uncen >= limit[i] ])
    B[i] = length(alldata[alldata <= limit[i]])
    C[i] = length(cen[ cen == limit[i] ])
    P[i] = A[i]/(A[i] + B[i])

    i = i - 1
    while (i > 0)
      {
        A[i] = length(uncen[ uncen >= limit[i] & uncen < limit[i + 1] ])
        B[i] = length(alldata[alldata <= limit[i]])
        C[i] = length(cen[ cen == limit[i] ])
        P[i] = P[i + 1] + ((A[i]/(A[i] + B[i])) * (1 - P[i + 1]))

        i = i - 1
      }
    return(list(A=A, B=B, C=C, P=P, limit=limit))
}

# hc.ppoints.uncen calculates plotting postions for uncensored data.
hc.ppoints.uncen =
function(obs, censored, cn)
{
    if (missing(cn)) { cn = .cohnN(obs, censored) }

    nonzero = (cn$A != 0)
    A     = cn$A[nonzero]
    B     = cn$B[nonzero]
    P     = cn$P[nonzero]
    limit = cn$limit[nonzero]

    pp = numeric()
    n = length(limit)
    for (i in 1:n)
      {
        R = 1:A[i] 

        k = P[i+1]
        if (is.na(k)) k = 0

        for (r in 1:length(R))
          {
            pp = c(pp, (1 - P[i]) + ((P[i] - k) * R[r])/(A[i] + 1))
          }
      }
    return(pp)
}

# hc.ppoints.cen calculates plotting postions for censored data.
hc.ppoints.cen =
function(obs, censored, cn)
{    
    if (missing(cn)) { cn = .cohnN(obs, censored) }

    C     = cn$C
    P     = cn$P
    limit = cn$limit

    if (P[1] == 1) 
      {
        C     = C[-1]
        P     = P[-1]
        limit = limit[-1]
      }
    
    pp = numeric()
    for (i in 1:length(limit)) 
      {
        c.i = C[i]
        for (r in 1:c.i) 
          {
            pp = c(pp, (1 - P[i]) * r/(c.i + 1))
          }
      }
    return(pp)
}

#-->> END Regression on Order Statistics (ROS) section
