[連載] フリーソフトによるデータ解析・マイニング 第45回

Rとカテゴリカルデータの操作と統計量

1.カテゴリカルデータ

  本欄で説明したデータ解析の方法のほとんどは、量的なデータを対象にしたものである。今月号から数回の誌面を用いてカテゴリカルデータ (categorical data) について説明を行うことにする。
  我々の周りには多くの物事の性質を数え上げる計数 (count) データがある。例えば、性別(男、女)、血液型(A,B,O,AB)やアンケート調査票の中の質問における選択肢などを集計したデータは、いずれも計数データである。このような計数データをカテゴリカルデータと呼び、性別や血液型などのカテゴリーをカテゴリカル変数 (categorical variable, 略して変数) と呼ぶ。
  カテゴリカル変数は、データの尺度によって2種類に分けることができる。1つはカテゴリカル変数が順序関係を持たない性別や血液型などのような吊義 (nominal) 尺度データであり、もう1つはアンケート調査質問票などでよく用いられている満足度に関する選択肢 (非常に満足する、満足する、なんともいえない、満足しない、非常に満足しない) のようなカテゴリーの順序に意味を持っている順序 (ordinal) 尺度データである。
  順序尺度は、カテゴリーの並べの順序に意味を持っているので、順序尺度データのために開発されたデータ解析方法は、吊義尺度には用いられないことに注意することが必要である。
  量的データは、ある一定の区間を区切りカテゴリカルデータに変換して解析を行ったり、カテゴリカルデータのカテゴリーに順序のスコアを付与(あるいはダミー変数0,1を導入)して量的のデータとして解析したりする場合もある。
  本稿では、Rによるカテゴリカルデータの解析のための操作方法や分割表の統計量や独立性の検定方法について説明を行う。カテゴリカルデータに関連する関数は幾つかのパッケージに用意されているが、本稿ではパッケージ stats、vcd の中の関数を用いる。

2.カテゴリカルデータの操作

  まずカテゴリカルデータの記録形式と基本操作について説明する。パッケージ vcd の中には86個体5変数の薬物実験に関するカテゴリカルデータフレーム Arthritis がある。その変数を次に示す。

 ID: 患者のID
 Treatment:
 実験グループの因子 (Placebo, Treated)
 Sex:
 性別 (Female, Male)
 Age:
 患者の年齢
 mproved:
 観察結果の因子 (None, Some, Marked)

  次のコマンドでデータを確認して見よう。関数 head を用いてデータ先頭の一部分を返すことができる。デフォルトでは先頭6行を返すように指定されている。

> library(vcd);data(Arthritis)
> head(Arthritis,n=4)

  ID Treatment Sex Age Improved
1 57 Treated Male 27 Some
2 46 Treated Male 29 None
3 77 Treated Male 30 None
4 17 Treated Male 32 Marked

  このようなデータの記録は1行が1つの個体 (ここでは1人の被実験者) であり、1列が1つの変数である。通常アンケートの調査票などの集計にも基本的にはこのような形式を用いる。
  このようなデータ形式はしばしば分割表 (contingency table) に集計して分析を行う。上記の形式から分割表に集計する関数としてパッケージ stats には関数 xtabs、パッケージ vcd には関数 structable がある。後者はフラットな分割表を返す。それぞれの書き式を次に示す。

xtabs(formula = ~., data …)
structable(formula, data,…)

  引数 formula には、応答変数と説明変数や分割表の表側、表頭の変数を指定する。データ Arthritis の変数 Treatment を表側、変数 Improved を表頭とした分割表を作成する例を示す。

> (ArtTI <- xtabs(~ Treatment + Improved, data = Arthritis))

Improved
Treatment None Some Marked
Placebo 29 7 7
Treated 13 7 21

  このよう分割表を2元分割表 (two-way table) と呼ぶ。
  引数 subset に条件を付けることにより、その条件に基づいた分割表を作成することができる。例えば、女性に限定した上記のような分割表を返すのには引数 subset = Sex = = "Female" を追加する。
  データ Arthritis の中の変数 Age は、量的な変数として扱うことも可能であり、カテゴリカルデータとして扱うことも可能である。
  仮に量的のデータとして扱うときには次のコマンドのように引数 formula を指定すると、分割表の各セルには年齢の合計が返される。

