Julia で統計解析　第 7 章　検定と推定¶

In [1]:
Version()

最新バージョン 2024-10-29 14:07

この章に示す統計検定を行う関数の詳細は以下を参照のこと。
外部リンク：HypothesisTests package

この章に示す統計検定を行う前に，以下を行っておく。

In [2]:
using HypothesisTests

検定関数関数の呼び出し方¶

検定関数は以下のように定義されている（引数の型宣言を省略）。

NameOfTest(a, b, c, d=2, e="abc"; f=100, g=20.0)

; の左にある「名前=値」の形式をした引数はオプション引数と呼ばれ，省略可能である。省略された場合は指定された値がデフォルトとして使われる。; がなければ，「名前=値」の形式をした引数はすべてオプション引数である。

たとえば，上の例のようにオプション引き数が 2 個ある場合は 3 通りの呼び出し方がある。

  • 2 個の引数すべてを指定する。
    NameOfTest(a, b, c, 3, "xyz")
    d=3, e="xyz" として呼び出される

  • 1 個の引数だけを指定する。
    NameOfTest(a, b, c, 4)
    d=4 で，省略された e はデフォルトの "abc" として呼び出される

  • 全く指定しない。
    NameOfTest(a, b, c)
    省略された d, e はデフォルトの 2, "abc" として呼び出される

もし，引数の記述において，; があれば，その後にあるすべての「名前=値」の形式をした引数はキーワード引数と呼ばれ，省略可能である。省略された場合は指定された値がデフォルトとして使われる。デフォルト以外の値を指定する必要がある場合は「名前=値」の形式で指定しなければならない。

たとえば，上の例のようにキーワード引き数が f, g のような名前を持つなら，f=123, g=25.6 のように指定する。

なお，キーワード変数を指定して呼び出す場合には，通常は区切り文字 ; は , にしてもかまわない。

検定関数関数により得られる結果の利用法¶

例えば，検定関数が

NameOfTest(a, b, c, d=2, e="abc"; f=100, g=20.0)

のように呼ばれると，REPL では検定結果をまとめて，見やすい形式で表示してくれる。 REPL でない場合は

result = NameOfTest(a, b, c, d=2, e="abc"; f=100, g=20.0)

のように，関数の戻り値を変数に代入しておき，println(result) とすれば同じように表示される。

また，戻り値はいくつかあり，個々の検定結果の数値を選択して表示することもできる。たとえば，戻り値の要素が a，b の場合には println(result.a)，println(result.b) のようにすればよい。

なお，HypothesisTests に含まれる関数はデフォルトで両側検定の場合の $p$ 値と，信頼率が95% の場合の信頼区間の表示するが，片側検定や信頼率が 95% 以外の信頼区間を表示する場合には，pvalue メソッド（関数）と confint メソッド（関数）を呼ばなければならない。

関数の戻り値を result に代入しているとして，

  • pvalue() の使用法

pvalue(result) # デフォルトで両側検定
pvalue(result, tail=:both) # 両側検定
pvalue(result, tail=:left) # 片側検定
pvalue(result, tail=:right) # 片側検定

  • confint() の使用法

confint(result, tail=:both) # デフォルトで両側に棄却域を持つ 95％ 信頼区間
confint(result, level=0.99) # デフォルトで両側に棄却域を持つ 99％ 信頼区間
confint(result, tail=:both, level=0.99) # 両側に棄却域を持つ 99％ 信頼区間
confint(result, tail=:left, level=0.99) # 右側に棄却域を持つ 99% 信頼区間
confint(result, tail=:right, level=0.99) # 左側に棄却域を持つ 99% 信頼区間

検定と推定は表裏一体である。HypothesisTest に含まれる検定関数は $p$ 値だけではなく信頼区間を求めるものもある。

$p$ 値が有意水準 (1 - level) より小さいことと，信頼区間に検定統計量が含まれないことは同じこと「帰無仮説が棄却される」である。

分布の検定¶

引数に観察値の度数分布ベクトルが与えられた場合，観察値が特定の確率に従っているかどうかを検定する。

ここでは，ピアソンの $\chi^2$ 検定と，対数尤度比検定（$G^2$ 検定）を提示する。

一般的に，$G^2$ 検定のほうが検定力は高い（$p$ 値が小さい）。

観察度数が一様かどうかの検定¶

たとえば，16 回サイコロを振ったところ，1~6 のそれぞれの目が，4, 2, 1, 3, 4, 2 回ずつ出たとする。「本来サイコロはどの目が出る確率も同じはずなので，観察された目の出方はおかしいのではないか」という疑問に答えるための検定である。

ピアソンの $\chi^2$ 検定¶

関数定義
PowerDivergenceTest(x; lambda = 1.0)

第 1 引数に観察度数ベクトル，第 2 引数に lambda=1.0 を指定する。注意すべき点は lambda=1 ではないということである（MethodError が出るので気づく）。

■ 使用例

In [3]:
result = PowerDivergenceTest([4, 2, 1, 3, 4, 2], lambda=1.0) # ピアソンの$\chi^2$検定
Out[3]:
Pearson's Chi-square Test
-------------------------
Population details:
    parameter of interest:   Multinomial Probabilities
    value under h_0:         [0.166667, 0.166667, 0.166667, 0.166667, 0.166667, 0.166667]
    point estimate:          [0.25, 0.125, 0.0625, 0.1875, 0.25, 0.125]
    95% confidence interval: [(0.0625, 0.5238), (0.0, 0.3988), (0.0, 0.3363), (0.0, 0.4613), (0.0625, 0.5238), (0.0, 0.3988)]

Test summary:
    outcome with 95% confidence: fail to reject h_0
    one-sided p-value:           0.7385

Details:
    Sample size:        16
    statistic:          2.75
    degrees of freedom: 5
    residuals:          [0.816497, -0.408248, -1.02062, 0.204124, 0.816497, -0.408248]
    std. residuals:     [0.894427, -0.447214, -1.11803, 0.223607, 0.894427, -0.447214]

■ 個別に表示できる統計量

  • 概要
In [4]:
result = PowerDivergenceTest([4, 2, 1, 3, 4, 2], lambda=1.0);
p = pvalue(result)
println("Chisq. = $(result.stat),  df = $(result.df),  p value = $p")

Chisq. = 2.75,  df = 5,  p value = 0.7384611787603711

以下のような関数を定義しておくと，chisq_test(観察ベクトル) とするだけで結果が表示できる。

In [5]:
function chisq_test(x; lambda=1.0, theta0=ones(length(x))/length(x))
    result = PowerDivergenceTest(x; lambda, theta0)
    p = pvalue(result)
    println("χ-sq. = $(result.stat),  df = $(result.df),  p value = $p")
end;
In [6]:
chisq_test([4, 2, 1, 3, 4, 2])

χ-sq. = 2.75,  df = 5,  p value = 0.7384611787603711
In [7]:
chisq_test([21, 12, 5, 4])

χ-sq. = 17.61904761904762,  df = 3,  p value = 0.0005270259708713794
  • $p$ 値　この検定では tail=:both，tail=:left は使用してはいけない
    上に示した結果の一覧表示で "one-sided p-value: 0.7385" とあるのは，この検定が "one-sided" はないので誤った記述である。棄却域は $\chi^2$ 分布の右側にしかないが，本来この検定は意味的に両側検定である（両側検定，片側検定という区分分け自体が不要である）。
In [8]:
pvalue(result)
Out[8]:
0.7384611787603711
In [9]:
pvalue(result, tail=:right)
Out[9]:
0.7384611787603711
  • 信頼区間　信頼区間は度数分布の確率の信頼区間である。
    この検定では，tail=:left，tail=:right は使用してはいけない（pvalue() の場合と違うので注意が必要である）。
In [10]:
confint(result)
Out[10]:
6-element Vector{Tuple{Float64, Float64}}:
 (0.0625, 0.5237538769784804)
 (0.0, 0.3987538769784804)
 (0.0, 0.3362538769784804)
 (0.0, 0.4612538769784804)
 (0.0625, 0.5237538769784804)
 (0.0, 0.3987538769784804)
In [11]:
confint(result, tail=:both)
Out[11]:
6-element Vector{Tuple{Float64, Float64}}:
 (0.0625, 0.5237538769784804)
 (0.0, 0.3987538769784804)
 (0.0, 0.3362538769784804)
 (0.0, 0.4612538769784804)
 (0.0625, 0.5237538769784804)
 (0.0, 0.3987538769784804)
  • lambda　$\lambda$
In [12]:
result.lambda
Out[12]:
1.0
  • n　サンプルサイズ
In [13]:
result.n
Out[13]:
16
  • observed　観察度数分布
In [14]:
result.observed
Out[14]:
6×1 Matrix{Int64}:
 4
 2
 1
 3
 4
 2
  • thetahat　観察度数の確率
    result.observed / result.n
In [15]:
result.thetahat
Out[15]:
6-element Vector{Float64}:
 0.25
 0.125
 0.0625
 0.1875
 0.25
 0.125
  • theta0　帰無仮説の下での観察度数の確率（等確率）
    ones(length(result.observed)) / length(result.observed)
In [16]:
round.(result.theta0, digits=5)
Out[16]:
6-element Vector{Float64}:
 0.16667
 0.16667
 0.16667
 0.16667
 0.16667
 0.16667
  • expected　帰無仮説の下での観察度数の期待値（同じ値）
    result.n * result.theta0
In [17]:
round.(result.expected, digits=5)
Out[17]:
6×1 Matrix{Float64}:
 2.66667
 2.66667
 2.66667
 2.66667
 2.66667
 2.66667
  • stat　検定統計量（$\chi^2$ 分布にしたがう）
    sum((result.observed - result.expected).^2 ./ result.expected)
In [18]:
result.stat
Out[18]:
2.75
  • df　$\chi^2$ 分布の自由度
    length(result.observed) - 1
In [19]:
result.df
Out[19]:
5
  • residuals　残差
    (result.observed - result.expected) ./ sqrt.(result.expected)
In [20]:
round.(result.residuals, digits=5)
Out[20]:
6×1 Matrix{Float64}:
  0.8165
 -0.40825
 -1.02062
  0.20412
  0.8165
 -0.40825
  • stdresiduals　標準化残差
    V = result.n .* result.theta0 .* (1 .- result.theta0)
    (result.observed - result.expected) ./ sqrt.(V)
In [21]:
round.(result.stdresiduals, digits=5)
Out[21]:
6×1 Matrix{Float64}:
  0.89443
 -0.44721
 -1.11803
  0.22361
  0.89443
 -0.44721

対数尤度比検定（$G^2$ 検定）¶

関数定義
PowerDivergenceTest(x; lambda = 0.0)

対数尤度比に基づく検定は，「3.1.1.　ピアソンの $\chi^2$ 検定」と同じ PowerDivergenceTest() で lambda=0.0 を指定することで行うことができる。

■ 使用例

In [22]:
result = PowerDivergenceTest([4, 2, 1, 3, 4, 2], lambda=0.0)
Out[22]:
Multinomial Likelihood Ratio Test
---------------------------------
Population details:
    parameter of interest:   Multinomial Probabilities
    value under h_0:         [0.166667, 0.166667, 0.166667, 0.166667, 0.166667, 0.166667]
    point estimate:          [0.25, 0.125, 0.0625, 0.1875, 0.25, 0.125]
    95% confidence interval: [(0.0625, 0.5238), (0.0, 0.3988), (0.0, 0.3363), (0.0, 0.4613), (0.0625, 0.5238), (0.0, 0.3988)]

Test summary:
    outcome with 95% confidence: fail to reject h_0
    one-sided p-value:           0.7106

Details:
    Sample size:        16
    statistic:          2.931024858031232
    degrees of freedom: 5
    residuals:          [0.816497, -0.408248, -1.02062, 0.204124, 0.816497, -0.408248]
    std. residuals:     [0.894427, -0.447214, -1.11803, 0.223607, 0.894427, -0.447214]

■ 個別に表示できる統計量

  • 概要

前項の chisq_test() と同様に以下の関数を定義して簡潔な結果を表示する。

In [23]:
function G2_test(x; theta0=ones(length(x))/length(x), lambda=0.0)
    result = PowerDivergenceTest(x; theta0, lambda);
    p = pvalue(result)
    println("G-sq. = $(result.stat),  df = $(result.df),  p value = $p")
end

G2_test([4, 2, 1, 3, 4, 2])

G-sq. = 2.931024858031232,  df = 5,  p value = 0.710619024390339
In [24]:
G2_test([21, 12, 5, 4])

G-sq. = 17.176914390863786,  df = 3,  p value = 0.0006499306837478541
  • $p$ 値　この検定では tail=:both，tail=:left は使用してはいけない
    上に示した結果の一覧表示で "one-sided p-value: 0.7106" とあるのは，この検定が "one-sided" はないので誤った記述である。棄却域は $\chi^2$ 分布の右側にしかないが，本来この検定は意味的に両側検定である（両側検定，片側検定という区分分け自体が不要である）。
In [25]:
pvalue(result)
Out[25]:
0.710619024390339
In [26]:
pvalue(result, tail=:right)
Out[26]:
0.710619024390339
  • 信頼区間　信頼区間は度数分布の確率の信頼区間である。
    この検定では，tail=:left，tail=:right は使用してはいけない（pvalue() の場合と違うので注意が必要である）。
