# R/nomogram.R In UncertainInterval: Uncertain Interval Methods for Three-Way Cut-Point Determination in Test Results

#### Documented in nomogram

# prob.pre.test=.3; LR=c(44, .022); probs.post.test=c(pos=NULL, neg=NULL)
# prob.pre.test=.3; probs.post.test=c(0.9778, 1-0.9785)); LR=c(NULL, NULL);

#' Fagan's nomogram to show the relationships between the prior probability, the
#' likelihood ratios, sensitivity and specificity, and the posterior
#' probability.
#'
#' Next to plotting Fagan's nomogram, this function also calculates the
#' minimally needed values for specificity and sensitivity to reach desired
#' posttest probabilities (or likelihood ratios) for a grey zone (Coste et al.,
#' 2003, 2006).
#'
#' @param prob.pre.test The prior test probability, with a default value of .5.
#'   Often, (local) prevalence is used.
#' @param probs.post.test A vector of two values that give the desired posttest
#'   probabilities of observing the event in the case of a positive test result
#'   (positive posttest probability: pos), and the posttest probability of
#'   observing the event in the case of a negative test result (negative
#'   posttest probability: neg). When not given, these probabilities are
#'   calculated using the likelihood ratios (LR).
#' @param SeSp A vector of two values that give the desired sensitivity and
#'   specificity. When not given, the Se and Sp values are calculated from the
#'   desired posttest probabilities.
#' @param LR A vector of two values that give the positive likelihood ratio
#'   (sensitivity / (1- specificity)): PLR of observing the event, and the
#'   negative likelihood ratio ((1 - sensitivity) / specificity): NLR of not
#'   observing the event. PLR is a value > 1, NLR is a value between 0 and 1.
#'   When not given, the LR values are calculated from the desired posttest
#'   probabilities.
#' @param plot A Boolean that indicates whether a plot is desired.
#'
#' @return Vector of values: \describe{ \item{$pre: }{The given pre-test #' probability.} \item{$min.LRpos: }{The given or calculated minimally
#'   required positive likelihood ratio. If no value is provided, it is
#'   calculated.} \item{$max.LRneg: }{The given or calculated maximally required #' negative likelihood ratio. If no value is provided, it is calculated.} #' \item{$post.pos: }{The given or calculated positive posttest probability.}
#'   \item{$minSp: }{The minimum value for the specificity, needed to reach the #' desired posttest probabilities.} \item{$minSe: }{The minimum value for the
#'   sensitivity, needed to reach the desired posttest probabilities.} }
#' @details Parameter probs.post.test or SeSp or LR must be supplied, the other
#'   two values are calculated. When more than one parameter is given the other
#'   two are ignored. The basis of this function is adapted from package
#'   TeachingDemos.
#' @references { Fagan, T. J. (1975). Nomogram for Bayes theorem. The New
#'   England Journal of Medicine, 293(5), 257-257.
#'
#'   Coste, J., Jourdain, P., & Pouchot, J. (2006). A gray zone assigned to
#'   inconclusive results of quantitative diagnostic tests: application to the
#'   use of brain natriuretic peptide for diagnosis of heart failure in acute
#'   dyspneic patients. Clinical Chemistry, 52(12), 2229-2235.
#'
#'   Coste, J., & Pouchot, J. (2003). A grey zone for quantitative diagnostic
#'   and screening tests. International Journal of Epidemiology, 32(2),
#'   304-313. }
#' @importFrom graphics axis segments text title

#' @export
#'
#' @examples
#' # Show calculated results (first 3 times about the same)
#' (nomogram(prob.pre.test = .10, probs.post.test=c(pos=.70, neg=.001), plot=FALSE))
#' (nomogram(prob.pre.test = .10, SeSp=c(Se=0.991416309, Sp=0.952789700), plot=FALSE))
#' (nomogram(prob.pre.test = .10, LR=c(pos=21, neg=0.0090090091), plot=FALSE))
#' (nomogram(prob.pre.test = .10, SeSp=c(Se=0.99, Sp=0.95), plot=FALSE))
#' # plot only
#' nomogram(prob.pre.test = .10, LR=c(pos=21, neg=0.0090090091))
#' # plot and display precise results
#' (nomogram(prob.pre.test = .10, probs.post.test=c(pos=.70, neg=.001)))
#'
#' # check the influence of different values of prevalence
#' i=1
#' out=matrix(0,nrow = 9, ncol= 7)
#' for (prev in (seq(.1, .9, by=.1))) {
#'   out[i,]=nomogram(prob.pre.test=prev, probs.post.test=c(.95, .05), plot=FALSE)
#'   i=i+1
#' }
#' colnames(out) = names(nomogram(prob.pre.test=prev, probs.post.test=c(.95, .05), plot=FALSE))
#' out
#'

# probs.post.test=c(pos=NULL, neg=.05)
# SeSp = c(.5, .5); prob.pre.test=.5; probs.post.test=c(pos=NULL, neg=NULL);
# LR=c(PLR=NULL, NLR=NULL); plot=FALSE

