相関係数の比較[Stata]

Stataでも相関係数の計算をしてみた。こちらのエントリ尺度水準に適した相関係数とシミュレーションをStataで行ったバージョンである。

データはこちらからダウンロード、もしくは下部にRでの作成方法を掲載しているので、そのままRで走らせると、dta形式のデータが作成される。

ピアソンの積率相関係数

cor Math Science

0.5605

ピアソンの積率相関係数相関係数を無理やり出してみた。本来は他の技法で導くべきものである。

. cor Math_4n Science_4n // 順序x順序
0.4942
. cor Math Science_4n // 連続x順序
0.5032
. cor Math_2n Science_2n //  2値x2値
0.3763
. cor Math_4 Science_2 //  順序x2値
0.4301
. cor Math Science_2n // 連続x2値
0.4246

適した方法は以下の方法である。

  • ピアソンの積率相関係数: 連続変数と連続変数
  • ポリコリック相関係数: 順序変数と順序変数
  • ポリシリアル相関係数: 順序変数と連続変数
  • 順位双列相関係数: 2値データと順序変数
  • テトラコリック相関係数: 2値データと2値データ
  • 点双列相関係数: 2値データと連続変数

ポリコリック相関係数

順序変数と順序変数の場合。

Stata13以下の場合

Kolenikovさんが作った下記のパッケージをインストールが必要。
polychoric, by any other 'namelist'

10未満の数を自動的にカテゴリカル変数、10以上を連続変数として捉えるように設計されている。

findit polychoric

こちらからインストール。

/* 順序x順序: ポリコリック相関係数 */
. polychoric Math_4 Science_4

Variables :  Math_4 Science_4
Type :       polychoric
Rho        = .5888424
S.e.       = .02450705
Goodness of fit tests:
Pearson G2 = 4.2263446, Prob( >chi2(8)) = .83614614
LR X2      = 4.2295511, Prob( >chi2(8)) = .83584131

Stata16ではパッケージのインストールなしで同じコマンドで出力が可能。(追記:2019/10/16)

ポリシリアル相関係数

順序変数と連続変数の場合。

. polychoric Science_4 Math

Variables :  Science_4 Math
Type :       polyserial
Rho        = .57925495
S.e.       = .02113936

Stata16ではパッケージのインストールなしで同じコマンドで出力が可能。(追記:2019/10/16)

順位双列相関係数

