Bing Image Creator で生成: Close up of yellow Chimonanthus praecox flowers, background is white clouds in the sky, photo
の続きです。
これまでの結果、沖縄県が一番、老齢化指数が低く、秋田県が一番、老齢化指数が高いことがわかりました。今回は、R の infer パッケージを利用して、二つの県の老年化指数が統計的に有意に違うのかを確認します。
まず infer パッケージを library() 関数を使って読み込みます。
秋田県と沖縄県の elder_index: 老齢化指数の平均値を調べるためのデータフレームを作成します。
２つの県の平均値を確認します。
秋田県の平均値は、328で、沖縄県の平均値は、116なので差は、212です。
infer パッケージのワークフローに従って差の平均値を求めます。
を参考にしています。
212となっています。
次に null distribution を求めます。これは、pref と elder_index をバラバラに組み合わせて、秋田県と沖縄県の平均値の差を計算することを何回も繰り返すことです。
1回目のバラバラの組み合わせでは、差は-3でした。2回目は-55.6、3回目は6.58となっています。
実際の差の212と、この null distibution を視覚化して比較します。
ヒストグラムが null distribution です。秋田県と沖縄県の elder_index の値をバラバラにシャッフルして差を計算していますので、0を中心とした分布になっています。赤い垂線が実際の差の212の線です。elder_index の値は秋田県、沖縄県と全く関係がなければ212という大きな差は起こりえない感じです。
p値を確認します。
p値は0です。
以上が infer パッケージのワークフローを使ったシミュレーションベースの統計的推定です。
昔からある理論ベースの t検定もやってみます。t.test() 関数を使います。
akita, okinawa とうオブジェクト名で、それぞれの elder_index のベクトルを作成して、t.test() 関数で処理しています。
p値は 1.966e-07 とほぼ0です。差の 95%信頼区間は、181 ~ 243 です。
infer パッケージのワークフローでも信頼区間は求めることができます。
やってみましょう。
を参考にしています。
d_hat は既に求めていますので、ブートストラップ法による 差の distribution を求めます。
get_ci() 関数で信頼区間を求めます。
グラフにしましょう。
緑の垂線が、infer パッケージのワークフローで求めたシミュレーションベースの95%信頼区間で、赤い垂線が、t.test() 関数で求めた理論ベースの95%信頼区間です。
シミュレーションベースの信頼区間のほうが幅が狭いですね。
数値を確認しておきます。
186 ~ 238 です。
今回は以上です。
次回は、
です。
初めから読むには、
です。
今回のコードは以下になります。
#
# inferパッケージを読み込む
library(infer)
#
# 沖縄県と秋田県のelder_index: 老齢化指数の比較
df2 <- df |>
filter(pref == "秋田県" | pref == "沖縄県") |>
mutate(pref = factor(pref, levels = c("秋田県", "沖縄県")))
df2
#
# 秋田県と沖縄県の平均値
df2 |>
group_by(pref) |>
summarize(avg = mean(elder_index))
#
# 平均値の差
d_hat <- df2 |>
specify(elder_index ~ pref) |>
calculate(stat = "diff in means",
order = c("秋田県", "沖縄県"))
d_hat
#
# null distribution
set.seed(123)
null_dist <- df2 |>
specify(elder_index ~ pref) |>
hypothesize(null = "independence") |>
generate(reps = 1000, type = "permute") |>
calculate(stat = "diff in means",
order = c("秋田県", "沖縄県"))
null_dist
#
# null distribution と d_hat
null_dist |>
visualize() +
shade_p_value(obs_stat = d_hat, direction = "two-sided")
#
# p値
null_dist |>
get_p_value(obs_stat = d_hat, direction = "two-sided")
#
# 理論ベースの t検定
akita <- df2 |>
filter(pref == "秋田県") |>
pull(elder_index)
akita
#
okinawa <- df2 |>
filter(pref == "沖縄県") |>
pull(elder_index)
okinawa
#
t.test(akita, okinawa)
#
# ブートストラップ法
set.seed(123)
boot_dist <- df2 |>
specify(elder_index ~ pref) |>
generate(reps = 1000, type = "bootstrap") |>
calculate(stat = "diff in means",
order = c("秋田県", "沖縄県"))
boot_dist
#
# 信頼区間
percentile_ci <- boot_dist |>
get_ci(type = "percentile")
#
# boot_dist と信頼区間のグラフ
visualize(boot_dist) +
shade_confidence_interval(endpoints = percentile_ci) +
geom_vline(xintercept = 181, color = "red", linewidth = 2) +
geom_vline(xintercept = 243, color = "red", linewidth = 2)
#
# 信頼区間の数値
percentile_ci
#
