The trouble with p-values

Irreproducibility in Science

Eduardo Serna, Wolfram Research

Join the Conversation #WolframTechConf

Talk structure

Motivation

Introduction

p-values are uniform

Statistical power

Publication bias

Example with LocationTest

Motivation

History: Computations of p-values date back to the 1770s, when they were calculated by Laplace.

“In the 1770s Laplace considered the statistics of almost half a million births. The statistics showed an excess of boys compared to girls. He concluded the excess was significant.”

p-values have acquired a bad reputation, they can be confusing and misleading.

A reproducibility crisis in science:

Reproducibility Project 2015: only 36% of studies in psychology meta-study were reproducible.

Ionnadis 2005: “Why Most Published Research Findings Are False”

Motivation - 2

Two frameworks for data analysis: Bayesian vs Frequentist.

Tricky topic: this is a very shallow opinionated explainer:

Grid[{{"Bayesian","Frequentist"},{HoldForm@Probability[model\[Conditioned]data],HoldForm@Probability[data\[Conditioned]model]}},Frame->All]

Bayesian is conceptually easier but often harder in practice, makes more sense with less data.

Frequentist is easier to compute, makes sense with more data, or with simulations.

There is no magic recipe in data analysis, you need good data and good assumptions for robust conclusions.

Introduction

We want to draw inferences from our data.

a simple question: does the data support a hypothesis?

First step: summarise data with a number (statistic)

e.g.

t == (n - 1) * Variance[{(Subscript[x, 1]),...,(Subscript[x, n])}]

Second step: use this statistic to try and determine whether the data supports our hypothesis.

A definition:

a p-value is the probability the sample statistic is smaller than the observed statistic ASSUMING(\[Conditioned]) our model (Subscript[ℋ, 0]).

pvalue == Probability[t ≤ x, t \[Distributed] (Subscript[ℋ, 0])].

To compute this probability, we must know the distribution of our statistic.

e.g. if our data is normally distributed, the statistic t will have a ChiSquareDistribution.

Example: deriving distribution of statistic:

n=10;
vars=Array[a,10];

TransformedDistribution[(x-μ)/σ,x\[Distributed]NormalDistribution[μ,σ]]

TransformedDistribution[x^2,x\[Distributed]NormalDistribution[0,1]]

TransformedDistribution[Total[vars],vars\[Distributed]ProductDistribution[{ChiSquareDistribution[1],n}]]

WL provides many built in tests that automatically compute a statistic and know about its underlying distribution to compute the p-value.

a modified version of the test statistic t is used in PearsonChiSquareTest.

Usage example

Let’s compute a p-value using the Wolfram Language:

data=RandomVariate[MaxwellDistribution[3],10^4];

Histogram[data,{.5}]

HypothesisTestData is an object that allows us to compute different properties from a test:

htd=DistributionFitTest[data,MaxwellDistribution[a],"HypothesisTestData"]

Is the data distributed as MaxwellDistribution?

htd[{"PValueTable","TestConclusion"}]

p-values are uniformly distributed

Intuitively, Probability[t ≤ x, t \[Distributed] (Subscript[ℋ, 0])] gives the position of x amongst a sorted sample from (Subscript[ℋ, 0]), normalized by the total number of samples. If x \[Distributed] (Subscript[ℋ, 0])  this position will be distributed uniformly:

ℋ0=EmpiricalDistribution[Array[a,5]]

Probability[t ≤ x, t \[Distributed] ℋ0]

The CDF is an efficient way of computing Probability[t ≤ x, t \[Distributed] (Subscript[ℋ, 0])]:

CDF[ℋ0,x]==Probability[t≤x,t\[Distributed]ℋ0]//Simplify

𝒟=EmpiricalDistribution[sample=RandomVariate[ParetoDistribution[1,1],10^6]]

Subsampled data has a uniformly distributed empirical p-value:

subsample=RandomSample[sample,10^4];

Histogram[CDF[𝒟,subsample]]

p-values are uniformly distributed - 2

We reject the Null hypothesis if the p-value is smaller than a SignificanceLevel:

this SignificanceLevel has traditionally been set at α = 5%

SignificanceLevel/.Options[DistributionFitTest]
UnitConvert[%,"Percent"]

