サイトトップ | 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