Rでまとめて無相関検定（相関分析）を行って、その結果を一覧したいと思った。
相関係数だけを一括して見たいならcor()という関数もあるが、無相関検定や95％信頼区間なども一緒に確認したいと思った。
cor(iris[1:4])
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length 1.0000000 -0.1175698 0.8717538 0.8179411
## Sepal.Width -0.1175698 1.0000000 -0.4284401 -0.3661259
## Petal.Length 0.8717538 -0.4284401 1.0000000 0.9628654
## Petal.Width 0.8179411 -0.3661259 0.9628654 1.0000000
データフレーム中の複数の列を組み合わせてまとめて無相関検定を行う関数を作る（答えはlist形式で返す）。
multi.cor.test <- function(df, vars, ...){ if(missing(vars)) vars <- colnames(df) vars_length <- length(vars) vars_list_x <- vector() vars_list_y <- vector() n <- vars_length - 1 for(i in vars){ if(n == 0) break vars_list_x <- c(vars_list_x, rep(i, n)) n <- n - 1 } for(i in 1:vars_length){ vars_list_y <- c(vars_list_y, vars[-c(1:i)]) } mapply(function(x, y) cor.test(as.formula(paste("~", x, "+", y)), data = df, ...), vars_list_x, vars_list_y, SIMPLIFY = FALSE) }
この関数だと、データフレームの指定された列同士（vars）のすべての組み合わせて無相関検定を行い結果をリストに保存する処理となる。
varsが指定されていない場合、データフレームのすべての列を用いて計算する。
このとき、変数にaとbとあったとき、aとaの組み合わせは相関係数1になることが明らかなので、その組み合わせは計算しない。また、aとbの組み合わせとbとaの組み合わせは相関分析においては同義なので、aとbのみの組み合わせだけを計算するようにしている。
たとえばirisデータのSpeciesを除く全ての組み合わせで無相関検定を行う場合は次のようにする。
# Species以外の行名を取り出しておく。 vars <- colnames(iris[1:4]) multi.cor.test(df = iris, vars = vars)
## $Sepal.Length
##
## Pearson's product-moment correlation
##
## data: Sepal.Length and Sepal.Width
## t = -1.4403, df = 148, p-value = 0.1519
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.27269325 0.04351158
## sample estimates:
## cor
## -0.1175698
##
##
## $Sepal.Length
##
## Pearson's product-moment correlation
##
## data: Sepal.Length and Petal.Length
## t = 21.646, df = 148, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.8270363 0.9055080
## sample estimates:
## cor
## 0.8717538
##
##
## $Sepal.Length
##
## Pearson's product-moment correlation
##
## data: Sepal.Length and Petal.Width
## t = 17.296, df = 148, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.7568971 0.8648361
## sample estimates:
## cor
## 0.8179411
##
##
## $Sepal.Width
##
## Pearson's product-moment correlation
##
## data: Sepal.Width and Petal.Length
## t = -5.7684, df = 148, p-value = 4.513e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5508771 -0.2879499
## sample estimates:
## cor
## -0.4284401
##
##
## $Sepal.Width
##
## Pearson's product-moment correlation
##
## data: Sepal.Width and Petal.Width
## t = -4.7865, df = 148, p-value = 4.073e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.4972130 -0.2186966
## sample estimates:
## cor
## -0.3661259
##
##
## $Petal.Length
##
## Pearson's product-moment correlation
##
## data: Petal.Length and Petal.Width
## t = 43.387, df = 148, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9490525 0.9729853
## sample estimates:
## cor
## 0.9628654
## $Sepal.Length
##
## Pearson's product-moment correlation
##
## data: Sepal.Length and Sepal.Width
## t = -1.4403, df = 148, p-value = 0.1519
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.27269325 0.04351158
## sample estimates:
## cor
## -0.1175698
##
##
## $Sepal.Length
##
## Pearson's product-moment correlation
##
## data: Sepal.Length and Petal.Length
## t = 21.646, df = 148, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.8270363 0.9055080
## sample estimates:
## cor
## 0.8717538
##
##
## $Sepal.Length
##
## Pearson's product-moment correlation
##
## data: Sepal.Length and Petal.Width
## t = 17.296, df = 148, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.7568971 0.8648361
## sample estimates:
## cor
## 0.8179411
##
##
## $Sepal.Width
##
## Pearson's product-moment correlation
##
## data: Sepal.Width and Petal.Length
## t = -5.7684, df = 148, p-value = 4.513e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5508771 -0.2879499
## sample estimates:
## cor
## -0.4284401
##
##
## $Sepal.Width
##
## Pearson's product-moment correlation
##
## data: Sepal.Width and Petal.Width
## t = -4.7865, df = 148, p-value = 4.073e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.4972130 -0.2186966
## sample estimates:
## cor
## -0.3661259
##
##
## $Petal.Length
##
## Pearson's product-moment correlation
##
## data: Petal.Length and Petal.Width
## t = 43.387, df = 148, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9490525 0.9729853
## sample estimates:
## cor
## 0.9628654
このままだと見た目の量が多くやや見にくいので、統計量や相関係数、p値や95％信頼区間をデータフレーム形式で表示できると嬉しい。
そこで、{broom}パッケージのtidy()関数をつかって、計算された値などをデータフレーム形式（tibble形式）で取り出す。
答えがlist形式なので{purrr}パッケージのmap_df()を使う。{broom}のtidy()はhtestクラスだと何と何の統計結果なのか（print()した時のdataのところ、list上のdata.nameのところ）を取り出さないので、それも取り出すようにすると次にようになる。
multi.cor.test(df = iris[1:4]) |> purrr::map_df(\(x) c(data = x$data.name, broom::tidy(x)))
## # A tibble: 6 × 9
## data estimate statistic p.value parameter conf.low conf.high method
## <chr> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <chr>
## 1 Sepal.Length … -0.118 -1.44 1.52e- 1 148 -0.273 0.0435 Pears…
## 2 Sepal.Length … 0.872 21.6 1.04e-47 148 0.827 0.906 Pears…
## 3 Sepal.Length … 0.818 17.3 2.33e-37 148 0.757 0.865 Pears…
## 4 Sepal.Width a… -0.428 -5.77 4.51e- 8 148 -0.551 -0.288 Pears…
## 5 Sepal.Width a… -0.366 -4.79 4.07e- 6 148 -0.497 -0.219 Pears…
## 6 Petal.Length … 0.963 43.4 4.68e-86 148 0.949 0.973 Pears…
## # ℹ 1 more variable: alternative <chr>
2つの変数はandで文章接続されているので次のようにすると変数名で取り出せる。
multi.cor.test(df = iris, vars = vars) |> purrr::map_df(\(x){ xy_name <- strsplit(x$data.name, " and ") c(x = xy_name[[1]][1], y = xy_name[[1]][2], broom::tidy(x))})
## # A tibble: 6 × 10
## x y estimate statistic p.value parameter conf.low conf.high method
## <chr> <chr> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <chr>
## 1 Sepal.L… Sepa… -0.118 -1.44 1.52e- 1 148 -0.273 0.0435 Pears…
## 2 Sepal.L… Peta… 0.872 21.6 1.04e-47 148 0.827 0.906 Pears…
## 3 Sepal.L… Peta… 0.818 17.3 2.33e-37 148 0.757 0.865 Pears…
## 4 Sepal.W… Peta… -0.428 -5.77 4.51e- 8 148 -0.551 -0.288 Pears…
## 5 Sepal.W… Peta… -0.366 -4.79 4.07e- 6 148 -0.497 -0.219 Pears…
## 6 Petal.L… Peta… 0.963 43.4 4.68e-86 148 0.949 0.973 Pears…
## # ℹ 1 more variable: alternative <chr>
