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

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

企業行動に関するアンケート調査のデータ分析4 - Rによる重回帰分析(理論ベースによる方法)。不均一分散があるときの信頼区間の算出方法

Bing Image Creatorで生成: Pictures of beautiful fall foliage landscapes, with waterfalls in the background.

www.crosshyou.info

の続きです。

前回はRで相関係数の信頼区間を調べました。理論ベース(公式ベース)による方法とシミュレーションベースの方法をしました。シミュレーションベースのほうが信頼区間の幅は広かったことがわかりました。

今回はRで重回帰分析をして、係数の信頼区間を調べます。前回と同じく、理論ベース(公式ベース)とシミュレーションベースの信頼区間を調べてみます。

まず、sectorが全産業だけのデータフレームを作成しておきます。

forecast, hensa, obsの推移のグラフを描いてみます。

forecastは2008年、リーマンショックのあった年は大きく落ち込んでいます。

その反対にhensaは2008年は大きくなっています。

obsはV字型の推移です。

それぞれの散布図を描いてみます。

一番上の行の散布図が、どれもY軸がforecastのものなので、注目します。

hensaは負の相関がありそうですね。

それでは、理論ベース(公式ベース)の回帰分析をしてみます。forecastを被説明変数にして、y, hensa, obsを説明変数にしてみます。

obsの係数のp値は0.6556と大きな値なので、obsは除外したほうがよさそうです。

yの係数は0.07351でp値は0.01454なので統計学的に有意です。1年経過するごとに、hensaが同じならば、forecastは0.07351増えるということですね。

hensaの係数は-0.94984でp値は0.00515なので統計学的に有意です。hensaが1増えるとforecastは0.94984低下する、ということですね。

forecastの平均値は1.11なので、1年経過するごとに0.07351の上昇というのはあまり大きな影響ではないかもしれません。

理論ベースでの回帰分析では、残差の分散が均一でないと正しい信頼区間がわかりませんので、lmtestパッケージのbptest()関数で確認します。

p値が0.0384と0.05よりも小さいので不均一分散があることが疑われます。

不均一分散があるときの信頼区間は、lmtestパッケージとsandwichパッケージを使って計算します。

ロバストな標準誤差を計算します。

回帰モデルの係数を取り出しておきます。

信頼区間を算出します。

不均一分散を考慮しなかった場合の信頼区間を計算してみます。

yはそれほど変わりませんが、hensaは結構変わりますね。不均一分散を考慮したロバストな信頼区間では、hensaの上限は-0.034とぎりぎり0よりも小さいですが、不均一分散を考慮していない信頼区間では、上限が-0.32と余裕で0よりも大きいです。

今回は以上です。

次回はシミュレーションベースの信頼区間を求めてみます。

www.crosshyou.info

 

初めから読むには、

www.crosshyou.info

です。

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

#
# sectorが全産業だけのデータフレームを作成する
df_all <- df |> 
  filter(sector == "全産業")
df_all
#
# forecastの推移
df_all |> 
  ggplot(aes(x = y, y = forecast)) +
  geom_line() +
  geom_point() +
  labs(x = "year")
#
# hensaの推移
df_all |> 
  ggplot(aes(x = y, y = hensa)) +
  geom_line() +
  geom_point() +
  labs(x = "year")
#
# obsの推移
df_all |> 
  ggplot(aes(x = y, y = obs)) +
  geom_line() +
  geom_point() +
  labs(x = "year")
#
# 散布図マトリックス
df_all |> 
  select(forecast, hensa, obs, y) |> 
  plot()
#
# forecastをy, hensa, obsで回帰分析
lm_mod <- lm(forecast ~ y + hensa + obs, data = df_all)
summary(lm_mod)
#
# obsを削除して回帰分析
lm_mod2 <- lm(forecast ~ y + hensa, data = df_all)
summary(lm_mod2)
#
# forecast, hensaの平均値
df_all |> 
  summarize(
    forecast_average = mean(forecast),
    hensa_average = mean(hensa)
  )
#
# Heteroskedasticity(不均一分散)の検定
lmtest::bptest(lm_mod2)
#
# lmtest, sandwichパッケージを読み込む
library(lmtest)
library(sandwich)
#
# ロバストな標準誤差を計算
robust_se <- sqrt(diag(vcovHC(lm_mod2, type = "HC1")))
robust_se
#
# モデルの係数
coef_est <- coef(lm_mod2)
coef_est
#
# ロバストな信頼区間を算出
robust_conf_int <- coef_est + cbind(-1.96 * robust_se, 1.96 * robust_se)
colnames(robust_conf_int) <- c("2.5 %", "97.5 %")
robust_conf_int
#
#
# ロバストでない95%信頼区間
non_robust_conf_int <- confint(lm_mod2)
non_robust_conf_int
#