Bayesian hypothesis testing with the Full Bayesian Evidence Test

The Bayesian evidence interval

First, the Bayesian evidence interval is required for the FBET. Let $s(\theta):=\frac{p(\theta|y)}{r(\theta)}$ be the surprise function of the Pereira-Stern FBST with a given reference function r(θ) and tangential set $\overline{T}(\nu):=\{\theta \in \Theta|s(\theta)> \nu \}$ to the null hypothesis H₀ : θ = θ₀. The expanded tangential set is defined as $\tilde{\overline{T}}(\nu)=\{\theta \in \Theta|s(\theta)\geq \nu \}$, which in contrast to $\overline{T}(\nu)$ also includes the parameter values θ for which s(θ) = ν.

The Bayesian evidence interval EI_r(ν) with reference function r(θ) to level ν includes all parameter values θ for which $\theta \in \tilde{\overline{T}}(\nu)$.

It can be shown that the Bayesian evidence interval includes standard HPD intervals, support intervals and credible intervals as special cases.

To unify Bayesian interval estimation and hypothesis testing based on the EI_r(ν) via an interval hypothesis (or the region of practical equivalence (ROPE)), the Bayesian evidence value is introduced, which is associated with a Bayesian evidence interval: For a given evidence interval EI_r(ν) with reference function r(θ) to level ν, the Bayesian evidence value Ev_EIr(ν)^R(H₀) to a ROPE R for the null hypothesis H₀ is given as Ev_EIr(ν)^R(H₀) := ∫_{EIr(ν) ∩ R}p(θ|y)dθ The corresponding Bayesian evidence value Ev_EIr(ν)^R(H₁) to a ROPE R for the alternative hypothesis H₁ is given as Ev_EIr(ν)^R(H₁) := ∫_{EIr(ν) ∩ R^c}p(θ|y)dθ where R^c is the set complement of the ROPE R. The Bayesian evidence value Ev_EIr(ν)^R(H₀) is the integral of the posterior density p(θ|y) over the intersection of the evidence interval EI_r(ν) with the ROPE (or interval hypothesis) R. This means that the evidence value Ev_EIr(ν)^R(H₀) is the integral over the evidence interval EI_r(ν) restricted to the ROPE (or interval hypothesis) R around the null value of H₀ : θ = θ₀. The Bayesian evidence value can also be called generalised e-value, as it can be shown to recover the e-value of the FBST as a special case.

The Bayesian evidence value depends on three quantities:

• The ROPE R := [a, b], a < b, which includes the null value θ₀ of H₀ : θ = θ₀, that is θ₀ ∈ R. Thus, the ROPE R can be interpreted as an interval hypothesis H₀ := [a, b] around the null value θ₀.
• the reference function r(θ) which is used for construction of the Bayesian evidence interval EI_r(ν)
• and the evidence-threshold ν used for deciding which values are included in the Bayesian evidence interval EI_r(ν)

The evidence value of course also depends on the prior density p(θ), because the Bayesian evidence interval depends on the posterior density, which itself depends on the prior density p(θ). Of course, the Bayesian evidence value also depends on the data y observed.

From a hypothesis testing perspective a ROPE R is interpreted as an interval hypothesis which includes the precise null value of the hypothesis H₀. Thus, from this point of view one could also omit the superscript R in Ev_EIr(ν)^R(H₀) and define the Bayesian evidence value directly for interval hypotheses H₀ ∈ Θ as Ev_EIr(ν)(H₀) := ∫_{EIr(ν) ∩ H0}p(θ|y)dθ The general approach to consider a (small) interval hypothesis instead of a point-null hypothesis goes back to Hodges & Lehmann (1954).

Example 1: Paired-sample t-test with flat reference function

We use Student’s sleep data as illustration, which is available in the data frame sleep.

sleep
#>    extra group ID
#> 1    0.7     1  1
#> 2   -1.6     1  2
#> 3   -0.2     1  3
#> 4   -1.2     1  4
#> 5   -0.1     1  5
#> 6    3.4     1  6
#> 7    3.7     1  7
#> 8    0.8     1  8
#> 9    0.0     1  9
#> 10   2.0     1 10
#> 11   1.9     2  1
#> 12   0.8     2  2
#> 13   1.1     2  3
#> 14   0.1     2  4
#> 15  -0.1     2  5
#> 16   4.4     2  6
#> 17   5.5     2  7
#> 18   1.6     2  8
#> 19   4.6     2  9
#> 20   3.4     2 10