In [27]:
confint(result)
Out[27]:
6-element Vector{Tuple{Float64, Float64}}:
 (0.0625, 0.5237538769784804)
 (0.0, 0.3987538769784804)
 (0.0, 0.3362538769784804)
 (0.0, 0.4612538769784804)
 (0.0625, 0.5237538769784804)
 (0.0, 0.3987538769784804)
In [28]:
confint(result, tail=:both)
Out[28]:
6-element Vector{Tuple{Float64, Float64}}:
 (0.0625, 0.5237538769784804)
 (0.0, 0.3987538769784804)
 (0.0, 0.3362538769784804)
 (0.0, 0.4612538769784804)
 (0.0625, 0.5237538769784804)
 (0.0, 0.3987538769784804)
  • lambda　$\lambda$
In [29]:
result.lambda
Out[29]:
0.0
  • n　サンプルサイズ
In [30]:
result.n
Out[30]:
16
  • observed　観察度数分布
In [31]:
result.observed
Out[31]:
6×1 Matrix{Int64}:
 4
 2
 1
 3
 4
 2
  • thetahat　観察度数の確率
    result.observed / result.n
In [32]:
result.thetahat
Out[32]:
6-element Vector{Float64}:
 0.25
 0.125
 0.0625
 0.1875
 0.25
 0.125
  • theta0　帰無仮説の下での観察度数の確率（等確率）
    ones(length(result.observed)) / length(result.observed)
In [33]:
round.(result.theta0, digits=5)
Out[33]:
6-element Vector{Float64}:
 0.16667
 0.16667
 0.16667
 0.16667
 0.16667
 0.16667
  • expected　帰無仮説の下での観察度数の期待値（同じ値）
    result.n * result.theta0
In [34]:
round.(result.expected, digits=5)
Out[34]:
6×1 Matrix{Float64}:
 2.66667
 2.66667
 2.66667
 2.66667
 2.66667
 2.66667
  • stat　検定統計量（$\chi^2$ 分布にしたがう）
    sum((result.observed - result.expected).^2 ./ result.expected)
In [35]:
result.stat
Out[35]:
2.931024858031232
  • df　$\chi^2$ 分布の自由度
    length(result.observed) - 1
In [36]:
result.df
Out[36]:
5
  • residuals　残差
    (result.observed - result.expected) ./ sqrt.(result.expected)
In [37]:
round.(result.residuals, digits=5)
Out[37]:
6×1 Matrix{Float64}:
  0.8165
 -0.40825
 -1.02062
  0.20412
  0.8165
 -0.40825
  • stdresiduals　標準化残差
    V = result.n .* result.theta0 .* (1 .- result.theta0)
    (result.observed - result.expected) ./ sqrt.(V)
In [38]:
round.(result.stdresiduals, digits=5)
Out[38]:
6×1 Matrix{Float64}:
  0.89443
 -0.44721
 -1.11803
  0.22361
  0.89443
 -0.44721

観察度数が理論比に從うかどうかの検定¶

観察値が理論比に從うかどうかの検定は PowerDivergenceTest() で，観察値の度数分布ベクトルと，理論比のベクトル theta0 を与える。

理論比は，合計して 1 にならなければならないが，分数で与えることもできる。

ピアソンの $\chi^2$ 検定¶

関数定義
PowerDivergenceTest(x; theta0 = ones(length(x))/length(x), lambda = 1.0)

PowerDivergenceTest() で，lambda=1.0 を指定する。

■ 使用例

In [39]:
result = PowerDivergenceTest([21, 12, 5, 4], theta0=[9/16, 3/16, 3/16, 1/16], lambda=1.0)
Out[39]:
Pearson's Chi-square Test
-------------------------
Population details:
    parameter of interest:   Multinomial Probabilities
    value under h_0:         [0.5625, 0.1875, 0.1875, 0.0625]
    point estimate:          [0.5, 0.285714, 0.119048, 0.0952381]
    95% confidence interval: [(0.3571, 0.6576), (0.1429, 0.4433), (0.0, 0.2766), (0.0, 0.2528)]

Test summary:
    outcome with 95% confidence: fail to reject h_0
    one-sided p-value:           0.2384

Details:
    Sample size:        42
    statistic:          4.22222222222222
    degrees of freedom: 3
    residuals:          [-0.540062, 1.46994, -1.0245, 0.848668]
    std. residuals:     [-0.816497, 1.63075, -1.13658, 0.876501]

■ 個別に表示できる統計量

  • 概要

前節で定義した chisq_test() を使うと以下のようになる。

In [40]:
chisq_test([21, 12, 5, 4], theta0=[9/16, 3/16, 3/16, 1/16]) # ピアソンのχ二乗検定

χ-sq. = 4.22222222222222,  df = 3,  p value = 0.23844641563476865
  • $p$ 値　この検定では tail=:both，tail=:left は使用してはいけない
    上に示した結果の一覧表示で "one-sided p-value: 0.2384" とあるのは，この検定が "one-sided" はないので誤った記述である。棄却域は $\chi^2$ 分布の右側にしかないが，本来この検定は意味的に両側検定である（両側検定，片側検定という区分分け自体が不要である）。
In [41]:
pvalue(result)
Out[41]:
0.23844641563476865
In [42]:
pvalue(result, tail=:right)
Out[42]:
0.23844641563476865
  • 信頼区間　信頼区間は度数分布の確率の信頼区間である。
    この検定では，tail=:left，tail=:right は使用してはいけない（pvalue() の場合と違うので注意が必要である）。
In [43]:
confint(result)
Out[43]:
4-element Vector{Tuple{Float64, Float64}}:
 (0.35714285714285715, 0.6575992779891421)
 (0.14285714285714285, 0.4433135637034278)
 (0.0, 0.2766468970367611)
 (0.0, 0.2528373732272373)
In [44]:
confint(result, tail=:both)
Out[44]:
4-element Vector{Tuple{Float64, Float64}}:
 (0.35714285714285715, 0.6575992779891421)
 (0.14285714285714285, 0.4433135637034278)
 (0.0, 0.2766468970367611)
 (0.0, 0.2528373732272373)
  • lambda　$\lambda$
In [45]:
result.lambda
Out[45]:
1.0
  • n　サンプルサイズ
In [46]:
result.n
Out[46]:
42
  • observed　観察度数分布
In [47]:
result.observed
Out[47]:
4×1 Matrix{Int64}:
 21
 12
  5
  4
  • thetahat　観察度数の確率
    result.observed / result.n
In [48]:
result.thetahat
Out[48]:
4-element Vector{Float64}:
 0.5
 0.2857142857142857
 0.11904761904761904
 0.09523809523809523
  • theta0　帰無仮説の下での観察度数の確率（等確率）
    ones(length(result.observed)) / length(result.observed)
In [49]:
round.(result.theta0, digits=5)
Out[49]:
4-element Vector{Float64}:
 0.5625
 0.1875
 0.1875
 0.0625
  • expected　帰無仮説の下での観察度数の期待値（同じ値）
    result.n * result.theta0
In [50]:
round.(result.expected, digits=5)
Out[50]:
4×1 Matrix{Float64}:
 23.625
  7.875
  7.875
  2.625
  • stat　検定統計量（$\chi^2$ 分布にしたがう）
    sum((result.observed - result.expected).^2 ./ result.expected)
In [51]:
result.stat
Out[51]:
4.22222222222222
  • df　$\chi^2$ 分布の自由度
    length(result.observed) - 1
In [52]:
result.df
Out[52]:
3
  • residuals　残差
    (result.observed - result.expected) ./ sqrt.(result.expected)
In [53]:
round.(result.residuals, digits=5)
Out[53]:
4×1 Matrix{Float64}:
 -0.54006
  1.46994
 -1.0245
  0.84867
  • stdresiduals　標準化残差
    V = result.n .* result.theta0 .* (1 .- result.theta0)
    (result.observed - result.expected) ./ sqrt.(V)
In [54]:
round.(result.stdresiduals, digits=5)
Out[54]:
4×1 Matrix{Float64}:
 -0.8165
  1.63075
 -1.13658
  0.8765

対数尤度比検定（$G^2$ 検定）¶

関数定義
PowerDivergenceTest(x; theta0 = ones(length(x))/length(x), lambda = 0.0)

PowerDivergenceTest() で，lambda=1.0 を指定する。

■ 使用例

In [55]:
result = PowerDivergenceTest([21, 12, 5, 4], theta0=[9/16, 3/16, 3/16, 1/16], lambda=0.0)
Out[55]:
Multinomial Likelihood Ratio Test
---------------------------------
Population details:
    parameter of interest:   Multinomial Probabilities
    value under h_0:         [0.5625, 0.1875, 0.1875, 0.0625]
    point estimate:          [0.5, 0.285714, 0.119048, 0.0952381]
    95% confidence interval: [(0.3571, 0.6576), (0.1429, 0.4433), (0.0, 0.2766), (0.0, 0.2528)]

Test summary:
    outcome with 95% confidence: fail to reject h_0
    one-sided p-value:           0.2626

Details:
    Sample size:        42
    statistic:          3.9893906620976454
    degrees of freedom: 3
    residuals:          [-0.540062, 1.46994, -1.0245, 0.848668]
    std. residuals:     [-0.816497, 1.63075, -1.13658, 0.876501]

■ 個別に表示できる統計量

  • 概要

前節で定義した G2_test() を使うと以下のようになる。

In [56]:
G2_test([21, 12, 5, 4], theta0=[9/16, 3/16, 3/16, 1/16]) # 対数尤度比検定（G2検定）

G-sq. = 3.9893906620976454,  df = 3,  p value = 0.26261202803389555
  • $p$ 値　この検定では tail=:both，tail=:left は使用してはいけない
    上に示した結果の一覧表示で "one-sided p-value: 0.2626" とあるのは，この検定が "one-sided" はないので誤った記述である。棄却域は $\chi^2$ 分布の右側にしかないが，本来この検定は意味的に両側検定である（両側検定，片側検定という区分分け自体が不要である）。
In [57]:
pvalue(result)
Out[57]:
0.26261202803389555
In [58]:
pvalue(result, tail=:right)
Out[58]:
0.26261202803389555
  • 信頼区間　信頼区間は度数分布の確率の信頼区間である。
    この検定では，tail=:left，tail=:right は使用してはいけない（pvalue() の場合と違うので注意が必要である）。
In [59]:
confint(result)
Out[59]:
4-element Vector{Tuple{Float64, Float64}}:
 (0.35714285714285715, 0.6575992779891421)
 (0.14285714285714285, 0.4433135637034278)
 (0.0, 0.2766468970367611)
 (0.0, 0.2528373732272373)
In [60]:
confint(result, tail=:both)
Out[60]:
4-element Vector{Tuple{Float64, Float64}}:
 (0.35714285714285715, 0.6575992779891421)
 (0.14285714285714285, 0.4433135637034278)
 (0.0, 0.2766468970367611)
 (0.0, 0.2528373732272373)
  • lambda　$\lambda$
In [61]:
result.lambda
Out[61]:
0.0
  • n　サンプルサイズ
In [62]:
result.n
Out[62]:
42
  • observed　観察度数分布
In [63]:
result.observed
Out[63]:
4×1 Matrix{Int64}:
 21
 12
  5
  4
  • thetahat　観察度数の確率
    result.observed / result.n
In [64]:
result.thetahat
Out[64]:
4-element Vector{Float64}:
 0.5
 0.2857142857142857
 0.11904761904761904
 0.09523809523809523
  • theta0　帰無仮説の下での観察度数の確率（等確率）
    ones(length(result.observed)) / length(result.observed)
In [65]:
round.(result.theta0, digits=5)
Out[65]:
4-element Vector{Float64}:
 0.5625
 0.1875
 0.1875
 0.0625
  • expected　帰無仮説の下での観察度数の期待値（同じ値）
    result.n * result.theta0
In [66]:
round.(result.expected, digits=5)
Out[66]:
4×1 Matrix{Float64}:
 23.625
  7.875
  7.875
  2.625
  • stat　検定統計量（$\chi^2$ 分布にしたがう）
    sum((result.observed - result.expected).^2 ./ result.expected)
In [67]:
result.stat
Out[67]:
3.9893906620976454
  • df　$\chi^2$ 分布の自由度
    length(result.observed) - 1
In [68]:
result.df
Out[68]:
3
  • residuals　残差
    (result.observed - result.expected) ./ sqrt.(result.expected)
In [69]:
round.(result.residuals, digits=5)
Out[69]:
4×1 Matrix{Float64}:
 -0.54006
  1.46994
 -1.0245
  0.84867
  • stdresiduals　標準化残差
    V = result.n .* result.theta0 .* (1 .- result.theta0)
    (result.observed - result.expected) ./ sqrt.(V)
In [70]:
round.(result.stdresiduals, digits=5)
Out[70]:
4×1 Matrix{Float64}:
 -0.8165
  1.63075
 -1.13658
  0.8765

独立性の検定¶

独立性の検定は，2 つのカテゴリー変数が独立であるか（関連があるか）を検定する。

ここでは，前節と同じピアソンの $\chi^2$ 検定 と，対数尤度比検定（$G^2$ 検定）を提示する。

