の続きです。前回は日銀の保有する国債の残高についてグラフでどのようなデータかを確認しました。その結果、国債1つ1つの保有残高の平均値は2025年3月10日と2026年3月10日で大きな変化は無いこと、国債の種類別の平均値は大きな違いがあることがグラフで見て取れました。
今回は統計学的な検定をして、年の違いは無いこと、種類別の違いはあることを確認します。
まずは、2025年と2026年で平均値に違いが無いことを確認します。これはt検定になりますので、t.test()関数でを使います。
p-valueが0.3899となっています。つまり2025年の平均値(17582.53)と2026年の平均値(16158.28)に統計学的に有意な違いがあることは言えない、つまり違いは無いということです。
これだけだと面白くないので、inferパッケージのワークフローでシミュレーションベースの検定もしてみましょう。
このウェブサイトを参考にしています。はじめにinferパッケージを読み込みます。
次に、2025年のamtの平均値と2026年のamtの平均値の差を計算します。
次はNull Distributionを生成します。これは、dateとamtをシャッフルしてdate別のamtの平均値の差を計算する、ということを何回も繰り返します。
このNull distributionをヒストグラムにして分布をみます。
赤い垂線が実際の差の位置です。こうしてみると、実際の差はシャッフルして計算された差の分布の中に入っていますね。
この分布を基にp値を計算します。上のグラフの赤いシャドーの部分の面積になります。
p値は0.398とt.test()関数での0.3899と同じくらいです。シミュレーションベースでも2025年と2026年の残高の平均値に差が無いことが確認できました。
次は、種類別の残高の平均値に差があるかどうかです。複数のカテゴリー別の平均値に差があるかどうかはANOVA(Analysis Of Variance)です。aov()関数でモデルを作り、summary()関数で結果を見る、という流れです。
F値は212.6でp値は <2e-16です。つまり種類別に違いが無いとは言えない、種類別に違いがある、ということです。
どの種類と、どの種類に違いがあるかは、TukeyHSD()関数でわかります。
y5とy2などが違いがあるとわかります。
ANOVAもinferパッケージのワークフローでシミュレーションベースでやってみましょう。
のサイトを参考にしました。
はじめにF値を計算します。
F値のNull distributionを生成します。
Null distributionを視覚化してみます。
赤い垂線が実際のF値の位置ですが、Null distributionの分布から大きく離れていて、Null distributionの形がよくわからないです。F値を書かないで分布をみてみます。
典型的なF値の分布です。
p値を計算します。
p値は0です。シミュレーションベースでも種類別の残高の平均値には違いがあることが確認できました。
今回は以上です。
次回は
です。
はじめから読むには、
です。
今回のコードは以下になります。
#
# dateでamtの平均値は違うのか？
# t検定(理論ベース)
t.test(df$amt[df$date == "2025-03-10"],
df$amt[df$date == "2026-03-10"])
#
# t検定(シミュレーションベース)
# step 1: inferパッケージの読み込み
library(infer)
#
# step 2: 両者の差を計算して保存
obs_stat <- df |>
specify(amt ~ date) |>
calculate(stat = "diff in means", order = c("2025-03-10", "2026-03-10"))
obs_stat
#
# step 3: null distributionの生成
set.seed(9)
null_dist <- df |>
specify(amt ~ date) |>
hypothesize(null = "independence") |>
generate(reps = 1000, type = "permute") |>
calculate(stat = "diff in means", order = c("2025-03-10", "2026-03-10"))
#
# step 3: null distributionのグラフ
null_dist |>
visualize() +
shade_p_value(obs_stat, direction = "two-sided")
#
# step 4: p値の計算
p_value_sample <- null_dist |>
get_p_value(obs_stat, direction = "two-sided")
p_value_sample
#
# typeでamtの平均値は違うのか?
# ANOVA検定(理論ベース)
aov_fit <- aov(amt ~ type, data = df)
summary(aov_fit)
#
# どのtype同士が違うか
TukeyHSD(aov_fit)$type |>
as.data.frame() |>
filter(`p adj` < 0.05)
#
# ANOVA(シミュレーションベース)
# Step 1: F値を計算して保存
observed_F <- df |>
specify(amt ~ type) |>
calculate(stat = "F")
observed_F
#
# Step 2: Null Distributionの生成
set.seed(15)
null_dist_F <- df |>
specify(amt ~ type) |>
hypothesize(null = "independence") |>
generate(reps = 1000, type = "permute") |>
calculate(stat = "F")
#
# Step 3: Null Distributionの視覚化
null_dist_F |>
visualize() +
shade_p_value(obs_stat = observed_F, direction = "greater")
#
# Step 3.1: p-value無しで視覚化
null_dist_F |> visualize()
#
# Step 4: p値の計算
p_value_F <- null_dist_F |>
get_p_value(obs_stat = observed_F, direction = "greater")
p_value_F
#
(冒頭の画像は、Bing Image Creator で生成しました。プロンプトは Long and wide view of natural wild field, close up of red Chaenomeles Japonica tree. Photo です。)
