四分表の検定についてシミュレーションを行う。四分表は2要因２水準のデータをtable形式で集計したものである。
総数Nのデータを生成し集計することで、1つの四分表が得られる。従って四分表のシミュレーションを行うには、
「データの生成と集計」というプロセスを無限回行う必要がある。
1回のプロセスの方法は次の通りである。

分割表の　縦 X=0,1 、横 Y=0,1　とする。

1 例4-5(p125)

library(dplyr)

## 
##  次のパッケージを付け加えます: 'dplyr'

##  以下のオブジェクトは 'package:stats' からマスクされています: 
## 
##      filter, lag

##  以下のオブジェクトは 'package:base' からマスクされています: 
## 
##      intersect, setdiff, setequal, union

data <- data.frame(X=c(0,0,1,1),Y=c(0,1,0,1))
data

X	Y
0	0
0	1
1	0
1	1

2 多項分布

p1=72/140;p2=80/140
s <- sample_n(tbl = data, size = 140,replace=T,weight=c(p1*p2,p1*(1-p2),p2*(1-p1),(1-p1)*(1-p2)))
table(s)

##    Y
## X    0  1
##   0 40 35
##   1 38 27

table(s)[1,1]

## [1] 40

for構文を用いてさらに２セット出力する。帰無仮説のもとでも出力は変動することが分かる。

p1=72/140;p2=80/140
for (i in 1:2){
  s <- sample_n(tbl = data, size = 140,replace=T,weight=c(p1*p2,p1*(1-p2),p2*(1-p1),(1-p1)*(1-p2)))
  print(table(s))
}

##    Y
## X    0  1
##   0 40 32
##   1 35 33
##    Y
## X    0  1
##   0 42 27
##   1 47 24

3 シミュレーション

3.1 データ生成

d <- rmultinom(n=1000, size = 140, prob = c(p1*p2,p1*(1-p2),p2*(1-p1),(1-p1)*(1-p2)))  # 乱数を発生
data <- as.data.frame(t(d))
colnames(data) <- c("a", "b", "c", "d")
head(data)

a	b	c	d
35	38	38	29
38	32	39	31
31	31	43	35
45	29	37	29
40	28	41	31
46	30	38	26

3.2 カイ二乗値を計算

X = c()
for (i in 1:1000) {
  a=data[i,1]
  b=data[i,2]
  c=data[i,3]
  d=data[i,4]
  kai=(a*d-b*c)^2*140/((a+b)*(c+d)*(a+c)*(b+d))
  X <-c(X,kai)
}
mean(X)

## [1] 0.9881003

平均はカイ二乗分布の適合度１に等しい。セルが4つあるのに、何故、自由度は１なのか？
4-1=3　ではないのかという疑問が生じると思う。
これは、公式 kai の導出（後述）による。
周辺分布の比率に母比率を使わずに、データごとの比率
　　　\(p_1=\dfrac{a+b}{N}~,~~p_2=\dfrac{a+c}{N}\)
を推定値として用いているからである。
（一般に母パラメータの代わりにr個の推定値を用いると自由度はr減ずる）

3.2.1 推定母数を用いないとき

X = c()
y <- c(c(80,60)*(72/140),c(80,60)*(68/140))

for (i in 1:1000) {
  s <- sum((data[i,]-y)^2/y)
  X <-c(X,s)
}
mean(X)

## [1] 2.978093

#hist(X)

平均は３になる。
これは、期待度数の計算に母数の推定値を用いていないためで、通常の適合度検定と同様、自由度は 4-1=3 となっている。

3.2.2 別の方法1

X = c()
Y = c()
for (i in 1:1000) {
  a1 <- rmultinom(n=1, size = 72, prob = c(80,60)/140)
  a2 <- rmultinom(n=1, size = 68, prob = c(80,60)/140)
  x <-c(a1,a2)
  y <- c(c(80,60)*(72/140),c(80,60)*68/140)
  s <- sum((abs(x-y))^2/y)
  X <-c(X,s)
  Y <-c(Y,x)
}
mean(X)

## [1] 1.980969

hist(X,prob=T)
curve(dchisq(x,2),add=T)

　　　\(p_1=\dfrac{72}{140}~,~~p_2=\dfrac{68}{140}\)
を推定値として用いている。実質1つの母パラメーターを推定値として用いているので、自由度は 4-1-1=2 である。
（一般に母パラメータの代わりにr個の推定値を用いると自由度はr減ずる）

3.2.3 別の方法2

周辺分布を固定して出力してみる。

X = c()
Y = c()
for (i in 1:1000) {
  a1 <- rbinom(n=1, size = 72, prob = c(80/140))
  a2 <- 80-a1
  x <-c(a1,72-a1,a2,68-a2)
  #print(x)
  y <- c(c(80,60)*(72/140),c(80,60)*68/140)
  s <- sum((abs(x-y))^2/y)
  X <-c(X,s)
  Y <-c(Y,x)
}
mean(X)

## [1] 2.221094

やはり自由度は２となるようである。

　　　\(\chi^2=\dfrac{(ad-bc)^2N}{(a+b)(c+d)(a+c)(b+d)}\)
で計算したカイ二乗値の自由度が１となるのは、
公式の計算過程で2つの比率を推定値として用いているからだと思う。（確証はないが）

X = c()
for (i in 1:1000) {
  a1 <- rmultinom(n=1, size = 72, prob = c(80,60)/140)
  a2 <- rmultinom(n=1, size = 68, prob = c(80,60)/140)
  data <-c(a1,a2)
  a=data[1]
  b=data[2]
  c=data[3]
  d=data[4]
  kai=(a*d-b*c)^2*140/((a+b)*(c+d)*(a+c)*(b+d))
  X <-c(X,kai)
}
mean(X)

## [1] 1.016991

hist(X,prob=T)
curve(dchisq(x,1),add=T)

4 2要因の検定

\(2\times2\)分割表の検定を行う。

4.1 研究

「2.3_研究」で行ったように、各セルの出現頻度を指定して抽出することもできる。
この場合“dplyr”というパッケージが必要になる。

library(dplyr)
data <- data.frame(X=c(0,0,0,1,1,1),Y=c(0,1,2,0,1,2))
data

X	Y
0	0
0	1
0	2
1	0
1	1
1	2

s <- sample_n(tbl = data, size = 1000,replace=T,weight=c(21,14,1,7,2,0))
table(s)

##    Y
## X     0   1   2
##   0 458 300  18
##   1 169  55   0

比率をweightで指定している。
点をずらして散布図をかくと

計数値データの検定