We select both groups and perform a traditional Bayes factor paired-sample t-test based on the model of Rouder et al. (2009). There, the null hypothesis states that the standardized effect size δ := (μ − μ₀)/σ is assumed to be zero under H₀, and Cauchy-distributed under the alternative H₁, δ ∼ C(0, γ) for γ > 0.

require(BayesFactor)
#> Loading required package: BayesFactor
#> Loading required package: coda
#> Loading required package: Matrix
#> ************
#> Welcome to BayesFactor 0.9.12-4.7. If you have questions, please contact Richard Morey ([email protected]).
#> 
#> Type BFManual() to open the manual.
#> ************
#> 
#> Attaching package: 'BayesFactor'
#> The following object is masked from 'package:ks':
#> 
#>     compare
grp1 = sleep[1:10,]$extra
grp2 = sleep[11:20,]$extra
diff = grp1 - grp2
ttestBF(x = grp1, y = grp2, rscale = "medium", paired = TRUE)
#> Bayes factor analysis
#> --------------
#> [1] Alt., r=0.707 : 17.25888 ±0%
#> 
#> Against denominator:
#>   Null, mu = 0 
#> ---
#> Bayes factor type: BFoneSample, JZS

Based on a traditional Bayes factor test, there is strong evidence for the alternative H₁ : δ ≠ 0: BF₁₀ = 17.25 according to Jeffreys (1961). However, note that the precise Bayes factor tests the point-null hypothesis H₀ : δ = 0 of exactly no effect, that is, no difference between both population means. Without observing any data y, we can be reasonable sure that the effect size δ will almost surely not be exactly zero, so an interval hypothesis which separates scientifically meaningful from scientifically irrelevant parameter values constitutes a more adequate test in such a situation (interval Bayes factors would be more adequate when sticking to Bayes factors).

Therefore, we now obtain posterior MCMC samples and conduct the FBET based on a flat reference function r(θ) := 1. We use the evidence threshold ν = 0, which implies that the surprise function s(θ) = p(θ|y)/r(θ) reduces to the posterior density p(θ|y) and the Bayesian evidence interval EI_r(ν) becomes EI_r(0), which includes all parameter values θ ∈ Θ with positive posterior density p(θ|y) ≥ 0.

set.seed(42)
posteriorDraws = ttestBF(x = grp1, y = grp2, rscale = "medium", paired = TRUE,
                         posterior = TRUE, iterations = 100000)[,3]

set.seed(42)
result = fbet(posteriorDensityDraws = posteriorDraws, interval = c(-0.2,0.2), nu=0)
summary(result)
#> Full Bayesian Evidence Test:
#> Reference function: Flat 
#> Testing Hypothesis H_0:=[ -0.2 , 0.2 ]
#> Bayesian e-value in favour of H_0: 0.01416656 
#> Bayesian e-value against H_0: 0.9858334

The Bayesian evidence value Ev_EIr(ν)^R(H₀) in favour of H₀ : θ ∈ [−0.1, 0.1] is equal to 0.014, while the Bayesian evidence value Ev_EIr(ν)^R(H₁) in favour of H₁ is equal to 0.986. Thus, the evidence for the alternative H₁ : θ ∉ [−0.1, 0.1] is much stronger than the evidence in favour of H₀. In fact, 98.6% of the posterior probability mass are located outside the interval [−0.1, 0.1], which is expressed by the Bayesian evidence value in favour of H₁. The following plot illustrates this result:

plot(result, type = "posterior")

The blue area amounts to the 0.5% of posterior probability mass which are located inside [−0.1, 0.1] and corroborate the null hypothesis H₀ based on the chosen ROPE R := [−0.1, 0.1] around θ₀ := 0. However, the evidence threshold ν = 0 requires no form of corroboration for the parameter values θ after observing the data y, which is shown in the above plot: The evidence threshold is not visible, as it is a horizontal dashed line at the ordinate y = 0.

If we choose a different hypothesis H₀ : δ ∈ (−0.25, 0.25), we obtain:

set.seed(42)
result2 = fbet(posteriorDensityDraws = posteriorDraws, interval = c(-0.25,0.25), nu=0)
summary(result2)
#> Full Bayesian Evidence Test:
#> Reference function: Flat 
#> Testing Hypothesis H_0:=[ -0.25 , 0.25 ]
#> Bayesian e-value in favour of H_0: 0.02020536 
#> Bayesian e-value against H_0: 0.9797946

plot(result2, type = "posterior")

Example 2: Paired-sample t-test with user-specified reference function

