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

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

都道府県別の電子レンジ所有数量のデータの分析5 - Rで線形重回帰分析をする

www.crosshyou.info

の続きです。今回は回帰分析をしてみます。

lm()関数でモデルを作ります。

このモデルは、

電子レンジ所有数量 = β0 + β1 x 年ダミー + β2 x 人口密度の対数 + β3 x コンビニエンスストアの数 + β4 x 単身世帯の割合 + β5 x 1人当たり県民所得 + u

です。β0が切片項で、uが誤差項です。

summary()関数で結果をみます。

p-value: 2.543e-09と表示されていますので、全体としてこのモデルは有意なモデルです。個々の係数のp値もどれも0.05以下なので、どの係数も有意です。

係数が1.161e+03のようにわかりずらいので、broomパッケージのtidy()関数で見やすく出力します。

これで係数がわかりやすくなりました。整理すると、

電子レンジ所有数量 = 1161 + 22.2 x 2014年度ダミー - 14.4 x 人口密度の対数 - 1.29 x コンビニエンスストアの数 - 3.19 x 単身世帯の割合 + 0.0373 x 1人当たり県民所得 + u 

です。

人口密度の対数やコンビニエンスストアの数など他の条件が変わらないとして、2014年度のほうが2004年度よりも、電子レンジ所有数量は22.2台多いということですね。

人口密度の対数が-14.4ということは、他の条件が変わらずに人口密度だけが1%上昇すると、電子レンジ所有数量が14.4台減る、ということです。

他の条件が変わらずに、コンビニエンスストアの数が1所増えると、電子レンジ所有数量は1.29台減る、ということです。

他の条件が変わらずに、単身世帯の割合が1%上昇すると、電子レンジ所有数量は3.19台減る、ということです。

他の条件が変わらずに、1人当たり県民所得が1千円増えると、電子レンジ所有数量は0.0373台増える、ということです。

残差が均一分散しているかどうかのチェックをlmtestパッケージのbptest()関数で行います。

p-value = 0.3508と0.05よりも大きな値なので、残差は不均一分散ではないようです。

このモデル式から逆算される電子レンジ所有数量をpredict()関数で求めます。

lm_mod1の列が今回のモデルから逆算した電子レンジ所有数量の推定値です。

グラフにして、実際の値と推定値の様子をみてみます。

直線の下にある都道府県はモデルの推計値よりも実際の電子レンジ所有数量が少ないところです。もう少し見やすく、2014年度だけにしてグラフにしてみます。

北日本が実際の所有が多いとか、九州は少ないとかなどの地理的な傾向は無いようですね。

最後にRMSE(Root Mean Square Error)を計算します。

RMSEは21.6です。だいだい平均すると、21.6台ぐらいの誤差があるということですね。電子レンジ所有数量の平均値は1048台なので、約2%ほどの誤差です。結構、いい感じな推定値ではないでしょうか?

今回は以上です。

次回は、

www.crosshyou.info

です。

 

はじめから読むには、

www.crosshyou.info

です。

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

#
# 最初のモデル
lm_mod1 <- lm(microwave ~ year + l_mitsudo + conve + single + income,
   data = df)
#
# モデルのサマリー
summary(lm_mod1)
#
# tidy()で出力
broom::tidy(lm_mod1)
#
# Hetroskedasticityのチェック
library(lmtest)
bptest(lm_mod1)
#
# dfにモデルの予測値を追加
df <- df |> 
  mutate(lm_mod1 = predict(lm_mod1, newdata = df))
df
#
# 実際の値と予測値の散布図
df |> 
  ggplot(aes(x = lm_mod1, y = microwave)) +
  geom_point(aes(color = year)) +
  geom_text(aes(label = pref), hjust = 1.1) +
  geom_abline() +
  theme_bw()
#
# 2014年度だけのグラフ
df |> 
  filter(year == "2014年度") |> 
  ggplot(aes(x = lm_mod1, y = microwave)) +
  geom_point() +
  geom_text(aes(label = pref), hjust = 1.1) +
  geom_abline() +
  theme_bw()
#
# RMSEの計算
df |> 
  summarize(
    RMSE = sqrt(mean*1,
    mean = mean(microwave)
  ) |> 
  mutate(RMSE_mean = RMSE / mean * 100)
#

 

(冒頭の画像は、Bing Image Creator で生成しました。プロンプトは、Long wide view of landscape, from the birds view, small river runs through green grass fields and wild colorful flowers, photo. です。)

 

*1:microwave - lm_mod1)^2