(冒頭の画像は、Bing Image Creator で生成しました。プロンプトは Close up photograph of yellow daisy flowers, flowering under the high blue sky, photo です。)
の続きです。今回は昼間の人口が多い都市と夜の人口が多い都市の違いを見てみます。
group_by()関数を利用して、day == 1 の都市と day == 0 の都市の平均値を計算してみます。
人口密度は、day == 0, つまり夜間人口のほうが多い都市のほうが高いですね。商品販売額と課税所得は昼間人口のほうが多い都市が多いです。
この２つのグループで、人口密度に違いはあるのか、調べてみます。
まず、グラフで分布を見てみます。
昼間の人口が多い都市(day == 1)のほうが左側(人口密度が低い)に集中してますね。
t検定をt.test()関数で実行してみます。
p-valueが0.00053となっていますので、二つの都市グループでmitsudoの平均値が同じとは言えませんね。2つのグループの差の95%信頼区間は、-4433 ~ -1251です。
ヒストグラムの分布形状を見ると、mitsudoは正規分布していませんので、シミュレーションベースでの検定もやってみます。
inferパッケージの読み込みをします。
まず。day == 0 の都市と day == 1 の都市の平均値の差を計算して保存します。
差はマイナス2842でした。5188 - 8031 = -2843 と一致します。(小数点以下の四捨五入のため1の位の値が違っています。)
次は null distribution, (これは、day と mitsudo をバラバラにシャッフルして、day == 0, dat == 1 の都市の平均値を計算して、その差を計算、という作業を何回も繰り返す作業です。)を生成します。
reps = 1000 として、1000回繰り返しました。
グラフにしてみます。
赤い垂線がマイナス2842の線です。マイナス2842という差が発生するのはとても稀なことだとわかります。get_p_value()関数でp値を計算します。
p値は0.002となりました。t.test()関数での値よりも大きいですが、やっぱりday == 0 の都市とday == 1の都市でmitsudoの平均値は同じとは言えませんね。
95%信頼区間も算出しましょう。
まず、Bootstrap distribution を生成します。
get_confidence_interval()関数で95%信頼区間を取得します。
-4443 ~ -1202 でした。t.test()関数では、-4433 ~ -1251でしたので、同じようなものですね。
グラフで表示します。
夜間人口のほうが多い都市のほうが人口密度は高いということが統計学的に確認できました。
今回は以上です。
次回は
です。
はじめから読むには、
です。
今回のコードは以下になります。
#
# day == 1, day == 0別の平均値
df |>
group_by(day) |>
summarize(
mitsudo_avg = mean(mitsudo),
hanbai_avg = mean(hanbai),
tax_mean = mean(tax)
)
#
# dayによって、mitsudo(の平均)に違いはあるか？
#
# グラフで確認
df |>
ggplot(aes(x = mitsudo)) +
geom_histogram(aes(fill = day)) +
facet_wrap(~ day, nrow = 2) +
theme_minimal()
#
# t.test()で確認
t.test(
df$mitsudo[df$day == 1],
df$mitsudo[df$day == 0]
)
#
# inferで確認
# パッケージの読み込み
library(infer)
#
# obs_statの計算
obs_stat <- df |>
mutate(day = as.factor(day)) |>
specify(mitsudo ~ day) |>
calculate(stat = "diff in means", order = c("1", "0"))
obs_stat
#
# null distributionの生成
null_dist <- df |>
mutate(day = as.factor(day)) |>
specify(mitsudo ~ day) |>
hypothesize(null = "independence") |>
generate(reps = 1000, type = "permute") |>
calculate(stat = "diff in means", order = c("1", "0"))
#
# null_distをグラフにする
null_dist |>
visualize() +
shade_p_value(obs_stat,
direction = "two-sided")
#
# p値の計算
null_dist |>
get_p_value(obs_stat,
direction = "two-sided")
#
# 信頼区間の計算
# Bootstrap distributionの生成
boot_dist <- df |>
mutate(day = as.factor(day)) |>
specify(mitsudo ~ day) |>
generate(reps = 1000, type = "bootstrap") |>
calculate(stat = "diff in means", order = c("1", "0"))
#
# 信頼区間の取得
percentile_ci <- boot_dist |>
get_confidence_interval(level = 0.95)
percentile_ci
#
# Bootdistributionと信頼区間のグラフ
boot_dist |>
visualize() +
shade_confidence_interval(endpoints = percentile_ci)
#
