サイトトップ | Software Tips | RについてのTips

ポリコリック相関係数について補足

最終更新:2014年5月14日

講義テキストの第8節「2つのカテゴリ変数の独立性」に記載した通り,2つのカテゴリ変数が独立でない場合,どのくらいの関連があるのかを示したいことがある。

この目的には,疫学研究でよく使われるオッズ比やリスク比の他,量的な変数でいえば相関係数に当たるものとして,2つのカテゴリ変数の水準が2つずつの場合はファイ係数やコンティンジェンシー係数(vcdパッケージのassocstats()関数にクロス集計表を与える),どちらかが3つ以上の場合はポリコリック相関係数(polycorパッケージのpolychor()関数に2つのカテゴリ変数を与える。標準誤差が欲しいときはstd.err=TRUEオプションを付ける)を使うのが普通である。実行例を以下に示す。

> library(MASS) # これは基本パッケージなので,R本体と同時にインストールされるはず。
> library(polycor) # 未インストールの際は一度だけinstall.packages("polycor", dep=TRUE)
 要求されたパッケージ mvtnorm をロード中です 
 要求されたパッケージ sfsmisc をロード中です 
> library(vcd) # 未インストールの際は一度だけinstall.packages("vcd", dep=TRUE)
 要求されたパッケージ grid をロード中です
> str(survey) # MASSパッケージに入っているデータ
'data.frame':   237 obs. of  12 variables:
 $ Sex   : Factor w/ 2 levels "Female","Male": 1 2 2 2 2 1 2 1 2 2 ...
 $ Wr.Hnd: num  18.5 19.5 18 18.8 20 18 17.7 17 20 18.5 ...
 $ NW.Hnd: num  18 20.5 13.3 18.9 20 17.7 17.7 17.3 19.5 18.5 ...
 $ W.Hnd : Factor w/ 2 levels "Left","Right": 2 1 2 2 2 2 2 2 2 2 ...
 $ Fold  : Factor w/ 3 levels "L on R","Neither",..: 3 3 1 3 2 1 1 3 3 3 ...
 $ Pulse : int  92 104 87 NA 35 64 83 74 72 90 ...
 $ Clap  : Factor w/ 3 levels "Left","Neither",..: 1 1 2 2 3 3 3 3 3 3 ...
 $ Exer  : Factor w/ 3 levels "Freq","None",..: 3 2 2 2 3 3 1 1 3 3 ...
 $ Smoke : Factor w/ 4 levels "Heavy","Never",..: 2 4 3 2 2 2 2 2 2 2 ...
 $ Height: num  173 178 NA 160 165 ...
 $ M.I   : Factor w/ 2 levels "Imperial","Metric": 2 1 NA 2 2 1 1 2 2 2 ...
 $ Age   : num  18.2 17.6 16.9 20.3 23.7 ...
> # 以下は,EZRでは「統計解析」「名義変数の解析」
> # 「分割表の作成と群間の比率の比較(Fisherの正確検定)」と選ぶ
> xtabs(~Sex+W.Hnd, data=survey)
        W.Hnd
Sex      Left Right
  Female    7   110
  Male     10   108
> fisher.test(xtabs(~Sex+W.Hnd, data=survey))

        Fisher's Exact Test for Count Data

data:  xtabs(~Sex + W.Hnd, data = survey)
p-value = 0.6158
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.2140244 2.0876656
sample estimates:
odds ratio 
 0.6883665 

> assocstats(xtabs(~Sex+W.Hnd, data=survey))
                     X^2 df P(> X^2)
Likelihood Ratio 0.54629  1  0.45984
Pearson          0.54351  1  0.46098

Phi-Coefficient   : 0.048 
Contingency Coeff.: 0.048 
Cramer's V        : 0.048 
> xtabs(~Smoke+Exer, data=survey)
       Exer
Smoke   Freq None Some
  Heavy    7    1    3
  Never   87   18   84
  Occas   12    3    4
  Regul    9    1    7
> fisher.test(xtabs(~Smoke+Exer, data=survey))

        Fisher's Exact Test for Count Data

data:  xtabs(~Smoke + Exer, data = survey)
p-value = 0.4138
alternative hypothesis: two.sided

> polychor(survey$Smoke, survey$Exer, std.err=TRUE)

Polychoric Correlation, 2-step est. = -0.042 (0.1006)
Test of bivariate normality: Chisquare = 5.628, df = 5, p = 0.3441
> # カテゴリ順序を付け直す。
> # EZRでは「アクティブデータセット」「変数の操作」「因子水準を再順序化する」
> survey$Smoke <- factor(survey$Smoke, levels=c("Never", "Occas", "Regul", "Heavy"), ordered=TRUE)
> survey$Exer <- factor(survey$Exer, levels=c("None","Some","Freq"), ordered=TRUE)
> xtabs(~Smoke+Exer, data=survey)
        Exer
Smoke   None Some Freq
  Never   18   84   87
  Occas    3    4   12
  Regul    1    7    9
  Heavy    1    3    7
> polychor(survey$Smoke, survey$Exer, std.err=TRUE)

Polychoric Correlation, 2-step est. = 0.1269 (0.1011)
Test of bivariate normality: Chisquare = 4.274, df = 5, p = 0.5107

> fisher.test(xtabs(~Smoke+Exer, data=survey))

        Fisher's Exact Test for Count Data

data:  xtabs(~Smoke + Exer, data = survey)
p-value = 0.4138
alternative hypothesis: two.sided

正しい順序の順序変数にすると,ポリコリック相関係数は,順序を無視したときとは結果が変わるけれども,フィッシャーの正確検定の結果は変わらない。順序変数の間では,量的変数としての扱いを強制し,スピアマンの順位相関係数を求めることもできるが,同順位が多いためp値が正しく得られないという問題を生じるので,ポリコリック相関係数を使う方が良い。

> cor.test(as.numeric(survey$Smoke), as.numeric(survey$Exer), method="spearman")

       Spearman's rank correlation rho

data:  as.numeric(survey$Smoke) and as.numeric(survey$Exer)
S = 1994287, p-value = 0.1699
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho 
0.08964512 

リンクと引用について