[連載]フリーソフトによるデータ解析・マイニング第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:   患者の年齢
Improved: 観察結果の因子(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元分割表の分析に用いる統計量の中で最も広く知られているのはピアソンのX2統計量(Pearson's chi-squared statistics)である。略してカイ2乗統計量と呼ぶことにする。カイ2乗統計量は次の式で定義されている。
カイ2乗統計量
 この式で得られた統計量は自由度(r-1)(c-1)のX2分布に従うことが知られている。式の中の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