HypothesisTests の OneWayANOVATest, LeveneTest，BrownForsytheTest，FlignerKilleenTest

# HypothesisTests に含まれる検定関数の不具合

OneWayANOVATest, LeveneTest，BrownForsytheTest，FlignerKilleenTest それぞれ，あまり関連がなさそうに見えるが，Julia の HypothesisTests ではそれらを統合して扱っている。

のは，良かったのだけど，それぞれに共通して使用している関数 `anova()` にバグがあるために，総崩れになってしまっている。

2021/08/02 に指摘があったのだが，2022/05/15 至っても，バグの修正はない。無責任極まりないと思うのは，私だけだろうか?

なので，とりあえず，修正版をここにアップしておく。

まずは，テストデータ。


```julia
x1 = [301, 311, 325, 291, 388, 412, 325, 361, 287];
x2 = [197, 180, 247, 260, 247, 199, 179, 134, 163, 200];
x3 = [209, 302, 187, 166, 234, 290, 175, 116];
x4 = [342, 216, 316, 386, 324, 145, 254, 228];

# R での分析のために
vx = vcat(x1, x2, x3, x4...);
fx = vcat([repeat([cat], inner=times) for (cat, times) in zip(1:4, [9, 10, 8, 8])]...);
```

## 1. バグのある状態での各関数の動作状況

### 1.1. OneWayANOVATest


```julia
using HypothesisTests
OneWayANOVATest(x1, x2, x3, x4)
```




    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: reject h_0
        p-value:                     <1e-04
    
    Details:
        number of observations: [9, 10, 8, 8]
        F statistic:            10.6266
        degrees of freedom:     (3, 31)





```julia
using RCall
R"""
oneway.test($vx ~ $fx, var.equal=TRUE)
"""
```




    RObject{VecSxp}
    
    	One-way analysis of means
    
    data:  `#JL`$vx and `#JL`$fx
    F = 10.624, num df = 3, denom df = 31, p-value = 5.791e-05
    




R では正しく 10.6266F = 10.624 であるが，OneWayANOVATest の 10.6266 は誤りである。

後に示す onewayANOVA(x1, x2, x3, x4) # W = 10.623885358647904 が正しい。

### 1.2. LeveneTest


```julia
LeveneTest(x1, x2, x3, x4)
```




    Levene's test
    -------------
    Population details:
        parameter of interest:   Variances
        value under h_0:         "all equal"
        point estimate:          NaN
    
    Test summary:
        outcome with 95% confidence: fail to reject h_0
        p-value:                     0.0850
    
    Details:
        number of observations: [9, 10, 8, 8]
        W statistic:            2.41877
        degrees of freedom:     (3, 31)





```julia
using RCall
R"""
library(lawstat)
#####
a = lawstat::levene.test($vx, $fx, location="mean")
print(a) # Test Statistic = 2.4024, p-value = 0.0865
"""
```
    	Classical Levene's test based on the absolute deviations from the mean
    	( none not applied because the location is not set to median )
    
    data:  `#JL`$vx
    Test Statistic = 2.4024, p-value = 0.0865

    RObject{VecSxp}
    
    	Classical Levene's test based on the absolute deviations from the mean
    	( none not applied because the location is not set to median )
    
    data:  `#JL`$vx
    Test Statistic = 2.4024, p-value = 0.0865
    




R では 2.4024 と正しいが，LeveneTest の 2.41877 は誤りである。

後に示す levene の 2.4024185774979356 はただしい

### 1.3. BrownForsytheTest


```julia
BrownForsytheTest(x1, x2, x3, x4)
```




    Levene's test
    -------------
    Population details:
        parameter of interest:   Variances
        value under h_0:         "all equal"
        point estimate:          NaN
    
    Test summary:
        outcome with 95% confidence: fail to reject h_0
        p-value:                     0.1014
    
    Details:
        number of observations: [9, 10, 8, 8]
        W statistic:            2.25728
        degrees of freedom:     (3, 31)





```julia
vx = vcat(x1, x2, x3, x4...);
fx = vcat([repeat([cat], inner=times) for (cat, times) in zip(1:4, [9, 10, 8, 8])]...);
```


```julia
println(fx)
```

    [1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4]



```julia
#####
R"""
b = lawstat::levene.test($vx, $fx, location="median")
print(b) # Test Statistic = 2.2425, p-value = 0.103
""";
```

    
    	Modified robust Brown-Forsythe Levene-type test based on the absolute
    	deviations from the median
    
    data:  `#JL`$vx
    Test Statistic = 2.2425, p-value = 0.103
    



```julia
R"""
library(car)
vx = c(301, 311, 325, 291, 388, 412, 325, 361, 287, 197, 180, 247, 260, 247, 199, 179, 134, 163, 200, 209, 302, 187, 166, 234, 290, 175, 116, 342, 216, 316, 386, 324, 145, 254, 228)
fx = factor(rep(c("A", "B", "C", "D"), c(9, 10, 8, 8)))
c = car:::leveneTest($vx, factor($fx), center=mean)
print(c) # group  3  2.4024 0.0865
#####
d = car:::leveneTest($vx, factor($fx)) # center=median
print(d) # group  3  2.2425  0.103
""";
```

    Levene's Test for Homogeneity of Variance (center = mean)
          Df F value Pr(>F)
    group  3  2.4024 0.0865
          31               
    Levene's Test for Homogeneity of Variance (center = median)
          Df F value Pr(>F)
    group  3  2.2425  0.103
          31               

