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

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

企業行動に関するアンケート調査のデータ分析5 - Rによる重回帰分析(シミュレーションベース)の信頼区間を算出

Bing Inage Creatorで生成: Full moon night landscape photo, silver grass fluttering in the wind.

www.crosshyou.info

の続きです。

前回は、理論ベース(公式ベース)で重回帰分析の係数の信頼区間を算出しました。

今回はシミュレーションベースの信頼区間を算出します。

Full infer Pipeline Examples • infer

に書かれている手順に従って実行していきます。

まずは、普通の回帰分析をして、データから係数を取得します。

これは、前回のlm_mod2と同じです。

続いて、ブートストラップ法で何回も仮想サンプルを作成して、そのサンプルでモデルの係数を推定するのを繰り返します。

1回目は、hensaの係数は-0.977で、2回目は-2.10で3回目は-1.66となっていることが上の図でわかります。

そして、get_confidence_interval()関数で信頼区間を算出するだけです。

簡単ですね。シミュレーションベースでは、いちいち残差が不均一分散なのかどうかを考慮する必要はありません。

visualization()関数で視覚化してみましょう。

hensaの信頼区間は0を含んでいませんので、hensaは統計学的に有意な変数であることがわかります。

前回の理論ベース(公式ベース)の信頼区間とも合わせて視覚化してみます。

こうしてみると、シミュレーションベースのほうが左側によっていますね。

y(year)も視覚化してみます。

y(year)の信頼区間はシミュレーションベースの信頼区間のほうが幅が広いですね。

今回は以上です。

次回は

www.crosshyou.info

です。

 

最初から読むには、

www.crosshyou.info

です。

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

# シミュレーションベースの回帰分析
# step1: 観察された係数を確認
obs_fit <- df_all |> 
  specify(forecast ~ y + hensa) |> 
  fit()
obs_fit
#
# lm_mod2
coef(lm_mod2)
#
# step2: ブートストラップ
set.seed(987)
boot_dist <- df_all |> 
  specify(forecast ~ y + hensa) |> 
  generate(reps = 1000, type = "bootstrap") |> 
  fit()
boot_dist
#
# step3: 信頼区間を算出
conf_ints <- boot_dist |> 
  get_confidence_interval(level = 0.95,
                          point_estimate = obs_fit)
conf_ints
#
# step4: 視覚化
boot_dist |>
  visualize() +
  shade_confidence_interval(endpoints = conf_ints)
#
# hensaの信頼区間の比較
boot_dist |> 
  filter(term == "hensa") |> 
  ggplot(aes(x = estimate)) +
  geom_histogram(color = "black", fill = "white") +
  geom_vline(aes(xintercept = -2.00), color = "red", linewidth = 2) +
  geom_vline(aes(xintercept = -0.0942), color = "red", linewidth = 2) +
  geom_vline(aes(xintercept = -1.865), color = "blue", linewidth = 2) +
  geom_vline(aes(xintercept = -0.0346), color = "blue", linewidth = 2) +
  geom_vline(aes(xintercept = 0), color = "black") +
  labs(title = "hensaの95%信頼区間の比較",
       subtitle = "赤はシミュレーションベース、青は理論ベース",
       x = "hensa")
#
# y(year)の信頼区間の比較
boot_dist |> 
  filter(term == "y") |> 
  ggplot(aes(x = estimate)) +
  geom_histogram(color = "black", fill = "white") +
  geom_vline(aes(xintercept = 0.00202), color = "red", linewidth = 2) +
  geom_vline(aes(xintercept = 0.141), color = "red", linewidth = 2) +
  geom_vline(aes(xintercept = 0.01492), color = "blue", linewidth = 2) +
  geom_vline(aes(xintercept = 0.1321), color = "blue", linewidth = 2) +
  geom_vline(aes(xintercept = 0), color = "black") +
  labs(title = "y(year)の95%信頼区間の比較",
       subtitle = "赤はシミュレーションベース、青は理論ベース",
       x = "y(year")