This vignette is based on Section 4.2 of Mulder, J., & Pfadt, J. M. (2026), “Going in the right direction: A tutorial to directional hypothesis testing using the BFpack module in JASP,” Advances in Methods and Practices in Psychological Science. The empirical case uses the verbal-memory data of Ichinose et al. (2019), distributed as the memory data set in BFpack.

1 Introduction and learning objectives

Many correlational questions are directional and joint: theory predicts not that a single association differs between two groups, but that a whole set of associations is uniformly stronger in one group than in another. The standard workflow answers such a question indirectly — one correlation-difference test per pair, followed by a multiplicity correction — and then asks the reader to reassemble the pieces. This tutorial shows a direct alternative: testing the single joint directional hypothesis with a Bayes factor using BFpack.

1.1 The running example

The concept of schizophrenia as a disorder of abnormal neural coordination — “dysconnectivity” — predicts that communication between brain regions is disrupted, which should have measurable behavioural consequences: the associations among cognitive-performance measures that are typically strong in healthy individuals should be weaker in patients. Ichinose et al. (2019) provide a data set to examine this. Twenty individuals diagnosed with schizophrenia (SZ) and twenty healthy controls (HC) each completed six assessments of verbal / working memory. This yields $\binom{6}{2} = 15$ pairwise correlations per group.

The substantive expectation is that every one of the 15 correlations is larger in the HC group than in the SZ group. We test it as a competition between three hypotheses:

\[ \begin{aligned} {H}_1 &:\ \rho_{\text{HC},pq} > \rho_{\text{SZ},pq} \quad\text{for all } 1\le p<q\le 6 \quad(\text{directional}),\\ {H}_2 &:\ \rho_{\text{HC},pq} = \rho_{\text{SZ},pq} \quad\text{for all } 1\le p<q\le 6 \quad(\text{equality}),\\ {H}_3 &:\ \text{complement (neither }{H}_1\text{ nor }{H}_2\text{).} \end{aligned} \]

1.2 What you will be able to do

By the end you will be able to:

estimate correlation matrices for two independent groups with cor_test();
read the exact correlation parameter names with get_estimates();
encode a 15-constraint directional expectation as a BFpack hypothesis string using >, =, and &;
compute Bayes factors and posterior probabilities with BF() and interpret them; and
see why the equivalent classical route — 15 Fisher $r$-to-$z$ tests with a multiplicity correction — gives only fragmented, lower-powered evidence.

2 Setup: packages, seed, and data

# install.packages("BFpack")   # once, if needed (correlation tests need v0.3+)
library("BFpack")

set.seed(1)   # cor_test() uses MCMC; fix the seed for reproducibility.

The memory data set ships with BFpack. It contains the six memory measures (Im, Del, Wmn, Cat, Fas, Rat) in the first six columns and a Group indicator (HC / SZ) in the last column. We split it into the two groups, dropping the grouping column:

memoryHC <- subset(memory, Group == "HC")[, -7]   # 20 healthy controls
memorySZ <- subset(memory, Group == "SZ")[, -7]   # 20 patients (schizophrenia)

colnames(memoryHC)   #> "Im" "Del" "Wmn" "Cat" "Fas" "Rat"

Measurement levels and the type of correlation. cor_test() chooses the appropriate correlation automatically from each variable’s measurement level: a product-moment correlation between two continuous variables, or a polyserial correlation when one variable is ordinal. If a variable such as Rat should be treated as continuous rather than ordinal, set its class accordingly (e.g. memoryHC$Rat <- as.numeric(memoryHC$Rat)) before calling cor_test().

3 The joint directional test in BFpack

3.1 Step 1 — Estimate the correlations per group

cor_test() performs Bayesian estimation of the correlation matrices. Passing two data frames fits the two groups jointly. The order matters: the first argument becomes group g1 and the second g2. We put HC first, so g1 = HC and g2 = SZ.

cor6 <- cor_test(memoryHC, memorySZ)

A joint uniform prior is placed over all admissible correlation matrices: every combination of correlations that yields a positive-definite matrix is equally likely a priori. This proper, “neutral” default is what makes the Bayes factor well defined for correlations without any manual prior tuning.

3.2 Step 2 — Read the exact parameter names

This step is essential. The names used in a hypothesis string are exactly the names printed by get_estimates(). They are built as [var]_with_[var]_in_[group], so the HC / SZ correlation between Del and Im appears as Del_with_Im_in_g1 and Del_with_Im_in_g2.

get_estimates(cor6)
#> Correlation names include, among others:
#>   Del_with_Im_in_g1,  Del_with_Im_in_g2,
#>   Del_with_Wmn_in_g1, Del_with_Wmn_in_g2,  ...  (30 in total: 15 per group)

3.3 Step 3 — Formulate the hypotheses

