using HypothesisTests
PowerDivergenceTest() が返すもの
x = [16 12 36
12 5 20
15 11 24
9 2 1];
obj = PowerDivergenceTest(x);
REPL の場合や println(obj) では，検定結果のまとめが出力される。
指定された lambda の値により，検定の名前が 1 行目に表示される。
#=
Pearson's Chi-square Test
-------------------------
Population details:
parameter of interest: Multinomial Probabilities
value under h_0: [0.125259, 0.0724152, 0.0978584, 0.023486, 0.0722647, 0.041778, 0.0564568, 0.0135496, 0.195115, 0.112801, 0.152433, 0.036584]
point estimate: [0.0981595, 0.0736196, 0.0920245, 0.0552147, 0.0736196, 0.0306748, 0.0674847, 0.0122699, 0.220859, 0.122699, 0.147239, 0.00613497]
95% confidence interval: [(0.0307, 0.1668), (0.0061, 0.1423), (0.0245, 0.1607), (0.0, 0.1239), (0.0061, 0.1423), (0.0, 0.0994), (0.0, 0.1362), (0.0, 0.0809), (0.1534, 0.2895), (0.0552, 0.1914), (0.0798, 0.2159), (0.0, 0.0748)]
Test summary:
outcome with 95% confidence: reject h_0
one-sided p-value: 0.0330
Details:
Sample size: 163
statistic: 13.713423881146099
degrees of freedom: 6
residuals: [-0.977568, 0.0571418, -0.238096, 2.64327, 0.0643514, -0.693532, 0.592554, -0.140357, 0.744097, 0.376286, -0.169847, -2.03246]
std. residuals: [-1.52004, 0.078758, -0.346528, 3.32797, 0.0914118, -0.873259, 0.787863, -0.161439, 1.34615, 0.603412, -0.287607, -2.97724]
=#
PowerDivergenceTest(x) の結果は，PowerDivergenceTest という型を持つ。
typeof(obj) # PowerDivergenceTest
'? obj' で，その内容を見ることができる。
obj は以下の名前の要素を持つ。
#=
lambda :: Float64
theta0 :: Vector{Float64}
stat :: Float64
df :: Int64
observed :: Matrix{Int64}
n :: Int64
thetahat :: Vector{Float64}
expected :: Matrix{Float64}
residuals :: Matrix{Float64}
stdresiduals :: Matrix{Float64}
=#
obj.xxx で，その要素を表示することができる。
使用された λ の値; obj.lambda
obj.lambda # 1.0
帰無仮説の下での各セルの確率; 'value under h_0'; obj.theta0
obj.theta0 # 要素数 12 のベクトル
#=
12-element Vector{Float64}:
0.12525876020926643
0.07241522074598215
0.0978584064134894
0.023486017539237458
0.07226466935149986
0.04177801196883586
0.056456772930859274
0.013549625503406226
0.19511460724904964
0.11280063231585684
0.15243328691332003
0.03658398885919681
=#
行列に直して表示すると以下のようになる。
pij = reshape(obj.theta0, size(x))
#=
4×3 Matrix{Float64}:
0.125259 0.0722647 0.195115
0.0724152 0.041778 0.112801
0.0978584 0.0564568 0.152433
0.023486 0.0135496 0.036584
=#
これは，行和，列和，サンプルサイズから計算される。
行和
rowsum = sum(x, dims=2)
#=
4×1 Matrix{Int64}:
64
37
50
12
=#
列和
colsum = sum(x, dims=1)
#=
1×3 Matrix{Int64}:
52 30 81
=#
サンプルサイズ; obj.n
obj.n # 163
n = sum(x) # 163
帰無仮説（二変数は独立）により，各セルの確率は，周辺確率の外積である。
pij = (rowsum ./ n) * (colsum ./ n)
#=
4×3 Matrix{Float64}:
0.125259 0.0722647 0.195115
0.0724152 0.041778 0.112801
0.0978584 0.0564568 0.152433
0.023486 0.0135496 0.036584
=#
期待値; obj.expected
obj.expected
#=
4×3 Matrix{Float64}:
20.4172 11.7791 31.8037
11.8037 6.80982 18.3865
15.9509 9.20245 24.8466
3.82822 2.20859 5.96319
=#
期待値はサンプルサイズと各セルの確率を掛けたものである。
eij = n .* pij
#=
4×3 Matrix{Float64}:
20.4172 11.7791 31.8037
11.8037 6.80982 18.3865
15.9509 9.20245 24.8466
3.82822 2.20859 5.96319
=#
各セルの確率; 'point estimate'; obj.thetahat
観察数をサンプルサイズで割ったものである。
reshape(obj.thetahat, size(x))
#=
4×3 Matrix{Float64}:
0.0981595 0.0736196 0.220859
0.0736196 0.0306748 0.122699
0.0920245 0.0674847 0.147239
0.0552147 0.0122699 0.00613497
=#
x ./ n
#=
4×3 Matrix{Float64}:
0.0981595 0.0736196 0.220859
0.0736196 0.0306748 0.122699
0.0920245 0.0674847 0.147239
0.0552147 0.0122699 0.00613497
=#
検定統計量; 'statistic'; obj.stat
obj.stat # 13.713423881146099
sum((x .- eij) .^ 2 ./ eij) # 13.713423881146104
自由度; 'degrees of freedom'; obj.df
obj.df # 6
nrow, ncol = size(x) # (4, 3)
(nrow - 1) * (ncol - 1) # 6
残差; 'residuals'; obj.residuals
obj.residuals
#=
4×3 Matrix{Float64}:
-0.977568 0.0643514 0.744097
0.0571418 -0.693532 0.376286
-0.238096 0.592554 -0.169847
2.64327 -0.140357 -2.03246
=#
これを，「標準化残差」とする場合もある。
residuals = (x .- eij) ./ sqrt.(eij)
#=
4×3 Matrix{Float64}:
-0.977568 0.0643514 0.744097
0.0571418 -0.693532 0.376286
-0.238096 0.592554 -0.169847
2.64327 -0.140357 -2.03246
=#
残差の二乗和は検定統計量である。sum((x .- eij) .^ 2 ./ eij) と同じ
sum(residuals .^ 2) # 13.713423881146102
標準化残差（調整された残差）; 'std. residuals'; obj.stdresiduals
obj.stdresiduals
#=
4×3 Matrix{Float64}:
-1.52004 0.0914118 1.34615
0.078758 -0.873259 0.603412
-0.346528 0.787863 -0.287607
3.32797 -0.161439 -2.97724
=#
PowerDivergenceTest では，以下のようにして標準化残差を計算する。
V = [(colsum[j]/n) * (rowsum[i]/n) * (1 - rowsum[i]/n) * (n - colsum[j]) for i in 1:nrow, j in 1:ncol]
#=
4×3 Matrix{Float64}:
8.44459 5.83748 9.71743
6.21349 4.29519 7.15004
7.53029 5.20545 8.66532
2.41503 1.66943 2.77904
=#
stdresiduals = (x .- eij) ./ sqrt.(V)
#=
4×3 Matrix{Float64}:
-1.52004 0.0914118 1.34615
0.078758 -0.873259 0.603412
-0.346528 0.787863 -0.287607
3.32797 -0.161439 -2.97724
=#
残差のことを「標準化残差」と呼ぶ場合には，PowerDivergenceTest での標準化残差は，「調整された残差」と呼ぶ。
調整された残差; adjusted residuals
この「調整された残差」と PowerDivergenceTest での標準化残差は同じものである。
adjresiduals = residuals ./ sqrt.((1 .- rowsum/n) * (1 .- colsum/n))
#=
4×3 Matrix{Float64}:
-1.52004 0.0914118 1.34615
0.078758 -0.873259 0.603412
-0.346528 0.787863 -0.287607
3.32797 -0.161439 -2.97724
=#
それぞれは標準正規分布に従うので，両側確率を求める。
using Rmath
2pnorm.(abs.(adjresiduals), false)
#=
4×3 Matrix{Float64}:
0.1285 0.927165 0.178255
0.937225 0.382522 0.546235
0.728946 0.430777 0.773648
0.00087482 0.871748 0.00290857
=#
p 値は obj に含まれないが，pvalue(obj) によって求める。
p = pvalue(obj) # 0.03300644487225656
なお，検定結果の表示においては 'one-sided p-value' とされているが，これは単に 'p-value' とすべきものである。
他に confint(obj) により各セルの確率の信頼限界を求めることができる。'point estimate' の '95% confidence interval'
reshape(confint(obj), size(x))
#=
4×3 Matrix{Tuple{Float64, Float64}}:
(0.0306748, 0.166838) (0.00613497, 0.142298) (0.153374, 0.289537)
(0.00613497, 0.142298) (0.0, 0.0993534) (0.0552147, 0.191378)
(0.0245399, 0.160703) (0.0, 0.136163) (0.0797546, 0.215918)
(0.0, 0.123893) (0.0, 0.0809485) (0.0, 0.0748135)
=#