利用する関数は，前節と同じ PowerDivergenceTest() である。引数として二重集計表を与える場合と，2 つのベクトルを与える場合に対応している。

しかし，2 つのベクトルを与える場合には完全対応していないので，二重集計表を与えて使用する。

ピアソンの $\chi^2$ 検定¶

関数定義
PowerDivergenceTest(x; lambda=1.0)

■ 使用例　2 × 2 分割表の場合

In [71]:
x = [12 35
     43 56]
Out[71]:
2×2 Matrix{Int64}:
 12  35
 43  56
In [72]:
result = PowerDivergenceTest(x, lambda=1.0)
Out[72]:
Pearson's Chi-square Test
-------------------------
Population details:
    parameter of interest:   Multinomial Probabilities
    value under h_0:         [0.12127, 0.255442, 0.200647, 0.42264]
    point estimate:          [0.0821918, 0.294521, 0.239726, 0.383562]
    95% confidence interval: [(0.0, 0.1716), (0.2123, 0.384), (0.1575, 0.3292), (0.3014, 0.473)]

Test summary:
    outcome with 95% confidence: reject h_0
    one-sided p-value:           0.0370

Details:
    Sample size:        146
    statistic:          4.350164943588549
    degrees of freedom: 1
    residuals:          [-1.35593, 0.934264, 1.05414, -0.726324]
    std. residuals:     [-2.0857, 2.0857, 2.0857, -2.0857]

■ 個別に表示できる統計量

  • 概要

前節で定義した chisq_test() もそのまま使える。

In [73]:
chisq_test(x)

χ-sq. = 4.350164943588549,  df = 1,  p value = 0.037005362019413034
  • $p$ 値　この検定では tail=:both，tail=:left は使用してはいけない
    上に示した結果の一覧表示で "one-sided p-value: 0.0370" とあるのは，この検定が "one-sided" はないので誤った記述である。棄却域は $\chi^2$ 分布の右側にしかないが，本来この検定は意味的に両側検定である（両側検定，片側検定という区分分け自体が不要である）。
In [74]:
pvalue(result)
Out[74]:
0.037005362019413034
In [75]:
pvalue(result, tail=:right)
Out[75]:
0.037005362019413034
  • イエーツの連続性の補正による $p$ 値

2×2 分割表の場合，PowerDivergenceTest() で得られる検定統計量はイエーツの連続性の補正をしないものである。

イエーツの連続性の補正をする場合には，以下の関数を定義しておくとよい。

In [76]:
using Distributions
function yates(x)
    result = PowerDivergenceTest(x, lambda=1.0)
    n = sum(x)
    a, c, b, d = x
    stat = n*max(0, abs(a*d - b*c) - 0.5n)^2 / ((a+b)*(c+d)*(a+c)*(b+d))
    p = ccdf(Chisq(1), stat)
    println("corrected χ-sq. = $stat,  df = 1,  p value = $p")
end;

yates(x)

corrected χ-sq. = 3.621119907386832,  df = 1,  p value = 0.05705045741699577
  • 信頼区間　信頼区間は分割表の各セルの確率の信頼区間である。
    この検定では，tail=:left，tail=:right は使用してはいけない（pvalue() の場合と違うので注意が必要である）。
In [77]:
confint(result)
Out[77]:
4-element Vector{Tuple{Float64, Float64}}:
 (0.0, 0.17163425153174042)
 (0.2123287671232877, 0.38396301865502813)
 (0.15753424657534246, 0.3291684981070829)
 (0.3013698630136986, 0.47300411454543906)
In [78]:
confint(result, tail=:both)
Out[78]:
4-element Vector{Tuple{Float64, Float64}}:
 (0.0, 0.17163425153174042)
 (0.2123287671232877, 0.38396301865502813)
 (0.15753424657534246, 0.3291684981070829)
 (0.3013698630136986, 0.47300411454543906)
  • lambda　$\lambda$
In [79]:
result.lambda
Out[79]:
1.0
  • n　サンプルサイズ
In [80]:
result.n
Out[80]:
146
  • observed　観察度数分布
In [81]:
result.observed
Out[81]:
2×2 Matrix{Int64}:
 12  35
 43  56
  • thetahat　観察度数の確率
    result.observed / result.n
In [82]:
result.thetahat
Out[82]:
4-element Vector{Float64}:
 0.0821917808219178
 0.2945205479452055
 0.23972602739726026
 0.3835616438356164
  • theta0　帰無仮説の下での観察度数の確率（等確率）
    ones(length(result.observed)) / length(result.observed)
In [83]:
round.(result.theta0, digits=5)
Out[83]:
4-element Vector{Float64}:
 0.12127
 0.25544
 0.20065
 0.42264
  • expected　帰無仮説の下での観察度数の期待値（同じ値）
    result.n * result.theta0
In [84]:
round.(result.expected, digits=5)
Out[84]:
2×2 Matrix{Float64}:
 17.7055  29.2945
 37.2945  61.7055
  • stat　検定統計量（$\chi^2$ 分布にしたがう）
    sum((result.observed - result.expected).^2 ./ result.expected)
In [85]:
result.stat
Out[85]:
4.350164943588549
  • df　$\chi^2$ 分布の自由度
    length(result.observed) - 1
In [86]:
result.df
Out[86]:
1
  • residuals　残差
    (result.observed - result.expected) ./ sqrt.(result.expected)
In [87]:
Base.print_matrix(stdout, round.(result.residuals, digits=7), "", "   ", "\n")

-1.3559332    1.0541416
 0.934264    -0.7263238
  • stdresiduals　標準化残差
    V = result.n .* result.theta0 .* (1 .- result.theta0)
    (result.observed - result.expected) ./ sqrt.(V)
In [88]:
Base.print_matrix(stdout, round.(result.stdresiduals, digits=7), "", "   ", "\n")

-2.0857049    2.0857049
 2.0857049   -2.0857049
  • 残差分析

残差分析の結果は PowerDivergenceTest(x, lambda=1.0) が返す std.residuals を使用する。

対応する $p$ 値を表示するために，以下の関数を用意しておくとよい。

In [89]:
using Distributions
function residual_analysis(x)
    result = PowerDivergenceTest(x, lambda=1.0);
    nrows, ncols = size(x)
    p = 2ccdf.(Normal(), abs.(result.stdresiduals))
    Base.print_matrix(stdout, round.(p, digits=7), "", "   ", "\n")
end

residual_analysis(x)

0.0370054   0.0370054
0.0370054   0.0370054

■ 使用例2　文字列ベクトルを与えて検定する

In [90]:
answer = ["Yes", "Yes", "No", "No", "No", "Yes", "No", "Yes", "Yes", "No",
          "No", "Yes", "Yes", "No", "Yes", "Yes", "No", "No", "No", "No"]
status = ["Hi", "Lo", "Med", "Med", "Lo", "Hi", "Lo", "Hi", "Lo", "Lo",
          "Med", "Med", "Lo", "Med", "Med", "Lo", "Lo", "Hi", "Lo", "Med"]

using FreqTables
y = freqtable(answer, status)
Out[90]:
2×3 Named Matrix{Int64}
Dim1 ╲ Dim2 │  Hi   Lo  Med
────────────┼──────────────
No          │   1    5    5
Yes         │   3    4    2
  • 概要
In [91]:
chisq_test(y)

χ-sq. = 2.2190155523488855,  df = 2,  p value = 0.32972121777778757
  • 残差分析の $p$ 値
In [92]:
residual_analysis(y)

0.1775299   0.9639693   0.2785029
0.1775299   0.9639693   0.2785029

対数尤度比検定（$G^2$ 検定）¶

関数定義
PowerDivergenceTest(x, lambda=0.0)

一般的には，$\chi^2$ 検定より $G^2$ 検定のほうが検定力は高い（$p$ 値が小さい）。

■ 使用例

In [93]:
x = [10 3; 4 12]
Out[93]:
2×2 Matrix{Int64}:
 10   3
  4  12
In [94]:
result = PowerDivergenceTest(x, lambda=0.0)
Out[94]:
Multinomial Likelihood Ratio Test
---------------------------------
Population details:
    parameter of interest:   Multinomial Probabilities
    value under h_0:         [0.216409, 0.26635, 0.231867, 0.285375]
    point estimate:          [0.344828, 0.137931, 0.103448, 0.413793]
    95% confidence interval: [(0.1724, 0.537), (0.0, 0.3301), (0.0, 0.2957), (0.2414, 0.606)]

Test summary:
    outcome with 95% confidence: reject h_0
    one-sided p-value:           0.0044

Details:
    Sample size:        29
    statistic:          8.1280145470097
    degrees of freedom: 1
    residuals:          [1.48658, -1.33999, -1.43618, 1.29455]
    std. residuals:     [2.7828, -2.7828, -2.7828, 2.7828]

■ 個別に表示できる統計量

  • 概要

前節で定義した G2_test() もそのまま使える。

In [95]:
G2_test(x) # 対数尤度比に基づく検定（G2検定）

G-sq. = 8.1280145470097,  df = 1,  p value = 0.00435864516471826
  • $p$ 値　この検定では tail=:both，tail=:left は使用してはいけない
    上に示した結果の一覧表示で "one-sided p-value: 0.0044" とあるのは，この検定が "one-sided" はないので誤った記述である。棄却域は $\chi^2$ 分布の右側にしかないが，本来この検定は意味的に両側検定である（両側検定，片側検定という区分分け自体が不要である）。
In [96]:
pvalue(result)
Out[96]:
0.00435864516471826
In [97]:
pvalue(result, tail=:right)
Out[97]:
0.00435864516471826
  • 信頼区間　信頼区間は分割表の各セルの確率の信頼区間である。
    この検定では，tail=:left，tail=:right は使用してはいけない（pvalue() の場合と違うので注意が必要である）。
In [98]:
confint(result)
Out[98]:
4-element Vector{Tuple{Float64, Float64}}:
 (0.1724137931034483, 0.5370409523965961)
 (0.0, 0.3301444006724581)
 (0.0, 0.2956616420517684)
 (0.24137931034482757, 0.6060064696379754)
In [99]:
confint(result, tail=:both)
Out[99]:
4-element Vector{Tuple{Float64, Float64}}:
 (0.1724137931034483, 0.5370409523965961)
 (0.0, 0.3301444006724581)
 (0.0, 0.2956616420517684)
 (0.24137931034482757, 0.6060064696379754)
  • lambda　$\lambda$
In [100]:
result.lambda
Out[100]:
0.0
  • n　サンプルサイズ
In [101]:
result.n
Out[101]:
29
  • observed　観察度数分布
In [102]:
result.observed
Out[102]:
2×2 Matrix{Int64}:
 10   3
  4  12
  • thetahat　観察度数の確率
    result.observed / result.n
In [103]:
result.thetahat
Out[103]:
4-element Vector{Float64}:
 0.3448275862068966
 0.13793103448275862
 0.10344827586206896
 0.41379310344827586
  • theta0　帰無仮説の下での観察度数の確率（等確率）
    ones(length(result.observed)) / length(result.observed)
In [104]:
round.(result.theta0, digits=5)
Out[104]:
4-element Vector{Float64}:
 0.21641
 0.26635
 0.23187
 0.28537
  • expected　帰無仮説の下での観察度数の期待値（同じ値）
    result.n * result.theta0
In [105]:
round.(result.expected, digits=5)
Out[105]:
2×2 Matrix{Float64}:
 6.27586  6.72414
 7.72414  8.27586
  • stat　検定統計量（$\chi^2$ 分布にしたがう）
    sum((result.observed - result.expected).^2 ./ result.expected)
In [106]:
result.stat
Out[106]:
8.1280145470097
  • df　$\chi^2$ 分布の自由度
    length(result.observed) - 1
In [107]:
result.df
Out[107]:
1
  • residuals　残差。ピアソンの $\chi^2$ 検定と同じである
    (result.observed - result.expected) ./ sqrt.(result.expected)
In [108]:
Base.print_matrix(stdout, round.(result.residuals, digits=7), "", "   ", "\n")

 1.4865827   -1.4361753
-1.3399875    1.2945509
  • stdresiduals　標準化残差。ピアソンの $\chi^2$ 検定と同じである

    V = result.n .* result.theta0 .* (1 .- result.theta0)
    (result.observed - result.expected) ./ sqrt.(V)
In [109]:
Base.print_matrix(stdout, round.(result.stdresiduals, digits=7), "", "   ", "\n")

 2.7827964   -2.7827964
-2.7827964    2.7827964

■ 分割表のマス目に 0 がある場合の $G^2$ 検定

PowerDivergenceTest() は　lambda=0.0 を指定することで対数尤度比に基づく検定（$G^2$ 検定）ができるが，分割表のマス目のうちに 0 であるものが 1 つでもある場合には正しい結果が得られない。

In [110]:
z = [4 5 2 0
     0 7 6 1
     1 0 3 1]
result = PowerDivergenceTest(z, lambda=0.0);
p = pvalue(result)
println("Chisq. = $(result.stat),  df = $(result.df),  p value = $p")

Chisq. = NaN,  df = 6,  p value = NaN

実は lambda を 0.0 に近づけていくと $G^2$ 検定になるということなので，たとえば lambda=0.0000001 では以下のようになる。

In [111]:
result = PowerDivergenceTest(z, lambda=0.000000001);
p = pvalue(result)
println("Chisq. = $(result.stat),  df = $(result.df),  p value = $p")