R の Levene 検定，Brown-Forsythe 検定それぞれを指定した検定結果はともに正しい。HypothesisTests の結果は誤りである。

## 2. Levene 検定の数式的説明

* $k$ グループ数
* $N_i$ 第 $i$ グループのサンプルサイズ
* $N$ 全体のサンプルサイズ $N = \displaystyle \sum_{i=1}^{k} N_i$
* $Y_{ij}$ 第 $i$ グループの $j$ 番目の測定値
* $\bar{Y}_{i\cdot}$ 第 $i$ グループの平均値
* $Z_{ij} = |Y_{ij} - \bar{Y}_{i\cdot}|$ 第 $i$ グループの測定値から第 $i$ グループの平均値を引いたもの
* $Z_{i\cdot} = \displaystyle \frac{1}{N_i} \sum_{j=1}^{N_i} Z_{ij}$ 第 $i$ グループの $Z_{ij}$ の平均値
* $Z_{\cdot\cdot} = \displaystyle \frac{1}{N} \sum_{i=1}^k \sum_{j=1}^{N_i} Z_{ij}$ $Z_{ij}$ の平均値

のとき，検定統計量 $W$ は以下のように定義される。

$\hspace{2cm}W = \displaystyle \frac{(N-k)}{(k-1)} \cdot \frac{\displaystyle \sum_{i=1}^k N_i (Z_{i\cdot}-Z_{\cdot\cdot})^2}
    {\displaystyle \sum_{i=1}^k \sum_{j=1}^{N_i} (Z_{ij}-Z_{i\cdot})^2}$



検定統計量 $W$ は，近似的に，自由度 ($k-1$ and $N-k$) の $F$ 分布にしたがう。
よって，$p$ 値は以下のように求められる。

$\hspace{2cm}p = \Pr\{F_{k-1, N-k}(x) \geqq W\}$



### 2.1. 関数定義

以下に，一応正しく動く関数を提示しておく。

元の関数定義より，簡略に表現できていると思うが。


```julia
using Distributions, StatsBase, FreqTables

function onewayANOVA(x...)
    k, Ni, SSti, SSei = anova(x...)
    printout(Ni, SSti, SSei,
        "One-way ANOVA test (not Welch)")
end

function levene(x...; statfunc=mean)
    name = statfunc == mean ? "Levene's test" : "Brown-Forsythe's test"
    Zij = [abs.(y .- statfunc(y)) for y in x]
    k, Ni, SSti, SSei = anova(Zij...)
    printout(Ni, SSti, SSei, name)
end

brownforsythe(x...) = levene(x...; statfunc=median)

function  flignerkilleen(x...)
    Zij = [abs.(y .- median(y)) for y in x]
    ranks = tiedrank(vcat(Zij...))
    ties = freqtable(ranks)
    tieadj = sum(ties.^3 .- ties) 
    qij = quantile.(Normal(), 0.5 .+ ranks ./ 2(length(ranks)+1))
    Ni = pushfirst!(cumsum([length(y) for y in x]), 0)
    Qij = [qij[(Ni[i]+1):(Ni[i+1])] for i in 1:length(Ni)-1]
    k, Ni, SSti, SSei = anova(Qij...)
    printout(Ni, SSti, SSei, "Fligner-Killeen test")
end

function anova(x...)
    Ni = [length(y) for y in x]
    Z̄i = mean.(x)
    Z̄ = sum(Z̄i .* Ni) / sum(Ni)
    SSti = Ni .* (Z̄i .- Z̄).^2
    SSei = sum.( (z .- z̄).^2 for (z, z̄) in zip(x, Z̄i) )
    (length(Ni), Ni, SSti, SSei)
end

function printout(Ni, SSti, SSei, name)
    println(name)
    DFt = length(Ni) - 1
    DFe = sum(Ni) - DFt - 1
    y = sum(SSti) / sum(SSei)
    if occursin("Fligner", name)
        chisq = y * (DFe + DFt) / (1 + y)
        p = ccdf(Chisq(DFt), chisq)
        println("Chi-sq. = $chisq,  df = $DFt,  p value = $p")
    else
        F = y * (DFe / DFt)
        p = ccdf(FDist(DFt, DFe), F)
        println("W = $F,  df = ($DFt, $DFe),  p value = $p")
    end
end;

```


### 2.2. 実行例



```julia
levene(x1, x2, x3, x4) # W = 2.4024185774979356,  df = (3, 31),  p value = 0.08650303516190276
```

    Levene's test
    W = 2.4024185774979356,  df = (3, 31),  p value = 0.08650303516190276



```julia
brownforsythe(x1, x2, x3, x4) # W = 2.242478379513558,  df = (3, 31),  p value = 0.10300664426579445
```

    Brown-Forsythe's test
    W = 2.242478379513558,  df = (3, 31),  p value = 0.10300664426579445



```julia
onewayANOVA(x1, x2, x3, x4) # W = 10.623885358647904,  df = (3, 31),  p value = 5.7907091792266794e-5
```

    One-way ANOVA test (not Welch)
    W = 10.623885358647904,  df = (3, 31),  p value = 5.7907091792266794e-5



```julia
flignerkilleen(x1, x2, x3, x4) # Chi-sq. = 5.380295369379746,  df = 3,  p value = 0.14597618474787735
```

    Fligner-Killeen test
    Chi-sq. = 5.380295369379746,  df = 3,  p value = 0.14597618474787735