In the above, a flat reference function r(θ) := 1 was used in the surprise function s(θ) and the evidence threshold ν was set to zero. However, the prior distribution on the parameter θ was a medium Cauchy prior $C(0,\sqrt{2}/2)$. As the prior beliefs are reflected via this prior, it is reasonable to choose the Cauchy prior as the reference function instead (then, the surprise function is equal to what Good (1968) motivated as statistical information, providing an independent justification of the expanded tangential set; the latter includes parameter values which are informative with regard to the hypothesis under consideration). Below, it is shown how to do this. Also, the evidence threshold ν is set to ν = 1, so that only parameter values θ are included in the Bayesian evidence interval EI_r(1) which attain a larger posterior density value than the corresponding prior density value. This time, we test H₀ : δ ∈ (−0.5, 0.5) against H₁ : δ ∉ (−0.5, 0.5):

set.seed(42)
result3 = fbet(posteriorDensityDraws = posteriorDraws, interval = c(-0.5,0.5), 
               nu=1, FUN = dcauchy, par=list(location = 0, scale = sqrt(2)/2))
summary(result3)
#> Full Bayesian Evidence Test:
#> Reference function: User-defined 
#> Testing Hypothesis H_0:=[ -0.5 , 0.5 ]
#> Bayesian e-value in favour of H_0: 0.01750199 
#> Bayesian e-value against H_0: 0.9108208

Now, based on the different reference function and evidence threshold, we see that the resulting Bayesian evidence values change accordingly. Below, the situation is illustrated:

plot(result3, type="surprise")

The horizontal dashed line shows the evidence-threshold ν = 1, and only parameter values θ which fulfill p(θ|y)/p(θ) ≥ 1 are included in the Bayesian evidence interval EI_r(1). The resulting Bayesian evidence interval EI_r(1) is shown as the vertical blue lines in the plot. The interval hypothesis H₀ := (−0.5, 0.5) is shown as the dashed vertical blue lines, and the resulting Bayesian evidence value Ev_EIr(ν)^R(H₀) in favour of H₀ is visualized as the blue area, which amounts to 1.75% of the posterior probability mass. The Bayesian evidence value Ev_EIr(ν)^R(H₁) in favour of H₁ is given as Ev_EIr(ν)^R(H₁) = 0.9117, which shows that 91.17% of the posterior probability mass confirm the alternative. It should be clear by now that omitting the ROPE R and instead switching towards an interval hypothesis H₀ from the start only simplifies the notation. Henceforth, we therefore only speak of H₀.

Now, if we choose a different interval hypothesis H₀ := (−0.2, 0.2), we obtain:

result4 = fbet(posteriorDensityDraws = posteriorDraws, interval = c(-0.2,0.2), 
               nu=1, FUN = dcauchy, par=list(location = 0, scale = sqrt(2)/2))
summary(result4)
#> Full Bayesian Evidence Test:
#> Reference function: User-defined 
#> Testing Hypothesis H_0:=[ -0.2 , 0.2 ]
#> Bayesian e-value in favour of H_0: 0 
#> Bayesian e-value against H_0: 0.9283223

Now, based on the different reference function and evidence threshold, we see that the resulting Bayesian evidence values change accordingly: The evidence for the null hypothesis reduces, while the evidence for the alternative has increased. Below, the situation is illustrated:

plot(result4, type="surprise")

Next to plotting the surprise function we can also plot the resulting posterior density and visualize the evidence values. Therefore, we just need to omit the type argument.

plot(result4)

Also, we could require that values need to be corroborated in a stronger form by choosing for example ν = 2. Then, every value inside the Bayesian evidence interval has at least twice the value of the corresponding reference function r(θ). In the above example, we obtain:

set.seed(42)
result5 = fbet(posteriorDensityDraws = posteriorDraws, interval = c(-0.2,0.2), nu=2, 
               FUN = dcauchy, par=list(location = 0, scale = sqrt(2)/2))
summary(result5)
#> Full Bayesian Evidence Test:
#> Reference function: User-defined 
#> Testing Hypothesis H_0:=[ -0.2 , 0.2 ]
#> Bayesian e-value in favour of H_0: 0 
#> Bayesian e-value against H_0: 0.8546793

Evidently, as the evidence for H₀ was already zero, the only thing which changes in this different setting now is that the evidence in favour of H₁ decreases to about 85.57%. Thus, compared to the previous example we are now less than 90% sure that the alternative hypothesis holds. To arrive at this conclusion, we must, however, accept that the reference function is sensible, and that we regard parameter values as evidential in favour of the null or alternative only if the ratio of posterior to prior density is at least ν = 2. The following plot visualizes the posterior in this setting:

plot(result5)

In contrast, note that the surprise function plot differs. Importantly, the evidence threshold ν is now equal to 2:

plot(result5, type = "surprise")