2値データと順序変数の場合。
まずはソマーズのDから。
Roger Newsonさん(http://www.rogernewsonresources.org.uk/stata.htm)のパッケージをインストールする必要がある。

ssc inst somersd, replace

書き方はsomersdと書いて先に2値の名義尺度、2つ目に順位尺度を書き入れる。

.somersd Science_2 Math_4

------------------------------------------------------------------------------
             |              Jackknife
   Science_2 |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
      Math_4 |   .5183179   .0294415    17.60   0.000     .4606136    .5760223
------------------------------------------------------------------------------

Stata16ではパッケージのインストールなしで同じコマンドで出力が可能。(追記:2019/10/16)

2値データと順序変数はポリコリック相関係数でもあるので、polychoricパッケージでも計算が可能である。

. polychoric Science_2 Math_4

Variables :  Science_2 Math_4
Type :       polychoric
Rho        = .59908999
S.e.       = .0300931
Goodness of fit tests:
Pearson G2 = 2.8323497, Prob( >chi2(2)) = .24264038
LR X2      = 2.83805, Prob( >chi2(2)) = .2419498

方法が異なるので値が異なる。 Stata16ではパッケージのインストールなしで同じコマンドで出力が可能。(追記:2019/10/16)

テトラコリック相関係数

2値データと2値データの場合。テトラコリックはポリコリックの2値の場合なので、特に何も考えず、polychoricで計算すればよい。

. polychoric Math_2 Science_2

Variables :  Math_2 Science_2
Type :       polychoric
Rho        = .60875007
S.e.       = .03596437
Goodness of fit tests:
Pearson G2 = ., Prob( >chi2(.)) = .
LR X2      = ., Prob( >chi2(.)) = .

Stata16ではパッケージのインストールなしで同じコマンドで出力が可能。(追記:2019/10/16)

Stata13ではテトラコリック相関係数はpolychoricパッケージを使わず計算できた。ただし2値の方が0/1でないと動作しないので、作成したサンプルデータは1/2になっていたので、リコードをする必要がある。

. recode Math_2 Science_2 (1=0) (2=1) // 上書きのリコード

結果は以下。

. tetrachoric Math_2 Science_2

   Number of obs =    1000
 Tetrachoric rho =       0.6087
       Std error =       0.0360

Test of Ho: Math_2 and Science_2 are independent
 2-sided exact P =       0.0000

値はpolychoricパッケージで出力すると同じである。

双列相関係数

2値データと連続変数の場合。

. polychoric Science_2 Math

Variables :  Science_2 Math
Type :       polyserial
Rho        = .58812229
S.e.       = .02720657

イシリアルはポリシリアルの2値の場合なので、ポリシリアルとして計算されている。 Stata16ではパッケージのインストールなしで同じコマンドで出力が可能。(追記:2019/10/16)

点双列相関係数

点双列相関係数はポリシリアル相関係数の一つのバリエーションと捉え方とピアソンの積率相関係数のバリエーションと捉え方があるようだ。ポリシリアル相関係数の一つのバリエーションと考えた方が実用的だが、ピアソンのバリエーションの方も紹介しておこう。ちなみに、ここで実用的といっているのは、連続変数でピアソンの積率相関係数を行った値に近くなるという意味である。

調べてみたところStataでは2つ方法があるようだ。 John A. Andersonによるパッケージ。Stata8のころからあるパッケージらしい。
https://www.stata.com/stb/stb17/sg20/pbis.hlp

help pbis

2値データは0/1でないと作動しないので、リコードが必要。

recode Math_2 Science_2 (1=0) (2=1) // 0/1に変更
pbis Science_2 Math

結果は以下。

------------------+------------------+------------------+------------------+
Coef.= 0.4244          t= 14.8072        P>|t| = 0.0001        df=    998

もう一つはEffect sizesを利用する方法。オプションにpbcorrをつけると、点双列相関係数が出力される。Stata13では標準機能である。

esize twosample Math, by(Science_2) pbcorr
                                   Obs per group:
---------------------------------------------------------
        Effect Size |   Estimate     [95% Conf. Interval]
--------------------+------------------------------------
   Point-Biserial r |  -.4246202    -.4712231   -.3742119
---------------------------------------------------------

ピアソンの積率相関係数と同じ値になるはずなのだが、pbisパッケージの値は微妙に違っており、Effect sizesを利用した方法は数字はあっているが符号が逆になる。2つとも微妙な仕上がりである。

結果

変数 相関係数 相関係数 ピアソンで計算した場合
連続x連続 ピアソンの積率相関係数 0.5605 -
順序x順序 ポリコリック相関係数 .5888424 0.4942
連続x順序 ポリシリアル相関係数 .57925495 0.5032
2値x順序 順位双列相関係数 .59908999 0.3763
2値x順序 順位双列相関係数(ソマーズのD) .5075835 0.3763
2値x2値 テトラコリック相関係数 .60875007 0.4301
連続x2値 双列相関係数(ポリシリアル) .58812229 0.4246
連続x2値 点双列相関係数(ピアソン) 0.4246 0.4246

データの作成

サンプルデータはRで作成している。こちらで作成したものより、点数の分散を大きくしてある。理由は点数が10点以上離れないと、ポリシリアル相関係数の計算ができないからである。

library(mvtnorm)
set.seed(1234) # 乱数シードの固定
N <-1000 # サンプルサイズの指定
Sigma = matrix(c(1,0.54,0.54,1), byrow=TRUE, ncol=2) # c()の中に相関係数を入れる。
mu <- c(10, 10) # 平均値
d1 <- data.frame(rmvnorm(n=N, mean=mu, sigma=Sigma)) # 乱数の作成
d1 <- d1*5
d1 <- round(d1) # まるめ(四捨五入で整数値に)
colnames(d1) <- c("Math","Science") # 列名を挿入
summary(d1)

4分位。順序変数に。

library(OneR)
d2 <- bin(d1, nbins = 4, labels = c(1,2,3,4),  method = "content")
library(dplyr)
d2 <- d2 %>% mutate_if(is.factor, as.ordered)
colnames(d2) <- c("Math.4","Science.4") # 四分位と変数名を変更
d2_num <- d2 %>% mutate_if(is.factor, as.numeric)
colnames(d2_num) <- c("Math.4n","Science.4n") # 四分位の連続と変数名を変更

2分位(中央値)。Factor型に。

d3 <- bin(d1, nbins = 2, labels = c(1,2),  method = "content")
colnames(d3) <- c("Math.2","Science.2") # 二分位と変数名を変更
d3_num <- d3 %>% mutate_if(is.factor, as.numeric)
colnames(d3_num) <- c("Math.2n","Science.2n") # 二分位の連続と変数名を変更

データの書き出し

data <- cbind(d1, d2, d2_num, d3, d3_num)
library(foreign)
write.dta(data, file = "TestScore.dta")