読者です 読者をやめる 読者になる 読者になる

カテゴリカルデータ解析読書メモ(第4章)

オッズとオッズ比

2つの割合p_1p_2に対して、オッズ\frac{p}{1-p}の比\frac{p_2(1-p_1)}{p_1(1-p_2)}を考えると、このオッズ比が1より大ならp_1\lt p_2、1より小ならp_1\gt p_2である。
オッズ比は、ロジスティック回帰において回帰係数がオッズ比の対数に一致することや、患者対照研究のようにオッズ比のみが推定できる研究方法があるためにしばしば用いられる。pが小さい場合には、オッズ比は割合の比に近似するので、リスクが何倍になる、といった文脈で使われる場合もある。

独立性のカイ二乗検定

2つの二項割合p_1p_2について、帰無仮説をp_1=p_2とする。このとき、割合の推定量の差を分散のルートで補正した量の2乗をピアソンのカイ二乗統計量\chi ^2と呼び、データ数が大きい時に\chi^2が自由度1のカイ二乗分布で近似されるため、この近似によって棄却域やp値が計算できる。
このとき、より近似を良くするためにイェーツの補正と呼ばれる補正を行う場合がある。Rのchisq.test()ではデフォルトである。

> Ctable <- matrix(c(512, 1385, 567, 2419), nr = 2)
> chisq.test(Ctable)

	Pearson's Chi-squared test with Yates' continuity correction

data:  Ctable
X-squared = 42.679, df = 1, p-value = 6.449e-11

この補正を行ってもデータ数が少ないと近似が不正確になり、警告メッセージが表示される場合がある。

> Ctable <- matrix(c(2, 10, 3, 15), nr = 2)
> chisq.test(Ctable)

	Pearson's Chi-squared test

data:  Ctable
X-squared = 0, df = 1, p-value = 1

 警告メッセージ: 
 chisq.test(Ctable):   カイ自乗近似は不正確かもしれません 

フィッシャーの直接確率法

紅茶が先かミルクが先かという話は、「ミルクが先」が美味しいということで決着が付いているらしい*1が、ともかく、サンプル数が少ない場合に超幾何分布モデルによって結果が偶然に起こる確率を直接計算する方法をフィッシャーの直接確率法と呼ぶ。
Rではfisher.test()で行う。なお、信頼区間や推定値はオッズ比として表示される。

> table <- matrix(c(3, 1, 1, 3), nr = 2)
> fisher.test(table)

	Fisher's Exact Test for Count Data

data:  table
p-value = 0.4857
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
   0.2117329 621.9337505
sample estimates:
odds ratio 
  6.408309 

かつてはデータ数が大きくなると計算量が増えるという点が問題であったが、2x2程度の表であれば多少データ数が大きくても実行に支障はない。

> Ctable <- matrix(c(512, 1385, 567, 2419), nr = 2)
> fisher.test(Ctable)

	Fisher's Exact Test for Count Data

data:  Ctable
p-value = 8.63e-11
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 1.372396 1.811955
sample estimates:
odds ratio 
  1.576932 

関連性の指標

2x2表が独立ではない場合にどの程度の関連性があるかを調べる指標をまとめて計算できる関数assocstats()がパッケージvcdに入っている。また、oddsratio()でオッズ比も計算できる。

> library(MASS)
> library(vcd)
> Ctable <- xtabs(~ Sex + W.Hnd, data = survey)
> assocstats(Ctable)
                     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 
> oddsratio(Ctable)
log odds ratios for Sex and W.Hnd 

[1] -0.3750241
> oddsratio(Ctable, log = FALSE) # デフォルトでは対数表示
 odds ratios for Sex and W.Hnd 

[1] 0.6872727

2xJ表の解析

カテゴリカル変数が2つよりも大きい場合もchisq.test()は実行できる。

> Ctable <- xtabs(~ Treatment + Improved, data = Arthritis)
> Ctable
         Improved
Treatment None Some Marked
  Placebo   29    7      7
  Treated   13    7     21
> chisq.test(Ctable)

	Pearson's Chi-squared test

data:  Ctable
X-squared = 13.055, df = 2, p-value = 0.00146

ただし、このときは単にPlaceboとTreatedの分布が違う、程度のことしか言えない。どの値が帰無仮説からのずれが大きいかを調べる方法として、調整済み残差を用いる方法がある。

> Ctable <- xtabs(~ Treatment + Improved, data = Arthritis)
> result <- chisq.test(Ctable)
> v <- result$expected * (1 - rowSums(Ctable)/sum(Ctable)) %o%
+ (1 - colSums(Ctable)/sum(Ctable))
> (Ctable - result$expected) / sqrt(v)
         Improved
Treatment        None        Some      Marked
  Placebo  3.27419655 -0.09761768 -3.39563632
  Treated -3.27419655  0.09761768  3.39563632

