PowerDivergenceTest--(5) まとめ

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)
=#