${H}_1$ is written as 15 one-sided constraints joined by &; ${H}_2$ as the matching 15 equality constraints. The two are separated by a semicolon. We do not write ${H}_3$ explicitly: BF() adds the complement automatically.

constraints_full <-
 "Del_with_Im_in_g1  > Del_with_Im_in_g2  &
  Del_with_Wmn_in_g1 > Del_with_Wmn_in_g2 &
  Del_with_Cat_in_g1 > Del_with_Cat_in_g2 &
  Del_with_Fas_in_g1 > Del_with_Fas_in_g2 &
  Del_with_Rat_in_g1 > Del_with_Rat_in_g2 &
  Im_with_Wmn_in_g1  > Im_with_Wmn_in_g2  &
  Im_with_Cat_in_g1  > Im_with_Cat_in_g2  &
  Im_with_Fas_in_g1  > Im_with_Fas_in_g2  &
  Im_with_Rat_in_g1  > Im_with_Rat_in_g2  &
  Wmn_with_Cat_in_g1 > Wmn_with_Cat_in_g2 &
  Wmn_with_Fas_in_g1 > Wmn_with_Fas_in_g2 &
  Wmn_with_Rat_in_g1 > Wmn_with_Rat_in_g2 &
  Cat_with_Fas_in_g1 > Cat_with_Fas_in_g2 &
  Cat_with_Rat_in_g1 > Cat_with_Rat_in_g2 &
  Fas_with_Rat_in_g1 > Fas_with_Rat_in_g2;
  Del_with_Im_in_g1  = Del_with_Im_in_g2  &
  Del_with_Wmn_in_g1 = Del_with_Wmn_in_g2 &
  Del_with_Cat_in_g1 = Del_with_Cat_in_g2 &
  Del_with_Fas_in_g1 = Del_with_Fas_in_g2 &
  Del_with_Rat_in_g1 = Del_with_Rat_in_g2 &
  Im_with_Wmn_in_g1  = Im_with_Wmn_in_g2  &
  Im_with_Cat_in_g1  = Im_with_Cat_in_g2  &
  Im_with_Fas_in_g1  = Im_with_Fas_in_g2  &
  Im_with_Rat_in_g1  = Im_with_Rat_in_g2  &
  Wmn_with_Cat_in_g1 = Wmn_with_Cat_in_g2 &
  Wmn_with_Fas_in_g1 = Wmn_with_Fas_in_g2 &
  Wmn_with_Rat_in_g1 = Wmn_with_Rat_in_g2 &
  Cat_with_Fas_in_g1 = Cat_with_Fas_in_g2 &
  Cat_with_Rat_in_g1 = Cat_with_Rat_in_g2 &
  Fas_with_Rat_in_g1 = Fas_with_Rat_in_g2"

3.4 Step 4 — Compute the Bayes factors

BF_full <- BF(cor6, hypothesis = constraints_full)
print(BF_full)
summary(BF_full)

By default the three hypotheses receive equal prior probabilities of $1/3$. To weight them differently, pass prior.hyp.conf (one weight per hypothesis, in order, including the complement) — this shifts the posterior probabilities but leaves the Bayes factors, which measure the evidence in the data, unchanged.

3.5 Results

Table 1. Bayes factors and posterior probabilities for the three hypotheses (equal prior probabilities).

Hypothesis	vs ${H}_1$	vs ${H}_2$	vs ${H}_3$	$P({H}_t\mid\text{Data})$
${H}_1$: all $\rho_{\text{HC},pq} > \rho_{\text{SZ},pq}$	1.000	$1.149\times10^{6}$	5779	1.000
${H}_2$: all $\rho_{\text{HC},pq} = \rho_{\text{SZ},pq}$	$8.70\times10^{-7}$	1.000	0.005	$\approx 0$
${H}_3$: complement	$1.73\times10^{-4}$	198.9	1.000	$\approx 0$

The joint one-sided hypothesis ${H}_1$ is favoured over the equality hypothesis by a Bayes factor of about $1.1\times10^{6}$ and over the complement by about $5.8\times10^{3}$. With equal priors this maps to a posterior probability of essentially 1 for ${H}_1$. A possible write-up:

“The joint one-sided hypothesis ${H}_1$ received overwhelming evidence against its equality-constrained alternative ${H}_2$ and the complement ${H}_3$, with Bayes factors of $1.1\times10^{6}$ and $5.8\times10^{3}$, respectively (equal prior probabilities; posterior probabilities 1.000, 0.000, 0.000). There is thus overwhelming evidence that the correlational structure among the memory measures is stronger in the healthy-control group than in the schizophrenia group across all pairs of variables.”

Note how a single number expresses support for the whole predicted pattern — exactly what the substantive theory claims, and exactly what a collection of separate tests cannot deliver.

4 Why not just test each correlation separately?

