A clairvoyant mouse?

Let’s think about it another way

Let’s think about it another way

But while you’ve been away studying, UCSF has announced a new policy requiring all graduate students to be tested for Covid-19 daily. You must provide a negative test result before entering any research building on-campus!

SARS-CoV-2 mRNA:

The scenario

N.B.: SF no longer tracks Covid case rates, so I’ll be using last year’s numbers.

The scenario

• To get to prevalence, let’s assume everyone infected has Covid for 14 days.

Recall

The \(p\)-value is defined by thinking about what our observations would be by chance if we presuppose that the null hypothesis is in fact true:

Usually, diagnostic tests are evaluated kind of like \(p\)-values: we use measures that ask how the test behaves when we pre-suppose the truth.

N.B.: The \(p\)-value can be thought of as

\[ p = \mathrm{Pr}(\textrm{test }+ \mid \textrm{actually }-) \]

Because there are only two test outcomes, \(+\) and \(-\), this means that

\[\begin{eqnarray*} p & = & 1 - \mathrm{Pr}(\textrm{test }- \mid \textrm{actually }-) \\ & = & 1 - \mathrm{specificity} \end{eqnarray*}\]

So, the \(p\)-value is in the same “family” of measures about a test as sensitivity and specificity.

Back to the code

# September 2021
# rtpcr_sensitivity = 0.777
# rtpcr_specificity = 0.988

# August 2022
# https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7350782/
rtpcr_sensitivity = 0.733
# https://www.uptodate.com/contents/covid-19-diagnosis
rtpcr_specificity = 0.97

# September 2023: These characteristics haven't appreciably changed

# September 2024: "

You already know that you’ve tested positive. What you would like to infer is whether or not you actually have Covid.

Positive predictive value is the probability one actually has the disease given a positive test result:

\[ \mathrm{PPV} := \mathrm{Pr}(\textrm{actually }+ \mid \textrm{test }+) \]

N.B.: Always write docstrings for your code:

help( sample_results )

Help on function sample_results in module __main__:

sample_results(n, prevalence, sensitivity, specificity)
    Sample disease presence and test results with the given characteristics

    Arguments:
    n -- The number of individiuals to sample
    prevalence -- The probability of finding the actual disease in the population
    sensitivity -- The probability of testing positive given one has the disease
    specificity -- The probability of testing negative given one does not have the disease

    Reutrns:
    has_disease -- Shape (n,), whether each individual has the disease
    test_results -- Shape (n,), whether each individual tests positive for the disease

Let’s simulate one day of Covid tests at UCSF:

# Approximately the number of employees at UCSF
n_employees = 24000
# Run our simulation and grab the results
# (See code for the slides online)
has_covid, rtpcr_positive = sample_results( n_employees,
                                            covid_prevalence,
                                            rtpcr_sensitivity,
                                            rtpcr_specificity )
# If I put an f in front of my string, I can use curly braces to
# put Python code inside of it
print( f'> Cases today: **{np.sum( has_covid )}**\\' )
print( f'> Positive tests today: **{np.sum( rtpcr_positive )}**' )

As is always good practice, let’s take a look at a subsample of the raw data:

# Subplots allows us to lay out some axes (`ax`) inside of a
# figure (`fig`) that we can manipulate later
fig, ax = plt.subplots( figsize = (12, 1.2) )

# Use the plotting function we just made on these axes
# (See code for the slides online)
sample_plot( ax, has_covid, rtpcr_positive )

# We're all done building the plots, we want to displaly them now
plt.show()

What are the PPV and NPV?

# (See code for the slides online)
ppv_naive, npv_naive = predictive_value( has_covid, rtpcr_positive )

# We can do math in our curly braces; here I convert the ratios into
# percents. The `:0.2f` at the end tells Python that I want 2 digits
# after the decimal point (not at all intuitive; blame the C people
# who wrote the original `printf` function).
print( f'> Naive PPV: **{ppv_naive * 100:0.2f}%**\\' )
print( f'> Naive NPV: **{npv_naive * 100:0.2f}%**' )

Let’s talk about the good news first:

Do you think this will scale well?

# 10 days is the policy, but a couple of those will be a weekend
quarantine_workdays = 8
# False positives *do not* have Covid *and* test positive
quarantine_fp_per_day = np.sum( ~has_covid & rtpcr_positive )
quarantine_fp_simultaneous = quarantine_fp_per_day * quarantine_workdays
quarantine_fraction = quarantine_fp_simultaneous / n_employees

print( f'> Fraction of workforce quarantined: **{quarantine_fraction * 100:0.2f}%**' )

But, n.b.! This test also has an extremely high accuracy—that is, the probability of the test being correct:

# == operates element by element on arrays
# `a == b` is an array where each element is `True` if the
# corresponding elements of `a` and `b` are equal
accuracy_naive = np.sum( has_covid == rtpcr_positive ) / len( has_covid )

print( f'> Naive accuracy: **{accuracy_naive * 100:0.2f}%**' )

A better way to use Covid tests

In our simulation, the prevalence of Covid was thought of as the probability that someone in our sample had Covid. We took this to be the prevalence of Covid in San Francisco, under the assumption that we were testing people essentially randomly, and had no criteria for how we were selecting people to test.

A better way to use Covid tests

After briefly quarantining a substantial portion of its workforce, UCSF has decided to try a new Covid testing schema.

A better way to use Covid tests

Before, we had no clue about the people being tested. This new criterion means we have more knowledge (of the exposure).

A better way to use Covid tests

All of the characteristics of the test itself are still exactly the same—the same reagents, the same technique, everything.

The only thing we’ll change to simulate this new scenario is the proportion of tested people who have Covid.

covid_belief_sus = 0.05

