(Bing Image Creator で生成: プロンプト: Close up of white Trifolium repens flowers, background is blue lake, blue sky, photo)
の続きです。
前回は、2011年度と2020年度の２つの年度のデータを使って、year のダミー変数を加えて回帰分析してみました。
その結果は、
income = -839.1 + 60.81 * popratio + 2122 * tokyo + 1.44 * energy + u
というモデルを推定しました。popratio: 15歳~64歳人口割合、energy: １人当たり最終エネルギー消費、ともに income: １人当たり県民総所得にはプラスの効果を持っていいる、という結果でした。
今回は、
https://www.econometrics-with-r.org/10.2-PDWTTP.html
こちらのサイトを参照して、"Before and After" 比較をしてみます。
このモデルは、
income = beta0 + beta1 * popratio + beta2 * energy + beta3 * 都道府県固有の効果 + u
というモデルです。都道府県固有の効果は、2011年度と2020年度で変わらずに同じだと仮定しています。beta0 は切片ですね。
すると、それぞれの年度で、
income(2011) = beta0 + beta1 * popratio(2011) + beta2 * energy(2011) + beta3 * 都道府県固有の効果 + u(2011)
という式と、
income(2020) = beta0 + beta1 * popratio(2020) + beta2 * energy(2020) + beta3 * 都道府県固有の効果 + u(2020)
という式が考えられます。
この２つの式を引き算すると、beta1の切片と『beta3 * 都道府県固有の効果』が消えてなくなります。
income(2020) - income(2011) = beta1 * (popratio(2020) - popratio(2011)) + beta2 * (energy(2020) - energy(2011) + (u(2020) - u(2011))
ですので、このモデル式を回帰分析して、beta1, beta2 の係数を求めることが可能です。
まずは、2011年だけのデータフレームを作成します。
2020年度だけのデータフレームは既に、df_2020 として作成済みですが、念のため、pcodeで並び替えて、df_2011と同じ都道府県の並びにしておきます。
そうしたら、income, popratio, energy の差分を生成します。
３つの差分の相関係数マトリックスを確認しましょう。
diff_income は diff_popratio とはマイナスの相関、diff_energy とはプラスの相関です。
では、lm() 関数でモデルの係数を推計します。
なんと！前回の推定結果とは逆に、popratio の係数がマイナスになりました。
summary() 関数でモデルのサマリーをみてみます。
diff_popratio の p値はかなり０に近いです。diff_energy の p値は 0.0532 なので、5%の有意水準ではわずかに有意ではありませんが、10%の有意水準では有意になります。
plot() 関数で model5 をプロットしてみます。
Breush-Pagan検定でHeteroskedasticity を確認します。
p-value は 0.6362 と 0.05 よりもうんと大きいですので、Heteroskedasticity ではないようです。
残差の2乗を回帰分析してみます。
こちらも p-value が 0.05 よりもかなり大きいです。
diff_income と diff_popratio の散布図を描いてみます。
diff_popratio が大きいほど、diff_income は小さいようです。
diff_income と diff_energy の散布図も描いてみます。
diff_energy が大きいほど、diff_income も大きいようです。
前回の結果とは popratio の効果の方向が逆になりました。どう判断したらいいでしょうか？
前回のデータで散布図を作成してみます。
まずは、popratio と income です。
どの件も2011年度が popratio が大きな値で、income は2020年度のほうが大きいようです。個々の都道府県だけに着目すると、popratio は2011年から2020年で減少、income は2011年から2020年で上昇、となっていますので、今回のモデルの結果のほうが正しいように見えます。
income と energy の散布図も描いてみます。
これはちょっとごちゃごちゃしていてわかりにくいですね。
沖縄県、愛知県、山口県、岡山県などは、energy は減少、income は増加のようですが、愛媛県、奈良県などは energyは減少、income も減少となっています。
この2つの散布図から判断すると、前回の結果よりも、今回のモデルのほうが、都道府県固有の効果を考慮しているので、あっているのかな～？と思いますが、よくわかりません。popratio: 15歳~64歳人口割合が増えると、income が減る、というちょっと直感とは違う結果なので、判断がつきません。
今回は以上です。
次回は
です。
初めから読むには、
です。
今回のコードは以下になります。
#
# https://www.econometrics-with-r.org/10.2-PDWTTP.html を参考にする
#
# 2011だけのデータフレームを作る
df_2011 <- df |>
filter(year == "2011年度") |>
arrange(pcode)
df_2011
#
# 2020だけのデータフレーム
df_2020 <- df_2020 |>
arrange(pcode)
df_2020
#
# income の差分
diff_income <- df_2020$income - df_2011$income
#
# popratio の差分
diff_popratio <- df_2020$popratio - df_2011$popratio
#
# energy の差分
diff_energy <- df_2020$energy - df_2011$energy
#
# diff_income, diff_popratio, diff_energy の相関係数マトリックス
tibble(diff_income, diff_popratio, diff_energy) |>
cor() |>
ggcorrplot(lab = TRUE)
#
# 差分のモデル
model5 <- lm(diff_income ~ diff_popratio + diff_energy - 1)
model5
#
# model5 のサマリー
summary(model5)
#
# model5 のプロット
par(mfrow = c(2, 2))
plot(model5)
par(mfrow = c(1, 1))
#
# Preush-Pagan Test
lmtest::bptest(model5)
#
# 残差の2乗を回帰分析
lm(I(model5$residuals^2) ~ diff_popratio + diff_energy) |>
summary()
#
# diff_income と diff_popratio
tibble(diff_income = diff_income,
diff_popratio = diff_popratio,
pref = df_2020$pref) |>
ggplot(aes(x = diff_popratio, y = diff_income)) +
geom_point(color = "red") +
geom_text(aes(label = pref), vjust = 1.2, alpha = 0.5) +
geom_smooth(method = "lm") +
theme_bw()
#
# diff_income と diff_energy
tibble(diff_income = diff_income,
diff_energy = diff_energy,
pref = df_2020$pref) |>
ggplot(aes(x = diff_energy, y = diff_income)) +
geom_point(color = "red") +
geom_text(aes(label = pref), vjust = 1.2, alpha = 0.5) +
geom_smooth(method = "lm") +
theme_bw()
#
# income と popratio
df_1120 |>
ggplot(aes(x = popratio, y = income)) +
geom_point(aes(color = year)) +
geom_text(aes(label = pref), vjust = 1.2, alpha = 0.5) +
geom_smooth(method = "lm") +
theme_bw()
#
# income と energy
df_1120 |>
ggplot(aes(x = energy, y = income)) +
geom_point(aes(color = year)) +
geom_text(aes(label = pref), vjust = 1.2, alpha = 0.5) +
geom_smooth(method = "lm") +
theme_bw()
#