Crucially, this means that we will have a false rejection (or false POSITIVE depending on the convention) rate of α.

Let’s compute the p-value of a collection of samples:

dist=NormalDistribution[];
sample=RandomVariate[dist,{10^4,10^3}];

It is uniformly distributed:

pvals=DistributionFitTest[#,dist]&/@sample;

Histogram[pvals,Automatic,PDF]

p-values are uniformly distributed - 3

p-values can only falsify hypotheses, not confirm them, if the p-value is very small, it is “unlikely” (more on this later) that our hypothesis is true(e.g. that the data comes from the model).

To lower the false rejection rate, set a lower SignificanceLevel.

A recent paper in nature recommends 0.5%.

dist=NormalDistribution[];
sample=RandomVariate[dist,{10^4,10^3}];

conclusions=DistributionFitTest[#,dist,"ShortTestConclusion"]&/@sample;

By default we falsely reject 5%:

Counts[conclusions]/10.^4

SignificanceLevel is the false rejection rate:

conclusions=DistributionFitTest[#,dist,"ShortTestConclusion",SignificanceLevel->0.005]&/@sample;

Counts[conclusions]/10.^4

p-values are uniformly distributed - 4

Now we can appreciate this joke

Reveal joke to maximise comedic effect

Reproducibility crisis in Science

Let’s go back to this scary fact:

Reproducibility Project 2015: only 36% of studies in psychology meta-study were reproducible.

To understand this, let’s look at statistical power.

Statistical Power

By lowering SignificanceLevel we can lower the false rejection rate.

Lowering the false acceptance rate is trickier.

Suppose we have some data distributed as (Subscript[ℋ, 1]) and we want to test whether it is distributed as (Subscript[ℋ, 0]).

Here we will abuse notation and use  (Subscript[ℋ, 1]) and (Subscript[ℋ, 0]) to denote the distribution of the statistic:

ℋ0=NormalDistribution[0,1];
ℋ1=NormalDistribution[1,1];

Let’s assume a 2-sided test:

we reject the Null Hypothesis if our test statistic t is:

t ≤ q1 || q2 ≤ t

α=.05;
{q1,q2}=InverseCDF[ℋ0,{α/2,1-α/2}]

Plot[{PDF[ℋ0,x],PDF[ℋ1,x]},{x,-2,3},Filling->Axis,Epilog->{Black, HalfLine[{q1,0},{0,1}],HalfLine[{q2,0},{0,1}]},Ticks->{{{q1,"q1"}, {q2,"q2"}}, Automatic},PlotLegends->{"(Subscript[ℋ, 0])","(Subscript[ℋ, 1])"},PlotLabel->"PDF of (Subscript[ℋ, 0]) and (Subscript[ℋ, 1])"]

What’s the probability β that we (falsely) accept the Null hypothesis (false negative)?

β=Probability[q1≤x≤q2,x\[Distributed]ℋ1]

The power is defined as 1 - β:

power=1-β

The power (1-β) is the probability that we don’t mistakenly accept the Null hypothesis, which obviously depends on the data (which comes from  (Subscript[ℋ, 1]) in this case).

If the power of the test is low we cannot accurately distinguish between (Subscript[ℋ, 0]) and (Subscript[ℋ, 1]):

statisticSample=RandomVariate[ℋ1,n=10^6];

empiricalℋ1=EmpiricalDistribution[statisticSample]

The empirical power is close to the theoretical one:

1-Probability[q1≤x≤q2,x\[Distributed]empiricalℋ1]

Publication Bias

Publication bias is an important contributor to the irreproducibility of published results.

Publication of negative results is suppressed, which results in most of the published ones being wrong.

Two necessary conditions for publication bias:

sought result ((Subscript[ℋ, 1])) is rare ⟹ most results will be negative

we only publish positive results  ⟹ most published results will be wrong

Publication Bias - 2

Sought result ((Subscript[ℋ, 1])) is rare (10%), e.g. DNA studies:

n=10^7;
effectFreq=.1;

ℋ1sample=RandomVariate[ℋ1,n1=n*effectFreq];
ℋ0sample=RandomVariate[ℋ0,n0=n*(1-effectFreq)];

statisticSample=Join[ℋ1sample,ℋ0sample];

Histogram[{ℋ0sample,ℋ1sample},100,ChartLegends->{"(Subscript[ℋ, 0]) sample","(Subscript[ℋ, 1]) sample"},ChartStyle->{ColorData[97,2],ColorData[97,1]}]

pvals=CDF[ℋ0,statisticSample];

Usually we seek to reject the Null Hypothesis:

rejectFn=#<α/2||#>1-α/2&;

rejectionVec=rejectFn/@pvals;

Only keep rejections (positive results, i.e. the ones journals want to see):

Counts[rejectionVec]
allRejections=%[True]

Only 6% are positive (remember, 5% for (Subscript[ℋ, 0])):

allRejections/n//N

Only 27% of the results are correct:

Counts[Take[rejectionVec,n1]]
%[True]/allRejections//N

ℋ1rejections=Pick[ℋ1sample,Take[rejectionVec,n1]];
ℋ0rejections=Pick[ℋ0sample,Take[rejectionVec,-n0]];
Histogram[{ℋ0rejections,ℋ1rejections},100,ChartLegends->{"(Subscript[ℋ, 0]) rejections","(Subscript[ℋ, 1]) rejections"}, PlotLabel->Column[{"α = 5%, power = 17%"," Probability[rejected\[Conditioned](Subscript[ℋ, 1])] = 27%"}], ChartStyle->{ColorData[97,2],ColorData[97,1]}]

Publication Bias - Direct computation

We can also compute this using Bayes’ Theorem and MixtureDistribution:

stat𝒟=MixtureDistribution[{effectFreq,1-effectFreq},{ℋ1,ℋ0}];

Probability[ℋ1\[Conditioned]rejected]
==
((Probability[rejected0⋂ℋ1])/(Probability[rejected])) 
==
(Probability[rejected\[Conditioned]ℋ1]*Probability[ℋ1])/Probability[rejected]
==
(((trueRejections*effectFreq))/allRejections)

{q1,q2}=InverseCDF[ℋ0,{α/2,1-α/2}]

trueRejections=Probability[x≤q1||q2≤x,x\[Distributed]ℋ1]

falseRejections == α

falseRejections=Probability[x≤q1||q2≤x,x\[Distributed]ℋ0]

allRejections=Probability[x≤q1||q2≤x,x\[Distributed]stat𝒟]

We get 27%:

trueRejections * effectFreq/allRejections

In Bayesian terms, because the prior distribution for (Subscript[ℋ, 1]) is so low, the posterior is low too (we can only update our information so much).

Publication Bias - Improving Power

This is obvious and ideal but can be hard:

ℋ2=NormalDistribution[5,1];

Plot[{PDF[ℋ0,x],PDF[ℋ2,x]},{x,-2,7},Filling->Axis,Epilog->{Green, HalfLine[{q,0},{0,1}],HalfLine[{-q,0},{0,1}]},Ticks->{{{q,α}, {-q,1-α}}, Automatic},PlotLegends->{"(Subscript[ℋ, 0])","(Subscript[ℋ, 2])"},PlotLabel->"PDF of (Subscript[ℋ, 0]) and (Subscript[ℋ, 2])"]

the power when (Subscript[ℋ, 2]) is the alternative hypothesis is very high:

1-Probability[q1≤x≤q2,x\[Distributed]ℋ2]

let’s look at the distribution of sample statistics:

ℋ2sample=RandomVariate[ℋ2,n2=n*effectFreq];

statisticSample2=Join[ℋ2sample,ℋ0sample];

Histogram[{ℋ0sample,ℋ2sample},100,ChartLegends->{"(Subscript[ℋ, 0]) sample","(Subscript[ℋ, 2]) sample"},PlotRange->Full,ChartStyle->{ColorData[97,2],ColorData[97,3]}]

Histogram[{ℋ0sample,ℋ1sample,ℋ2sample},100,ChartLegends->{"(Subscript[ℋ, 0]) sample","(Subscript[ℋ, 1]) sample","(Subscript[ℋ, 2]) sample"},PlotRange->Full,ChartStyle->{ColorData[97,2],ColorData[97,1], ColorData[97,3]}]

pvals2=CDF[ℋ0,statisticSample2];

rejectionVec2=rejectFn/@pvals2;

Again, only keep positive results:

Counts[rejectionVec2]
allRejections2=%[True]

We have more rejections this time (previously 6%):

allRejections2/n//N

And 69% of them are correct:

Counts[Take[rejectionVec2,n2]]
%[True]/allRejections2//N

ℋ2rejections=Pick[ℋ2sample,Take[rejectionVec2,n2]];
Histogram[{ℋ0rejections,ℋ2rejections},100,ChartLegends->{"(Subscript[ℋ, 0]) rejections","(Subscript[ℋ, 2]) rejections"}, PlotLabel->Column[{"α = 5%, power = 99.8%","Probability[rejected\[Conditioned](Subscript[ℋ, 2])] = 69%"}],ChartStyle->{ColorData[97,2],ColorData[97,3]}]

Note that because the prior on (Subscript[ℋ, 2]) is low, we only get 69% correct results, even with high power

it is CRUCIAL to publish negative results.

Publication Bias - SignificanceLevel

Let’s lower the SignificanceLevel:

α2=.0001;
rejectFn3=#<α2/2||#>1-α2/2&;

rejectionVec3=rejectFn3/@pvals;

Select positive results:

Counts[rejectionVec3]
allRejections=%[True]

Now 67% of the results are correct, at the cost of  failing to reject some samples from (Subscript[ℋ, 1]).

Counts[Take[rejectionVec3,n1]]
%[True]/allRejections//N

ℋ1rejections=Pick[ℋ1sample,Take[rejectionVec3,n1]];
ℋ0rejections=Pick[ℋ0sample,Take[rejectionVec3,-n0]];
Histogram[{ℋ0rejections,ℋ1rejections},100,ChartLegends->{"(Subscript[ℋ, 0]) rejections","(Subscript[ℋ, 1]) rejections"},PlotLabel->Column[{"α = .01%, power = 17%"," Probability[rejected\[Conditioned](Subscript[ℋ, 1])] = 68%"}],PlotRange->Full,ChartStyle->{ColorData[97,2],ColorData[97,1]}]

LocationTest

Let’s look at an example where we are trying test whether two datasets have the same mean.

This is a typical case in pharmaceutical tests.

b=3;
n=10^2;
SeedRandom[2]

data1=RandomVariate[CauchyDistribution[0,b],n];
data2=RandomVariate[CauchyDistribution[1,b],n];

CauchyDistribution is a heavy-tailed distribution:

Variance/@{data1,data2}

LocationTest against (Subscript[ℋ, 0]), (Subscript[μ, 1])==(Subscript[μ, 2]).

we can use HypothesisTestData to generate an object we can query:

htd=LocationTest[{data1,data2},Automatic,"HypothesisTestData"]

htd["Properties"]

The power is very low due to the low number of points and high variance data, so some tests fail to reject the Null hypothesis:

htd[{"TestDataTable","ShortTestConclusion"},All];
MapAt[Column,%,2]

DistributionFitTest has different tests at its disposal:

htd=DistributionFitTest[data1,data2,"HypothesisTestData"]

htd[{"TestDataTable","ShortTestConclusion"},All];
MapAt[Column,%,2]

A bigger sample means we have more information to make a decision, which raises the statistical power:

n=10^3;

data1=RandomVariate[CauchyDistribution[0,b],n];
data2=RandomVariate[CauchyDistribution[1,b],n];

LocationTest[{data1,data2},Automatic,"HypothesisTestData"][{"TestDataTable","ShortTestConclusion"},All];
MapAt[Column,%,2]

Closing Remarks

‘The plural of Anecdote is not ‘Data’.’

Small datasets don’t have a lot of information, the inferences one can draw are limited.

‘Data analysis is simply a dialogue with the data.’— Stephen F. Gull, Cambridge 1994

There is no magic recipe in data analysis, you need good data and good assumptions for robust conclusions.

It is easy to delude oneself into thinking that one is making informed decisions as a result of being quantitative and sloppy.

Can only falsify hypothesis, not confirm it, some suggest that p-values only be used as exploratory qualitative numbers to decide whether one should carry on with an analysis.

p-values are easy to compute but hard to interpret, there is no free lunch.

Use 0.5% SignificanceLevel

“Redefine Statistical Significance”  2017