Nothing
#-->> BEGIN Regression on Order Statistics (ROS) section
## Generics
setGeneric("ros",
function(obs, censored, forwardT = "log", reverseT = "exp",
na.action=getOption("na.action"))
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.
setMethod("ros", signature(obs="numeric", censored="logical"),
function(obs, censored, forwardT = "log", reverseT = "exp",
na.action=getOption("na.action"))
{
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 (anyNA(obs)||anyNA(censored)) {
if (is.null(na.action)) na.action = "na.exclude"
obs = do.call(na.action, list(obs))
censored = do.call(na.action, list(censored))
}
if ( (length(obs[censored])/length(obs)) > 0.8 ) {
warning("Input > 80% censored -- Results are tenuous.\n")
}
ix = which(obs > max(obs[!censored]))
if (length(ix)) {
warning("Dropped censored values that exceed max of uncensored values.")
obs = obs[-ix]
censored = censored[-ix]
}
ix = order(obs)
obs = obs[ix]
censored = censored[ix]
pp = hc.ppoints(obs, censored, na.action)
pp.nq = qnorm(pp[!censored])
obs.transformed = get(forwardT)(obs[!censored])
hc = lm(obs.transformed ~ pp.nq, na.action=na.action)
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]), na.action=na.action)
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
# Also, since the rest of the package is evolving to use functions
# names like cen*
cenros = 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("show", signature(object="ros"), function(object)
{
print(object)
})
setMethod("print", signature(x="ros"), function(x, ...)
{
n = length(x$modeled)
n.cen = length(x$modeled[x$censored])
median = median(x)
mean = mean(x)
sd = sd(x)
ret = c(n, n.cen, median, mean, sd)
names(ret) = c("n", "n.cen", "median", "mean", "sd")
print(ret, ...)
invisible(ret)
})
setMethod("summary", "ros", function(object, plot=FALSE, ...)
{
ret = summary.lm(object, ...)
if (plot)
{
oldClass(object) = "lm"
plot(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))
### min and max methods don't work yet cuz I don't know how to eval
## the expanded dots correctly yet.
#setMethod("min", "ros", function(..., na.rm=FALSE)
#{
# x = match.call(expand.dots=TRUE)[[2]]
# min(eval.parent(x$modeled))
#})
#setMethod("max", "ros", function(..., na.rm=FALSE)
#{
# x = match.call(expand.dots=TRUE)[[2]]
# max(eval.parent(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))
})
# pexceed method -- given new values or concentrations,
# returns the probability of exceedance.
setMethod("pexceed", "ros",
function(object, newdata, conf.int=FALSE, conf.level=0.95, ...)
{
cen = object$censored
forwardT = get(object$forwardT)
pp.nq = qnorm(object$pp[!cen])
obs.transformed = forwardT(object$obs[!cen])
n.lm = lm(pp.nq~obs.transformed)
nd = data.frame(obs.transformed=forwardT(newdata))
ret = NULL
if (!conf.int) ret = 1 - pnorm(as.vector(predict.lm(n.lm, nd)))
else
{
ret = 1 - pnorm(predict.lm(n.lm, nd,
interval="confidence", level=conf.level))
lcl = paste(as.character(conf.level), "LCL", sep='')
ucl = paste(as.character(conf.level), "UCL", sep='')
ret[,c(2,3)] = ret[,c(3,2)]
colnames(ret) = c("prob", lcl, ucl)
rownames(ret) = as.character(newdata)
}
return(ret)
})
## 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), ...)
})
## Broken for the time being -- use lines
#setMethod("abline", signature(a="ros"),
# function(a, b, h, v, reg, coef, untf, col, lty, lwd, ...)
#{
# minmax.nq = qnorm(c(min(a$pp), max(a$pp)))
# lines(x=minmax.nq, y=predict(a, 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")
}
})
#setMethod("boxplot", signature(x="ros"), function(x, log="y", range=0,...)
setMethod("boxplot", signature(x="ros"), function(x, ...)
{
# use boxplot.default instead
bx = boxplot(x$modeled, plot=FALSE, ...)
# The vector stats is the famous five -- with outlier limts.
# We relpace the quantiles with the ros-modeled ones
bx$stats[2:4] = as.vector(quantile(x, c(0.25, 0.5, 0.75)))
bx$stats = as.matrix(bx$stats)
bxp(bx, ...)
invisible(bx)
})
## 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, na.action = getOption("na.action"))
{
if (!is.logical(censored)) stop("censored indicator must be logical vector!\n")
if (anyNA(obs)||anyNA(censored)) {
if (is.null(na.action)) na.action = "na.omit"
obs = do.call(na.action, list(obs))
censored = do.call(na.action, list(censored))
}
pp = numeric(length(obs))
if (!any(censored)) pp = ppoints(obs)
else
{
cn = cohn(obs, censored)
pp[!censored] = hc.ppoints.uncen(obs, censored, cn, na.action)
pp[censored] = hc.ppoints.cen(obs, censored, cn, na.action)
}
return(pp)
}
# cohn Calculates "Cohn" Numbers -- quantities described by
# Helsel and Cohn's (1988) reformulation of a prob-plotting formula
# described by Hirsch and Stedinger (1987).
#
# The Cohn Numbers are:
# A_j = the number of uncensored obs above the jth threshold.
# B_j = the number of observations (cen & uncen) below the jth threshold.
# C_j = the number of censored observations at the jth threshold
# P_j = the probability of exceeding the jth threshold
cohn =
function(obs, censored, na.action = getOption("na.action"))
{
if (!is.logical(censored)) stop("censored indicator must be logical vector!\n")
if (anyNA(obs)||anyNA(censored)) {
if (is.null(na.action)) na.action = "na.omit"
obs = do.call(na.action, list(obs))
censored = do.call(na.action, list(censored))
}
uncen = obs[!censored]
cen = obs[censored]
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(obs[obs <= limit[i]])-length(uncen[uncen==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(obs[obs <= limit[i]])-length(uncen[uncen==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, na.action = getOption("na.action"))
{
if (!is.logical(censored)) stop("censored indicator must be logical vector!\n")
if (anyNA(obs)||anyNA(censored)) {
if (is.null(na.action)) na.action = "na.omit"
obs = do.call(na.action, list(obs))
censored = do.call(na.action, list(censored))
}
if (missing(cn)) { cn = cohn(obs, censored) }
#cn = cohn(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, na.action = getOption("na.action"))
{
if (!is.logical(censored)) stop("censored indicator must be logical vector!\n")
if (anyNA(obs)||anyNA(censored)) {
if (is.null(na.action)) na.action = "na.omit"
obs = do.call(na.action, list(obs))
censored = do.call(na.action, list(censored))
}
if (missing(cn)) { cn = cohn(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
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.