相関係数行列の計算と検定     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

Made with Macintosh