nomogram <- function(prob.pre.test=.5, probs.post.test=c(pos=NULL, neg=NULL),
SeSp=c(Se=NULL, Sp=NULL), LR=c(PLR=NULL, NLR=NULL), plot=T) {

if (!is.null(probs.post.test[1])) if (probs.post.test[1] < prob.pre.test)
stop('Invalid assumption. Probs.post.test[1] (pos) must be equal or larger than prob.pre.test.')
if (!is.null(probs.post.test[2])) if (probs.post.test[2] > prob.pre.test)
stop('Invalid assumption. Probs.post.test[2] (neg) must be equal or smaller than prob.pre.test.')

logits <- function(p) log(p/(1 - p))

# LR = c(LR+, LR-)
# calc.LR2 = function(prob.pre.test, probs.post.test){
#   exp(logits(probs.post.test)-logits(prob.pre.test))
# }
calc.LR = function(prob.pre.test, probs.post.test){
(probs.post.test * (1-prob.pre.test)) /
((1-probs.post.test) *prob.pre.test)
}
# prob.pre.test=.1; probs.post.test=.001;
# calc.LR(prob.pre.test, probs.post.test);
# calc.LR2(prob.pre.test, probs.post.test)

calc.Sp <- function(LR)  (1-LR[1])/(LR[2]-LR[1])
calc.Se <- function(LR) LR[1]*(1-(1-LR[1])/(LR[2]-LR[1]))

logits.pre <- logits(prob.pre.test)
Se = SeSp[1]; Sp=SeSp[2]
if (!is.null(probs.post.test)){
LR = calc.LR(prob.pre.test, probs.post.test) # LR positive
logits.post <- log(LR) + logits.pre
Se = calc.Se(LR); Sp = calc.Sp(LR)
} else if (!is.null(SeSp)){
LR = c(PLR=SeSp[1] / (1- SeSp[2]), NLR=(1 - SeSp[1]) / SeSp[2])
logits.post <- log(LR) + logits.pre
probs.post.test <- exp(logits.post)/(1 + exp(logits.post))
} else  if (!is.null(LR)){
logits.post <- log(LR) + logits.pre
probs.post.test <- exp(logits.post)/(1 + exp(logits.post))
Se = calc.Se(LR); Sp = calc.Sp(LR)
} else stop("Parameter probs.post.test or SeSP or LR must be supplied.")

if (any(prob.pre.test > 1) | any(prob.pre.test < 0) |
any(probs.post.test > 1) | any(probs.post.test < 0) |
any(LR < 0) | any(is.infinite(LR)) | any(is.nan(LR))) {
stop("Function call has wrong parameter values.")
}

# if (LR[1]!=LR[2]) {Se = calc.Se(LR); Sp = calc.Sp(LR) }

compl.logit.pre <- logits(1 - prob.pre.test)

if (plot) {
on.exit(par(opar))
par(mar = c(1.5, 6, 2, 6))
LR.vec <- c(0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2,
0.5, 1, 2, 5, 10, 20, 50, 100, 200, 500, 1000)
prob.vec <- c(0.001, 0.002, 0.003, 0.005, 0.007, 0.01, 0.02,
0.03, 0.05, 0.07, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7,
0.8, 0.9, 0.93, 0.95, 0.97, 0.98, 0.99, 0.993, 0.995,
0.997, 0.998, 0.999)
plot(0, 0, type = "n", ylim = range(logits(prob.vec)), axes = FALSE,
xlab = "", ylab = "")
axis(2, rev(logits(prob.vec)), 100 * prob.vec, pos = -1,
las = 1, cex.axis = 0.7)
axis(2, rev(logits(prob.vec)), 100 * prob.vec, pos = -1,
tck = 0.03, labels = FALSE)
axis(4, logits(prob.vec), 100 * prob.vec, pos = 1, las = 1,
cex.axis = 0.7)
axis(4, logits(prob.vec), 100 * prob.vec, pos = 1, tck = 0.03,
labels = FALSE)
axis(2, log(LR.vec[1:10])/2, LR.vec[1:10], pos = 0, las = 1,
cex.axis = 0.7)
axis(2, log(LR.vec[1:10])/2, LR.vec[1:10], pos = 0, tck = 0.03,
labels = FALSE)
axis(4, log(LR.vec[10:19])/2, LR.vec[10:19], pos = 0, las = 1,
cex.axis = 0.7)
axis(4, log(LR.vec[10:19])/2, LR.vec[10:19], pos = 0, tck = 0.03,
labels = FALSE)
text(0, 4.5, "Likelihood ratio", cex = 1.2)
segments(-1, compl.logit.pre, 1, logits.post, lwd = 1.5, col = c(3, 2))
mtext(side = 2, text = "Pretest %", line = 2,
cex = 1.2)
mtext(side = 4, text = "Posttest %", line = 2,
cex = 1.2, las = 3)
title(main = "Fagan's nomogram")
text(0, -5.3, paste("Pre test % of disease =",
round(100 * prob.pre.test, 2), "% \n",
"Likelihood ratios (+;-)=",
round(LR[1], 4),';', round(LR[2], 4),"\n", "posttest %. of disease (+;-) =",
round(100 * probs.post.test[1], 2), ';', round(100 * probs.post.test[2], 4),"%\n",
'Min Se =', round(Se, 2), 'Min Sp =', round(Sp, 2)),
cex = 0.7) }
invisible(c(pre = unname(prob.pre.test),
post.pos= unname(probs.post.test[1]), post.neg = unname(probs.post.test[2]),
min.LRpos = unname(LR[1]), max.LRneg = unname(LR[2]),
minSe = unname(Se), minSp = unname(Sp)))

}


## Try the UncertainInterval package in your browser

Any scripts or data that you put into this service are public.

UncertainInterval documentation built on March 3, 2021, 1:10 a.m.