Let’s run a simulation of a group of employees with suspected Covid exposure:

n_sus = 1000
has_covid_sus, rtpcr_positive_sus = sample_results( n_sus,
                                                    covid_belief_sus,
                                                    rtpcr_sensitivity,  # Same as before!
                                                    rtpcr_specificity ) # Same as before!

# Python is fully Unicode compatible for strings, so there's no
# reason not to include emoji 🤪
print( f'> Cases among the suspicious 🤔: **{np.sum( has_covid_sus )}**\\' )
print( f'> Positive tests among the suspicious 🤔: **{np.sum( rtpcr_positive_sus )}**' )

# You can't use emoji in variable names, though 😢
# Got to switch to Haskell for that ...

Dummy check, round 2

fig, ax = plt.subplots( figsize = (12, 1.2) )
sample_plot( ax, has_covid_sus, rtpcr_positive_sus )
plt.show()

Already this is looking much more reasonable. What are the PPV and the NPV?

ppv_sus, npv_sus = predictive_value( has_covid_sus, rtpcr_positive_sus )

print( f'> Suspicious PPV: **{ppv_sus * 100:0.2f}%**\\' )
print( f'> Suspicious NPV: **{npv_sus * 100:0.2f}%**' )

Aha!

In this example, we saw that the PPV of our Covid test depended strongly on the scheme we used to select which people to test, even while the test itself remained entirely the same.

What does a test do?

Before we ran the test, we had some prior knowledge about whether we had Covid:

\[ \mathrm{Pr}(\textrm{actually }+) \]

What does a test do?

When we get new knowledge—the result of the test—we update our belief of whether we have Covid. Our knowledge about whether we have Covid after the update caused by the test is precisely the PPV:

\[\begin{eqnarray*} &\mathrm{Pr}(\textrm{actually }+) & \\ & \downarrow & \\ & \mathrm{Pr}(\textrm{actually }+ \mid \textrm{test }+) & \end{eqnarray*}\]

This is also known as the posterior probability (as in, the probability that we have after) given a positive test.

Tests in science

For Covid testing, we obtained some bit of noisy data (a positive test result), and sought to infer something about what was actually happening (whether we truly have Covid).

We made this inference quantitative using the Covid RT-PCR test’s PPV:

\[ \mathrm{PPV} = \mathrm{Pr}(\textrm{actually }+ \mid \textrm{test }+) \]

Tests in science

In science, we use statistical tests.

Just like the Covid RT-PCR, statistical tests are noisy bits of data. What we would like to do is to infer something about what is actually happening in the world on the basis of these test results.

Example—Correlation

Let’s say that you’re analyzing some data. You plug two variables you’re working on, x and y, into a magic black box like scipy.stats.pearsonr, and it pops out:

>>> scipy.stats.pearsonr( x, y )
(0.21365304326850618, 0.03281172835181021)

Example—Correlation

Given this:

>>> scipy.stats.pearsonr( x, y )
(0.21365304326850618, 0.03281172835181021)

… what do you now know about whether there is actually a relationship between x and y?

Example—Correlation

>>> scipy.stats.pearsonr( x, y )
(0.21365304326850618, 0.03281172835181021)

This information, by itself, is not a statement about whether the null hypothesis is actually true or false—just as a positive Covid test is not, by itself, a statement about whether you actually have Covid.

Example—Correlation

>>> scipy.stats.pearsonr( x, y )
(0.21365304326850618, 0.03281172835181021)

With our shiny positive test result, what we really care about is the posterior probability, given we saw this test result, of whether there is a relationship between x and y. This is precisely the PPV:

\[ \mathrm{PPV} := \mathrm{Pr}(\textrm{actually a relationship} \mid \textrm{test rejects null hyp.}) \]

Example—Correlation

>>> scipy.stats.pearsonr( x, y )
(0.21365304326850618, 0.03281172835181021)

\[ \mathrm{PPV} := \mathrm{Pr}(\textrm{actually a relationship} \mid \textrm{test rejects null hyp.}) \]

But, just as before with Covid RT-PCR,

PPV depends on our prior knowledge of whether there is a real relationship between x and y—

not just the characteristics of pearsonr!

Example—Correlation

>>> scipy.stats.pearsonr( x, y )
(0.21365304326850618, 0.03281172835181021)

Suppose I got this result by going out into the world and running the pearsonr function between every pair of two datasets I could get my hands on.

Example—Correlation

>>> scipy.stats.pearsonr( x, y )
(0.21365304326850618, 0.03281172835181021)

What is the probability that any two of those randomly chosen datasets are actually causally related to one another? That is, what is the prior,

\[ \mathrm{Pr}(\textrm{actually a relationship}) \]

It is astronomically small.

Example—Correlation

>>> scipy.stats.pearsonr( x, y )
(0.21365304326850618, 0.03281172835181021)

In this experimental setup of randomly testing correlations, what is the PPV?

Example—Correlation

>>> scipy.stats.pearsonr( x, y )
(0.21365304326850618, 0.03281172835181021)

So, how do I interpret this result I just saw, that pearsonr rejected the null hypothesis of no relationship between x and y?

Recall

When looking at Covid testing, increasing our prior belief of Covid infection to just 5% raised the PPV from 5% to 60%, for the exact same test. To make the test’s results more meaningful, we used knowledge of prior exposure to Covid to choose who to test.

A clairvoyant mouse?

You look at your data, and you see that you reject the null hypothesis that manipulated mice do not predict the future with a \(p\)-value of 0.014.

A clairvoyant mouse?

You look at your data, and you see that you reject the null hypothesis that manipulated mice do not predict the future with a \(p\)-value of 0.014.

What do you do with this information?

Statistical testing is not just a matter of selecting the “correct” test for a given question.