Chisq. = 15.364590917511304,  df = 6,  p value = 0.017602889013051494

次に示す trueG2_test() は前節で定義した関数と同じように $G^2$ 検定を実施するが，分割表のマス目のうちに 0 があるような場合でも，正確に lambda=0.0 にした場合の，正しい答えを返す。

In [112]:
using Distributions
function trueG2_test(x; correct=false)
    # セルに 0 があっても正しい答えを出す
    # correct=true で，連続性の補正も行うことができる
    ln(n) = sum(n .== 0 ? 0 : n .* log(n) for n in vcat(n...))
    nrows, ncols = size(x)
    n = sum(x) # 全サンプルサイズ
    n1 = sum(x, dims=2) # 行和
    n2 = sum(x, dims=1) # 列和
    G2 = 2*(ln(x) - ln(n1) - ln(n2) + ln(n)) # G 統計量
    correct && (G2 /= 1 + (n * sum(1 ./ n1) - 1) * (n * sum(1 ./ n2) - 1) / (6n * nrows * ncols)) # 連続性の補正
    df = (nrows - 1) * (ncols - 1) # G の自由度
    p = ccdf(Chisq(df), G2)
    name = correct ? "corrected G-sq." : "G-sq."
    println("$name = $G2,  df = $df,  p value = $p")
end

trueG2_test(z)

G-sq. = 15.364591286599591,  df = 6,  p value = 0.01760288650305146

また，ピアソンの $\chi^2$ 検定では 2×2 分割表の場合にしか連続性の補正は行えないが，$G^2$ 検定は分割表の大きさに関わらず，連続性の補正を行うことができる。

上に示した関数 trueG2_test() で correct=true を指定すればよい。

In [113]:
trueG2_test(z)

G-sq. = 15.364591286599591,  df = 6,  p value = 0.01760288650305146
In [114]:
trueG2_test(z, correct=true) # 連続性の補正

corrected G-sq. = 13.776490649307341,  df = 6,  p value = 0.03223501559558125

パワーダイバージェンス検定¶

関数定義
PowerDivergenceTest(x; lambda = 1.0, theta0 = ones(length(x))/length(x))

パワーダイバージェンス検定（Power divergence test）は前節の ピアソンの $\chi^2$，対数尤度比検定（$G^2$ 検定）などを包含する，一般化した検定である。

$$ \frac{2}{\lambda(\lambda+1)} \sum_{i=1}^I \sum_{j=1}^J n_{ij} \left [\left ( \frac{n_{ij}}{\hat{n}_{ij}}\right )^\lambda - 1\right ] $$

ここで，lambda($\lambda$) と適用される検定統計量の対応は以下のようになる。

lambda 適用される検定統計量
λ = 1 ピアソンの $\chi^2$ 検定統計量
λ → 0 対数尤度比検定統計量 $G^2$ に収束
λ → −1 最小判別情報量統計量に収束(Gokhale and Kullback, 1978)
λ = −2 修正された $\chi^2$ 検統計量(Neyman, 1949)
λ = −1/2 Freeman-Tukey 統計量(Freeman and Tukey, 1950)
In [115]:
using Plots
x = [2 4 3; 1 5 3; 4 2 3];
lambdas = collect(-2:0.1:1)
stats = Float64[]
for lambda in lambdas
    result = PowerDivergenceTest(x, lambda=lambda)
    append!(stats, result.stat)
end
plot(lambdas, stats, xlabel="λ", ylabel="test statistic",
     grid=false, tick_direction=:out, tickfont=("Times", 6),
        size=(340, 230), label="")
Out[115]:
−2 −1 0 1 3.4 3.6 3.8 4.0 4.2

$\chi^2$ 分布で近似する場合には λ ≒ 2/3 のときに最も近似効率が高いとされる。

