[Feature] cor_diff() : test for differences between correlations - #338
[Feature] cor_diff() : test for differences between correlations#338DominiqueMakowski wants to merge 8 commits into
Conversation
|
@bwiernik @mattansb do you know of any way to compute a Bayes factor for that? |
|
|
||
| # If pairs are passed | ||
| if(length(x) == 2 & length(y) == 2) { | ||
| if(length(x) == 2 && length(y) == 2) { |
|
I'm unfamiliar with a Bayes factor for such a comparison. |
|
|
|
Didn't look into code yet. The straightforward approach to these tests is to do a Fisher z transform and then do a z test for the difference (potentially accounting for covariance). So BF methods analogous to a z test would work |
|
Test: library(ggplot2)
library(correlation)
rez <- data.frame()
for(r in seq(0, 0.6, length.out=200)) {
for(n in seq(20, 100, by=30)) {
data <- bayestestR::simulate_correlation(
n = n,
r = matrix(c(
1.0, r, 0.0,
r, 1.0, 0.0,
0.0, 0.0, 1.0), nrow = 3))
rez <- cor_diff(data, x="V1", y="V2", x2="V1", y2="V3", method = "parametric") |>
datawizard::data_rename("t", "z") |>
rbind(cor_diff(data, x="V1", y="V2", x2="V1", y2="V3", method = "bootstrapping", iterations=2000)) |>
datawizard::data_modify(n = n, r = r) |>
rbind(rez)
}
}
rez |>
datawizard::reshape_longer(select=-c("r", "n", "Method")) |>
datawizard::data_modify(type = ifelse(name %in% c("p"), "p-value", "z-value")) |>
ggplot(aes(x=r, y=value)) +
geom_smooth(aes(group=interaction(Method, n), color=n), method="loess", se=FALSE, formula = 'y ~ x') +
geom_point(aes(color=n), alpha=0.1) +
facet_grid(type~Method, scales="free_y", switch="y") +
theme_minimal() |
|
bump - what's the status of this PR? |
There was a problem hiding this comment.
Pull request overview
This PR introduces a new user-facing API, cor_diff(), to test whether two correlations differ, with both a parametric approach (via psych::r.test()) and a bootstrap-based alternative. It also adds an S3 print method plus documentation and a basic test.
Changes:
- Added
cor_diff()with parametric and bootstrapped implementations, plus an S3print.cor_diff()method. - Exported the new function and registered its S3 print method.
- Added generated Rd docs and a new test file.
Reviewed changes
Copilot reviewed 3 out of 4 changed files in this pull request and generated 6 comments.
| File | Description |
|---|---|
| tests/testthat/test-cor_diff.R | Adds a unit test for cor_diff() (currently only parametric path). |
| R/cor_diff.R | Implements cor_diff(), internal helpers for parametric/bootstrapped testing, and an S3 print method. |
| NAMESPACE | Exports cor_diff and registers print.cor_diff. |
| man/cor_diff.Rd | Adds generated documentation for the new API. |
Files not reviewed (1)
- man/cor_diff.Rd: Generated file
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| test_that("cor_diff", { | ||
| expect_equal( | ||
| cor_diff(iris, "Sepal.Length", "Sepal.Width", "Sepal.Length", "Petal.Width")$t, | ||
| -10, | ||
| tolerance = 0.001 | ||
| ) | ||
| }) |
| # Compute | ||
| if (method %in% c("bootstrapping")) { | ||
| out <- .cor_diff_bootstrapping(data, x, y, x2, y2, ...) | ||
| } else { | ||
| out <- .cor_diff_parametric(data, x, y, x2, y2, ...) | ||
| } |
| args <- list(n = nrow(data), r12 = stats::cor(data[[x]], data[[y]])) | ||
| if (x == x2 && y != y2) { | ||
| args$r13 <- stats::cor(data[[x]], data[[y2]]) | ||
| args$r23 <- stats::cor(data[[y]], data[[y2]]) | ||
| } else if (y == y2 && x != x2) { | ||
| args$r13 <- stats::cor(data[[y]], data[[x2]]) | ||
| args$r23 <- stats::cor(data[[x]], data[[x2]]) | ||
| } else { | ||
| args$r34 <- stats::cor(data[[x2]], data[[y2]]) | ||
| } |
| .cor_diff_bootstrapping <- function(data, x, y, x2, y2, iterations = 1000, robust = FALSE, ...) { | ||
| diff <- rep(NA, iterations) # Initialize vector | ||
|
|
||
| # Bootstrap | ||
| for (i in 1:iterations) { | ||
| # Take random sample of data | ||
| dat <- data[sample(nrow(data), nrow(data), replace = TRUE), ] | ||
| # Compute diff | ||
| diff[i] <- stats::cor(dat[[x]], dat[[y]]) - stats::cor(dat[[x2]], dat[[y2]]) | ||
| } |
| # Summarize | ||
| if (robust == FALSE) { | ||
| out <- data.frame( | ||
| Method = "bootstrapping", | ||
| z = mean(diff) / stats::sd(diff), | ||
| p = bayestestR::pd_to_p(as.numeric(bayestestR::p_direction(diff))) | ||
| ) | ||
| } else { | ||
| out <- data.frame( | ||
| Method = "bootstrapping_robust", | ||
| z = stats::median(diff) / stats::mad(diff), | ||
| p = bayestestR::pd_to_p(as.numeric(bayestestR::p_direction(diff))) | ||
| ) | ||
| } |
| #' @param x,y,x2,y2 The variable names in `data` to be used. `x` and `y` can also | ||
| #' be pairs of variables, in which case the second variable is used as `x2` and `y2`. | ||
| #' @param method Can be `"parametric"` or `"bootstrapping"`. If `"parametric"`, | ||
| #' the [psych::r.test()] function is used. If `"bootstrapping"`, a bootstrapping | ||
| #' procedure is used. | ||
| #' @param ... Other arguments to be passed, for instance `iterations` (default: 1000) | ||
| #' if method is bootstrapping. |


Simple function with parametric vs. bootstrapped approaches. I ran some simulations and it seems to nicely work.
Created on 2025-01-29 with reprex v2.1.1