---
title: "Jointly testing Directional Expectations on Correlations with BFpack"
author: "Joris Mulder"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
    number_sections: true
vignette: >
  %\VignetteIndexEntry{Testing Directional Expectations on Correlations with BFpack}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
header-includes:
  - \usepackage{mathrsfs}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE,
  eval = FALSE,        # Code is shown but not run when knitting (keeps the
  message = FALSE,     # vignette fast for CRAN). The expected results are
  warning = FALSE,     # reported in the text and tables. Set eval = TRUE, or
  collapse = TRUE,     # paste the chunks into R, to reproduce them live.
  comment = "#>"
)
```

> **How to use this document.** Every code block below can be pasted into an R
> session and run in order. To keep the vignette light (and fast to check on
> CRAN) the chunks are *not* executed while knitting (`eval = FALSE` in the
> setup chunk); instead the expected results are reported and discussed in the
> text and tables. To run the analysis live, change `eval = FALSE` to
> `eval = TRUE`, or simply copy each block into R.

*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`.*

---

# 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`.

## 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}
$$

## 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.

---

# Setup: packages, seed, and data

```{r packages}
# 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:

```{r data}
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()`.

---

# The joint directional test in BFpack

## 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`.

```{r cortest}
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.

## 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`.

```{r estimates}
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)
```

## 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.

```{r hypotheses}
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"
```

## Step 4 — Compute the Bayes factors

```{r BFtest}
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.

## 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.

---

# 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).

```{r posthoc}
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.

---

# References

Ichinose, M. C., Han, G., Polyn, S., Park, S., & Tomarken, A. J. (2019). Verbal
memory performance discordance in schizophrenia: A reflection of cognitive
dysconnectivity? *(Data distributed as the `memory` set in `BFpack`.)*

Mulder, J. (2016). Bayes factors for testing order-constrained hypotheses on
correlations. *Journal of Mathematical Psychology.*

Mulder, J., & Gelissen, J. P. (2019). Bayes factor testing of equality and
order constraints on measures of association in social research.
*Journal of Applied Statistics.* <https://doi.org/10.1080/02664763.2021.1992360>

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.*

Mulder, J., et al. (2021). BFpack: Flexible Bayes factor testing of scientific
expectations in R. *Journal of Statistical Software, 100*(18), 1–63.

---

```{r session-info, eval=FALSE, echo=FALSE}
# sessionInfo()   # uncomment when running live to record the environment
```
