【Rで統計】カイ二乗検定

カイ二乗検定とは

カイχ二乗検定というと、一般にピアソンのカイ二乗検定を指す事が多い。

仮説によって、適合度検定や独立性の検定などがある。

基本はｔ検定と似ているけれど、ｔ検定が「２つの群の連続変数（numerical variable）の平均値の比較検定」（帰無仮説（null hypothesis）は、「平均値に差がない」）を検証するのに対して、

カイ二乗検定は「２群間のカテゴリカル変数（categorical variables）の比較検定」（帰無仮説（null hypothesis）は、「相関がない」）を検証する。

カイ二乗検定は、Yes/No で表されるようなカテゴリカル変数同士の相関をみる場合などによく使われる。

２群間のカテゴリカル変数の比較検定だと、

• カイ二乗検定
• フィッシャーの正確確率検定

の２つが挙げられるが、カイ二乗検定は近似したP値を算出するのに対して、フィッシャーの正確確率検定は直接計算したP値を算出する。

カイ二乗検定とフィッシャーの正確確率検定の使い分け

基本的にはほぼ同じ値にあるけれども、サンプル数の少ない研究では、カイ二乗検定とフィッシャーの正確確率検定との誤差が大きくなる。

一般的に、２ｘ２表のとき、

• 症例数が２０例以下の場合
• または症例数が４０以下で隠せるの期待値が５以下のマスが存在する場合

には、フィッシャーの正確確率検定を選択した方がよいとされている。

カイ二乗検定を行う

必要なパッケージのインストール

カイ二乗検定は、R標準パッケージの

chisq.test()

で求める事ができる。

それ以外にも、データ整理で下記パッケージがあると便利。

{dityverse} パッケージ

データ整理ではほぼ必須。

install.packages("tidyverse")
library(tidyverse)

{epitools} パッケージ

上記のようなテーブル型のデータをデータフレームに展開する時、便利。

install.packages("epitools")
library(epitools)

{flextable} パッケージ

検定結果を表として出力するため。

install.packages("flextable")
library(flextable)

{webshot} パッケージ

最終的に出てきた表を画像として出力する時に必要。

install.packages("webshot")
library(webshot)

独立性の検定

2つの観察（カテゴリ変数）に関連性があるか（or 互いに独立か）どうかを検定する。

計算値はエクセルでも簡単に計算できるけど、

今回は、Rで行う方法。

データ準備

今回は練習を兼ねて、上記例のクロス集計表作成から。

１．２ｘ２のマトリックスを作る。

# 2x2のマトリックスを作る
x <- matrix(data=c(15, 620, 485, 655), byrow=TRUE, nc=2)

２．それぞれの項目に "dimnames" を指定する。

# dimnamesを指定する
dimnames(x) <- list(c("Yes", "No" ), c("Yes", "No" ))
names(dimnames(x)) <- c("Disease", "Treatment" )

カイ二乗検定：独立性の検定

独立性の検定は、

chisq.test

の因数ｘに分割表を表す行列またはデータフレームを代入する

chisq.test(x)

か、引数ｘとｙにそれぞれのベクトルを代入する。

chisq.test(x, y)

マトリックスを代入する場合

# カイ二乗検定：独立性の検定、マトリックスを使う場合
chisq.test(x)

データフレームとして代入する場合

# カイ二乗検定：独立性の検定、データフレームを使う場合
chisq.test(data.frame(x))

x, y のベクトルを代入する場合

クロス集計表を、

epitools::expand()

でｘ，ｙのデータフレームに展開する。

# ２ｘ２のマトリックスを作る
x <- matrix(data=c(15, 620, 485, 655), byrow=TRUE, nc=2)

# dimnmaes etc. を指定
dimnames(x) <- list(c("Yes", "No"), c("Yes", "No"))
names(dimnames(x)) <- c("Disease", "Treatment")

# データフレームに展開
df <- expand.table(x)
df

# 確認
head(df)

２列のデータフレームに変換された。

※ epitools::expand() でデータフレームに展開するためには、それぞれのセルに dimnames を指定しておく必要がある。

で、これを

chisq.test()

に代入。

# カイ二乗検定、データフレーム
chisq.test(df$Treatment, df$Disease)

上記の結果は、いずれも P < 2.22-16 となっており、「Treatment と Disease の分布は独立ではない（つまり関係がある）」といえる。

P値以外のパラメーターにアクセスする方法

# カイ二乗検定、結果をlistに代入
Res <- chisq.test(df$Treatment, df$Disease)            

# 検定統計量
Res$statistic          # もしくは Res[["statistic"]]

