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