diagstats(): calculate confidence intervals for PPV, NPV, and PD (#2)

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: diagmeta
 Title: Meta-Analysis of Diagnostic Accuracy Studies with Several Cutpoints
 Version: 0.6-0
-Date: 2024-09-10
+Date: 2024-09-17
 Depends: meta (>= 5.0-0)
 Imports: lme4, grDevices
 Authors@R: c(person("Gerta", "Rücker",

R/diagstats.R

@@ -8,39 +8,43 @@
 #' predictive values (NPV), and probabilities of disease (PD) are
 #' calculated if the prevalence is provided.
 #' 
-#' @param x An object of class \code{diagmeta}
-#' @param cutoff A numeric or vector with cutoff value(s)
-#' @param sens A numeric or vector with sensitivity value(s)
-#' @param spec A numeric or vector with specificity value(s)
-#' @param prevalence A numeric or vector with the prevalence(s)
-#' @param level The level used to calculate confidence intervals
+#' @param x An object of class \code{diagmeta}.
+#' @param cutoff A numeric or vector with cutoff value(s).
+#' @param sens A numeric or vector with sensitivity value(s).
+#' @param spec A numeric or vector with specificity value(s).
+#' @param prevalence A numeric or vector with the prevalence(s).
+#' @param level The level used to calculate confidence intervals.
 #' 
 #' @return
 #' A data frame of class "diagstats" with the following variables:
-#' \item{cutoff}{Cutoffs provided in argument "cutoff" and / or
+#' \item{cutoff}{Cutoffs provided in argument \code{cutoff} and / or
 #'   model-based cutoff values for given sensitivities /
 #'   specificities.}
-#' \item{Sens}{Sensitivities provided in argument "sens" and / or
+#' \item{Sens}{Sensitivities provided in argument \code{sens} and / or
 #'   model-based estimates of the sensitivity for given cutoffs /
-#'   specificities}
-#' \item{seSens}{Standard error of sensitivity}
+#'   specificities.}
 #' \item{lower.Sens, upper.Sens}{Lower and upper confidence limits of
-#'   the sensitivity}
-#' \item{Spec}{Specificities provided in argument "spec" and / or
+#'   the sensitivities.}
+#' \item{Spec}{Specificities provided in argument \code{spec} and / or
 #'   model-based estimates of the specificity for given cutoffs /
-#'   sensitivities}
-#' \item{seSpec}{Standard error of specificity}
+#'   sensitivities.}
 #' \item{lower.Spec, upper.Spec}{Lower and upper confidence limits of
-#'   the specificity}
+#'   the specificities.}
 #' \item{prevalence}{As defined above.}
-#' \item{PPV}{Positive predictive value (based on the cutoff)}
-#' \item{NPV}{Negative predictive value (based on the cutoff)}
+#' \item{PPV}{Positive predictive value (based on the prevalence).}
+#' \item{lower.PPV, upper.PPV}{Lower and upper confidence limits of
+#'   positive predictive values.}
+#' \item{NPV}{Negative predictive value (based on the prevalence)}
+#' \item{lower.NPV, upper.NPV}{Lower and upper confidence limits of
+#'   negative predictive values.}
 #' \item{PD}{Probability of disease if the given cutoff value was
-#'   observed as the measurement for an individual}
+#'   observed as the measurement for an individual.}
+#' \item{lower.PD, upper.PD}{Lower and upper confidence limits of
+#'   probabilities of disease.}
 #' \item{dens.nondiseased}{Value of the model-based density function at the
-#'    cutoff(s) for non-diseased individuals}
+#'    cutoff(s) for non-diseased individuals.}
 #' \item{dens.diseased}{Value of the model-based density function at the
-#'    cutoff(s) for diseased individuals}
+#'    cutoff(s) for diseased individuals.}
 #' 
 #' @author
 #' Gerta Rücker \email{gerta.ruecker@@uniklinik-freiburg.de},
@@ -115,7 +119,18 @@ diagstats <- function(x,
   alpha1 <- regr$alpha1
   beta0 <- regr$beta0
   beta1 <- regr$beta1
-  ##
+  #
+  var.alpha0 <- regr$var.alpha0
+  var.alpha1 <- regr$var.alpha1
+  var.beta0 <- regr$var.beta0
+  var.beta1 <- regr$var.beta1
+  cov.alpha0.alpha1 <- regr$cov.alpha0.alpha1
+  cov.alpha0.beta0 <- regr$cov.alpha0.beta0
+  cov.alpha0.beta1 <- regr$cov.alpha0.beta1
+  cov.alpha1.beta0 <- regr$cov.alpha1.beta0 
+  cov.alpha1.beta1 <- regr$cov.alpha1.beta1
+  cov.beta0.beta1 <- regr$cov.beta0.beta1  
+  #
   distr <- x$distr
 
 
@@ -206,27 +221,83 @@ diagstats <- function(x,
   ##
   PPV <- Sens * prevalence /
     (prevalence * Sens + (1 - prevalence) * (1 - Spec))
-  ##
+  varPPV <-
+    Spec^2 * (var.alpha0 + cutoff.tr^2 * var.beta0 +
+                2 * cutoff.tr * cov.alpha0.beta0 + x$var.nondiseased) +
+    (1 - Sens)^2 * (var.alpha1 + cutoff.tr^2 * var.beta1 +
+                      2 * cutoff.tr * cov.alpha1.beta1 + x$var.diseased) -
+    2 * Spec * (1 - Sens) *
+    (cov.alpha0.alpha1 + cutoff.tr * cov.alpha0.beta1 +
+       cutoff.tr * cov.alpha1.beta0 + cutoff.tr^2 * cov.beta0.beta1)
+  ciPPV <- ci(qdiag(PPV, distr), sqrt(varPPV), level = level)
+  lower.PPV <- pdiag(ciPPV$lower, distr)
+  upper.PPV <- pdiag(ciPPV$upper, distr)  
+  #
   NPV <- Spec * (1 - prevalence) /
     (prevalence * (1 - Sens) + (1 - prevalence) * Spec)
-  ##
+  varNPV <-
+    (1 - Spec)^2 * (var.alpha0 + cutoff.tr^2 * var.beta0 +
+                      2 * cutoff.tr * cov.alpha0.beta0 + x$var.nondiseased) +
+    Sens^2 * (var.alpha1 + cutoff.tr^2 * var.beta1 +
+                2 * cutoff.tr * cov.alpha1.beta1 + x$var.diseased) -
+    2 * Sens * (1 - Spec) *
+    (cov.alpha0.alpha1 + cutoff.tr * cov.alpha0.beta1 +
+       cutoff.tr * cov.alpha1.beta0 + cutoff.tr^2 * cov.beta0.beta1)
+  ciNPV <- ci(qdiag(NPV, distr), sqrt(varNPV), level = level)
+  lower.NPV <- pdiag(ciNPV$lower, distr)
+  upper.NPV <- pdiag(ciNPV$upper, distr)  
+  #
   PD <- prevalence * dens.diseased /
     (prevalence * dens.diseased + (1 - prevalence) * dens.nondiseased)
+  varPD <- 
+    (1 - 2 * Spec)^2 *
+    (var.alpha0 + cutoff.tr^2 * var.beta0 +
+       2 * cutoff.tr * cov.alpha0.beta0 + x$var.nondiseased) +
+    (2 * Sens - 1)^2 *
+    (var.alpha1 + cutoff.tr^2 * var.beta1 +
+       2 * cutoff.tr * cov.alpha1.beta1 + x$var.diseased) -
+    2 * (2 * Sens - 1) * (1 - 2 * Spec) *
+    (cov.alpha0.alpha1 + cutoff.tr * cov.alpha0.beta1 +
+       cutoff.tr * cov.alpha1.beta0 + cutoff.tr^2 * cov.beta0.beta1)
+  ciPD <- ci(qdiag(PD, distr), sqrt(varPD), level = level)
+  lower.PD <- pdiag(ciPD$lower, distr)
+  upper.PD <- pdiag(ciPD$upper, distr)  
+  #
+  DOR <- Sens * Spec / (1 - Sens) / (1 - Spec)
+  varDOR <- 
+    var.alpha0 + cutoff.tr^2 * var.beta0 +
+    2 * cutoff.tr * cov.alpha0.beta0 + x$var.nondiseased +
+    var.alpha1 + cutoff.tr^2 * var.beta1 +
+    2 * cutoff.tr * cov.alpha1.beta1 + x$var.diseased - 
+    2 * cov.alpha0.alpha1 -
+    2 * cutoff.tr * cov.alpha0.beta1 - 
+    2 * cutoff.tr * cov.alpha1.beta0 -
+    2 * cutoff.tr^2 * cov.beta0.beta1
+  ciDOR <- ci(log(DOR), sqrt(varDOR), level = level)
+  lower.DOR <- exp(ciDOR$lower)
+  upper.DOR <- exp(ciDOR$upper)
+  #
+  LRpos <- Sens / (1 - Spec)
+  ci.LRpos <- ci(log(LRpos), sqrt(varPPV), level = level)
+  lower.LRpos <- exp(ci.LRpos$lower)
+  upper.LRpos <- exp(ci.LRpos$upper)
+  #
+  LRneg <- (1 - Sens) / Spec
+  ci.LRneg <- ci(log(LRneg), sqrt(varNPV), level = level)
+  lower.LRneg <- exp(ci.LRneg$lower)
+  upper.LRneg <- exp(ci.LRneg$upper)
 
 
-  res <- data.frame(cutoff = cutoff,
-                    Sens = Sens, seSens = seSens,
-                    lower.Sens = lower.Sens, upper.Sens = upper.Sens,
-                    Spec = Spec, seSpec = seSpec,
-                    lower.Spec = lower.Spec, upper.Spec = upper.Spec,
-                    LRpos = Sens / (1 - Spec),
-                    LRneg = (1 - Sens) / Spec,
-                    prevalence = prevalence,
-                    PPV = PPV, NPV = NPV, PD = PD,
-                    dens.nondiseased = dens.nondiseased,
-                    dens.diseased = dens.diseased)
-  #
-  #res <- res[order(res$cutoff), ]
+  res <- data.frame(cutoff,
+                    Sens, lower.Sens, upper.Sens,
+                    Spec, lower.Spec, upper.Spec,
+                    LRpos, lower.LRpos, upper.LRpos,
+                    LRneg, lower.LRneg, upper.LRneg,
+                    prevalence,
+                    PPV,lower.PPV, upper.PPV,
+                    NPV, lower.NPV, upper.NPV,
+                    PD, lower.PD, upper.PD,
+                    dens.nondiseased, dens.diseased)
   #
   class(res) <- c("diagstats", "data.frame")
   #

R/ipd2diag.R

@@ -4,12 +4,12 @@
 #' Function to transform individual patient data (IPD) to enter them
 #' into \code{diagmeta}
 #' 
-#' @param studlab A vector with study labels
+#' @param studlab A vector with study labels.
 #' @param value A vector with individual patients' measurements of a
-#'   discrete or continuous variable
+#'   discrete or continuous variable.
 #' @param status A vector with information of the individual's status
-#'   (0 = non-diseased, 1 = diseased)
-#' @param data An optional data frame containing the study information
+#'   (0 = non-diseased, 1 = diseased).
+#' @param data An optional data frame containing the study information.
 #' @param direction A character string specifying whether the probability of
 #'   the target condition (e.g., a disease) is \code{"increasing"} or
 #'   \code{"decreasing"} with higher values of the biomarker, can be

R/plot.diagmeta.R

@@ -5,74 +5,74 @@
 #' Provides several plots for meta-analysis of diagnostic test
 #' accuracy studies with the multiple cutoffs model
 #' 
-#' @param x An object of class \code{diagmeta}
+#' @param x An object of class \code{diagmeta}.
 #' @param which A character vector indicating the type of plot, either
 #'   \code{"regression"} or \code{"cdf"} or \code{"survival"} or
 #'   \code{"Youden"} or \code{"ROC"} or \code{"SROC"} or
-#'   \code{"density"} or \code{"sensspec"}, can be abbreviated
-#' @param xlab An x axis label
-#' @param main A logical indicating title to the plot
+#'   \code{"density"} or \code{"sensspec"}, can be abbreviated.
+#' @param xlab An x axis label.
+#' @param main A logical indicating title to the plot.
 #' @param ci A logical indicating whether confidence intervals should
 #'   be plotted for \code{"regression"}, \code{"cdf"},
-#'   \code{"survival"}, \code{"Youden"}, and \code{"sensspec"}
+#'   \code{"survival"}, \code{"Youden"}, and \code{"sensspec"}.
 #' @param ciSens A logical indicating whether confidence intervals
 #'   should be plotted for sensitivity, given the specificity in
-#'   \code{"SROC"} plot
+#'   \code{"SROC"} plot.
 #' @param ciSpec A logical indicating whether confidence intervals
 #'   should be plotted for specificity, given the sensitivity in
-#'   \code{"SROC"} plot
+#'   \code{"SROC"} plot.
 #' @param mark.optcut A logical indicating whether the optimal cutoff
-#'   should be marked on \code{"SROC"} plot
+#'   should be marked on \code{"SROC"} plot.
 #' @param mark.cutpoints A logical indicating whether the given
-#'   cutoffs should be marked on \code{"SROC"} plot
+#'   cutoffs should be marked on \code{"SROC"} plot.
 #' @param points A logical indicating whether points should be plotted
 #'   in plots \code{"regression"}, \code{"cdf"}, \code{"survival"},
-#'   \code{"Youden"}, \code{"ROC"}, and \code{"sensspec"}
+#'   \code{"Youden"}, \code{"ROC"}, and \code{"sensspec"}.
 #' @param lines A logical indicating whether polygonal lines
 #'   connecting points belonging to the same study should be printed
 #'   in plots \code{"regression"}, \code{"cdf"}, \code{"survival"},
-#'   \code{"Youden"}, and \code{"sensspec"}
+#'   \code{"Youden"}, and \code{"sensspec"}.
 #' @param rlines A logical indicating whether regression lines or
 #'   curves should be plotted for plots \code{"regression"},
 #'   \code{"cdf"}, \code{"survival"}, \code{"Youden"}, and
-#'   \code{"sensspec"}
+#'   \code{"sensspec"}.
 #' @param line.optcut A logical indicating whether a vertical line
 #'   should be plotted at the optimal cutoff line for plots
 #'   \code{"cdf"}, \code{"survival"}, \code{"Youden"}, and
-#'   \code{"density"}
+#'   \code{"density"}.
 #' @param col.points A character string indicating color of points,
 #'   either \code{"rainbow"}, \code{"topo"}, \code{"heat"},
 #'   \code{"terrain"}, \code{"cm"}, \code{"grayscale"}, or any color
-#'   defined in \code{\link[grDevices]{colours}}
+#'   defined in \code{\link[grDevices]{colours}}.
 #' @param cex A numeric indicating magnification to be used for
-#'   plotting text and symbols
-#' @param pch.points A numeric indicating plot symbol(s) for points
-#' @param col A character string indicating color of lines
+#'   plotting text and symbols.
+#' @param pch.points A numeric indicating plot symbol(s) for points.
+#' @param col A character string indicating color of lines.
 #' @param col.ci A character string indicating color of confidence
-#'   lines
+#'   lines.
 #' @param col.optcut A character string indicating color of optimal
-#'   cutoff line
+#'   cutoff line.
 #' @param cex.marks A numeric indicating magnification(s) to be used
-#'   for marking cutoffs
-#' @param lwd A numeric indicating line width
-#' @param lwd.ci A numeric indicating line width of confidence lines
-#' @param lwd.optcut A numeric indicating line width of optimal cutoff
+#'   for marking cutoffs.
+#' @param lwd A numeric indicating line width.
+#' @param lwd.ci A numeric indicating line width of confidence lines.
+#' @param lwd.optcut A numeric indicating line width of optimal cutoff.
 #' @param lwd.study A numeric indicating line width of individual
-#'   studies
+#'   studies.
 #' @param shading A character indicating shading and hatching
 #'   confidence region in \code{"SROC"} plot, either \code{"none"} or
-#'   \code{"shade"} or \code{"hatch"}
+#'   \code{"shade"} or \code{"hatch"}.
 #' @param col.hatching A character string indicating color used in
-#'   hatching of the confidence region
+#'   hatching of the confidence region.
 #' @param lwd.hatching A numeric indicating line width used in hatching
-#'   of the confidence region
+#'   of the confidence region.
 #' @param ellipse A logical indicating whether a confidence ellipse
-#'   should be drawn around the optimal cutoff
+#'   should be drawn around the optimal cutoff.
 #' @param xlim A character or numerical vector indicating the minimum
-#'   and maximum value for the horizontal axes
+#'   and maximum value for the horizontal axes.
 #' @param ylim A numerical vector indicating the minimum and maximum value for
-#'   the vertical axes
-#' @param \dots Additional graphical arguments
+#'   the vertical axes.
+#' @param \dots Additional graphical arguments.
 #' 
 #' @details
 #' The first argument of the plot function is an object of class

R/print.diagstats.R

@@ -60,11 +60,15 @@ print.diagstats <- function(x,
   drop.names <- c()
   ##
   if (!sensspec)
-    drop.names <- c("Sens", "seSens", "lower.Sens", "upper.Sens",
-                    "Spec", "seSpec", "lower.Spec", "upper.Spec")
+    drop.names <- c("Sens", "lower.Sens", "upper.Sens",
+                    "Spec", "lower.Spec", "upper.Spec")
   ##
   if (all(is.na(x$prevalence)) | !predicted)
-    drop.names <- c(drop.names, c("prevalence", "PPV", "NPV", "PD"))
+    drop.names <-
+      c(drop.names, c("prevalence",
+                      "PPV", "lower.PPV", "upper.PPV",
+                      "NPV", "lower.NPV", "upper.NPV",
+                      "PD", "lower.PD", "upper.PD"))
   ##
   if (!density)
     drop.names <- c(drop.names, c("dens.nondiseased", "dens.diseased"))