Stataでも相関係数の計算をしてみた。こちらのエントリ尺度水準に適した相関係数とシミュレーションをStataで行ったバージョンである。
データはこちらからダウンロード、もしくは下部にRでの作成方法を掲載しているので、そのままRで走らせると、dta形式のデータが作成される。
ピアソンの積率相関係数
cor Math Science
0.5605
ピアソンの積率相関係数で相関係数を無理やり出してみた。本来は他の技法で導くべきものである。
. cor Math_4n Science_4n // 順序ｘ順序 0.4942 . cor Math Science_4n // 連続ｘ順序 0.5032 . cor Math_2n Science_2n // ２値ｘ２値 0.3763 . cor Math_4 Science_2 // 順序ｘ２値 0.4301 . cor Math Science_2n // 連続ｘ２値 0.4246
適した方法は以下の方法である。
- ピアソンの積率相関係数: 連続変数と連続変数
- ポリコリック相関係数: 順序変数と順序変数
- ポリシリアル相関係数: 順序変数と連続変数
- 順位双列相関係数: ２値データと順序変数
- テトラコリック相関係数: ２値データと2値データ
- 点双列相関係数: ２値データと連続変数
ポリコリック相関係数
順序変数と順序変数の場合。
Stata13以下の場合
Kolenikovさんが作った下記のパッケージをインストールが必要。
polychoric, by any other 'namelist'
10未満の数を自動的にカテゴリカル変数、10以上を連続変数として捉えるように設計されている。
findit polychoric
こちらからインストール。
/* 順序ｘ順序: ポリコリック相関係数 */ . 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)
順位双列相関係数
２値データと順序変数の場合。
まずはソマーズのＤから。
Roger Newsonさん(http://www.rogernewsonresources.org.uk/stata.htm)のパッケージをインストールする必要がある。
ssc inst somersd, replace
書き方はsomersdと書いて先に２値の名義尺度、２つ目に順位尺度を書き入れる。
.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)
２値データと順序変数はポリコリック相関係数でもあるので、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値データの場合。テトラコリックはポリコリックの２値の場合なので、特に何も考えず、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パッケージを使わず計算できた。ただし２値の方が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パッケージで出力すると同じである。
双列相関係数
２値データと連続変数の場合。
. polychoric Science_2 Math Variables : Science_2 Math Type : polyserial Rho = .58812229 S.e. = .02720657
バイシリアルはポリシリアルの２値の場合なので、ポリシリアルとして計算されている。 Stata16ではパッケージのインストールなしで同じコマンドで出力が可能。(追記：2019/10/16)
点双列相関係数
点双列相関係数はポリシリアル相関係数の一つのバリエーションと捉え方とピアソンの積率相関係数のバリエーションと捉え方があるようだ。ポリシリアル相関係数の一つのバリエーションと考えた方が実用的だが、ピアソンのバリエーションの方も紹介しておこう。ちなみに、ここで実用的といっているのは、連続変数でピアソンの積率相関係数を行った値に近くなるという意味である。
調べてみたところStataでは２つ方法があるようだ。
John A. Andersonによるパッケージ。Stata8のころからあるパッケージらしい。
https://www.stata.com/stb/stb17/sg20/pbis.hlp
help pbis
２値データは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つとも微妙な仕上がりである。
結果
| 変数 | 相関係数名 | 相関係数 | ピアソンで計算した場合 |
|---|---|---|---|
| 連続ｘ連続 | ピアソンの積率相関係数 | 0.5605 | - |
| 順序ｘ順序 | ポリコリック相関係数 | .5888424 | 0.4942 |
| 連続ｘ順序 | ポリシリアル相関係数 | .57925495 | 0.5032 |
| ２値ｘ順序 | 順位双列相関係数 | .59908999 | 0.3763 |
| ２値ｘ順序 | 順位双列相関係数(ソマーズのＤ) | .5075835 | 0.3763 |
| ２値ｘ２値 | テトラコリック相関係数 | .60875007 | 0.4301 |
| 連続ｘ２値 | 双列相関係数(ポリシリアル) | .58812229 | 0.4246 |
| 連続ｘ２値 | 点双列相関係数(ピアソン) | 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)
４分位。順序変数に。
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")
