裏 RjpWiki

文字通り,RjpWiki の裏を行きます
R プログラム コンピュータ・サイエンス 統計学

Julia で一元配置分散分析(データフレーム)

2021年01月15日 | ブログラミング

Julia でANOVA」 という記事がある。記事の日付は2019年3月15日とわずか2年足らず前のものである。Julia v1.1.1 がリリースされたのが 2019/05/16 であるから,少なくともそれより前のバージョンということになる。

しかしその後の Julia の発展が著しく,今では廃止された関数があったりして,そのままでは参考にならない。また,示された一元配置分散分析は各群の例数が等しく,また,各群の等分散を仮定する検定なので,等分散を仮定しない検定も示しておく。

廃止あるいは仕様変更のあった関数

dat = CSV.read("foo.csv", allowmissing=:none)

dat = CSV.read("lfoo.csv", DataFrame)

と書くようになった。今でも多くのページで単にファイル名を指定すればよいと書いてあるので,小生も最初は困り果てた。 
また,allowmissing は不要になった(当然のものとして組み込まれた)。

by() によるグループごとの統計処理は廃止された。

by(dat, :parasite, :growth => mean)

は groupby() と combine() で以下のように書く。 combine() は非常に強力である。

gd = groupby(dat, :parasite);
combine(gd, :growth => mean)
stat = combine(gd, nrow => :sample_size, :growth_rate => mean => :group_mean, :growth_rate => var => :group_var)

名前の変更は names!() ではなく rename!(),
列の削除は deletecols!() ではなく,select!()
列の並べ替えも permutecols!() ではなく,select!() で行える。

さて,以上の点に留意して,以下のような処理手順を示す。

まずは,データ入力から前処理まで。

using CSV, DataFrames, Statistics, Distributions

dat = CSV.read("Daphniagrowth.csv", DataFrame);
# '.' を使っている列名は,最初に '_' に変更しておく
rename!(dat, "growth.rate" => "growth_rate");
# 群ごとの統計は groupby(), combine() を使う
gd = groupby(dat, :parasite);
stat = combine(gd, nrow => :sample_size, :growth_rate => mean => :group_mean, :growth_rate => var => :group_var)
println(stat)

4×4 DataFrame

 Row │ parasite                   sample_size  group_mean  group_var 
     │ String                     Int64        Float64     Float64   
─────┼───────────────────────────────────────────────────────────────
   1 │ control                             10    1.21391   0.0130376
   2 │ Metschnikowia bicuspidata           10    0.801154  0.0465879
   3 │ Pansporella perplexa                10    1.07636   0.0322315
   4 │ Pasteuria ramosa                    10    0.482203  0.0375737

後は,計算処理を黙々と行う(ドット演算子に注意)
まずは,各群の分散が等しいことを仮定する検定。

k = length(gd)
sample_size = stat[!, :sample_size];
group_mean  = stat[!, :group_mean ];
group_var   = stat[!, :group_var  ];
n = sum(sample_size)
grand_mean = sum(sample_size .* group_mean) / n
St = var(dat[!, :growth_rate]) * (n - 1)
Sb = sum(sample_size .* (group_mean .- grand_mean) .^ 2)
Sw = St - Sb
dfb = k - 1
dfw = n - k
dft = n - 1
Vb = Sb / dfb
Vw = Sw / dfw
# 最終的に必要なのは,F 統計量と,P 値
F = Vb / Vw
p_value = ccdf(FDist(dfb, dfw), F)
# 見やすくするためにデータフレームにしておこう
ANOVA_table = DataFrame(SS=[Sb, Sw, St], df=[dfb, dfw, dft],
      MS=[Vb, Vw, ""], F=[F, "", ""], p_value=[p_value, "", ""])
println(ANOVA_table)

3×5 DataFrame

 Row │ SS       df     MS         F        p_value     
     │ Float64  Int64  Any        Any      Any         
─────┼─────────────────────────────────────────────────
   1 │ 3.13791      3  1.04597    32.3252  2.57128e-10
   2 │ 1.16488     36  0.0323577
   3 │ 4.30278     39


次は,各群の分散が等しいとは仮定しない,いわゆるウェルチの方法によるものである。

grand_mean = mean(dat[!, :growth_rate])
k = length(gd)
sample_size = stat[!, :sample_size];
group_mean = stat[!, :group_mean];
group_var = stat[!, :group_var];
wi = sample_size ./ group_var;
w = sum(wi)
m = sum(wi .* group_mean) / w
a = sum(((1 .- wi ./ w) .^ 2) ./ (sample_size .- 1)) / (k^2 - 1)
F = sum(wi .* (group_mean .- m) .^ 2) / (k - 1) / (1 + 2a*(k - 2))
dfb = k - 1
dfw = 1 / 3a
p_value = ccdf(FDist(dfb, dfw), F)
println("F = $F, p_value = $p_value")

コメント   この記事についてブログを書く
« 指定された平均値と標準偏差... | トップ | 数学の問題を Julia で解く »
最新の画像もっと見る

コメントを投稿

ブログ作成者から承認されるまでコメントは反映されません。

ブログラミング」カテゴリの最新記事