この残差が1.96より大きければ帰無仮説からのずれが大きいと考える。

カテゴリカル変数に順序関係のある場合の解析

マンテル検定

カテゴリカル変数間に順序関係が有るとき、これに1, 2, 3...のようなスコアを与える検定手法としてマンテル検定がある。

manteltrend.test <- function(x, score = c(1:dim(x)[2])){
  ctotal <- colSums(x) # 列和
  rtotal <- rowSums(x) # 行和
  n <- sum(x) # データ総数
  expected <- (rtotal %o% ctotal) / n # 期待度数
  tscore <- (x - expected) %*% score
  v <- n * (score^2 %*% ctotal) - (score %*% ctotal)^2
  statistics <- n^2 * (n-1) * tscore[2]^2 / rtotal[1] / rtotal[2] / v
  p.value <- pchisq(statistics, 1, lower.tail = FALSE)
  list(statistics = statistics, p_value = p.value)
}

スコアは等間隔である必要は無いが、どのようなスコア間隔を用いるかは大勢に影響を与えない場合が多い。

> manteltrend.test(Ctable, score=0:2)
$statistics
         [,1]
[1,] 12.85902

$p_value
             [,1]
[1,] 0.0003358568

> manteltrend.test(Ctable, score=c(1, 2, 5))
$statistics
         [,1]
[1,] 12.60627

$p_value
             [,1]
[1,] 0.0003844559
Wilcoxonの順位和検定

Wilcoxonの順位和検定では、スコアを使わず、各データを順位に変換して検定を行う。Rではwilcox.test()で行えるが、あらかじめデータを1人1人のデータに戻しておく必要がある。

> x <- rep(col(Ctable)[1,], Ctable[1,])
> y <- rep(col(Ctable)[2,], Ctable[2,])
> wilcox.test(x, y)

	Wilcoxon rank sum test with continuity correction

data:  x and y
W = 517.5, p-value = 0.0003666
alternative hypothesis: true location shift is not equal to 0

 警告メッセージ: 
 wilcox.test.default(x, y): 
   タイがあるため、正確な p 値を計算することができません 
コクラン・アーミテージ検定

薬の投与量に応じて効果のありなしが分かれるようなデータがあったとして、効果ありの割合p_jが投与量c_jに対して、次のような関係にあったとする。
p_j = \alpha + \beta c_j
ここで\beta = 0であるとすると、投与量と効果ありの割合の間に関係がないということになる。このとき、ロジスティック回帰モデルのようにp_jを何らかの関数で変換したとしても、検定方法は変わらない。この検定はRではprop.trend.test()で実行できる。

> gd <- c(2, 3, 9, 15, 18)
> n <- rep(20, 5)
> prop.trend.test(gd, n, score = 1:5)

	Chi-squared Test for Trend in Proportions

data:  gd out of n ,
 using scores: 1 2 3 4 5
X-squared = 38.86, df = 1, p-value = 4.553e-10

I x J表の解析

独立性の検定

chisq.test()がこれまでと同様に使用できる。

> Ctable <- xtabs(Freq ~ Status + Alcohol, data = DanishWelfare)
> chisq.test(Ctable)

	Pearson's Chi-squared test

data:  Ctable
X-squared = 87.972, df = 4, p-value < 2.2e-16
線形関連の検定

行と列が共に順序カテゴリカルデータの場合、これらにスコアを与えて相関係数を求めることができる。相関係数をr、標本数をnとしたとき、
M^2 = (n - 1)r^2
は、標本数nが大きくなると近似的に自由度1のカイ二乗分布に近づく。このことを利用した検定を線形関連の検定と呼ぶ。
Rではcoinパッケージのlbl_test()関数で実行できる。

> library(coin)
> Ctable <- apply(jobsatisfaction, 1:2, sum) # 性別を合計
> Ctable
             Job.Satisfaction
Income        Very Dissatisfied A Little Satisfied Moderately Satisfied
  <5000                       2                  4                   13
  5000-15000                  2                  6                   22
  15000-25000                 0                  1                   15
  >25000                      0                  3                   13
             Job.Satisfaction
Income        Very Satisfied
  <5000                    3
  5000-15000               4
  15000-25000              8
  >25000                   8
> lbl_test(as.table(Ctable))

	Asymptotic Linear-by-Linear Association Test

data:  Job.Satisfaction (ordered) by
	 Income (<5000 < 5000-15000 < 15000-25000 < >25000)
Z = 2.7623, p-value = 0.005739
alternative hypothesis: two.sided

スコアは指定しなければ自動的に1, 2, 3, ...と設定されるが、指定したければlistで与える。

lbl_test(as.table(Ctable), 
         scores = list(Job.Satisfaction = 1:4,
                       Income = c(2.5, 10, 20, 30)))

なお、I x J表においてI=2の場合はマンテル検定と同じ結果が得られる。