Bing Image Creatorで生成: Pictures of beautiful fall foliage landscapes, with waterfalls in the background.
の続きです。
前回は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よりも大きいです。
今回は以上です。
次回はシミュレーションベースの信頼区間を求めてみます。
初めから読むには、
です。
今回の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
#