Julia で統計解析 第 7 章 検定と推定¶
Version()
最新バージョン 2024-07-18 07:55
この章に示す統計検定を行う関数の詳細は以下を参照のこと。
外部リンク:HypothesisTests package
この章に示す統計検定を行う前に,以下を行っておく。
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
が出るので気づく)。
■ 使用例
result = PowerDivergenceTest([4, 2, 1, 3, 4, 2], lambda=1.0) # ピアソンの$\chi^2$検定
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]
■ 個別に表示できる統計量
- 概要
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(観察ベクトル)
とするだけで結果が表示できる。
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;
chisq_test([4, 2, 1, 3, 4, 2])
χ-sq. = 2.75, df = 5, p value = 0.7384611787603711
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$ 分布の右側にしかないが,本来この検定は意味的に両側検定である(両側検定,片側検定という区分分け自体が不要である)。
pvalue(result)
0.7384611787603711
pvalue(result, tail=:right)
0.7384611787603711
- 信頼区間 信頼区間は度数分布の確率の信頼区間である。
この検定では,tail=:left
,tail=:right
は使用してはいけない(pvalue()
の場合と違うので注意が必要である)。
confint(result)
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)
confint(result, tail=:both)
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$
result.lambda
1.0
n
サンプルサイズ
result.n
16
observed
観察度数分布
result.observed
6×1 Matrix{Int64}: 4 2 1 3 4 2
thetahat
観察度数の確率
result.observed / result.n
result.thetahat
6-element Vector{Float64}: 0.25 0.125 0.0625 0.1875 0.25 0.125
theta0
帰無仮説の下での観察度数の確率(等確率)
ones(length(result.observed)) / length(result.observed)
round.(result.theta0, digits=5)
6-element Vector{Float64}: 0.16667 0.16667 0.16667 0.16667 0.16667 0.16667
expected
帰無仮説の下での観察度数の期待値(同じ値)
result.n * result.theta0
round.(result.expected, digits=5)
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)
result.stat
2.75
df
$\chi^2$ 分布の自由度
length(result.observed) - 1
result.df
5
residuals
残差
(result.observed - result.expected) ./ sqrt.(result.expected)
round.(result.residuals, digits=5)
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)
round.(result.stdresiduals, digits=5)
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
を指定することで行うことができる。
■ 使用例
result = PowerDivergenceTest([4, 2, 1, 3, 4, 2], lambda=0.0)
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()
と同様に以下の関数を定義して簡潔な結果を表示する。
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
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$ 分布の右側にしかないが,本来この検定は意味的に両側検定である(両側検定,片側検定という区分分け自体が不要である)。
pvalue(result)
0.710619024390339
pvalue(result, tail=:right)
0.710619024390339
- 信頼区間 信頼区間は度数分布の確率の信頼区間である。
この検定では,tail=:left
,tail=:right
は使用してはいけない(pvalue()
の場合と違うので注意が必要である)。
confint(result)
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)
confint(result, tail=:both)
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$
result.lambda
0.0
n
サンプルサイズ
result.n
16
observed
観察度数分布
result.observed
6×1 Matrix{Int64}: 4 2 1 3 4 2
thetahat
観察度数の確率
result.observed / result.n
result.thetahat
6-element Vector{Float64}: 0.25 0.125 0.0625 0.1875 0.25 0.125
theta0
帰無仮説の下での観察度数の確率(等確率)
ones(length(result.observed)) / length(result.observed)
round.(result.theta0, digits=5)
6-element Vector{Float64}: 0.16667 0.16667 0.16667 0.16667 0.16667 0.16667
expected
帰無仮説の下での観察度数の期待値(同じ値)
result.n * result.theta0
round.(result.expected, digits=5)
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)
result.stat
2.931024858031232
df
$\chi^2$ 分布の自由度
length(result.observed) - 1
result.df
5
residuals
残差
(result.observed - result.expected) ./ sqrt.(result.expected)
round.(result.residuals, digits=5)
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)
round.(result.stdresiduals, digits=5)
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
を指定する。
■ 使用例
result = PowerDivergenceTest([21, 12, 5, 4], theta0=[9/16, 3/16, 3/16, 1/16], lambda=1.0)
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() を使うと以下のようになる。
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$ 分布の右側にしかないが,本来この検定は意味的に両側検定である(両側検定,片側検定という区分分け自体が不要である)。
pvalue(result)
0.23844641563476865
pvalue(result, tail=:right)
0.23844641563476865
- 信頼区間 信頼区間は度数分布の確率の信頼区間である。
この検定では,tail=:left
,tail=:right
は使用してはいけない(pvalue()
の場合と違うので注意が必要である)。
confint(result)
4-element Vector{Tuple{Float64, Float64}}: (0.35714285714285715, 0.6575992779891421) (0.14285714285714285, 0.4433135637034278) (0.0, 0.2766468970367611) (0.0, 0.2528373732272373)
confint(result, tail=:both)
4-element Vector{Tuple{Float64, Float64}}: (0.35714285714285715, 0.6575992779891421) (0.14285714285714285, 0.4433135637034278) (0.0, 0.2766468970367611) (0.0, 0.2528373732272373)
- lambda $\lambda$
result.lambda
1.0
n
サンプルサイズ
result.n
42
observed
観察度数分布
result.observed
4×1 Matrix{Int64}: 21 12 5 4
thetahat
観察度数の確率
result.observed / result.n
result.thetahat
4-element Vector{Float64}: 0.5 0.2857142857142857 0.11904761904761904 0.09523809523809523
theta0
帰無仮説の下での観察度数の確率(等確率)
ones(length(result.observed)) / length(result.observed)
round.(result.theta0, digits=5)
4-element Vector{Float64}: 0.5625 0.1875 0.1875 0.0625
expected
帰無仮説の下での観察度数の期待値(同じ値)
result.n * result.theta0
round.(result.expected, digits=5)
4×1 Matrix{Float64}: 23.625 7.875 7.875 2.625
stat
検定統計量($\chi^2$ 分布にしたがう)
sum((result.observed - result.expected).^2 ./ result.expected)
result.stat
4.22222222222222
df
$\chi^2$ 分布の自由度
length(result.observed) - 1
result.df
3
residuals
残差
(result.observed - result.expected) ./ sqrt.(result.expected)
round.(result.residuals, digits=5)
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)
round.(result.stdresiduals, digits=5)
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
を指定する。
■ 使用例
result = PowerDivergenceTest([21, 12, 5, 4], theta0=[9/16, 3/16, 3/16, 1/16], lambda=0.0)
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()
を使うと以下のようになる。
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$ 分布の右側にしかないが,本来この検定は意味的に両側検定である(両側検定,片側検定という区分分け自体が不要である)。
pvalue(result)
0.26261202803389555
pvalue(result, tail=:right)
0.26261202803389555
- 信頼区間 信頼区間は度数分布の確率の信頼区間である。
この検定では,tail=:left
,tail=:right
は使用してはいけない(pvalue()
の場合と違うので注意が必要である)。
confint(result)
4-element Vector{Tuple{Float64, Float64}}: (0.35714285714285715, 0.6575992779891421) (0.14285714285714285, 0.4433135637034278) (0.0, 0.2766468970367611) (0.0, 0.2528373732272373)
confint(result, tail=:both)
4-element Vector{Tuple{Float64, Float64}}: (0.35714285714285715, 0.6575992779891421) (0.14285714285714285, 0.4433135637034278) (0.0, 0.2766468970367611) (0.0, 0.2528373732272373)
- lambda $\lambda$
result.lambda
0.0
n
サンプルサイズ
result.n
42
observed
観察度数分布
result.observed
4×1 Matrix{Int64}: 21 12 5 4
thetahat
観察度数の確率
result.observed / result.n
result.thetahat
4-element Vector{Float64}: 0.5 0.2857142857142857 0.11904761904761904 0.09523809523809523
theta0
帰無仮説の下での観察度数の確率(等確率)
ones(length(result.observed)) / length(result.observed)
round.(result.theta0, digits=5)
4-element Vector{Float64}: 0.5625 0.1875 0.1875 0.0625
expected
帰無仮説の下での観察度数の期待値(同じ値)
result.n * result.theta0
round.(result.expected, digits=5)
4×1 Matrix{Float64}: 23.625 7.875 7.875 2.625
stat
検定統計量($\chi^2$ 分布にしたがう)
sum((result.observed - result.expected).^2 ./ result.expected)
result.stat
3.9893906620976454
df
$\chi^2$ 分布の自由度
length(result.observed) - 1
result.df
3
residuals
残差
(result.observed - result.expected) ./ sqrt.(result.expected)
round.(result.residuals, digits=5)
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)
round.(result.stdresiduals, digits=5)
4×1 Matrix{Float64}: -0.8165 1.63075 -1.13658 0.8765
独立性の検定¶
独立性の検定は,2 つのカテゴリー変数が独立であるか(関連があるか)を検定する。
ここでは,前節と同じピアソンの $\chi^2$ 検定 と,対数尤度比検定($G^2$ 検定)を提示する。
利用する関数は,前節と同じ PowerDivergenceTest()
である。引数として二重集計表を与える場合と,2 つのベクトルを与える場合に対応している。
しかし,2 つのベクトルを与える場合には完全対応していないので,二重集計表を与えて使用する。
x = [12 35
43 56]
2×2 Matrix{Int64}: 12 35 43 56
result = PowerDivergenceTest(x, lambda=1.0)
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()
もそのまま使える。
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$ 分布の右側にしかないが,本来この検定は意味的に両側検定である(両側検定,片側検定という区分分け自体が不要である)。
pvalue(result)
0.037005362019413034
pvalue(result, tail=:right)
0.037005362019413034
- イエーツの連続性の補正による $p$ 値
2×2 分割表の場合,PowerDivergenceTest()
で得られる検定統計量はイエーツの連続性の補正をしないものである。
イエーツの連続性の補正をする場合には,以下の関数を定義しておくとよい。
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()
の場合と違うので注意が必要である)。
confint(result)
4-element Vector{Tuple{Float64, Float64}}: (0.0, 0.17163425153174042) (0.2123287671232877, 0.38396301865502813) (0.15753424657534246, 0.3291684981070829) (0.3013698630136986, 0.47300411454543906)
confint(result, tail=:both)
4-element Vector{Tuple{Float64, Float64}}: (0.0, 0.17163425153174042) (0.2123287671232877, 0.38396301865502813) (0.15753424657534246, 0.3291684981070829) (0.3013698630136986, 0.47300411454543906)
- lambda $\lambda$
result.lambda
1.0
n
サンプルサイズ
result.n
146
observed
観察度数分布
result.observed
2×2 Matrix{Int64}: 12 35 43 56
thetahat
観察度数の確率
result.observed / result.n
result.thetahat
4-element Vector{Float64}: 0.0821917808219178 0.2945205479452055 0.23972602739726026 0.3835616438356164
theta0
帰無仮説の下での観察度数の確率(等確率)
ones(length(result.observed)) / length(result.observed)
round.(result.theta0, digits=5)
4-element Vector{Float64}: 0.12127 0.25544 0.20065 0.42264
expected
帰無仮説の下での観察度数の期待値(同じ値)
result.n * result.theta0
round.(result.expected, digits=5)
2×2 Matrix{Float64}: 17.7055 29.2945 37.2945 61.7055
stat
検定統計量($\chi^2$ 分布にしたがう)
sum((result.observed - result.expected).^2 ./ result.expected)
result.stat
4.350164943588549
df
$\chi^2$ 分布の自由度
length(result.observed) - 1
result.df
1
residuals
残差
(result.observed - result.expected) ./ sqrt.(result.expected)
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)
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$ 値を表示するために,以下の関数を用意しておくとよい。
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 文字列ベクトルを与えて検定する
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)
2×3 Named Matrix{Int64} Dim1 ╲ Dim2 │ Hi Lo Med ────────────┼────────────── No │ 1 5 5 Yes │ 3 4 2
- 概要
chisq_test(y)
χ-sq. = 2.2190155523488855, df = 2, p value = 0.32972121777778757
- 残差分析の $p$ 値
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$ 値が小さい)。
■ 使用例
x = [10 3; 4 12]
2×2 Matrix{Int64}: 10 3 4 12
result = PowerDivergenceTest(x, lambda=0.0)
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() もそのまま使える。
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$ 分布の右側にしかないが,本来この検定は意味的に両側検定である(両側検定,片側検定という区分分け自体が不要である)。
pvalue(result)
0.00435864516471826
pvalue(result, tail=:right)
0.00435864516471826
- 信頼区間 信頼区間は分割表の各セルの確率の信頼区間である。
この検定では,tail=:left
,tail=:right
は使用してはいけない(pvalue()
の場合と違うので注意が必要である)。
confint(result)
4-element Vector{Tuple{Float64, Float64}}: (0.1724137931034483, 0.5370409523965961) (0.0, 0.3301444006724581) (0.0, 0.2956616420517684) (0.24137931034482757, 0.6060064696379754)
confint(result, tail=:both)
4-element Vector{Tuple{Float64, Float64}}: (0.1724137931034483, 0.5370409523965961) (0.0, 0.3301444006724581) (0.0, 0.2956616420517684) (0.24137931034482757, 0.6060064696379754)
- lambda $\lambda$
result.lambda
0.0
n
サンプルサイズ
result.n
29
observed
観察度数分布
result.observed
2×2 Matrix{Int64}: 10 3 4 12
thetahat
観察度数の確率
result.observed / result.n
result.thetahat
4-element Vector{Float64}: 0.3448275862068966 0.13793103448275862 0.10344827586206896 0.41379310344827586
theta0
帰無仮説の下での観察度数の確率(等確率)
ones(length(result.observed)) / length(result.observed)
round.(result.theta0, digits=5)
4-element Vector{Float64}: 0.21641 0.26635 0.23187 0.28537
expected
帰無仮説の下での観察度数の期待値(同じ値)
result.n * result.theta0
round.(result.expected, digits=5)
2×2 Matrix{Float64}: 6.27586 6.72414 7.72414 8.27586
stat
検定統計量($\chi^2$ 分布にしたがう)
sum((result.observed - result.expected).^2 ./ result.expected)
result.stat
8.1280145470097
df
$\chi^2$ 分布の自由度
length(result.observed) - 1
result.df
1
residuals
残差。ピアソンの $\chi^2$ 検定と同じである
(result.observed - result.expected) ./ sqrt.(result.expected)
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)
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 つでもある場合には正しい結果が得られない。
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
では以下のようになる。
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
にした場合の,正しい答えを返す。
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
を指定すればよい。
trueG2_test(z)
G-sq. = 15.364591286599591, df = 6, p value = 0.01760288650305146
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) |
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="")
$\chi^2$ 分布で近似する場合には λ ≒ 2/3 のときに最も近似効率が高いとされる。
PowerDivergenceTest()
での lambda
は実数で指定しなければならない(lambda=1
ではなく,lambda=1.0
と指定する)。
lambda=1.0
でピアソンの $\chi^2$ 検定を行う。
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
chisq_test([10 3; 4 12]) # 上の例と同じである
χ-sq. = 7.743956043956044, df = 1, p value = 0.005389260671694263
lambda=0.0
で,対数尤度比に基づく $G^2$ 検定を行う。
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
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$
■ 使用例
result = FisherExactTest(10, 3, 4, 12)
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$ 値
pvalue(result) # デフォールトで両側検定の p 値
0.014589609220157725
pvalue(result, tail=:both) # 両側検定の p 値
0.014589609220157725
pvalue(result, tail=:left) # 片側検定の p 値
0.9994164940233701
pvalue(result, tail=:right) # 片側検定の p 値
0.007294804610078863
FisherExactTest()
の両側検定の $p$ 値は該当する片側確率を 2 倍するものである。
上の例では,0.014589609220157725 = 2 × 0.007294804610078863 になっている。
pvalue(result, method=:central) # 両側検定の p 値(デフォルト)
0.014589609220157725
このような計算方法は R や matlab などの他の統計ソフトと違う。同じ結果を得るためには,オプション設定 method=:minlike
としなければならない。
pvalue(result, method=:minlike)
0.009220570313398513
using RCall # R で計算すると以下のようになる
R"""
options(digits=15)
x = matrix(c(10, 3, 4, 12), 2)
fisher.test(x)$p.value
"""
RObject{RealSxp} [1] 0.00922057031339851
- ω オッズ比
result.ω
9.067840887689634
confint
オッズ比の信頼区間
confint(result)
(1.4222662898175877, 80.34790877482563)
信頼率は,level
で指定できる(0 < level < 1
)
confint(result, level=0.95)
(1.4222662898175877, 80.34790877482563)
■ オッズ比の $p$ 値と信頼区間について
オッズ比($OR$)の $p$ 値,信頼区間には何通りかの方法がある。
FisherExactTest()
は,R の fisher.test()
と同じく,条件付き最尤推定量を求める。 外部リンク:Fisherの正確検定
一方,広く使われているのは条件なしの最尤推定量(Wald のオッズ比)で,以下の関数で計算できる。 外部リンク:統計学入門-第3章
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;
oddsratio(10, 3, 4, 12)
オッズ比 = 10.0 95% 信頼区間 [1.7975962537192371, 55.62984446206951]
R での実行結果を添付しておく。
using RCall
R"fisher.test(matrix(c(10, 3, 4, 12), 2))"
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)` とも書ける。
■ 使用例
result = BinomialTest(12, 20)
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$ 値が表示される。
pvalue(result, tail=:left)
0.8684120178222656
pvalue(result, tail=:right)
0.2517223358154298
pvalue(result, tail=:both)
0.5034446716308596
なお,両側検定の $p$ 値は該当する片側検定の $p$ 値を 2 倍したものであり,R や matlab などの他の統計ソフトとは違う。
計算方法により結果が違うのは,「母比率が 0.5 以外の両側検定」の場合に限られる。
以下の例,x=18, n=24, p=0.68 の場合,$p$ 値は 0.6206975448201958 が表示されるが,R では 0.5205908910845283 となる。
result2 = BinomialTest(18, 24, 0.68); # p=0.68 の場合。数値だけ指定するので注意。
pvalue(result2, tail=:both)
0.6206975448201958
using RCall
R"""
options(digits=16)
binom.test(18, 24, 0.68)$p.value
"""
RObject{RealSxp} [1] 0.5205908910845283
BinomialTest()
では,$n = 18,..., 24$ の二項分布の確率の合計が 0.3103487724100979 で,
pvalue(result, tail=:right)
0.2517223358154298
その 2 倍が両側 $p$ 値になっている。
2 * pvalue(result, tail=:right)
0.5034446716308596
R などでは,まず $n = 18,...,24$ を右側の棄却域とするところまでは同じであるが,もう一方の側の棄却域を次のようにして決める。
$x = 0, 1, ...$ と順に見ていき,観察値 18 のときの確率より小さい最大の $n$ まで,図でいえば $x = 0, ..., 14$ までを棄却域とする。
それぞれの棄却域における確率の合計を両側確率とする。
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="")
sum(p_value[1:15]) + sum(p_value[19:25]) # Julia の添字は 1 から始まるので p_value[n+1] が n のときの確率である
0.5205908910845278
FisherTest()
の場合と異なり,同じ結果を得るためのオプション設定がないので,以下のようなプログラムを使用すればよい。
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
binomial_pvalue (generic function with 2 methods)
binomial_pvalue(12, 20)
0.5034446716308596
binomial_pvalue(4, 20, 0.6, tail=:left)
0.00031703112116863026
binomial_pvalue(16, 20, 0.6, tail=:right)
0.05095195319416655
- 信頼区間の表示
result = BinomialTest(5, 26)
confint(result)
(0.06554810873678252, 0.3935055279393218)
confint(result, level=0.99)
(0.04400829697421974, 0.4550043049581578)
confint(result, tail=:right)
(0.07898571028519472, 1.0)
confint(result, tail=:left)
(0.0, 0.36259486197815627)
$t$ 検定¶
いわゆる $t$ 検定は,一標本の場合の母平均の検定と独立二標本の平均値の差の検定に大別される。更に,二群の平均値の差の検定は,二群の分散が等しいと仮定する場合と仮定しない場合に分かれる。
なお,対応のある二標本の平均値の差の検定もあるが,これは対応のあるデータの差をとって一標本の場合の母平均の検定をするのと同じである。
これらの検定は OneSampleTTest()
,EqualvarianceTTest()
,UnequalvarianceTTest()
により行う。
一標本の検定(母平均の検定)¶
関数定義
OneSampleTTest(v::AbstractVector{T<:Real}, μ0::Real = 0)
対応のある二標本の場合,対応のあるデータの差をとって一標本の検定を行うのと同じである。
■ 使用例
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;
using HypothesisTests
result = OneSampleTTest(x, y)
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
result = OneSampleTTest(z)
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$ 値
pvalue(result)
0.4453534218554348
pvalue(result, tail=:both)
0.4453534218554348
pvalue(result, tail=:left)
0.7773232890722825
pvalue(result, tail=:right)
0.2226767109277174
- 信頼区間 信頼区間は母平均の信頼区間である。
confint(result)
(-1.6142612308053204, 3.3742612308053195)
confint(result, tail=:both)
(-1.6142612308053204, 3.3742612308053195)
confint(result, level=0.99)
(-2.7032783553640023, 4.463278355364001)
confint(result, tail=:right)
(-1.141195783743237, Inf)
confint(result, tail=:right, level=0.99)
(-2.2309258663578753, Inf)
confint(result, tail=:left)
(-Inf, 2.901195783743236)
n
サンプルサイズ
result.n
10
result.xbar
0.8799999999999997
t
検定統計量($t$ 分布にしたがう)
result.t
0.7981113921334072
- df $t$ 分布の自由度
result.df
9
- xbar 標本平均
xbar = mean(x) - mean(y) = mean(x - y) = mean(z)
result.xbar
0.8799999999999997
- μ0 母平均
result.μ0
0
- stderr 標準誤差
result.t = (result.xbar - result.μ0) / result.stderr
result.stderr
1.1026029808291529
等分散の場合の $t$ 検定¶
独立 2 標本の分散が等しいことを仮定する平均値の差の検定である。
一般的には,等分散を仮定しない $t$ 検定(Welch の方法)が推奨される。
関数定義
EqualVarianceTTest(x::AbstractVector{T<:Real}, y::AbstractVector{T<:Real})
■ 使用例
x = [1, 2, 3, 6, 3, 1, 2, 1, 1, 1, 3, 4];
y = [2, 1, 2, 3, 4, 2, 1, 2, 3, 5];
result = EqualVarianceTTest(x, y)
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$ 値
pvalue(result) # tail=:both 両側検定がデフォルト
0.788931157688499
pvalue(result, tail=:left) # 棄却域が左にある片側検定
0.3944655788442495
pvalue(result, tail=:right) # 棄却域が右にある片側検定
0.6055344211557505
- 信頼区間
confint(result) #デフォルトで両側に棄却域がある場合の信頼率 95% の信頼区間
(-1.448068275504171, 1.114734942170838)
confint(result, level=0.99) # 信頼率 99%
(-1.9145510251333067, 1.5812176917999736)
confint(result, tail=:left, level=0.99) # 片側と信頼率の指定
(-Inf, 1.386262653673622)
t
検定統計量 $t$
result.t
-0.271312734545191
df
$t$ 分布の自由度
result.df
20
n_x
,n_y
各群のサンプルサイズ
result.n_x, result.n_y
(12, 10)
xbar
,μ0
平均値の差と母平均
xbar = mean(x) - mean(y)
result.xbar, result.μ0
(-0.16666666666666652, 0)
- stderr 標準偏差
t = (result.xbar - result.μ0) / result.stderr
result.stderr
0.6142972497994038
等分散でない場合 Welch の方法¶
関数定義
UnequalVarianceTTest(x::AbstractVector{T<:Real}, y::AbstractVector{T<:Real})
■ 使用例
x = [1, 2, 3, 6, 3, 1, 2, 1, 1, 1, 3, 4];
y = [2, 1, 2, 3, 4, 2, 1, 2, 3, 5];
result = UnequalVarianceTTest(x, y)
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$ 値
pvalue(result) # tail=:both 両側検定がデフォルト
0.784943053185432
pvalue(result, tail=:left) # 棄却域が左にある片側検定
0.392471526592716
pvalue(result, tail=:right) # 棄却域が右にある片側検定
0.607528473407284
- 信頼区間
confint(result) #デフォルトで両側に棄却域がある場合の信頼率 95% の信頼区間
(-1.423689159406056, 1.0903558260727229)
confint(result, level=0.99) # 信頼率 99%
(-1.8813078766019782, 1.5479745432686451)
confint(result, tail=:left, level=0.99) # 片側と信頼率の指定
(-Inf, 1.3567230542137687)
t
検定統計量 $t$
result.t
-0.2765775336645321
df
$t$ 分布の自由度
result.df
19.99676443938765
n_x
,n_y
各群のサンプルサイズ
result.n_x, result.n_y
(12, 10)
xbar
,μ0
平均値の差と母平均
xbar = mean(x) - mean(y)
result.xbar, result.μ0
(-0.16666666666666652, 0)
- stderr 標準偏差
t = (result.xbar - result.μ0) / result.stderr
result.stderr
0.6026037778933294
result.t
-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()
が適用される。
- 同順位データ(tie)がなく,プールされたサンプルサイズが 50 以下の場合
- 同順位があるが,プールされたサンプルサイズが 10 以下の場合
それ以外の場合には漸近近似検定 ApproximateMannWhitneyUTest()
が適用される。
しかし,ExactMannWhitneyUTest()
はデータの並べ替えのすべてを評価するので,程々のサンプルサイズの場合でも計算時間がかかりすぎることがある。
そのため,条件を満たすかどうかに限らず ExactMannWhitneyUTest()
か ApproximateMannWhitneyUTest()
を直接選択することもできる。
■ 使用例1
以下のデータは,プールされたサンプルサイズが 22 であるが同順位があるので,条件 1 に反する。また,同順位があるがプールされたサンプルサイズが 11 以上なので,条件 2 にも反する。
したがって,何も指示しない場合は ApproximateMannWhitneyUTest(x, y)
が適用される。
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)
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() と同じ結果になる。
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/dDAVd/src/io.jl:172
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() を適用した結果を以下に示す。
exact_result = ExactMannWhitneyUTest(x, y)
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
パッケージを使った正確な検定とも違う。
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/dDAVd/src/io.jl:172
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 以下のサンプルサイズデータでも確認してみよう。
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)
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)
exact_result = ExactMannWhitneyUTest(x, y)
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
pvalue(approx_result) |> println
pvalue(exact_result) |> println
0.03211417930790749 0.029960751322980424
同順位がまったくないデータの場合には,R での結果と一致する。
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")
"""
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$ 値を求めることができる。
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
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
は帰無仮説下の中央値。
■ 使用例
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)
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
result2 = SignTest([1, 0, 1, 1, -1, 1, 1, -1, 1, 1])
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$ 値が同じになる)。
using HypothesisTests
result3 = BinomialTest(7, 9)
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
pvalue(result) |> println
pvalue(result2) |> println
pvalue(result3) |> println
0.17968750000000006 0.17968750000000006 0.17968750000000006
■ 個別に表示できる統計量
p
$p$ 値
pvalue(result) # デフォルトで両側検定
0.17968750000000006
pvalue(result, tail=:left) # 棄却域が左にある場合の片側検定
0.98046875
pvalue(result, tail=:right) # 棄却域が右にある場合の片側検定
0.08984375000000003
- confint 信頼区間
confint(result) # デフォルトで両側に棄却域がある場合で,信頼率が 95% の信頼区間
(0, 2)
confint(result, level=0.99) # 信頼率の指定
(-1, 2)
confint(result, tail=:left, level=0.99) # 棄却域が左にある場合の信頼区間
(-1, 0.0)
confint(result, tail=:right, level=0.99) # 棄却域が右にある場合の信頼区間
(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()
が実行される。
- 同順位がなくて,サンプルサイズが 50 以下である。
- 同順位があるが,サンプルサイズが 15 以下である。
どちらの条件も満たさないときは,ApproximateSignedRankTest()
が実行される。
条件を満たすかどうかに限らず,ApproximateSignedRankTest()
か ExactSignedRankTest()
を直接選択できる。
■ 使用例
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)
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
z = x .- y
result2 = SignedRankTest(z)
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$ 値
pvalue(result)
0.15625
pvalue(result, tail=:left)
0.9609375
pvalue(result, tail=:right)
0.078125
confint
信頼区間
confint(result)
(0.0, 2.0)
confint(result, level=0.99)
(-1.0, 2.0)
confint(result, tail=:left, level=0.99)
(-0.5, Inf)
confint(result, tail=:right, level=0.99)
(-Inf, 2.0)
W
検定統計量
result3.W
■ 漸近近似検定の場合
using Random; Random.seed!(123)
x = randn(100)
y = randn(100) .+ 0.1
result3 = ApproximateSignedRankTest(x, y)
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()
の場合のみ
result3.mu, result3.sigma
(-163.0, 290.8393027085576)
using Distributions
z = result3.mu / result3.sigma
-0.5604469495078456
2*cdf(Normal(), z)
0.5751746153622865
pvalue(result3) # 一致しない?
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" |
ケンドールの順位相関係数 |
スピアマンの順位相関係数,ケンドールの順位相関係数においては,同順位がある場合には近似検定になる。
ケンドールの順位相関係数においては,信頼区間は求めない。
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 変数をベクトルで与える場合
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]
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]
CorTest(x, y, method="kendall")
n = 15, r = 0.31429, z = 1.63308, p-value = 0.10245
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]
■ サンプルサイズと相関係数を与える場合
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]
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]
一元元配置分散分析¶
HypothesisTests の OneWayANOVATest()
には,2022/05/21 現在でバグがある。
関数定義
OneWayANOVATest(groups::AbstractVector{<:Real}...)
引数として,各標本の観測データベクトルを列挙する。
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)
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)
正しいプログラム例を挙げておく。
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
EqualVarianceOneWayANOVATest (generic function with 1 method)
■ 使用例
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}...)
の場合と同じく,引数として各標本の観測データベクトルを列挙する。
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;
■ 使用例
UnequalVarianceOneWayANOVATest(g1, g2, g3)
F = 0.51339, df1 = 2, df2 = 13.589, p-value = 0.60962
クラスカル・ウォリス検定¶
対応のない 2 標本においけるマン・ホイットニーの $U$ 検定の拡張である。
関数定義
KruskalWallisTest(groups::AbstractVector{<:Real}...)
引数として各標本の観測データベクトルを列挙する。
■ 使用例
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)
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$ 分布にしたがう)
result.H
0.6899999999999977
df
$p$ 値を求めるとき,$\chi^2$ 分布で近似するときの自由度
result.df
2
p
$p$ 値
"one-sided p-value:" と表示されているが,そもそもこの検定は両側/片側の区別はない。単に "p-value:" とすべきである。当然,tail=:left
もtail=:right
は言うに及ばずtail=:both
さえも指定できない。
pvalue(result)
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(μ,σ) 等がある。 外部リンク
■ 使用例
using HypothesisTests
using Distributions
using Random; Random.seed!(123)
x = randn(1000);
- 正確な検定
result = ExactOneSampleKSTest(x, Normal()) # 標準正規分布に從うか
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$ 値が表示できる。
pvalue(result) # 両側検定
0.32218030504370365
pvalue(result, tail=:left) # 左片側検定
0.985395153969924
pvalue(result, tail=:right) # 右片側検定
0.16177235356103092
- 漸近近似検定
result2 = ApproximateOneSampleKSTest(randn(100), Normal()) # 一様分布に從うか
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() により,両側検定の 𝑝 値が表示できる。
pvalue(result2) # 両側検定
0.3348891803417554
pvalue(result2, tail=:left) # 左片側検定
0.16824574754375227
pvalue(result2, tail=:right) # 右片側検定
0.980499551913982
2 標本の分布の差の検定¶
2 標本の分布に差があるかどうかの検定である。
関数定義
ApproximateTwoSampleKSTest(x::AbstractVector{<:Real}, y::AbstractVector{<:Real})
ApproximateTwoSampleKSTest()
は以下の戻り値を持つ。
■ 使用例
using Distributions
using Random; Random.seed!(456)
x = randn(100); # 標準正規乱数
y = randn(100) .+ 0.15; # 母平均が 0.15
result = ApproximateTwoSampleKSTest(x, y)
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$ 値
pvalue(result)
0.03663105270711932
pvalue(result, tail=:left)
0.9607894391523232
pvalue(result, tail=:right)
0.018315638888734147
アンダーソン・ダーリング検定¶
1 標本の場合指定した分布に従うかどうかの検定,$k$ 標本の場合はそれらすべてが同じ分布に従うかどうかの検定である。
using Distributions
using Random; Random.seed!(123)
x = randn(1000);
result = OneSampleADTest(x, Normal()) # 標準正規分布に從うか
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.05002771798850483 sample SD: 0.9862398192221139 A² statistic: 1.2358386455035488
■ 個別に表示できる統計量
p
$p$ 値
"one-sided p-value:" と表示されているが,そもそもこの検定は両側/片側の区別はない。単に "p-value:" とすべきである。当然,tail=:left
もtail=:right
は言うに及ばずtail=:both
さえも指定できない。
pvalue(result)
0.2542809445499802
- A² 検定統計量
result.A²
1.2358386455035488
$k$ 標本の場合¶
関数定義
KSampleADTest(xs::AbstractVector{<:Real}...; modified = true, nsim = 0)
引数は,各標本の測定値ベクトルを列挙する。
$p$ 値の計算を標準正規分布により漸近検定で求めるとき modified=true
にする。
シミュレーションにより $p$ 値を求めるときは,シミュレーションの回数を nsim
で指定する。
nsim=0
のときはシミュレーションを行わず,漸近検定で $p$ 値を求める。
■ 使用例1
using Distributions
using Random; Random.seed!(456)
x = randn(100); # 標準正規乱数
y = randn(100) .+ 0.15; # 母平均が 0.15
result = KSampleADTest(x, y)
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
さえも指定できない。
pvalue(result)
0.00403630676510518
A²k
検定統計量
result.A²k
4.655692755994917
使用例2
result = KSampleADTest(x, y, modified=false)
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
result = KSampleADTest(x, y, randn(100) .+ 1.0, nsim=10000)
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})
引数が実数ベクトルの場合は,平均値より大きいか小さいかの論理ベクトルに変換して検定に用いられる。
■ 使用例
using HypothesisTests
using Distributions
using Random; Random.seed!(456)
x = rand(1:5, 100);
result = WaldWolfowitzTest(x)
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$ 値
pvalue(result)
0.350214000986417
pvalue(result, tail=:left)
0.8248929995067915
pvalue(result, tail=:right)
0.1751070004932085
nabove
,nbelow
,nruns
,μ
,σ
,z
平均値より大きいデータの個数,平均値より小さいデータの個数,連の個数,標準正規分布で近似するときの連の平均と標準偏差,標準化得点
result.nabove, result.nbelow, result.nruns, result.μ, result.σ, result.z
(61, 39, 53, 48.58, 4.731451183625411, 0.9341742794043235)
(result.nruns - result.μ) / result.σ # = result.z
0.9341742794043235
using Distributions
2ccdf(Normal(), result.z) # = pvalur(result)
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 平均値の差を検定する
using Statistics
x = [2, 1, 3, 2, 5];
y = [4, 3, 2, 5, 4];
result = ExactPermutationTest(x, y, mean)
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$ 値
pvalue(result)
0.373015873015873
pvalue(result, tail=:left)
0.1865079365079365
pvalue(result, tail=:right)
0.9126984126984127
observation
観察値
result.observation
-1.0
samples
シミュレーションの値
first(result.samples, 5)
5-element Vector{Float64}: -1.0 -1.0 -1.0 -1.0 -1.0
result.samples
にシミュレーションの値が返されるので,ヒストグラムを描いてみる。
using Plots
histogram(result.samples)
度数分布を求めてみる。
using FreqTables
freq = freqtable(result.samples)
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
回行われる。
n = sum(freq)
3628800
観察されたデータでの平均値の差は -1 なので,シミュレーションの結果が -1 以下になる場合の数を求める。
leftn = sum(freq[names(freq)[1] .<= -1])
676800
(leftn / n)
が片側検定の $p$ 値
leftn/n
0.1865079365079365
2 倍すれば両側検定の $p$ 値
■ 使用例2 中央値の差を検定する
result = ExactPermutationTest(x, y, median)
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 回のシミュレーションを行う。
result = ApproximatePermutationTest(x, y, mean, 10000)
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.3749 Details: observation: -1.0 samples: [-1.4, -0.6, 0.6, -1.4, -0.2, 0.6, 1.0, 0.6, 0.2, -1.0 … 0.6, -1.0, -1.4, 1.0, 0.2, 1.0, -1.0, -1.4, -0.2, -1.0]
using Plots
histogram(result.samples)
freqtable(result.samples)
12-element Named Vector{Int64} Dim1 │ ──────┼───── -2.2 │ 82 -1.8 │ 178 -1.4 │ 616 -1.0 │ 1016 -0.6 │ 1327 -0.2 │ 1764 0.2 │ 1769 0.6 │ 1391 1.0 │ 1014 1.4 │ 567 1.8 │ 206 2.2 │ 70
mean(abs.(result.samples) .>= abs(result.observation)) # -1 以下,または 1 以上の割合が $p$ 値
0.3749