> (ArtTI.age<-xtabs(Age~ Treatment +Improved, data = Arthritis))

Improved
Treatment None Some Marked
Placebo 1439 402 403
Treated 648 397 1193

  分割表の各セルの度数は異なるので合計より、度数当たりの平均度数を求めた方がよい場合もある。そのコマンドを次に示す。

> ArtTI.age/ArtTI

Improved
Treatment None Some Marked
Placebo 49.62069 57.42857 57.57143
Treated 49.84615 56.71429 56.80952

  次のコマンドのように書くことで2つの変数 Sex と Treatment を表側 (あるいは表頭) に配置することができる。

> xtabs(~ I(Sex:Treatment) + Improved, data = Arthritis)

Improved
I(Sex:Treatment) None Some Marked
Female:Placebo 19 7 6
Female:Treated 6 5 16
Male:Placebo 10 0 1
Male:Treated 7 2 5

  引数 formula の説明変数を3つに指定すると次のような3元分割表が返される。3元分割表には3つの添え字で分割表を操作する。添え字の順序は引数 formula で指定した変数の順序と一致する。

> (ArtTIS<- xtabs(~ Treatment +Improved+Sex, data = Arthritis))

, , Sex = Female

Improved
Treatment None Some Marked
Placebo 19 7 6
Treated 6 5 16

, , Sex = Male

Improved
Treatment None Some Marked
Placebo 10 0 1
Treated 7 2 5

> dim(ArtTIS)
[1] 2 3 2

> ArtTIS[,,1]  #女の分割表
<結果は省略>

> ArtTIS[,,2]  #男の分割表
<結果は省略>

  多元分割表をフラット分割表に変換する関数としてパッケージ stats にはftable、パッケージ vcd には structable がある。その例を次に示す。

> ftable(ArtTIS)
                Sex Female Male
Treatment Improved
Placebo
None
19 10
 
Some
7 0
 
Marked
6 1
Treated
None
6 7
 
Some
5 2
 
Marked
16 5

  表側、表頭の組み合わせなどは引数 row.vars , col.vars を次のように指定することができる。

> ftable(ArtTIS,row.vars ="Improved")

Treatment  Placebo Treated  
Sex Female Male Female Male
Improved
     
None
19 10 6 7
Some
7 0 5 2
Marked
6 1 16 5

  関数 structable の使用例を次に示す。

> structable(ArtTIS)

  Improved None Some Marked
Treatment
Sex
     
Placebo
Female
19 7 6
 
Male
10 0 1
Treated
Female
6 5 16
 
Male
7 2 5

  場合によっては、多元分割表の変数を併合して分析すると都合がよい場合がある。Rには配列を併合する関数 margin.table がある。3元分割表 ArtTIS を2元に併合する例を次に示す。

> margin.table(ArtTIS,1:2)

Improved
Treatment None Some Marked
Placebo 29 7 7
Treated 13 7 21

  Rでは関数 prop.table を用いて分割表の各セルの比率を求めることができる。分割表の度数の総合計を基準とする各セルの比率は次のように求める。

> ArtTI.p<-prop.table(ArtTI)
> round(ArtTI.p,3)

Improved
Treatment None Some Marked
Placebo 0.345 0.083 0.083
Treated 0.155 0.083 0.250

  列の合計を基準とするときには引数 margin=1,行の合計を基準とするときには引数 margin=2 を追加する。

> prop.table(ArtTI, margin=1)

Improved
Treatment None Some Marked
Placebo 0.6744186 0.1627907 0.1627907
Treated 0.3170732 0.1707317 0.5121951

3.統計量と仮説検定

  説明の便利のため、まず2元分割表の一般形式を表1に示す。

表1 2元分割表の一般形式と周辺度数
  b1 b2 bj b1
