You tested positive. Now what?
Testing positive for a disease is not exactly good news, but how alarmed you should be ultimately depends on a few factors: the accuracy of the test and the overall rarity of the disease (the so-called base rate). Before you walk into the doctor’s office, the probability that you are sick is \(P(D=+)\), the prevalence. This is the unconditional (or marginal) probability that a random member of the population (ie you in that moment) is sick. After you test positive, the revised or conditional probability that you have the disease is given by Bayes’ theorem:
\[ \begin{aligned} P(D = +\,|\,T=+) &= \frac{{\color{blue}{P(T=+\,|\,D=+)}}{\color{red}{P(D=+)}}}{P(T=+)} \\ &= \frac{{\color{blue}{P(T=+\,|\,D=+)}}{\color{red}{P(D=+)}}}{{\color{blue}{P(T=+\,|\,D=+)}}{\color{red}{P(D=+)}}+[1-{\color{green}{P(T=-\,|\,D=-)}}][1-{\color{red}{P(D=+)}}]} \\ &= \frac{{\color{blue}{\text{TPR}}}\times{\color{red}{\text{PREV}}}}{{\color{blue}{\text{TPR}}}\times{\color{red}{\text{PREV}}}+(1-{\color{green}{\text{TNR}}})(1-{\color{red}{\text{PREV}}})}. \end{aligned} \]
\(P(T=+\,|\,D=+)\) is the true positive rate or sensitivity: the probability that a sick person is correctly diagnosed. \(P(T=-\,|\,D=-)\) is the true negative rate or specificity: the probability that a well person is correctly diagnosed. It is the nonlinear interaction of three numbers (prevalence, TPR, and TNR) that determines the revised likelihood that you have the disease. Two numbers interacting linearly is about all a regular human can handle off the dome, hence people are generally bad at reasoning about this stuff.
The application below visualizes the sample space for our simple model of disease testing. We have \(D\in\{-,,\,+\}\) and \(T\in\{-,\,+\}\), and so there are four individual outcomes:
- \(P(T = +\cap D=+)=P(T=+\,|\,D=+)P(D=+)=\text{TPR}\times P(D=+)\);
- \(P(T = +\cap D=-)=P(T=+\,|\,D=-)P(D=-)=\text{FPR}\times P(D=-)\);
- \(P(T = -\cap D=-)=P(T=-\,|\,D=-)P(D=-)=\text{TNR}\times P(D=-)\);
- \(P(T = -\cap D=+)=P(T=-\,|\,D=+)P(D=+)=\text{FNR}\times P(D=+)\).
The application plots the \([0,\,1]\times[0,\,1]\) unit square, partitioned into four rectangles representing the outcomes listed above. The area of each rectangle is equal to the probability of that outcome.
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 650
library(shiny)
ui <- fluidPage(
# Application title
titlePanel("The sample space of a simple diagnostic test"),
sidebarLayout(
sidebarPanel(
sliderInput("TPR",
"Sensitivity (TPR):",
min = 0,
max = 1,
value = 0.5),
sliderInput("TNR",
"Specificity (TNR):",
min = 0,
max = 1,
value = 0.5),
sliderInput("p",
"Prevalence:",
min = 0,
max = 1,
value = 0.5)
),
mainPanel(
plotOutput("distPlot"),
# ---- Legend goes here ----
tags$div(
style = "margin-top:20px;",
tags$h4("Legend"),
tags$div(style="display:flex; align-items:center; margin-bottom:5px;",
tags$div(style="width:20px; height:20px; background-color:rgba(255,0,0,0.25); margin-right:10px;"),
"True Positive (T=+, D=+)"
),
tags$div(style="display:flex; align-items:center; margin-bottom:5px;",
tags$div(style="width:20px; height:20px; background-color:rgba(0,255,0,0.25); margin-right:10px;"),
"False Negative (T=-, D=+)"
),
tags$div(style="display:flex; align-items:center; margin-bottom:5px;",
tags$div(style="width:20px; height:20px; background-color:rgba(0,0,255,0.25); margin-right:10px;"),
"True Negative (T=-, D=-)"
),
tags$div(style="display:flex; align-items:center; margin-bottom:5px;",
tags$div(style="width:20px; height:20px; background-color:rgba(255,165,0,0.25); margin-right:10px;"),
"False Positive (T=+, D=-)"
)
)
)
)
)
server <- function(input, output) {
output$distPlot <- renderPlot({
p <- input$p
TNR <-input$TNR
TPR <- input$TPR
FNR <- 1 - TPR
FPR <- 1 - TNR
par(mar = c(8, 8, 1, 8))
plot(0, col = "white", xlim = c(0, 1), ylim = c(0, 1),
xaxt = "n", yaxt = "n", bty = "n", xlab = "", ylab = "")
polygon(c(0, p, p, 0), c(FNR, FNR, 1.0, 1.0), col = rgb(1, 0, 0, alpha = 0.25), border = NA)
polygon(c(0, p, p, 0), c(0, 0, FNR, FNR), col = rgb(0, 1, 0, alpha = 0.25), border = NA)
polygon(c(p, 1, 1, p), c(0, 0, TNR, TNR), col = rgb(0, 0, 1, alpha = 0.25), border = NA)
polygon(c(p, 1, 1, p), c(TNR, TNR, 1, 1), col = rgb(1, 0.647, 0, alpha = 0.25), border = NA)
mtext("TNR", side = 4, at = TNR / 2, las = 1)
mtext("FPR", side = 4, at = (1 + TNR) / 2, las = 1)
mtext("TPR", side = 2, at = (1 + FNR) / 2, las = 1, line = 2)
mtext("FNR", side = 2, at = FNR / 2, las = 1, line = 2)
mtext("p(D = +)", side = 1, at = p / 2, las = 1, line = 2)
mtext("p(D = -)", side = 1, at = (1 + p) / 2, las = 1, line = 2)
mtext(c("0", "1"), side = 1, at = c(0, 1), las = 1)
mtext(c("0", "1"), side = 2, at = c(0, 1), las = 1)
})
}
shinyApp(ui = ui, server = server)
The application below displays how testing positive shifts the probability of disease from \(P(D=+)\) to \(P(D = +\,|\,T=+)\). The size of the revision ultimately depends on how big \(P(D=+)\) was to begin with, and how good the test is (TPR and TNR). Play around with different values and see what happens:
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 500
library(shiny)
# Define UI for application that draws a histogram
ui <- fluidPage(
# Application title
titlePanel("What can a diagnostic test teach you?"),
# Sidebar with a slider input for number of bins
sidebarLayout(
sidebarPanel(
sliderInput("TPR",
"Sensitivity (TPR):",
min = 0,
max = 1,
value = 0.5),
sliderInput("TNR",
"Specificity (TNR):",
min = 0,
max = 1,
value = 0.5)
),
# Show a plot of the generated distribution
mainPanel(
plotOutput("distPlot")
)
)
)
# Define server logic required to draw a histogram
server <- function(input, output) {
output$distPlot <- renderPlot({
TPR <- input$TPR
TNR <- input$TNR
P <- seq(0, 1, length.out = 1000)
B <- TPR * P / (TPR * P + (1 - TNR) * (1 - P))
# draw the histogram with the specified number of bins
par(mar = c(5, 5, 5, 0.5))
plot(P, B, type = "l", xlim = c(0, 1), ylim = c(0, 1),
xlab = "P(D = +)", ylab = "P(D = + | T = +)",
lwd = 2, cex.lab = 2,
panel.first = abline(a = 0, b = 1, lty = 2, lwd = 2, col = "grey"))
})
}
# Run the application
shinyApp(ui = ui, server = server)