# 自由度
Res$parameter          # もしくは Res[["parameter"]]

# P値
Res$p.value            # もしくは Res[["p.value"]]

• カイ二乗値: 323.4
• 自由度: 1
• P値: <0.001

という結果が得られた。

表として出力

出力結果を

as.data.frame

でデータフレームに変換し、

dplyr::mutate

で行名を追加し、

dplyr::select

で列を並び替える。

# データフレームに
Res2 <- as.data.frame(c(Res$statistic, Res$parameter, Res$p.value, Res$method, Res$data.name))                        
colnames(Res2) <- "Chi square test"

Res3 <- Res2 %>%
as.data.frame(row.names = c("statistic", "parameter", "p.value", "method", "data name")) %>%
mutate(" " = rownames(Res3)) %>%
select(" ", everything())

{flextable}

で表として出力する。

# 表にする
Res3 %>%                        
flextable::regulartable() %>%
flextable::theme_booktabs() %>%
flextable::font(fontname = "serif", part = "all") %>%
flextable::autofit()

docxで保存する場合は、

flextable::save_as_docx

を使用。

私は出力場所を "Export" という別フォルダにまとめたいので、下記のように出力。

# docx で保存
save_as_docx(Res3, path = "Export/Chi_square_test.docx", zoom = 3, expand = 10, webshot = "webshot")            

write.csv

や

readr::write_csv

を使って、csv で保存してもOK。

私は出力場所を "Export" という別フォルダにまとめたいので、下記のように出力。

# cvs で保存
write_csv(Res3, "Export/Chi_square_test.csv")            

カテゴリカルの要素が３つ以上の場合

上記例では２ｘ２の分割表で表されているけど、２つの群がカテゴリカルであれば、要素は３つ以上でも適応できる。

検定の方法は２ｘ２の時と同じ。

１．２ｘ４のマトリックスを作る。

# 2x4のマトリックスを作る
x <- matrix(data=c(0,5,100,250,150,200,15,0),byrow=TRUE,nc=4)

２．それぞれの項目に "dimnames" を指定する。

# dimnamesを指定する
dimnames(x) <- list(c("Yes","No"),c("Not","Low","Int","High"))
names(dimnames(x)) <- c("Dementia","ADNC")

３．カイ二乗検定

# Chi二乗検定
chisq.test(x)

結果は P < 2.22-16 となっており、「Dementia と ADNC の分布は独立ではない（つまり関係がある）」といえる。

適合度検定

「適合度」とは、想定している分布と、実際の測定結果の一致度の事。

例えばサイコロを60回振った場合、通常なら、1の目がでる回数は10回（1/6の確率）だけど、もしサイコロに仕掛けがしてあって、20回とか30回とかでてたら、「適合度が悪い」という事になる。

データ作り

上例のデータを作る。

# cohortの人種別人数
cohort <- c(205, 26, 25, 5, 15)                        

# アメリカ全体の人種比
USA <- c(0.615, 0.123, 0.176, 0.053, 0.033)

カイ二乗検定：適合度検定

chisq.test(x = 観測度数, p = 理論分布)

で求められる。

# 適合度検定
chisq.test(x = cohort, p = USA)        

結果は、

P = 3.4e-06 で、仮説は否定された。

とゆーことで、今回のコホートはアメリカ全体の人種の分布と異なる、という結果になった。

表として出力

chisq.test()

を res に代入し、それぞれのパラメーターを見てみる。

# 適合度検定
res &lt- chisq.test(x=cohort, p=USA)        

# それぞれのパラメーター
res$statistic
res$p.value
res$parameter
res$method
res$data.name
res$observed
res$expected
res$residuals

この結果を、

{flextable}

を使って表として出力してみる。

# データフレームにする
res2 &lt- as.data.frame(c(res$statistic, res$p.value, res$parameter, res$method ))    
colnames(res2) &lt- "chi-square goodness-of-fit test"

res3 &lt- res2 %>%
    as.data.frame(row.names = c("statistic",
                                "p.value",
                                "parameter",
                                "method"
                                )) %>%
    mutate(" " = rownames(res3)) %>%                            
    select(" ", everything())

# 表にする    
res3 %>%
regulartable() %>%
theme_booktabs() %>%
font(fontname = "serif", part = "all") %>%
autofit()

# docxで保存
save_as_docx(res3, path = "Export/chi_square_goodness_of_fit_test.docx", zoom = 3, expand = 10, webshot =
"webshot")

# csv で保存
write_csv(res3, "Export/chi_square_goodness_of_fit_test.csv")

References