The classical route tests each of the 15 correlation differences on its own with Fisher’s $r$-to-$z$ transformation, then corrects for multiplicity. With 15 tests at $\alpha = .05$ the probability of at least one false positive under the global null is $1 - (1-.05)^{15} \approx .54$, which forces a correction — and the correction further drains power from an already small sample ($n = 20$ per group).

vars  <- colnames(memoryHC)
n1    <- nrow(memoryHC)
n2    <- nrow(memorySZ)
pairs <- combn(vars, 2, simplify = FALSE)

# Fisher r-to-z comparison of two independent correlations.
compare_corrs <- function(x1, x2, y1, y2, n1, n2) {
  r1 <- cor(x1, y1, use = "pairwise.complete.obs")
  r2 <- cor(x2, y2, use = "pairwise.complete.obs")
  z1 <- atanh(r1); z2 <- atanh(r2)
  se <- sqrt(1 / (n1 - 3) + 1 / (n2 - 3))
  zstat <- (z1 - z2) / se
  data.frame(r_HC = r1, r_SZ = r2, z = zstat,
             p_equal        = 2 * pnorm(abs(zstat), lower.tail = FALSE),
             p_HC_greater   = pnorm(zstat, lower.tail = FALSE))
}

results <- do.call(rbind, lapply(pairs, function(pr) {
  res <- compare_corrs(memoryHC[[pr[1]]], memorySZ[[pr[1]]],
                       memoryHC[[pr[2]]], memorySZ[[pr[2]]], n1, n2)
  res$pair <- paste0(pr[1], "_with_", pr[2]); res
}))

# Holm correction across the 15 one-sided tests.
results$p_HC_greater_holm <- p.adjust(results$p_HC_greater, method = "holm")
results <- results[order(results$pair), ]
print(cbind(results$pair, round(results[, c("r_HC","r_SZ","z",
      "p_equal","p_HC_greater","p_HC_greater_holm")], 3)), row.names = FALSE)

Table 2. Classical post-hoc comparisons (Fisher $r$-to-$z$, Holm-adjusted).

Pair	$r_{\text{HC}}$	$r_{\text{SZ}}$	$z$	$p_{\text{equal}}$	$p_{\text{HC}>\text{SZ}}$	Adj. $p_{\text{HC}>\text{SZ}}$
Cat–Fas	0.734	0.219	2.08	0.037	0.019	0.119
Cat–Rat	0.769	−0.251	3.71	<0.001	<0.001	0.002
Del–Cat	0.395	0.164	0.74	0.461	0.231	0.461
Del–Fas	0.322	0.269	0.17	0.866	0.433	0.461
Del–Rat	0.466	0.088	1.21	0.225	0.112	0.337
Del–Wmn	0.499	−0.223	2.26	0.024	0.012	0.095
Fas–Rat	0.673	−0.144	2.80	0.005	0.003	0.033
Im–Cat	0.563	−0.280	2.70	0.007	0.003	0.042
Im–Del	0.833	0.346	2.44	0.015	0.007	0.069
Im–Fas	0.388	−0.173	1.70	0.089	0.044	0.222
Im–Rat	0.544	0.078	1.55	0.121	0.061	0.243
Im–Wmn	0.654	−0.066	2.47	0.013	0.007	0.069
Wmn–Cat	0.773	−0.053	3.15	0.002	0.001	0.011
Wmn–Fas	0.700	0.010	2.50	0.012	0.006	0.069
Wmn–Rat	0.611	−0.016	2.12	0.034	0.017	0.119

Every correlation is numerically larger in HC than in SZ — perfectly consistent with the directional expectation — yet after Holm correction only a handful of the 15 comparisons remain significant. This is the familiar bind:

Multiplicity and power. Guarding 15 tests against false positives requires a correction that removes exactly the power a 20-per-group study can least afford, so genuinely uniform effects fail to reach significance one by one.
No single statement. The output is a scatter of 15 $p$-values, not an answer to the actual question (“is the whole pattern stronger in HC?”). Reassembling them into support for the theory is subjective and invites selective emphasis on the significant rows.

The joint Bayes-factor test in the previous section sidesteps both problems: it aggregates the mild-to-moderate evidence in each pair into one very large evidential statement, with no multiplicity correction, because it is a single test of a single hypothesis.

Hypothesis	vs \({H}_1\)	vs \({H}_2\)	vs \({H}_3\)	\(P({H}_t\mid\text{Data})\)
\({H}_1\): all \(\rho_{\text{HC},pq} > \rho_{\text{SZ},pq}\)	1.000	\(1.149\times10^{6}\)	5779	1.000
\({H}_2\): all \(\rho_{\text{HC},pq} = \rho_{\text{SZ},pq}\)	\(8.70\times10^{-7}\)	1.000	0.005	\(\approx 0\)
\({H}_3\): complement	\(1.73\times10^{-4}\)	198.9	1.000	\(\approx 0\)

Jointly testing Directional Expectations on Correlations with BFpack

Joris Mulder

2026-07-10