相関係数行列の計算と検定 Last modified: Dec 27, 2013
目的
相関係数行列を求め，各相関係数について無相関検定を行う。
使用法
mycor(i, df, use=c("all.obs", "complete.obs", "pairwise.complete.obs"),
method=c("pearson", "kendall", "spearman"), latex=TRUE, caption=NULL, label=NULL,
digits=3, output="")
引数
i データフレーム上で，相関係数を計算する変数が入っている，列の番号または変数名（ベクトルとして複数指定可）
df 読み込んだデータフレームの名前
use 欠損値の処理法の指定
"all.obs", "complete.obs" i で指定した変数の内 1 つでも欠損値がある
ケースは除外する
"pairwise.complete.obs" 1 個の相関係数を算出するとき，
二変数の何れかまたは両方が欠損値のとき除外する
method "pearson", "kendall", "spearman" それぞれ，ピアソンの積率相関係数，
ケンドールの順位相関係数，スピアマンの順位相関係数を計算する
latex LaTeX ソースを出力する。Word 用なら，latex=FALSE にする
caption LaTeX 出力のときの表題
label LaTeX 出力のときのラベル
digits 相関係数を書き出すときに小数点以下何桁までにするかを指定
デフォルトでは 3
output 出力コネクション
デフルトはコンソールに出力
ソース
インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/mycor.R", encoding="euc-jp")
#####
#
# 相関係数行列を作成し，無相関検定の結果と共に表示する
#
#####
mycor <- function(i, # 相関係数行列を求める変数が入っているデータフレーム上の列番号または変数名のベクトル
df, # データフレーム
use=c("all.obs", "complete.obs", "pairwise.complete.obs"), # 欠損値の処理法
method=c("pearson", "kendall", "spearman"), # 計算する相関係数の種類
latex=TRUE, # LaTeX 形式で度数分布表を出力する（デフォルトは LaTeX 形式）
caption=NULL, # LaTeX 出力のときの表題
label=NULL, # LaTeX 出力のときのラベル
digits=3, # 相関係数を表示するときの小数点以下の桁数
output="", # ファイルに出力するときはファイル名（デフォルトはコンソールに出力）
encoding=getOption("encoding")) # ファイルに出力するときのエンコーディング（デフォルトは OS による）
{
# 下請け関数
getNum <- function(str, df) { # 変数名から列番号を得る
names <- colnames(df)
seq_along(names)[names %in% str]
}
outfunc <- function(r, n) { # 相関係数行列に基づき，若干の計算を行い，結果を出力する
nr <- nrow(r) # 相関係数行列のサイズ
cname <- colnames(r) # 列名（変数名）をベクトルに入れておく
rname <- rownames(r) # 行名（変数名）をベクトルに入れておく（正方行列を扱うので，colnames と同じでよい）
if (!is.array(n)) { # n が行列でない（欠損値をリストワイズ除去した場合）
nall <- sprintf(ifelse(latex, "，$n$=%i", "，n=%i"), n) # n を nall に文字列としてセットする（latex=TRUE のときは，$n$ とする）
n <- matrix(n, nr, nr) # n を 行列にする
}
else {
nall <- "" # nall を空にセットする
}
p <- matrix(0, nr, nr) # P 値を入れるための行列を作る
if (method.name == "ケンドールの順位相関係数") { # ケンドールの順位相関係数の場合の P 値の計算
z <- abs(r)/sqrt((4*n+10)/(9*n*(n-1))) # 検定統計量（対称行列だが，無駄に 2 倍の計算をしている）
p[upper.tri(p, diag=FALSE)] <- # P 値を計算（こちらは上三角行列だけ）
pnorm(z[upper.tri(z, diag=FALSE)], lower.tail=FALSE)*2
}
else { # ピアソンの積率相関係数，スピアマンの順位相関係数の場合の P 値の計算
t <- abs(r)*sqrt(n-2)/sqrt(1-r^2) # 検定統計量（対角成分は Inf になるが，P 値を計算するのには何の支障もない）
p[upper.tri(p, diag=FALSE)] <- # P 値を計算（こちらは上三角行列だけ）
pt(t[upper.tri(t, diag=FALSE)], n[upper.tri(n, diag=FALSE)]-2, lower.tail=FALSE)*2
}
p <- p+t(p) # 下側三角行列を補完し
diag(p) <- NA # 対角成分は NA を代入する
if (latex) { # LaTeX 形式で集計結果を出力する
cat("\n\\begin{table}[htbp]\n", file=output) # \begin{table}[htbp]
if (is.null(caption)) {
cat(sprintf("\\caption{%s行列（欠損値は%s%s）}\n", # \caption{スピアマンの順位相関係数行列（欠損値はリストワイズ除去，n = xx）}
method.name, use.name, nall), file=output) # --- のようになる
}
else {
cat(sprintf("\\caption{%s}\n", caption), file=output) # 指定した表題
}
if (!is.null(label)) {
cat(sprintf("\\label{%s}\n", label), file=output) # 指定したラベル
}
cat("\\centering\n", file=output) # \centering
cat("\\begin{tabular}{l", rep("c", nr), "} \\hline\n", # \begin{tabular}{cc…c} \hline
sep="", file=output)
cat("", sprintf("(%03i)", 1:nr), sep=" & ", file=output)# 列名は書けないので (001), (002) … のようにする
cat("\\\\ \\hline\n", file=output) # \\ \hline
for (i in 1:nr) { # 各行ごとに
cat(sprintf("(%03i)%s", # (00i)変数名 & 相関係数i1 & 相関係数i2 & …
i, rname[i]), sprintf(ifelse(r[i,] < 0, formatl, formatw), r[i,]), sep=" & ", file=output)
cat("\\\\\n", file=output) # \\
cat("$P$値", sprintf(formatw, p[i,]), # P 値 & 数値i1 & 数値i2 & …
sep=" & ", file=output)
if (use.name == "ペアワイズ除去") { # 欠損値をペアワイズ除去したのならば，
cat("\\\\\n", file=output) # \\
cat("$n$", n[i,], sep=" & ",file=output)# n & 数値i1 & 数値i2 & …
}
cat("\\\\ \\hline\n", file=output) # \\ \hline
}
cat("\\end{tabular}\n", file=output) # \end{tabular}
cat("\\end{table}\n", file=output) # \end{table}
}
else { # tab で区切って出力する
cat("表 ", method.name, # 表 ピアソンの積率相関係数行列（欠損値はペアワイズ除去）
"行列（欠損値は", use.name, nall, "）\n\n", # --- のようになる
sep="", file=output)
cat(" ", sprintf("(%03i)", 1:nr), # 列名は書けないので (001), (002) … のようにする
sep="\t", file=output, fill=TRUE)
for (i in 1:nr) { # 各行ごとに
cat(sprintf("(%03i)%s", # (00i)変数名 & 相関係数i1 & 相関係数i2 & …
i, rname[i]), sprintf(formatw, r[i,]),
sep="\t", file=output, fill=TRUE)
cat("P 値", sprintf(formatw, p[i,]), # P 値 & 数値i1 & 数値i2 & …
sep="\t", file=output, fill=TRUE)
if (use.name == "ペアワイズ除去") { # 欠損値をペアワイズ除去したのならば，
cat("n", n[i,], sep="\t", # n & 数値i1 & 数値i2 & …
file=output, fill=TRUE)
}
}
}
}
# 関数本体
if (output != "") { # 結果をファイルに出力する場合の処理
output <- file(output, open="w", encoding=encoding)
}
formatw <- paste("%.", digits, "f", sep="") # 相関係数の小数点以下の桁数を設定する(Word)
formatl <- paste("$%.", digits, "f$", sep="") # 相関係数の小数点以下の桁数を設定する(LaTeX)
method <- match.arg(method)
method.name <- switch(match.arg(method), # 引数 method と変数 method.name の対応付け
pearson = "ピアソンの積率相関係数",
kendall = "ケンドールの順位相関係数",
spearman = "スピアマンの順位相関係数")
if (method == "pearson" && !all(sapply(i, function(k) is.numeric(df[,k])))) {
warning("ピアソンの積率相関係数を計算するためには，全ての変数は間隔尺度・比尺度変数であることが必要です")
}
else if (!all(sapply(i, function(k) is.numeric(df[,k]) || is.ordered(df[,k])))) {
warning("順位相関係数を計算するためには，全ての変数は順序尺度・間隔尺度・比尺度変数であることが必要です")
}
if (is.character(i[1])) {
i <- getNum(i, df)
}
if (match.arg(use) != "pairwise.complete.obs") { # ペアワイズ除去でない場合（リストワイズ除去の場合）
use.name <- "リストワイズ除去" # 欠損値の処理法の名前
df2 <- df[, i] # データフレームから分析対象変数を抽出
df2 <- subset(df2, complete.cases(df2)) # 欠損値 NA を一つも含まないケースを抽出
n <- nrow(df2) # サンプルサイズ
if (n == 0) {
stop("欠損値をリストワイズ除去した結果，有効なデータがなくなりました。分析に使用する変数を再確認するか，ペアワイズ欠損値処理を選択してやり直してみてください。")
}
r <- cor(df2, method=method) # cor 関数により相関係数行列を求める
outfunc(r, n) # 出力関数 outfunc を呼び出す
}
else { # ペアワイズ除去の場合
use.name <- "ペアワイズ除去" # 欠損値の処理法の名前
ln <- length(i) # 分析対象とする変数の個数（相関係数行列のサイズ）
r <- n <- matrix(0, ln, ln)
for (i1 in 1:ln) {
for (j1 in i1:ln) { # 全ての変数の組み合わせについて，
df2 <- df[, c(i[i1], i[j1])] # データフレームから分析対象となる 2 変数を抽出
df2 <- subset(df2, complete.cases(df2)) # 2 変数ともに欠損値でないケースを抽出
n[i1, j1] <- nrow(df2) # サンプルサイズ
r[i1, j1] <- cor(df2[1], df2[2], method=method) # cor 関数により 2 変数間の相関係数を求める
}
}
r <- r+t(r) # もう半分の三角行列を完成させる
diag(r) <- 1 # 対角成分は 1 に決まっている
n <- n+t(n) # もう半分の三角行列を完成させる
diag(n) <- diag(n)/2 # このようにして作ると，対角成分は半分にしないといけない
rownames(r) <- colnames(r) <- colnames(df)[i] # 相関係数行列の行と列の名前を付ける
outfunc(r, n) # 出力関数 outfunc を呼び出す
}
if (output != "") { # 結果をファイルに出力した場合の後始末
close(output)
}
}
使用例
set.seed(12345)
df <- round(data.frame(matrix(rnorm(100), 20)), 3)
for (i in 1:20) df[sample(20, 1, replace=TRUE), sample(5, 1, replace=TRUE)] <- NA
df
mycor(1:5, df, latex=FALSE)
mycor(1:5, df, use="pairwise", method="kendall", latex=FALSE)
出力結果例
> set.seed(12345)
> df <- round(data.frame(matrix(rnorm(100), 20)), 3)
> for (i in 1:20) df[sample(20, 1, replace=TRUE), sample(5, 1, replace=TRUE)] <- NA
> df
X1 X2 X3 X4 X5
1 0.586 0.780 1.129 0.150 0.645
2 NA 1.456 -2.380 -1.343 1.043
3 NA -0.644 -1.060 0.553 -0.304
4 -0.453 -1.553 0.937 NA 2.477
5 0.606 -1.598 0.854 -0.587 0.971
6 -1.818 NA 1.461 -1.832 1.867
7 0.630 -0.482 -1.413 0.888 0.672
8 -0.276 0.620 NA 1.593 -0.308
9 -0.284 0.612 0.583 0.517 NA
10 -0.919 -0.162 -1.307 -1.296 0.825
11 NA NA -0.540 NA -0.964
12 1.817 NA NA -0.785 -0.855
13 0.371 2.049 0.054 -1.049 1.887
14 0.520 NA NA 2.331 -0.392
15 -0.751 0.254 -0.671 1.403 NA
16 0.817 0.491 0.278 NA 0.687
17 -0.886 -0.324 0.691 0.826 -0.505
18 -0.332 -1.662 NA -0.812 2.158
19 1.121 1.768 2.145 0.476 NA
20 0.299 0.026 -2.347 1.021 -0.695
> mycor(1:5, df, latex=FALSE)
表 ピアソンの積率相関係数行列（欠損値はリストワイズ除去，n=7）
(001) (002) (003) (004) (005)
(001)X1 1.000 0.076 0.076 0.166 0.324
P 値 NA 0.871 0.872 0.722 0.479
(002)X2 0.076 1.000 0.036 -0.235 0.389
P 値 0.871 NA 0.939 0.612 0.388
(003)X3 0.076 0.036 1.000 -0.228 0.304
P 値 0.872 0.939 NA 0.623 0.507
(004)X4 0.166 -0.235 -0.228 1.000 -0.778
P 値 0.722 0.612 0.623 NA 0.039
(005)X5 0.324 0.389 0.304 -0.778 1.000
P 値 0.479 0.388 0.507 0.039 NA
> mycor(1:5, df, use="pairwise", method="kendall", latex=FALSE)
表 ケンドールの順位相関係数行列（欠損値はペアワイズ除去）
(001) (002) (003) (004) (005)
(001)X1 1.000 0.253 0.051 0.124 -0.231
P 値 NA 0.208 0.807 0.520 0.250
n 17 14 13 15 14
(002)X2 0.253 1.000 0.033 -0.033 -0.103
P 値 0.208 NA 0.870 0.870 0.625
n 14 16 14 14 13
(003)X3 0.051 0.033 1.000 -0.154 0.205
P 値 0.807 0.870 NA 0.464 0.329
n 13 14 16 13 13
(004)X4 0.124 -0.033 -0.154 1.000 -0.495
P 値 0.870 0.464 0.625 NA 0.014
n 15 14 13 17 14
(005)X5 -0.231 -0.103 0.205 -0.495 1.000
P 値 0.520 0.250 0.329 0.014 NA
n 14 13 13 14 17
直前のページへ戻る E-mail to Shigenobu AOKI
