カイ二乗検定（独立性の検定）

データの読み込みと準備

サンプルデータは、統計教育推進委員会が提供してくださっている大学生の睡眠時間と起床・就寝時刻[1]を使います。このデータには男女別の睡眠時間や就寝時間、起床時間が収集されていますが、ここでは、男女別の睡眠時間を使います。

厚生労働省の平成27年「国民健康・栄養調査」の結果[2]を見ると睡眠時間6時間未満と以上で分類して分析しているので、ここでも6時間未満を睡眠時間不十分（less）、6時間以上を睡眠時間十分（enough）と便宜的に分類してカイ二乗（χ²）検定を行うことにします。

R Console

> x <- read.table(pipe("pbpaste"), header=TRUE) # クリップボードからデータを読み込み、xに代入。ちなみにこちらはMac用のコードでWindowsは、read.table("clipboard", header=TRUE)。
> by(x$sleep, x$sex, summary) # by関数を使ってsex（性別ごと）のsleep（睡眠時間）のsummary（要約）を計算します。男性（sex: 0）では2時間から13時間、女性（sex: 1）では3時間から13時間の範囲。
x$sex: 0
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  2.000   6.000   7.000   7.067   8.000  13.000 
------------------------------------------------------------------ 
x$sex: 1
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  3.000   6.000   7.000   6.901   8.000  13.000

睡眠時間6時間未満と6時間以上でデータを分割して、それぞれにラベルを付けるためにはcut関数を使います。cut関数の説明はこちらのブログ記事が分かりやすかったです。

R Console

> x$sleep2 <- cut(x$sleep, breaks = c(2, 6, 13), labels = c("less", "enough"), right=TRUE, include.lowest=TRUE)
# breaksで区切り位置を指定。最小値2、6、最大値13を指定。labelsで2と6の間、6と13の間に対してラベルを付ける。rightで2と6区間における6、6と13区間における13（区間の右端）を区間に含むかどうかを指定、今回は2と6区間をlessと設定するが6時間未満で6を含まないのでright=FALSEと設定する。しかし、そうすると6と13区間で13が含まれないので、include.lowest=TRUEを指定し、13を含める。

ついでにsexで男性が0、女性が1に設定されているのを、それぞれm、fに置換してfactor型に変えてみます。

R Console

> x$sex2 <- gsub("0", "m", x$sex) # gsub関数で0をmに置換。
> x$sex2 <- gsub("1", "f", x$sex2) # gsub関数で1をfに置換。
> x$sex2 <- as.factor(x$sex2) # そのままでは文字列関数になっているので、as.factorで因子に変更。

[1] 統計教育推進委員会 – 統計活用提案のための教材サイト

[2] 厚生労働省 – 平成27年「国民健康・栄養調査」の結果

カイ二乗検定の実行

R Console

> y <- table(x$sex2, x$sleep2) # カイ二乗検定はクロス集計表の形式のデータを指定するので、table関数を使って2変数をまとめる。
> chisq.test(y) # デフォルトでは、イェーツ補正ありになっているが、イェーツ補正をしない場合はcorrect=FALSEと設定する。

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

data:  y
X-squared = 0.77526, df = 1, p-value = 0.3786

p値が0.3786なので、睡眠時間6時間未満の者は男女間で差はないことになります。

MEMO

イェーツ補正をするかどうかについては、サンプルサイズが小さい場合（50未満？）はイェーツ補正を行うようですが、詳しいことは分かりません。クロス集計表に5未満の数値があるときはフィッシャーの直接確率検定を行うようです。

関数chisq.testは、漸近的にχ²分布に従うことを利用して計算されるそうですが、正確に求めるためには、以下のようにフィッシャーの正確確率検定を行います。

R Console

> fisher.test(y)

	Fisher's Exact Test for Count Data

data:  y
p-value = 0.3669
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.7241556 2.5420650
sample estimates:
odds ratio 
  1.358115

効果量の計算

パッケージvcdに含まれるassocstats関数を使ってφを計算します。2×2分割表はφ、2×2以外はCramer’s Vを効果量として用います[1]。

R Console

> library(vcd) # パッケージvcdを読み込む。
> assocstats(y) # 効果量を計算。
                    X^2 df P(> X^2)
Likelihood Ratio 1.0566  1  0.30398
Pearson          1.0627  1  0.30259

Phi-Coefficient   : 0.072 
Contingency Coeff.: 0.072 
Cramer's V        : 0.072

φが0.072なので効果量からみても差があるとは言えないようです。

[1] 英語教育研究 – 研究論文における効果量の報告のために―基礎的概念と注意点―（水本篤）