a1 n11 n12 n1j n1c n1+
a2 n21 n22 n2j n2c n2+
ai ni1 ni2 nij nic ni+
ar nr1 nr2 nrj nrc nr+
n1+ n2+ ni+ nr+ n++

 2元分割表の分析に用いる統計量の中で最も広く知られているのはピアソンの χ2 統計量 (Pearson's chi-squared statistics) である。略してカイ2乗統計量と呼ぶことにする。カイ2乗統計量は次の式で定義されている。

カイ2乗統計量

  この式で得られた統計量は自由度 ( r - 1 ) ( c - 1 ) の χ2 分布に従うことが知られている。式の中の nij は分割表のij 列セルの度数であり、Eijij 列セルの期待度数 (expected frequencies) である。
  期待度数 Eij は、i 行の和 ni+j 列和 n+j、分割表の度数の総合計 n++ で次のように求める。分割表の ni+ , n+j , n++ を周辺度数とも呼ぶ。

期待度数Eij

  パッケージ vcd には、2元分割表の周辺度数を求める関数 mar_table、期待度数を求める関数 independence_table がある。
  第2節で作成した2元分割表 ArtTI を用いた例を次に示す。

> mar_table(ArtTI)

Improved
Treatment None Some Marked TOTAL
Placebo 29 7 7 43
Treated 13 7 21 41
TOTAL 42 14 28 84

> independence_table(ArtTI)

Improved
Treatment None Some Marked
Placebo 21.5 7.166667 14.33333
Treated 20.5 6.833333 13.66667

  分割表の各セルの値

分割表の各セルの値

をピアソン残差 (Person residuals) と呼ぶ。ピアソン残差は近似的に標準正規分布に従うことが知られている。
  パッケージ stats には、適合検定や分割表の独立性の検定に必要となるカイ2乗統計量を求める関数 chisq.test がり、カイ2乗統計量だけではなく、期待度数、ピアソン残差も返す。2元分割表 ArtTI を用いた例を次に示す。

> (ArtTI.chi<-chisq.test(ArtTI))

    Pearson's Chi-squared test

data: ArtTI
X-squared = 13.055, df = 2, p-value = 0.001463

> summary(ArtTI.chi)

  Length Class Mode
statistic
1 -none- numeric
parameter
1 -none- numeric
p.value
1 -none- numeric
method
1 -none- character
data.name
1 -none- character
observed
6 xtabs numeric
expected
6 -none- numeric
residuals
6 xtabs numeric

> ArtTI.chi$expected

Improved
Treatment None Some Marked
Placebo 21.5 7.166667 14.33333
Treated 20.5 6.833333 13.66667

> round(ArtTI.chi$residuals,2)

Improved
Treatment None Some Marked
Placebo 1.62 -0.06 -1.94
Treated -1.66 0.06 1.98

  関数 chisq.test のデフォルトは、2 × 2 分割表のときには,連続性の補正を行うように引数が "correct=TRUE" に指定されている。連続修正は次の式を用いている。

連続修正

  分割表のカイ2自乗統計量の近似値として尤度比統計量 (likelihood ratio statistics)

尤度比統計量

がある。分割表の独立検定にはカイ2乗統計量が多く用いられているが、モデルの比較や選択には尤度比統計量が多く用いられている。
  パッケージ vcd の中の関数 assocstats は分割表のカイ2乗統計量、尤度比統計量、ファイ連関係数 φ、ピアソンの連関係数 C、グラメール (Gramer) の連関係数 V を返す。これらの連関係数は次のように定義されている。式の中の k は行と列の数の中の少ない方の数を取る。

各連関係数

  分割表 ArtTI を用いた例を次に示す。

> assocstats (ArtTI)

  X^2 df P(> X^2)
Likelihood Ratio
13.530 2 0.0011536
Pearson
13.055 2 0.0014626

Phi-Coefficient
: 0.394
Contingency Coeff.
: 0.367
Cramer's V
: 0.394

  表2のような2×2の分割表は、しばしば特殊なケースとして扱われている。

表2 2×2の分割表
  b1 b2
a1 n11 n12 n1+
a2 n21 n22 n2+
n+1 n+2 n++

  2×2の分割表のピアソンカイ2乗値は

ピアソンカイ2乗値

で計算されるが、イェーツ (Yates) はよりカイ2乗分布に近似するように次のように連続性の補正を施して用いることを提唱した。

連続性の補正ピアソンカイ2乗値

  セルの度数が非常に低く、期待値が5以下のセルが全セルの25%以上を含むときには上記の式によるカイ2乗独立性の検定は上適切である。そこでフィッシャー (Fisher) は超幾何分布を用いて、2×2の分割表の正確検定統計量の計算式を導出した。

2×2の分割表の正確検定統計量の計算式

  この確率を用いた検定をフィッシャーの正確検定 (Fisher's exact test) と呼ぶ。Rにはフィッシャーの正確検定値を求める関数 fisher.test がある。
  ここではフィッシャーが1935年に公表した有吊な紅茶の実験データを用いることにする。フィッシャーの紅茶の実験は「紅茶にミルクを注ぐ《と「ミルクに紅茶を注ぐ《という異なる2つの方法で入れた紅茶を、その味から紅茶の入れ方を推定するテストである。2つの入れ方によるそれぞれ4杯の紅茶を、味からその入れ方を識別することが可能であると吊乗る女性に識別テストを行ったデータを表3に示す。

表3 フィッシャーの紅茶実験データ
(表の中のミルク、紅茶は先に入れたもの)
実際 識別結果
ミルク 紅茶 合計
ミルク 3 1 4
紅茶 1 3 4
合計 4 4 8

  本例では、フィッシャーの両側検定は p(0) + p(1) + p(3) + p (4) 、片側検定は p(0) + p (1) (左側検定)、あるいは p(3) + p (4) (右側検定)値を用いる。
  表3のデータを用いたフィッシャーの正確検定値を求めるコマンドを次に示す。

> (TeaTasting<-matrix(c(3,1,1,3),nr = 2,dimnames = list(実際 = c("ミルク","紅茶"),識別 = c("ミルク", "紅茶"))))

識別
実際 ミルク 紅茶
ミルク 3 1
紅茶 1 3

> fisher.test(TeaTasting, alternative = "greater") #片側検定統計量を算出

  Fisher's Exact Test for Count Data
data: TeaTasting
p-value = 0.2429
<後略>

  この結果では、独立性の仮説が棄却されないため、残念ながらこの小規模のテストでは、紅茶の味からその入れ方を識別することが可能であるという女性の主張は統計的に認めることはできない。
  表側 a と表頭 b が独立であるならば、θ が1になるはずである。この比をオッズ比 (Odds Ratios) と呼ぶ。フィッシャーの正確検定はオッズ比 θ=1 である仮説検定である。関数 fisher.test では、両側検定がデフォルトに指定されている。
  2×2の分割表の場合はオッズ比の検定が多く用いられている。パッケージ vcd にはオッズ比に関連する関数 oddsratio, woolf_test などがある。誌面上の都合によりその使用例は省略する。
  カテゴリカルデータによく用いられている幾つか検定に関する関数を表4に示す。パッケージ Hmisc の中の関数 summary.formula は、分割表の作成や比率と検定統計量などを求めることができる。

表4 検定に関する幾つかの関数
関数 機能 パッケージ
 assocstats
 分割表の統計量と検定:カイ2乗、尤度比、連関係数
vcd
 binom.test
 二項分布の検定
stats
 coindep_test
 条件付き独立検定
vcd
 fisher.test
 Fisherの直接検定
stats
 goodfit
 適合検定
vcd
 kruskal.test
 ランク和の検定
stats
 mantelhaen.test
 Cochran-Mantel -Haenszelの3元分割表のカイ2乗検定
stats
 mcnemArtest
 正方形の分割表のMcnemarのカイ2乗検定
stats
 oddsratio
 オッズ比と検定統計量
vcd
 prop.test
 比率の検定
stats
 summary.formula
 分割表の統計量など
Hmisc
 woolf_test
 同種の2×2×k分割表のWoolfのオッズ比検定
vcd