PowerDivergenceTest() での lambda は実数で指定しなければならない(lambda=1 ではなく，lambda=1.0 と指定する）。

  • lambda=1.0 でピアソンの $\chi^2$ 検定を行う。
In [116]:
result1 = PowerDivergenceTest([10 3; 4 12], lambda=1.0)
p = pvalue(result1)
println("χ-sq. = $(result1.stat),  df = $(result1.df),  p value = $p")

χ-sq. = 7.743956043956044,  df = 1,  p value = 0.005389260671694263
In [117]:
chisq_test([10 3; 4 12]) # 上の例と同じである

χ-sq. = 7.743956043956044,  df = 1,  p value = 0.005389260671694263
  • lambda=0.0 で，対数尤度比に基づく $G^2$ 検定を行う。
In [118]:
result2 = PowerDivergenceTest([10 3; 4 12], lambda=0.0)
p = pvalue(result2)
println("G-sq. = $(result2.stat),  df = $(result2.df),  p value = $p")

G-sq. = 8.1280145470097,  df = 1,  p value = 0.00435864516471826
In [119]:
G2_test([10 3; 4 12]) # 上の例と同じである

G-sq. = 8.1280145470097,  df = 1,  p value = 0.00435864516471826

フィッシャーの正確検定¶

関数定義
FisherExactTest(a::Integer, b::Integer, c::Integer, d::Integer)

FisherExactTest() で用意されているのは 2×2 分割表の場合だけである。4 つのセルの度数を指定する。

帰無仮説は何通りかの表現法があるが，同じことを意味している。

  • 2 要因が独立である
  • 2 群の比率が等しい $\displaystyle \frac{a}{a+b} = \frac{c}{c+d}$
  • オッズ比が 1 である $\displaystyle \frac{a/c}{b/d} = \frac{a\cdot d}{b\cdot c}=1$

■ 使用例

In [120]:
result = FisherExactTest(10, 3, 4, 12)
Out[120]:
Fisher's exact test
-------------------
Population details:
    parameter of interest:   Odds ratio
    value under h_0:         1.0
    point estimate:          9.06784
    95% confidence interval: (1.422, 80.35)

Test summary:
    outcome with 95% confidence: reject h_0
    two-sided p-value:           0.0146

Details:
    contingency table:
        10   3
         4  12

■ 個別に表示できる統計量

  • $p$ 値
In [121]:
pvalue(result) # デフォールトで両側検定の p 値
Out[121]:
0.014589609220157725
In [122]:
pvalue(result, tail=:both) # 両側検定の p 値
Out[122]:
0.014589609220157725
In [123]:
pvalue(result, tail=:left) # 片側検定の p 値
Out[123]:
0.9994164940233701
In [124]:
pvalue(result, tail=:right) # 片側検定の p 値
Out[124]:
0.007294804610078863

FisherExactTest() の両側検定の $p$ 値は該当する片側確率を 2 倍するものである。

上の例では，0.014589609220157725 = 2 × 0.007294804610078863 になっている。

In [125]:
pvalue(result, method=:central) # 両側検定の p 値（デフォルト）
Out[125]:
0.014589609220157725

このような計算方法は R や matlab などの他の統計ソフトと違う。同じ結果を得るためには，オプション設定 method=:minlike としなければならない。

In [126]:
pvalue(result, method=:minlike)
Out[126]:
0.009220570313398513
In [127]:
using RCall # R で計算すると以下のようになる
R"""
options(digits=15)
x = matrix(c(10, 3, 4, 12), 2)
fisher.test(x)$p.value
"""
Out[127]:
RObject{RealSxp}
[1] 0.00922057031339851
  • ω　オッズ比
In [128]:
result.ω
Out[128]:
9.067840887689634
  • confint　オッズ比の信頼区間
In [129]:
confint(result)
Out[129]:
(1.4222662898175877, 80.34790877482563)

信頼率は，level で指定できる（0 < level < 1）

In [130]:
confint(result, level=0.95)
Out[130]:
(1.4222662898175877, 80.34790877482563)

■ オッズ比の $p$ 値と信頼区間について

オッズ比（$OR$）の $p$ 値，信頼区間には何通りかの方法がある。

FisherExactTest() は，R の fisher.test() と同じく，条件付き最尤推定量を求める。　外部リンク：Fisherの正確検定

一方，広く使われているのは条件なしの最尤推定量（Wald のオッズ比）で，以下の関数で計算できる。　外部リンク：統計学入門-第3章

In [131]:
using Distributions
function oddsratio(a, b, c, d; conf=0.95, correct=false)
    if correct || a*b*c*d == 0
        a, b, c, d = a+0.5, b+0.5, c+0.5, d+0.5
    end
    or = a*d / (b*c)
    LCL, UCL = or .* exp.([1, -1] .* quantile(Normal(), 0.5 - conf/2) * sqrt(1/a+1/b+1/c+1/d))
    println("オッズ比 = $or")
    println("$(Int(100conf))% 信頼区間  [$LCL, $UCL]")
end;
In [132]:
oddsratio(10, 3, 4, 12)

オッズ比 = 10.0
95% 信頼区間  [1.7975962537192371, 55.62984446206951]

R での実行結果を添付しておく。

In [133]:
using RCall
R"fisher.test(matrix(c(10, 3, 4, 12), 2))"
Out[133]:
RObject{VecSxp}

	Fisher's Exact Test for Count Data

data:  matrix(c(10, 3, 4, 12), 2)
p-value = 0.0092205703134
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
  1.42227515556236 80.34058365829780
sample estimates:
      odds ratio 
9.06776007785795

二項検定¶

二項検定は BinomialTest() で行う。

関数定義1
BinomialTest(x::Integer, n::Integer, p::Real = 0.5)
成功母比率 p の試行を n 回行い，x 回成功したとき，母比率が特定の値（デフォルトでは 0.5）であるかどうか検定する。

関数定義2
BinomialTest(x::AbstractVector{Bool}, p::Real = 0.5)
成功のとき 1，不成功のとき 0 の Bool 試行結果ベクトル x と母比率 p を与える場合。

たとえば，20回の実験結果が x = Bool[1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0] の場合，関数定義2 では ``BinomialTest(x)と書く。 成功が 12 回なので， 関数定義1 で書くとBinomialTest(12, 20)` とも書ける。

■ 使用例

In [134]:
result = BinomialTest(12, 20)
Out[134]:
Binomial test
-------------
Population details:
    parameter of interest:   Probability of success
    value under h_0:         0.5
    point estimate:          0.6
    95% confidence interval: (0.3605, 0.8088)

Test summary:
    outcome with 95% confidence: fail to reject h_0
    two-sided p-value:           0.5034

Details:
    number of observations: 20
    number of successes:    12

■ 個別に表示できる統計量

概要

  • $p$ 値

デフォルトでは両側検定の $p$ 値が表示される。

In [135]:
pvalue(result, tail=:left)
Out[135]:
0.8684120178222656
In [136]:
pvalue(result, tail=:right)
Out[136]:
0.2517223358154298
In [137]:
pvalue(result, tail=:both)
Out[137]:
0.5034446716308596

なお，両側検定の $p$ 値は該当する片側検定の $p$ 値を 2 倍したものであり，R や matlab などの他の統計ソフトとは違う。

計算方法により結果が違うのは，「母比率が 0.5 以外の両側検定」の場合に限られる。

以下の例，x=18, n=24, p=0.68 の場合，$p$ 値は 0.6206975448201958 が表示されるが，R では 0.5205908910845283 となる。

In [138]:
result2 = BinomialTest(18, 24, 0.68); # p=0.68 の場合。数値だけ指定するので注意。
pvalue(result2, tail=:both)
Out[138]:
0.6206975448201958
In [139]:
using RCall
R"""
options(digits=16)
binom.test(18, 24, 0.68)$p.value
"""
Out[139]:
RObject{RealSxp}
[1] 0.5205908910845283

BinomialTest() では，$n = 18,..., 24$ の二項分布の確率の合計が 0.3103487724100979 で，

In [140]:
pvalue(result, tail=:right)
Out[140]:
0.2517223358154298

その 2 倍が両側 $p$ 値になっている。

In [141]:
2 * pvalue(result, tail=:right)
Out[141]:
0.5034446716308596

R などでは，まず $n = 18,...,24$ を右側の棄却域とするところまでは同じであるが，もう一方の側の棄却域を次のようにして決める。

$x = 0, 1, ...$ と順に見ていき，観察値 18 のときの確率より小さい最大の $n$ まで，図でいえば $x = 0, ..., 14$ までを棄却域とする。

それぞれの棄却域における確率の合計を両側確率とする。

In [142]:
using Plots, Distributions
n = 24
p = 0.68
p_value = Float64[]
col = [i <= 14 || i >= 18 ? :red : :blue for i in 0:24]
for x = 0:n
    probability = pdf(Binomial(n, p), x)
    append!(p_value, probability)
end
bar(0:n, p_value, grid=false, tick_direction=:out,
    xlabel= "x", ylabel="probability", color=col, alpha=0.4,
    tickfont=("Times", 6), size=(340, 230), label="")
Out[142]:
0 5 10 15 20 25 0.00 0.05 0.10 0.15
In [143]:
sum(p_value[1:15]) + sum(p_value[19:25]) # Julia の添字は 1 から始まるので p_value[n+1] が n のときの確率である
Out[143]:
0.5205908910845278

FisherTest() の場合と異なり，同じ結果を得るためのオプション設定がないので，以下のようなプログラムを使用すればよい。

In [144]:
using Distributions
function binomial_pvalue(x, n, p=0.5; tail=:both)
    obj = Binomial(n, p)
    tail == :left  && return cdf(obj, x)
    tail == :right && return ccdf(obj, x - 1)
    p == 0 && return Float64(x == 0)
    p == 1 && return Float64(x == n)
    m = n * p
    x == m && return 1
    limit = 1.0000001pdf(obj, x)
    if x < m
        y = sum(pdf.(obj, ceil(Int, m):n) .<= limit)
        cdf(obj, x) + ccdf(obj, n - y)
    else
        y = sum(pdf.(obj, 0:floor(Int, m)) .<= limit)
        cdf(obj, y - 1) + ccdf(obj, x - 1)
    end
end
Out[144]:
binomial_pvalue (generic function with 2 methods)
In [145]:
binomial_pvalue(12, 20)
Out[145]:
0.5034446716308596
In [146]:
binomial_pvalue(4, 20, 0.6, tail=:left)
Out[146]:
0.00031703112116863026
In [147]:
binomial_pvalue(16, 20, 0.6, tail=:right)
Out[147]:
0.05095195319416655
  • 信頼区間の表示
In [148]:
result = BinomialTest(5, 26)
confint(result)
Out[148]:
(0.06554810873678252, 0.3935055279393218)
In [149]:
confint(result, level=0.99)
Out[149]:
(0.04400829697421974, 0.4550043049581578)
In [150]:
confint(result, tail=:right)
Out[150]:
(0.07898571028519472, 1.0)
In [151]:
confint(result, tail=:left)
Out[151]:
(0.0, 0.36259486197815627)

$t$ 検定¶

いわゆる $t$ 検定は，一標本の場合の母平均の検定と独立二標本の平均値の差の検定に大別される。更に，二群の平均値の差の検定は，二群の分散が等しいと仮定する場合と仮定しない場合に分かれる。

なお，対応のある二標本の平均値の差の検定もあるが，これは対応のあるデータの差をとって一標本の場合の母平均の検定をするのと同じである。

これらの検定は OneSampleTTest()，EqualvarianceTTest()，UnequalvarianceTTest() により行う。

一標本の検定（母平均の検定）¶

関数定義
OneSampleTTest(v::AbstractVector{T<:Real}, μ0::Real = 0)

対応のある二標本の場合，対応のあるデータの差をとって一標本の検定を行うのと同じである。

■ 使用例

In [152]:
x = [5.6, 8.1, 6.6, 4.7, 4.8, 9.6, 5.6, 1.3, 4.1, 5.4];
y = [2.0, 3.4, 5.2, 8.0, 9.0, 3.6, 2.7, 0.7, 6.7, 5.7];
z = x .- y;
In [153]:
using HypothesisTests
result = OneSampleTTest(x, y)
Out[153]:
One sample t-test
-----------------
Population details:
    parameter of interest:   Mean
    value under h_0:         0
    point estimate:          0.88
    95% confidence interval: (-1.614, 3.374)

Test summary:
    outcome with 95% confidence: fail to reject h_0
    two-sided p-value:           0.4454

Details:
    number of observations:   10
    t-statistic:              0.7981113921334072
    degrees of freedom:       9
    empirical standard error: 1.1026029808291529
In [154]:
result = OneSampleTTest(z)
Out[154]:
One sample t-test
-----------------
Population details:
    parameter of interest:   Mean
    value under h_0:         0
    point estimate:          0.88
    95% confidence interval: (-1.614, 3.374)

Test summary:
    outcome with 95% confidence: fail to reject h_0
    two-sided p-value:           0.4454

Details:
    number of observations:   10
    t-statistic:              0.7981113921334072
    degrees of freedom:       9
    empirical standard error: 1.1026029808291529

■ 個別に表示できる統計量

  • $p$ 値
In [155]:
pvalue(result)
Out[155]:
0.4453534218554348
In [156]:
pvalue(result, tail=:both)
Out[156]:
0.4453534218554348
In [157]:
pvalue(result, tail=:left)
Out[157]:
0.7773232890722825
In [158]:
pvalue(result, tail=:right)
Out[158]:
0.2226767109277174
  • 信頼区間　信頼区間は母平均の信頼区間である。
In [159]:
confint(result)
Out[159]:
(-1.6142612308053204, 3.3742612308053195)
In [160]:
confint(result, tail=:both)
Out[160]:
(-1.6142612308053204, 3.3742612308053195)
In [161]:
confint(result, level=0.99)
Out[161]:
(-2.7032783553640023, 4.463278355364001)
In [162]:
confint(result, tail=:right)
Out[162]:
(-1.141195783743237, Inf)
In [163]:
confint(result, tail=:right, level=0.99)
Out[163]:
(-2.2309258663578753, Inf)
In [164]:
confint(result, tail=:left)
Out[164]:
(-Inf, 2.901195783743236)
  • n　サンプルサイズ
In [165]:
result.n
Out[165]:
10
In [166]:
result.xbar
Out[166]:
0.8799999999999997
  • t　検定統計量（$t$ 分布にしたがう）
In [167]:
result.t
Out[167]:
0.7981113921334072
  • df　$t$ 分布の自由度
In [168]:
result.df
Out[168]:
9
  • xbar　標本平均
    xbar = mean(x) - mean(y) = mean(x - y) = mean(z)
In [169]:
result.xbar
Out[169]:
0.8799999999999997
  • μ0　母平均
In [170]:
result.μ0
Out[170]:
0
  • stderr　標準誤差
    result.t = (result.xbar - result.μ0) / result.stderr
In [171]:
result.stderr
Out[171]:
1.1026029808291529

等分散の場合の $t$ 検定¶

独立 2 標本の分散が等しいことを仮定する平均値の差の検定である。

一般的には，等分散を仮定しない $t$ 検定（Welch の方法）が推奨される。

関数定義
EqualVarianceTTest(x::AbstractVector{T<:Real}, y::AbstractVector{T<:Real})

■ 使用例

In [172]:
x = [1, 2, 3, 6, 3, 1, 2, 1, 1, 1, 3, 4];
y = [2, 1, 2, 3, 4, 2, 1, 2, 3, 5];
In [173]:
result = EqualVarianceTTest(x, y)
Out[173]:
Two sample t-test (equal variance)
----------------------------------
Population details:
    parameter of interest:   Mean difference
    value under h_0:         0
    point estimate:          -0.166667
    95% confidence interval: (-1.448, 1.115)

Test summary:
    outcome with 95% confidence: fail to reject h_0
    two-sided p-value:           0.7889

Details:
    number of observations:   [12,10]
    t-statistic:              -0.271312734545191
    degrees of freedom:       20
    empirical standard error: 0.6142972497994038

■ 個別に表示できる統計量

  • p　$p$ 値
In [174]:
pvalue(result) # tail=:both 両側検定がデフォルト
Out[174]:
0.788931157688499
In [175]:
pvalue(result, tail=:left) # 棄却域が左にある片側検定
Out[175]:
0.3944655788442495
In [176]:
pvalue(result, tail=:right) # 棄却域が右にある片側検定
Out[176]:
0.6055344211557505
  • 信頼区間
In [177]:
confint(result) #デフォルトで両側に棄却域がある場合の信頼率 95% の信頼区間
Out[177]:
(-1.448068275504171, 1.114734942170838)
In [178]:
confint(result, level=0.99) # 信頼率 99%
Out[178]:
(-1.9145510251333062, 1.5812176917999732)
In [179]:
confint(result, tail=:left, level=0.99) # 片側と信頼率の指定
Out[179]:
(-Inf, 1.386262653673622)
  • t　検定統計量 $t$
In [180]:
result.t
Out[180]:
-0.271312734545191
  • df　$t$ 分布の自由度
In [181]:
result.df
Out[181]:
20
  • n_x, n_y　各群のサンプルサイズ
In [182]:
result.n_x, result.n_y
Out[182]:
(12, 10)
  • xbar, μ0　平均値の差と母平均
    xbar = mean(x) - mean(y)
In [183]:
result.xbar, result.μ0
Out[183]:
(-0.16666666666666652, 0)
  • stderr　標準偏差
    t = (result.xbar - result.μ0) / result.stderr
In [184]:
result.stderr
Out[184]:
0.6142972497994038

等分散でない場合 Welch の方法¶

関数定義
UnequalVarianceTTest(x::AbstractVector{T<:Real}, y::AbstractVector{T<:Real})

■ 使用例

In [185]:
x = [1, 2, 3, 6, 3, 1, 2, 1, 1, 1, 3, 4];
y = [2, 1, 2, 3, 4, 2, 1, 2, 3, 5];
In [186]:
result = UnequalVarianceTTest(x, y)
Out[186]:
Two sample t-test (unequal variance)
------------------------------------
Population details:
    parameter of interest:   Mean difference
    value under h_0:         0
    point estimate:          -0.166667
    95% confidence interval: (-1.424, 1.09)

Test summary:
    outcome with 95% confidence: fail to reject h_0
    two-sided p-value:           0.7849

Details:
    number of observations:   [12,10]
    t-statistic:              -0.2765775336645321
    degrees of freedom:       19.99676443938765
    empirical standard error: 0.6026037778933294

■ 個別に表示できる統計量

  • p　$p$ 値
In [187]:
pvalue(result) # tail=:both 両側検定がデフォルト
Out[187]:
0.784943053185432
In [188]:
pvalue(result, tail=:left) # 棄却域が左にある片側検定
Out[188]:
0.392471526592716
In [189]:
pvalue(result, tail=:right) # 棄却域が右にある片側検定
Out[189]:
0.607528473407284
  • 信頼区間
In [190]:
confint(result) #デフォルトで両側に棄却域がある場合の信頼率 95% の信頼区間
Out[190]:
(-1.423689159406056, 1.0903558260727229)
In [191]:
confint(result, level=0.99) # 信頼率 99%
Out[191]:
(-1.8813078766019782, 1.5479745432686451)
In [192]:
confint(result, tail=:left, level=0.99) # 片側と信頼率の指定
Out[192]:
(-Inf, 1.3567230542137687)
  • t　検定統計量 $t$
In [193]:
result.t
Out[193]:
-0.2765775336645321
  • df　$t$ 分布の自由度
In [194]:
result.df
Out[194]:
19.99676443938765
  • n_x, n_y　各群のサンプルサイズ
In [195]:
result.n_x, result.n_y
Out[195]:
(12, 10)
  • xbar, μ0　平均値の差と母平均
    xbar = mean(x) - mean(y)
In [196]:
result.xbar, result.μ0
Out[196]:
(-0.16666666666666652, 0)
  • stderr　標準偏差
    t = (result.xbar - result.μ0) / result.stderr
In [197]:
result.stderr
Out[197]:
0.6026037778933294
In [198]:
result.t
Out[198]:
-0.2765775336645321

マン・ホイットニーの $U$ 検定¶

関数定義
MannWhitneyUTest(x::AbstractVector{<:Real}, y::AbstractVector{<:Real})

ExactMannWhitneyUTest(x::AbstractVector{<:Real}, y::AbstractVector{<:Real})
ApproximateMannWhitneyUTest(x::AbstractVector{<:Real}, y::AbstractVector{<:Real})

マン・ホイットニーの $U$ 検定は，ウィルコクソンの順位和検定とも呼ばれる。

順序尺度以上の独立二標本の代表値の差の検定である。

MannWhitneyUTest() はラッパーであり，実際に適用される関数は ExactMannWhitneyUTest() か ApproximateMannWhitneyUTest() である。

以下の 2 つの条件のいずれかを満たす場合には正確な検定 ExactMannWhitneyUTest() が適用される。

  1. 同順位データ（tie）がなく，プールされたサンプルサイズが 50 以下の場合
  2. 同順位があるが，プールされたサンプルサイズが 10 以下の場合

それ以外の場合には漸近近似検定 ApproximateMannWhitneyUTest() が適用される。

しかし，ExactMannWhitneyUTest() はデータの並べ替えのすべてを評価するので，程々のサンプルサイズの場合でも計算時間がかかりすぎることがある。

そのため，条件を満たすかどうかに限らず ExactMannWhitneyUTest() か ApproximateMannWhitneyUTest() を直接選択することもできる。

■ 使用例1

以下のデータは，プールされたサンプルサイズが 22 であるが同順位があるので，条件 1 に反する。また，同順位があるがプールされたサンプルサイズが 11 以上なので，条件 2 にも反する。 したがって，何も指示しない場合は ApproximateMannWhitneyUTest(x, y) が適用される。

In [199]:
x = [1, 2, 3, 6, 3, 1, 2, 1, 1, 1, 3, 4];
y = [2, 1, 2, 3, 4, 2, 1, 2, 3, 5];

approx_result = ApproximateMannWhitneyUTest(x, y)
Out[199]:
Approximate Mann-Whitney U test
-------------------------------
Population details:
    parameter of interest:   Location parameter (pseudomedian)
    value under h_0:         0
    point estimate:          0.0

Test summary:
    outcome with 95% confidence: fail to reject h_0
    two-sided p-value:           0.6334

Details:
    number of observations in each group: [12, 10]
    Mann-Whitney-U statistic:             52.5
    rank sums:                            [137.0, 116.0]
    adjustment for ties:                  672.0
    normal approximation (μ, σ):          (-7.5, 14.6784)

この結果は R の wilcox.test() と同じ結果になる。

In [200]:
using RCall
R"""
x = c(1, 2, 3, 6, 3, 1, 2, 1, 1, 1, 3, 4)
y = c(2, 1, 2, 3, 4, 2, 1, 2, 3, 5)
wilcox.test(x, y)
"""

┌ Warning: RCall.jl: wilcox.test.default(x, y) で警告がありました:
│   タイがあるため、正確な p 値を計算することができません
└ @ RCall ~/.julia/packages/RCall/0ggIQ/src/io.jl:172
Out[200]:
RObject{VecSxp}

	Wilcoxon rank sum test with continuity correction

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

あえて ExactMannWhitneyUTest() を適用した結果を以下に示す。

In [201]:
exact_result = ExactMannWhitneyUTest(x, y)
Out[201]:
Exact Mann-Whitney U test
-------------------------
Population details:
    parameter of interest:   Location parameter (pseudomedian)
    value under h_0:         0
    point estimate:          0.0

Test summary:
    outcome with 95% confidence: fail to reject h_0
    two-sided p-value:           0.6305

Details:
    number of observations in each group: [12, 10]
    Mann-Whitney-U statistic:             52.5
    rank sums:                            [137.0, 116.0]
    adjustment for ties:                  672.0

$p$ 値は 0.6305 となり，ApproximateMannWhitneyUTest() の場合の 0.6334388940675 とは異った値になった。

しかし，R の coin パッケージを使った正確な検定とも違う。

In [202]:
using RCall
R"""
x = c(1, 2, 3, 6, 3, 1, 2, 1, 1, 1, 3, 4)
y = c(2, 1, 2, 3, 4, 2, 1, 2, 3, 5)
val = c(x, y)
g = factor(rep(c("a", "b"), c(length(x), length(y))))
df = data.frame(val, g)
library(coin)
wilcox_test(val ~ g, data=df, distribution="exact")
"""

┌ Warning: RCall.jl: 要求されたパッケージ survival をロード中です
└ @ RCall ~/.julia/packages/RCall/0ggIQ/src/io.jl:172
Out[202]:
RObject{S4Sxp}

	Exact Wilcoxon-Mann-Whitney Test

data:  val by g (a, b)
Z = -0.51095591724442, p-value = 0.6326831063673
alternative hypothesis: true mu is not equal to 0

分割表形式のデータを入力して正確な $p$ 値を計算できるオンラインサイトで計算しても，coin ライブラリの結果と同じになる。

同順位のない 50 以下のサンプルサイズデータでも確認してみよう。

In [203]:
x = [1.99, 14.93, 7.51, 11.58, 4.75, 3.86, 2.74, 5.66, 11.32, 4.63, 8.89, 17.99]
y = [3.04, 19.18, 13.72, 19.31, 15.48, 17.92, 16.94, 7.99, 18.75, 5.4]
approx_result = ApproximateMannWhitneyUTest(x, y)
Out[203]:
Approximate Mann-Whitney U test
-------------------------------
Population details:
    parameter of interest:   Location parameter (pseudomedian)
    value under h_0:         0
    point estimate:          -9.625

Test summary:
    outcome with 95% confidence: reject h_0
    two-sided p-value:           0.0321

Details:
    number of observations in each group: [12, 10]
    Mann-Whitney-U statistic:             27.0
    rank sums:                            [164.0, 89.0]
    adjustment for ties:                  0.0
    normal approximation (μ, σ):          (-33.0, 15.1658)
In [204]:
exact_result = ExactMannWhitneyUTest(x, y)
Out[204]:
Exact Mann-Whitney U test
-------------------------
Population details:
    parameter of interest:   Location parameter (pseudomedian)
    value under h_0:         0
    point estimate:          -9.625

Test summary:
    outcome with 95% confidence: reject h_0
    two-sided p-value:           0.0300

Details:
    number of observations in each group: [12, 10]
    Mann-Whitney-U statistic:             27.0
    rank sums:                            [164.0, 89.0]
    adjustment for ties:                  0.0
In [205]:
pvalue(approx_result) |> println
pvalue(exact_result)  |> println

0.03211417930790749
0.029960751322980424

同順位がまったくないデータの場合には，R での結果と一致する。

In [206]:
using RCall
R"""
x = c(1.99, 14.93, 7.51, 11.58, 4.75, 3.86, 2.74, 5.66, 11.32, 4.63, 8.89, 17.99)
y = c(3.04, 19.18, 13.72, 19.31, 15.48, 17.92, 16.94, 7.99, 18.75, 5.4)
val = c(x, y)
g = factor(rep(c("a", "b"), c(length(x), length(y))))
df = data.frame(val, g)
library(coin)
wilcox_test(val ~ g, data=df, distribution="exact")
"""
Out[206]:
RObject{S4Sxp}

	Exact Wilcoxon-Mann-Whitney Test

data:  val by g (a, b)
Z = -2.1759555622061, p-value = 0.02996075132298
alternative hypothesis: true mu is not equal to 0

外部リンク:Mann-Whitney U testにおいて，

When there are no tied ranks, the exact p-value is computed using the pwilcox function from the Rmath package.

と述べているが，ここでいっている「exact p-value」は漸近分布ではなく（$U$ 検定統計量の平均値と分散により，標準正規分布から $p$ 値を求めるのではなく），pwilcox で計算したといっているだけで，「exact test による $p$ 値ではない」。

また，

In the presence of tied ranks, a p-value is computed by exhaustive enumeration of permutations, which can be very slow for even moderately sized data sets.

といってはいるが，同順位のあるデータの場合はデータの順列についてすべての評価をしているようには見えない。

なお，プールしたサンプルサイズが 10 以下であれば，「1.12. 並べ替え検定（無作為検定）」を用いて，同順位があっても正確な $p$ 値を求めることができる。

In [207]:
using StatsBase # tiedrank
function ExactU(x, y)
    z = tiedrank(vcat(x, y...))
    length(z) > 10 && throw(ArgumentError("sample size must be less than 11."))
    pvalue(ExactPermutationTest(z[1:length(x)], z[length(x)+1:end], mean))
end

x = [1, 2, 3, 6, 3];
y = [2, 1, 2, 3, 4];
ExactU(x, y) # p value = 0.7142857142857143
Out[207]:
0.7142857142857143

符号検定¶

対応のある二標本の母代表値（中央値）に差があるかどうかの検定である。2 変数の組で，いずれが優れているか，劣っているか，あるいは同等であるかしかわからないときに使われる。どの程度優れているかがわかるときには「11.　ウィルコクソンの符号付き順位和検定」を用いる。

関数定義
SignTest(x::AbstractVector{T<:Real}, median::Real = 0)
SignTest(x::AbstractVector{T<:Real}, y::AbstractVector{T<:Real}, median::Real = 0)

SignTest(x, y, median=0)
SignTest(z, median=0)
x, y は対応のあるデータベクトル。z は x .- y を事前に計算した場合。median は帰無仮説下の中央値。

■ 使用例

In [208]:
using HypothesisTests

x = [3, 1, 5, 4, 2, 2, 4, 1, 5, 4]
y = [1, 1, 3, 2, 4, 1, 3, 2, 4, 2]
result = SignTest(x, y)
Out[208]:
Sign Test
---------
Population details:
    parameter of interest:   Median
    value under h_0:         0.0
    point estimate:          1.0
    95% confidence interval: (0.0, 2.0)

Test summary:
    outcome with 95% confidence: fail to reject h_0
    two-sided p-value:           0.1797

Details:
    number of observations:       9
    observations larger than 0.0: 7
In [209]:
result2 = SignTest([1, 0, 1, 1, -1, 1, 1, -1, 1, 1])
Out[209]:
Sign Test
---------
Population details:
    parameter of interest:   Median
    value under h_0:         0.0
    point estimate:          1.0
    95% confidence interval: (0.0, 1.0)

Test summary:
    outcome with 95% confidence: fail to reject h_0
    two-sided p-value:           0.1797

Details:
    number of observations:       9
    observations larger than 0.0: 7

対応する観察値が同じではならないデータ対の数を n，median より大きいものの対の数を x としたときの BinomialTest(x, n) と同じである（$p$ 値が同じになる）。

In [210]:
using HypothesisTests
result3 = BinomialTest(7, 9)
Out[210]:
Binomial test
-------------
Population details:
    parameter of interest:   Probability of success
    value under h_0:         0.5
    point estimate:          0.777778
    95% confidence interval: (0.3999, 0.9719)

Test summary:
    outcome with 95% confidence: fail to reject h_0
    two-sided p-value:           0.1797

Details:
    number of observations: 9
    number of successes:    7
In [211]:
pvalue(result)  |> println
pvalue(result2) |> println
pvalue(result3) |> println

0.17968750000000006
0.17968750000000006
0.17968750000000006

■ 個別に表示できる統計量

  • p　$p$ 値
In [212]:
pvalue(result) # デフォルトで両側検定
Out[212]:
0.17968750000000006
In [213]:
pvalue(result, tail=:left) # 棄却域が左にある場合の片側検定
Out[213]:
0.98046875
In [214]:
pvalue(result, tail=:right) # 棄却域が右にある場合の片側検定
Out[214]:
0.08984375000000003
  • confint　信頼区間
In [215]:
confint(result) # デフォルトで両側に棄却域がある場合で，信頼率が 95% の信頼区間
Out[215]:
(0, 2)
In [216]:
confint(result, level=0.99) # 信頼率の指定
Out[216]:
(-1, 2)
In [217]:
confint(result, tail=:left, level=0.99) # 棄却域が左にある場合の信頼区間
Out[217]:
(-1, 0.0)
In [218]:
confint(result, tail=:right, level=0.99) # 棄却域が右にある場合の信頼区間
Out[218]:
(0.0, 2)

ウィルコクソンの符号付き順位和検定¶

対応のある二標本の母代表値（中央値）に差があるかどうかの検定である。2 変数の組で，どちらがどの程度優れているかがわかるときにはウィルコクソンの符号付き順位和検定を用いる。

関数定義
SignedRankTest(x::AbstractVector{<:Real})
SignedRankTest(x::AbstractVector{<:Real}, y::AbstractVector{<:Real})

ExactSignedRankTest(x::AbstractVector{<:Real}[, y::AbstractVector{<:Real}])
ApproximateSignedRankTest(x::AbstractVector{<:Real}[, y::AbstractVector{<:Real}])

SignedRankTest(z)
SignedRankTest(x，y)
x, y は対応のあるデータベクトル。z は x .- y を事前に計算した場合。

SignedRankTest() はラッパーであり，実際に適用される関数は ExactSignedRankTest() か ApproximateSignedRankTest() である。

以下の 2 つの条件を満たすときには，ExactSignedRankTest() が実行される。

  1. 同順位がなくて，サンプルサイズが 50 以下である。
  2. 同順位があるが，サンプルサイズが 15 以下である。

どちらの条件も満たさないときは，ApproximateSignedRankTest() が実行される。

条件を満たすかどうかに限らず，ApproximateSignedRankTest() か ExactSignedRankTest() を直接選択できる。

■ 使用例

In [219]:
using HypothesisTests

x = [3, 1, 5, 4, 2, 2, 4, 1, 5, 4]
y = [1, 1, 3, 2, 4, 1, 3, 2, 4, 2]
result = SignedRankTest(x, y)
Out[219]:
Exact Wilcoxon signed rank test
-------------------------------
Population details:
    parameter of interest:   Location parameter (pseudomedian)
    value under h_0:         0
    point estimate:          1.0
    95% confidence interval: (0.0, 2.0)

Test summary:
    outcome with 95% confidence: fail to reject h_0
    two-sided p-value:           0.1562

Details:
    number of observations:      10
    Wilcoxon rank-sum statistic: 35.5
    rank sums:                   [35.5, 9.5]
    adjustment for ties:         180.0
In [220]:
z = x .- y
result2 = SignedRankTest(z)
Out[220]:
Exact Wilcoxon signed rank test
-------------------------------
Population details:
    parameter of interest:   Location parameter (pseudomedian)
    value under h_0:         0
    point estimate:          1.0
    95% confidence interval: (0.0, 2.0)

Test summary:
    outcome with 95% confidence: fail to reject h_0
    two-sided p-value:           0.1562

Details:
    number of observations:      10
    Wilcoxon rank-sum statistic: 35.5
    rank sums:                   [35.5, 9.5]
    adjustment for ties:         180.0

■ 個別に表示できる統計量

  • p　$p$ 値
In [221]:
pvalue(result)
Out[221]:
0.15625
In [222]:
pvalue(result, tail=:left)
Out[222]:
0.9609375
In [223]:
pvalue(result, tail=:right)
Out[223]:
0.078125
  • confint　信頼区間
In [224]:
confint(result)
Out[224]:
(0.0, 2.0)
In [225]:
confint(result, level=0.99)
Out[225]:
(-1.0, 2.0)
In [226]:
confint(result, tail=:left, level=0.99)
Out[226]:
(-0.5, Inf)
In [227]:
confint(result, tail=:right, level=0.99)
Out[227]:
(-Inf, 2.0)
  • W　検定統計量

result3.W

■ 漸近近似検定の場合

In [228]:
using Random; Random.seed!(123)
x = randn(100)
y = randn(100) .+ 0.1
result3 = ApproximateSignedRankTest(x, y)
Out[228]:
Approximate Wilcoxon signed rank test
-------------------------------------
Population details:
    parameter of interest:   Location parameter (pseudomedian)
    value under h_0:         0
    point estimate:          0.0187289
    95% confidence interval: (-0.3522, 0.1781)

Test summary:
    outcome with 95% confidence: fail to reject h_0
    two-sided p-value:           0.5763

Details:
    number of observations:      100
    Wilcoxon rank-sum statistic: 2362.0
    rank sums:                   [2362.0, 2688.0]
    adjustment for ties:         0.0
    normal approximation (μ, σ): (-163.0, 290.839)
  • mu, sigma　ApproximateSignedRankTest() の場合のみ
In [229]:
result3.mu, result3.sigma
Out[229]:
(-163.0, 290.8393027085576)
In [230]:
using Distributions
z = result3.mu / result3.sigma
Out[230]:
-0.5604469495078456
In [231]:
2*cdf(Normal(), z)
Out[231]:
0.5751746153622865
In [232]:
pvalue(result3) # 一致しない？
Out[232]:
0.576347512365529

相関係数の検定¶

ここでは母相関係数が 0 である（無相関）かどうかの検定を取り上げる。HypothesTests にはないが，簡単な関数を作っておこう。

関数定義　2 変数をベクトル x, y で与える場合
CorTest(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}; level=0.95, method="pearson")

関数定義　サンプルサイズ n と 相関係数 r で与える場合 CorTest(n::Int, r::Float64; level=0.95, method="pearson")

method には以下のものが指定できる。

method 相関係数
"pearson" ピアソンの積率相関係数
"spearman" スピアマンの順位相関係数
"kendall" ケンドールの順位相関係数

スピアマンの順位相関係数，ケンドールの順位相関係数においては，同順位がある場合には近似検定になる。

ケンドールの順位相関係数においては，信頼区間は求めない。

In [233]:
using Statistics, StatsBase, Distributions
roundn(x, digits) = round(x, digits=digits)
formatp(p) = p < 0.00001 ? "< 0.00001" : "= " * string(roundn(p, 5))
index(method) = indexin(method[1], ['p', 's', 'k'])[1]
function CorTest(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}; level=0.95, method="pearson")
    func = [cor, corspearman, corkendall][index(method)]
    CorTest(length(x), func(x, y); level, method)
end
function CorTest(n::Int, r::Float64; level=0.95, method="pearson")
    if method != "kendall"
        t = r * sqrt((n-2)/(1-r^2))
        df = n-2
        p = 2ccdf(TDist(df), abs(t))
        q = quantile(Normal(), 0.5 - 0.5level)
        conf = tanh.(atanh(r) .+ [q, -q] ./ sqrt(n-3))
        println("n = $n, r = $(roundn(r, 5)), t = $(roundn(t, 5)), df = $df, p-value $(formatp(p))")
        println("$(Int(100level)) percent confidence interval = $(roundn.(conf, 5))")
    else
        z = r / sqrt((4n + 10) / (9n * (n - 1)))
        p = 2ccdf(Normal(), abs(z))
        println("n = $n, r = $(roundn(r, 5)), z = $(roundn(z, 5)), p-value $(formatp(p))")
    end
end;

■ 2 変数をベクトルで与える場合

In [234]:
x = [58.3, 43.6, 49.9, 64.1, 51.6, 36.6, 58.7, 66.7, 50.7, 53.1, 
     56.1, 42.0, 45.5, 30.9, 42.3]
y = [46.3, 51.5, 64.5, 55.6, 63.9, 35.6, 61.2, 56.9, 44.2, 57.3, 
     42.7, 47.1, 32.4, 39.7, 51.0]
CorTest(x, y)

n = 15, r = 0.52225, t = 2.20803, df = 13, p-value = 0.04582
95 percent confidence interval = [0.01363, 0.81616]
In [235]:
CorTest(x, y, method="spearman")

n = 15, r = 0.47143, t = 1.92737, df = 13, p-value = 0.07607
95 percent confidence interval = [-0.05384, 0.79234]
In [236]:
CorTest(x, y, method="kendall")

n = 15, r = 0.31429, z = 1.63308, p-value = 0.10245
In [237]:
CorTest(x, y, level=0.99) # 信頼率の設定は，level= で指定する

n = 15, r = 0.52225, t = 2.20803, df = 13, p-value = 0.04582
99 percent confidence interval = [-0.16269, 0.86753]

■ サンプルサイズと相関係数を与える場合

In [238]:
CorTest(length(x), cor(x, y))

n = 15, r = 0.52225, t = 2.20803, df = 13, p-value = 0.04582
95 percent confidence interval = [0.01363, 0.81616]
In [239]:
CorTest(24, 0.476, level=0.9)

n = 24, r = 0.476, t = 2.53869, df = 22, p-value = 0.01871
90 percent confidence interval = [0.15754, 0.70478]

対応のない $k$ 標本（独立 $k$ 標本）¶

対応のない $k$ 標本は，対応のない 2 標本の拡張である。

対応のない 2 標本に対する $t$ 検定は，一元配置分散分析に対応する。

一元元配置分散分析¶

HypothesisTests の OneWayANOVATest() には，2022/05/21 現在でバグがある。

関数定義
OneWayANOVATest(groups::AbstractVector{<:Real}...)

引数として，各標本の観測データベクトルを列挙する。

In [240]:
using HypothesisTests

g1 = [27.7, 45.2, 32.8, 49.5, 31.0, 55.5, 31.7, 53.9];
g2 = [38.0, 47.4, 55.3, 52.5, 39.2, 34.9];
g3 = [44.0, 29.8, 42.5, 23.4, 20.7, 57.3, 42.0, 56.8, 44.1, 32.7];

result = OneWayANOVATest(g1, g2, g3)
Out[240]:
One-way analysis of variance (ANOVA) test
-----------------------------------------
Population details:
    parameter of interest:   Means
    value under h_0:         "all equal"
    point estimate:          NaN

Test summary:
    outcome with 95% confidence: fail to reject h_0
    p-value:                     0.6612

Details:
    number of observations: [8, 6, 10]
    F statistic:            0.421994
    degrees of freedom:     (2, 21)

正しいプログラム例を挙げておく。

In [241]:
using Statistics, Distributions

function EqualVarianceOneWayANOVATest(x::AbstractVector{<:Real}...)
    roundn(x, digits) = round(x, digits=digits)
    formatp(p) = p < 0.00001 ? "< 0.00001" : "= " * string(roundn(p, 5))
    X = vcat(x...)
    n, grandmean = length(X), mean(X)
    St = var(X) * (n -1)
    k = length(x)
    ns = zeros(Int, k)
    means = zeros(k)
    for (i, d) in enumerate(x)
        ns[i] = length(d)
        means[i] = mean(d)
    end
    Sb = sum(ni * (meani - grandmean)^2 for (ni, meani) in zip(ns, means))
    Sw = St - Sb
    dfb, dfw = k - 1, n - k
    Vb, Vw = Sb / dfb, Sw / dfw
    F = Vb / Vw
    p = ccdf(FDist(dfb, dfw), F)
    println("F = $(roundn(F, 5)),  df1 = $dfb,  df2 = $dfw,  p-value $(formatp(p))")
end
Out[241]:
EqualVarianceOneWayANOVATest (generic function with 1 method)

■ 使用例

In [242]:
g1 = [27.7, 45.2, 32.8, 49.5, 31.0, 55.5, 31.7, 53.9];
g2 = [38.0, 47.4, 55.3, 52.5, 39.2, 34.9];
g3 = [44.0, 29.8, 42.5, 23.4, 20.7, 57.3, 42.0, 56.8, 44.1, 32.7];

EqualVarianceOneWayANOVATest(g1, g2, g3)

F = 0.40416,  df1 = 2,  df2 = 21,  p-value = 0.67262

一元元配置分散分析（ウェルチの方法）¶

対応のない 2 標本において，等分散を仮定しない平均値の差の検定（ウェルチの方法）の拡張である。

R ではウェルチの方法による一元元配置分散分析がデフォルトなのであるが，Julia にはないので以下の関数を定義しておく。

関数定義
UnequalVarianceOneWayANOVATestOneWayANOVATest(groups::AbstractVector{<:Real}...)

OneWayANOVATest(groups::AbstractVector{<:Real}...) の場合と同じく，引数として各標本の観測データベクトルを列挙する。

In [243]:
using Statistics, Distributions

function UnequalVarianceOneWayANOVATest(groups::AbstractVector{<:Real}...)
    roundn(x, digits) = round(x, digits=digits)
    formatp(p) = p < 0.00001 ? "< 0.00001" : "= " * string(roundn(p, 5))
    ni = Int[]
    mi = Float64[]
    vi = Float64[]
    for x in groups
        append!(ni, length(x))
        append!(mi, mean(x))
        append!(vi, var(x))
    end
    k = length(ni)
    wi = ni ./ vi
    sum_wi = sum(wi)
    tmp = sum((1 .- wi ./ sum_wi).^2 ./ (ni .- 1)) ./ (k^2 - 1)
    m = sum(wi .* mi) / sum_wi
    df1, df2 = k - 1, 1 / 3tmp
    F = sum(wi .* (mi .- m).^2) / (df1 * (1 + 2 * (k - 2) * tmp))
    p = ccdf(FDist(df1, df2), F)
    println("F = $(roundn(F, 5)),  df1 = $df1,  df2 = $(roundn(df2, 5)),  p-value $(formatp(p))")
end;

■ 使用例

In [244]:
UnequalVarianceOneWayANOVATest(g1, g2, g3)

F = 0.51339,  df1 = 2,  df2 = 13.589,  p-value = 0.60962

クラスカル・ウォリス検定¶

対応のない 2 標本においけるマン・ホイットニーの $U$ 検定の拡張である。

関数定義
KruskalWallisTest(groups::AbstractVector{<:Real}...)

引数として各標本の観測データベクトルを列挙する。

■ 使用例

In [245]:
g1 = [27.7, 45.2, 32.8, 49.5, 31.0, 55.5, 31.7, 53.9];
g2 = [38.0, 47.4, 55.3, 52.5, 39.2, 34.9];
g3 = [44.0, 29.8, 42.5, 23.4, 20.7, 57.3, 42.0, 56.8, 44.1, 32.7];

result = KruskalWallisTest(g1, g2, g3)
Out[245]:
Kruskal-Wallis rank sum test (chi-square approximation)
-------------------------------------------------------
Population details:
    parameter of interest:   Location parameters
    value under h_0:         "all equal"
    point estimate:          NaN

Test summary:
    outcome with 95% confidence: fail to reject h_0
    one-sided p-value:           0.7082

Details:
    number of observation in each group: [8, 6, 10]
    χ²-statistic:                        0.69
    rank sums:                           [98.0, 87.0, 115.0]
    degrees of freedom:                  2
    adjustment for ties:                 1.0

■ 個別に表示できる統計量

  • H　検定統計量（$\chi^2$ 分布にしたがう）
In [246]:
result.H
Out[246]:
0.6899999999999977
  • df　$p$ 値を求めるとき，$\chi^2$ 分布で近似するときの自由度
In [247]:
result.df
Out[247]:
2
  • p　$p$ 値
    "one-sided p-value:" と表示されているが，そもそもこの検定は両側/片側の区別はない。単に "p-value:" とすべきである。当然，tail=:left も tail=:right は言うに及ばず tail=:both さえも指定できない。
In [248]:
pvalue(result)
Out[248]:
0.7082203534678008

コルモゴロフ・スミルノフ検定¶

1 標本と母分布の差の検定，及び，2 標本の分布の差の検定がある。

1 標本の分布の検定¶

関数定義
ExactOneSampleKSTest(x::AbstractVector{<:Real}, d::UnivariateDistribution)
ApproximateOneSampleKSTest(x::AbstractVector{<:Real}, d::UnivariateDistribution)

  • x は観察ベクトル
  • d は 1 変数分布関数
    Normal(μ,σ), Uniform(a,b), Weibull(α,θ), Beta(α, β), Cauchy(μ, σ), Exponential(θ), Gamma(α,θ), Gumbel(μ, θ), Logistic(μ,θ), LogNormal(μ,σ) 等がある。　外部リンク

■ 使用例

In [249]:
using HypothesisTests
using Distributions
using Random; Random.seed!(123)
x = randn(1000);
  • 正確な検定
In [250]:
result = ExactOneSampleKSTest(x, Normal()) # 標準正規分布に從うか
Out[250]:
Exact one sample Kolmogorov-Smirnov test
----------------------------------------
Population details:
    parameter of interest:   Supremum of CDF differences
    value under h_0:         0.0
    point estimate:          0.0300133

Test summary:
    outcome with 95% confidence: fail to reject h_0
    two-sided p-value:           0.3222

Details:
    number of observations:   1000

pvalue() により，両側検定の $p$ 値が表示できる。

In [251]:
pvalue(result) # 両側検定
Out[251]:
0.32218030504370365
In [252]:
pvalue(result, tail=:left) # 左片側検定
Out[252]:
0.985395153969924
In [253]:
pvalue(result, tail=:right) # 右片側検定
Out[253]:
0.16177235356103092
  • 漸近近似検定
In [254]:
result2 = ApproximateOneSampleKSTest(randn(100), Normal()) # 一様分布に從うか
Out[254]:
Approximate one sample Kolmogorov-Smirnov test
----------------------------------------------
Population details:
    parameter of interest:   Supremum of CDF differences
    value under h_0:         0.0
    point estimate:          0.0944015

Test summary:
    outcome with 95% confidence: fail to reject h_0
    two-sided p-value:           0.3349

Details:
    number of observations:   100
    KS-statistic:             0.9440152502489629

pvalue() により，両側検定の 𝑝 値が表示できる。

In [255]:
pvalue(result2) # 両側検定
Out[255]:
0.3348891803417554
In [256]:
pvalue(result2, tail=:left) # 左片側検定
Out[256]:
0.16824574754375227
In [257]:
pvalue(result2, tail=:right) # 右片側検定
Out[257]:
0.9803720417635794

2 標本の分布の差の検定¶

2 標本の分布に差があるかどうかの検定である。

関数定義
ApproximateTwoSampleKSTest(x::AbstractVector{<:Real}, y::AbstractVector{<:Real})

ApproximateTwoSampleKSTest() は以下の戻り値を持つ。

■ 使用例

In [258]:
using Distributions
using Random; Random.seed!(456)
x = randn(100);         # 標準正規乱数
y = randn(100) .+ 0.15; # 母平均が 0.15
result = ApproximateTwoSampleKSTest(x, y)
Out[258]:
Approximate two sample Kolmogorov-Smirnov test
----------------------------------------------
Population details:
    parameter of interest:   Supremum of CDF differences
    value under h_0:         0.0
    point estimate:          0.2

Test summary:
    outcome with 95% confidence: reject h_0
    two-sided p-value:           0.0366

Details:
    number of observations:   [100,100]
    KS-statistic:              1.4142135623730954

■ 個別に表示できる統計量

  • p　$p$ 値
In [259]:
pvalue(result)
Out[259]:
0.03663105270711932
In [260]:
pvalue(result, tail=:left)
Out[260]:
0.9607894391523232
In [261]:
pvalue(result, tail=:right)
Out[261]:
0.018315638888734147

アンダーソン・ダーリング検定¶

1 標本の場合指定した分布に従うかどうかの検定，$k$ 標本の場合はそれらすべてが同じ分布に従うかどうかの検定である。

1 標本の場合¶

関数定義
OneSampleADTest(x::AbstractVector{<:Real}, d::UnivariateDistribution)

x は観察ベクトル
d は 1 変数分布関数
Normal(μ,σ), Uniform(a,b), Weibull(α,θ), Beta(α, β), Cauchy(μ, σ), Exponential(θ), Gamma(α,θ), Gumbel(μ, θ), Logistic(μ,θ), LogNormal(μ,σ) 等がある。　外部リンク

■

In [262]:
using Distributions
using Random; Random.seed!(123)
x = randn(1000);

result = OneSampleADTest(x, Normal()) # 標準正規分布に從うか
Out[262]:
One sample Anderson-Darling test
--------------------------------
Population details:
    parameter of interest:   not implemented yet
    value under h_0:         NaN
    point estimate:          NaN

Test summary:
    outcome with 95% confidence: fail to reject h_0
    one-sided p-value:           0.2543

Details:
    number of observations:   1000
    sample mean:              -0.050027717988504875
    sample SD:                0.9862398192221139
    A² statistic:             1.2358386455035488

■ 個別に表示できる統計量

  • p　$p$ 値
    "one-sided p-value:" と表示されているが，そもそもこの検定は両側/片側の区別はない。単に "p-value:" とすべきである。当然，tail=:left も tail=:right は言うに及ばず tail=:both さえも指定できない。
In [263]:
pvalue(result)
Out[263]:
0.2542809445499802
  • A²　検定統計量
In [264]:
result.A²
Out[264]:
1.2358386455035488

$k$ 標本の場合¶

関数定義
KSampleADTest(xs::AbstractVector{<:Real}...; modified = true, nsim = 0)

引数は，各標本の測定値ベクトルを列挙する。
$p$ 値の計算を標準正規分布により漸近検定で求めるとき modified=true にする。
シミュレーションにより $p$ 値を求めるときは，シミュレーションの回数を nsim で指定する。
nsim=0 のときはシミュレーションを行わず，漸近検定で $p$ 値を求める。

■ 使用例1

In [265]:
using Distributions
using Random; Random.seed!(456)
x = randn(100);         # 標準正規乱数
y = randn(100) .+ 0.15; # 母平均が 0.15
result = KSampleADTest(x, y)
Out[265]:
k-sample Anderson-Darling test
------------------------------
Population details:
    parameter of interest:   not implemented yet
    value under h_0:         NaN
    point estimate:          NaN

Test summary:
    outcome with 95% confidence: reject h_0
    one-sided p-value:           0.0040

Details:
    number of samples:        2
    number of observations:   200
    SD of A²k:                0.7541939849473112
    A²k statistic:            4.655692755994917
    standardized statistic:   4.847151832230945
    modified test:            true
    p-value calculation:      asymptotic

■ 個別に表示できる統計量

  • p　$p$ 値
    "one-sided p-value:" と表示されているが，そもそもこの検定は両側/片側の区別はない。単に "p-value:" とすべきである。当然，tail=:left も tail=:right は言うに及ばず tail=:both さえも指定できない。
In [266]:
pvalue(result)
Out[266]:
0.00403630676510518
  • A²k　検定統計量
In [267]:
result.A²k
Out[267]:
4.655692755994917

使用例2

In [268]:
result = KSampleADTest(x, y, modified=false)
Out[268]:
k-sample Anderson-Darling test
------------------------------
Population details:
    parameter of interest:   not implemented yet
    value under h_0:         NaN
    point estimate:          NaN

Test summary:
    outcome with 95% confidence: reject h_0
    one-sided p-value:           0.0041

Details:
    number of samples:        2
    number of observations:   200
    SD of A²k:                0.7541939849473112
    A²k statistic:            4.639424885412627
    standardized statistic:   4.825581956433769
    modified test:            false
    p-value calculation:      asymptotic

■ 使用例3

In [269]:
result = KSampleADTest(x, y, randn(100) .+ 1.0, nsim=10000)
Out[269]:
k-sample Anderson-Darling test
------------------------------
Population details:
    parameter of interest:   not implemented yet
    value under h_0:         NaN
    point estimate:          NaN

Test summary:
    outcome with 95% confidence: reject h_0
    one-sided p-value:           <1e-99

Details:
    number of samples:        3
    number of observations:   300
    SD of A²k:                1.0681421196124696
    A²k statistic:            38.555683139480266
    standardized statistic:   34.22361357002096
    modified test:            true
    p-value calculation:      simulation
    number of simulations:    10000

ワルド・ウォルフォビッツ連検定¶

関数定義
WaldWolfowitzTest(x::AbstractVector{Bool})
WaldWolfowitzTest(x::AbstractVector{<:Real})

引数が実数ベクトルの場合は，平均値より大きいか小さいかの論理ベクトルに変換して検定に用いられる。

■ 使用例

In [270]:
using HypothesisTests
using Distributions
using Random; Random.seed!(456)
x = rand(1:5, 100);
result = WaldWolfowitzTest(x)
Out[270]:
Wald-Wolfowitz Test
-------------------
Population details:
    parameter of interest:   Number of runs
    value under h_0:         48.58
    point estimate:          53

Test summary:
    outcome with 95% confidence: fail to reject h_0
    two-sided p-value:           0.3502

Details:
    number of runs:  53
    z-statistic:     0.9341742794043235

■ 個別に表示できる統計量

  • p　$p$ 値
In [271]:
pvalue(result)
Out[271]:
0.350214000986417
In [272]:
pvalue(result, tail=:left)
Out[272]:
0.8248929995067915
In [273]:
pvalue(result, tail=:right)
Out[273]:
0.1751070004932085
  • nabove, nbelow, nruns, μ, σ, z
    平均値より大きいデータの個数，平均値より小さいデータの個数，連の個数，標準正規分布で近似するときの連の平均と標準偏差，標準化得点
In [274]:
result.nabove, result.nbelow, result.nruns, result.μ, result.σ, result.z
Out[274]:
(61, 39, 53, 48.58, 4.731451183625411, 0.9341742794043235)
In [275]:
(result.nruns - result.μ) / result.σ # = result.z
Out[275]:
0.9341742794043235
In [276]:
using Distributions
2ccdf(Normal(), result.z) # = pvalur(result)
Out[276]:
0.350214000986417

並べ替え検定（無作為検定）¶

2 標本に対して，指定した関数を使って並べ替え検定（無作為検定）を行う。

関数定義

  • 正確検定
    ExactPermutationTest(x::Vector, y::Vector, f::Function)
    $n_x$, $n_y$ をサンプルサイズとしたとき，$(n_x + n_y)!$ 通りのデータの並べ替えに対して $f(x)$, $f(y)$ に差があるか検定する。

  • 近似検定
    ApproximatePermutationTest(x::Vector, y::Vector, f::Function, n::Int)
    $(n_x + n_y)!$ 通りのデータの並べ替えのなかから $n$ 通りだけ抽出して $f(x)$, $f(y)$ に差があるか検定する。

■ 使用例1　平均値の差を検定する

In [277]:
using Statistics
x = [2, 1, 3, 2, 5];
y = [4, 3, 2, 5, 4];
result = ExactPermutationTest(x, y, mean)
Out[277]:
Permutation Test
----------------
Population details:
    parameter of interest:   not implemented yet
    value under h_0:         NaN
    point estimate:          NaN

Test summary:
    outcome with 95% confidence: fail to reject h_0
    p-value:                     0.3730

Details:
    observation: -1.0
    samples: [-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]

■ 個別に表示できる統計量

  • p　$p$ 値
In [278]:
pvalue(result)
Out[278]:
0.373015873015873
In [279]:
pvalue(result, tail=:left)
Out[279]:
0.1865079365079365
In [280]:
pvalue(result, tail=:right)
Out[280]:
0.9126984126984127
  • observation　観察値
In [281]:
result.observation
Out[281]:
-1.0
  • samples　シミュレーションの値
In [282]:
first(result.samples, 5)
Out[282]:
5-element Vector{Float64}:
 -1.0
 -1.0
 -1.0
 -1.0
 -1.0

result.samples にシミュレーションの値が返されるので，ヒストグラムを描いてみる。

In [283]:
using Plots
histogram(result.samples)
Out[283]:

度数分布を求めてみる。

In [284]:
using FreqTables
freq = freqtable(result.samples)
Out[284]:
12-element Named Vector{Int64}
Dim1  │ 
──────┼───────
-2.2  │  28800
-1.8  │  72000
-1.4  │ 216000
-1.0  │ 360000
-0.6  │ 504000
-0.2  │ 633600
0.2   │ 633600
0.6   │ 504000
1.0   │ 360000
1.4   │ 216000
1.8   │  72000
2.2   │  28800

シミュレーションは (length(x) + length(y))! = 3628800 回行われる。

In [285]:
n = sum(freq)
Out[285]:
3628800

観察されたデータでの平均値の差は -1 なので，シミュレーションの結果が -1 以下になる場合の数を求める。

In [286]:
leftn = sum(freq[names(freq)[1] .<= -1])
Out[286]:
676800

(leftn / n) が片側検定の $p$ 値

In [287]:
leftn/n
Out[287]:
0.1865079365079365

2 倍すれば両側検定の $p$ 値

■ 使用例2　中央値の差を検定する

In [288]:
result = ExactPermutationTest(x, y, median)
Out[288]:
Permutation Test
----------------
Population details:
    parameter of interest:   not implemented yet
    value under h_0:         NaN
    point estimate:          NaN

Test summary:
    outcome with 95% confidence: fail to reject h_0
    p-value:                     0.3333

Details:
    observation: -2.0
    samples: [-2.0, -2.0, -2.0, -2.0, -2.0, -2.0, -2.0, -2.0, -2.0, -2.0  …  2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0]

■ 使用例3　全数探索は無理なので，10000 回のシミュレーションを行う。

In [289]:
result = ApproximatePermutationTest(x, y, mean, 10000)
Out[289]:
Permutation Test
----------------
Population details:
    parameter of interest:   not implemented yet
    value under h_0:         NaN
    point estimate:          NaN

Test summary:
    outcome with 95% confidence: fail to reject h_0
    p-value:                     0.3710

Details:
    observation: -1.0
    samples: [1.0, -0.2, 0.2, -0.6, 0.6, -1.0, -0.6, 0.6, 1.8, -1.0  …  -1.0, 0.6, -0.2, -0.2, 0.2, 0.2, -0.2, -1.4, -1.0, -0.2]
In [290]:
using Plots
histogram(result.samples)
Out[290]:
In [291]:
freqtable(result.samples)
Out[291]:
12-element Named Vector{Int64}
Dim1  │ 
──────┼─────
-2.2  │   79
-1.8  │  172
-1.4  │  609
-1.0  │ 1039
-0.6  │ 1432
-0.2  │ 1661
0.2   │ 1761
0.6   │ 1436
1.0   │  974
1.4   │  582
1.8   │  169
2.2   │   86
In [292]:
mean(abs.(result.samples) .>=  abs(result.observation)) # -1 以下，または 1 以上の割合が $p$ 値
Out[292]:
0.371