Rで何かをしたり、読書をするブログ

政府統計の総合窓口のデータや、OECDやUCIやのデータを使って、Rの練習をしています。ときどき、読書記録も載せています。

全国の主要都市の人口密度、昼夜人口比率、商品販売額、課税所得のデータの分析5 - 昼夜人口比率と商品販売額の関係

www.crosshyou.info

前回の続きです。前回は昼夜人口比率と人口密度の関係を調べてみました、その結果、昼間の人口のほうが多い都市のほうが人口密度が低いことがわかりました。

今回は、昼夜人口比率と商品販売額の関係を調べてみましょう。前回は、dayという変数(これは昼間の人口が多い都市は1で、そうでない都市は0というダミー変数)を使いました。今回はdaynight(比率そのもの)を使ってみます。

まず、散布図を描いてみます。

daynightの値が大きいほど、hanbaiの値が大きい傾向のようです。

2つの数値データなので、cor.test()関数で相関関係の検定をしてみます。

p-valueが2.2e-16とほとんど0ですので、相関関係があるのは間違いないようです。相関係数の95%信頼区間は0.74~0.84となっています。

今回も前回と同じように、inferパッケージのワークフローでシミュレーションベースでの相関関係を検定してみます。

相関係数を計算して保存します。

相関係数は0.797です。これは当然ですが、cor.test()関数と同じです。

続いて、null distributionを計算します。これは、daynightとhanbaiをシャッフルして相関係数を計算する作業を何回も繰り返すことです。今回は10000回繰り返してみましょう。

私のラップトップでは、3秒くらいかかりました。

こうして生成したnull distributionをグラフにします。

赤い垂線が実際の相関係数です。このグラフを見ると、実際の相関係数は、シャッフルした組み合わせでは実現できないほど大きな相関係数だとわかります。

p値を計算します。

p値は0でした。

95%信頼区間を求めます。まずはブートストラップ法で何回も相関係数を計算します。

2万回計算しましょう。

8秒くらいかかりました。

95%信頼区間を計算します。

0.697~0.894となりました。これは、cor.test()関数で算出された0.74~0.84よりも幅広いですね。なんでかというと、散布図をみるとわかりますが、daynight, hanbaiが極端に大きい都市があるからですね。

グラフで視覚化します。

今回は以上です。

次回は

www.crosshyou.info

です。

はじめから読むには、

www.crosshyou.info

です。

今回のコードは以下になります。

#
# daynightと販売の散布図
df |> 
  ggplot(aes(x = daynight, y = hanbai)) +
  geom_point()
#
# 相関関係の検定
cor.test(df$daynight, df$hanbai)
#
# 相関係数を計算して保存
corr_hat <- df |> 
  specify(hanbai ~ daynight) |> 
  calculate(stat = "correlation")
corr_hat
#
# null distributionの生成
set.seed(273)
null_dist <- df |> 
  specify(hanbai ~ daynight) |> 
  hypothesize(null = "independence") |> 
  generate(reps = 10000, type = "permute") |> 
  calculate(stat = "correlation")
#
# null distributionのグラフ
null_dist |> 
  visualize() +
  shade_p_value(obs_stat = corr_hat, direction = "two-sided")
#
# p値の計算
null_dist |> 
  get_p_value(obs_stat = corr_hat, direction = "two-sided")
#
# Bootstrap法で何回も相関係数を計算
set.seed(290)
boot_dist <- df |> 
  specify(hanbai ~ daynight) |> 
  generate(reps = 200000, type = "bootstrap") |> 
  calculate(stat = "correlation")
#
# 95%信頼区間を計算
percentile_ci <- boot_dist |> 
  get_confidence_interval(level = 0.95)
percentile_ci
#
# boot distと信頼区間のグラフ
boot_dist |> 
  visualize() +
  shade_confidence_interval(endpoints = percentile_ci,
                            color = "lightblue", fill = "pink") +
  theme_bw()
#

(冒頭の画像は、Bing Image Creator で生成しました。プロンプトは Close up of flowers of New Zealand's National Flower, photo です。)