(Bing Image Creator で生成, プロンプト: Close up of purple Muscari flowers, background is blue sky and white clouds, photo)
の続きです。今回は、R で単回帰分析をしてみます。
2020年度のデータだけを使って回帰分析してみようと思いますので、まずは、2020年度だけのデータフレームを作成します。
説明変数を popratio: 15歳~64歳人口割合、被説明変数を income: 1人当たり県民総所得にして回帰分析をしてみようと思います。
まずは、散布図を描いて傾向を把握します。
geom_smooth() 関数で追加した回帰直線が右肩上がりです。popratio が大きい値ほど、income も大きい値になる、ということですね。
lm() 関数で回帰モデルを推定します。
income = -2859.33 + 99.91 * popratio + u という推定式が得られました。
popratio が 1 増えると income は 99.91 増える、という式です。popratio はパーセント表示で、income は千円単位なので、約10万円増える、ということですね。結構大きな影響なのではないでしょうか？
summary() 関数で model1 のサマリーを確認します。
R-squared が 0.3795 ということですので、popratio が income を説明する割合は 37.95% ということですね。
plot() 関数で、残差が popratio を関係しているかどうかをみてみます。
左の 2 つのグラフを見ると、残差の分散が均一ではないようです。
lmtest パッケージの bptest() 関数で確認してみます。
p-value が 0.005 よりも小さい値なので、分散が均一である、という帰無仮説を棄却できます。
単純に、残差の 2乗を popratio で回帰分析しても確認できます。
p-value が 0.0001254 とやはり 0.05 よりも小さい値ですね。
冒頭の散布図を見ると、一番右上にポツンと点があります。これは東京都の点です。
この東京都を除いて回帰分析したらどうなるでしょうか？
まず、東京都を除いたデータフレームを作成します。
lm() 関数で回帰モデルを推定します。
popratio の係数が 47.87 と半分くらいになりました。
summary() 関数をつかってモデルの詳細を確認します。
R-squared が 0.1938 と東京都を含んだモデルよりも小さい値になっていますね。
plot() 関数でモデルのグラフを確認します。
左の２つを見ると、東京都を含むモデルよりは、残差の分散は均一のように見えます。
lmtest パッケージの bptest() 関数で確認します。
p-value が 0.2329 なので、このモデルの残差が均一分散である、という帰無仮説を棄却できません。
単純に残差の 2乗を popratio で回帰分析する方法でも確認します。
p-value が 0.2488 と 0.05 よりも大きいので、東京都を除いた線形単回帰分析モデルの分散が不均一とは言えませんね。
あらためて、model1_ex_tokyo の係数を確認しましょう。
popratio の係数は 47.86549 です。15歳~64歳人口割合が 1ポイント上昇すると、１人当たり県民総所得は約4万8千円上昇する、という関係です。
15歳~64歳はいわゆる生産年齢人口ですので、この割合が大きいほど、県民総所得も大きくなる、というのは直感と同じかと思います。
今回は以上です。
次回は
です。
初めから読むには、
です。
今回のコードは以下になります。
#
# 2020年度だけのデータフレームを作成
df_2020 <- df |>
filter(year == "2020年度")
summary(df_2020)
#
# income と popratio の散布図
df_2020 |>
ggplot(aes(x = popratio, y = income)) +
geom_point() +
geom_smooth(method = "lm")
#
# income を popratio で説明するモデル
model1 <- lm(income ~ popratio, data = df_2020)
model1
#
# model1 のサマリー
summary(model1)
#
# model1 のプロット
par(mfrow = c(2, 2))
plot(model1)
par(mfrow = c(1, 1))
#
# Breush-Pagan テスト
lmtest::bptest(model1)
#
# Heteroskedasticity のチェック
lm(I(model1$residuals^2) ~ df_2020$popratio) |>
summary()
#
# 東京都を除いたデータフレームを作成
df_2020_ex_tokyo <- df_2020 |>
filter(pref != "東京都")
#
# 東京都を除いて回帰分析
model1_ex_tokyo <- lm(income ~ popratio,
data = df_2020_ex_tokyo)
model1_ex_tokyo
#
# model1_ex_tokyo のサマリー
summary(model1_ex_tokyo)
#
# model1_ex_tokyo のグラフ
par(mfrow = c(2, 2))
plot(model1_ex_tokyo)
par(mfrow = c(1, 1))
#
# Breuch-Pagan Test
lmtest::bptest(model1_ex_tokyo)
#
# model1_ex_tokyo の残差の2乗と popratio
lm(I(model1_ex_tokyo$residuals^2) ~ df_2020_ex_tokyo$popratio) |>
summary()
#
# model1_ex_tokyo の係数
coef(model1_ex_tokyo)
confint(model1_ex_tokyo)
#
