の続きです。今回は、第1次産業、第2次産業、第3次産業に従事している人口比率、one, two, threeを変数に加えて回帰分析をしてみます。
はじめに、都市別のone, two, threeの平均値のデータフレームを作成します。
このデータフレームを前回作成したデータフレームと合体します。
crime: 刑法犯認知件数をone, two, three で回帰分析してみます。
one, two, three どれもp値は0.05よりも小さく有意です。どの係数の符号もマイナスなので、就業者比率が高いほど、刑法犯認知件数が多いという結果ですね。
不思議な結果です。
モデルの予測値と実際の値の散布図を描いてみます。
このモデルの決定係数(R2)は0.3175です。実際のcrimeが大きなところはモデルの予測値とは大きく離れています。
onw, two, three に加えて人口や、人口密度、年度、都道府県を加えて回帰分析をしてみます。
summary()関数で結果をみてみます。
(中略)
人口などその他の説明変数を加えても、one, two, three の係数の符号はマイナスで、p値も0.05よりも小さいので有意です。決定係数(R2)は0.8731になりました。
mitsudo:人口密度は有意ではないですね。
モデルの予測値と実際の値の散布図を描いてみましょう。
crimeが多いところの予測値が改善されていることがわかります。
人口や調査年度、都道府県をコントロールしてもone, two, threeは有意で、しかも符号の係数はマイナスでした。
数式どおりに解釈すると、第1次、第2次、第3次産業の比率が高いほうが刑法犯認知件数が多いということですね。
確認のため、コントロール変数と第1次産業、コントロール変数と第2次産業、コントロール変数と第3次産業のモデルを試してみます。
update()関数で、twoとthreeを除外したモデルを作ります。
summary()関数で結果をみます。
(後略)
oneは係数の符号はマイナスで変わらないですが、p値が0.543と0.05よりも大きくなり、有意な変数ではなくなりました。そのかわり、mitsudo:人口密度が有意な変数になりました。
twoだけのモデルはどうでしょうか？
(後略)
towもoneと同じように係数の符号はマイナスですが、p値は0.297と0.05よりも大きく有意な変数ではなくなりました。これもmitsudoが有意な変数になりました。
threeだけのモデルも確認します。
threeも有意な変数ではなくなりました。これもmitsudoが有意な変数になりました。
one, two, threeは単独では有意な変数ではなく、かわりにmitsudoが有意な変数になりました。
one, two, three, mitsudoの相関係数マトリックスをみてみます。
人口密度は、第3次産業の比率と正の相関で、第1次、第2次とは負の相関です。そして、第2次産業と第3次産業は-0.96とかなり強い逆相関ですね。
oneとtwoとコントロール変数のモデルを作ってみます。
one, twoともに有意ではないですね。
one + two + threeを合計した変数を調べてみます。
まずは分布からです。
allという名前でone + two + threeの変数を作りました。大半は95%以上ですが90%以下の都市もありますね。
このallとその他のコントロール変数で回帰分析してみます。
(中略)
allの係数は-3.5でp値は < 2e-16となっていますので有意な変数です。
産業従事比率が高いほど、刑法犯認知件数が多いという解釈ですね。これならなんとなくわかります。
今回は以上です。
はじめから読むには、
です。
今回のコードは以下になります。
#
# 都市別のone, two, threeの平均値のデータフレームを作成
df_123 <- df_raw |>
group_by(code) |>
summarize(
one = mean(one, na.rm = TRUE),
two = mean(two, na.rm = TRUE),
three = mean(three, na.rm = TRUE)
)
df_123
#
# df_mitsudoとdf_123を合体
df_123 <- df_mitsudo |>
inner_join(df_123, by = "code")
df_123
#
# crimeをone, two, threeで回帰分析してみます。
one23_mod <- lm(crime ~ one + two + three, data = df_123)
#
# 結果
summary(one23_mod)
#
# 実際のcrimeと予測値の散布図
df_123 |>
mutate(estimate = predict(one23_mod)) |>
ggplot(aes(x = crime, y = estimate)) +
geom_point(aes(color = year)) +
geom_abline(color = "red") +
theme_minimal()
#
# crimeをyear, mitsudo, pop, pop^2, l_pop, l_pop^2, pref, one, two, three
# で回帰分析
one23_mod2 <- lm(crime ~ one + two + three + pop + I(pop^2) + l_pop +
I(l_pop^2) + mitsudo + year + pref, data = df_123)
#
# 結果
summary(one23_mod2)
#
# 実際のcrimeと予測値の散布図
df_123 |>
mutate(estimate = predict(one23_mod2)) |>
ggplot(aes(x = crime, y = estimate)) +
geom_point(aes(color = year)) +
geom_abline(color = "red") +
theme_minimal()
#
# コントロール変数と第1産業
one_mod <- update(one23_mod2, . ~ . -two - three)
#
# 結果
summary(one_mod)
#
# コントロール変数と第2次産業
two_mod <- update(one23_mod2, . ~ . - one - three)
#
# 結果
summary(two_mod)
#
# コントロール変数と第3次産業
three_mod <- update(one23_mod2, . ~ . - one - two)
#
# 結果
summary(three_mod)
#
# one, two, three, mitsudoの相関係数マトリックス
df_123 |>
select(one, two, three, mitsudo) |>
cor()
#
# コントロール変数と第1次産業と第2次産業
one_two_mod <- update(one23_mod2, . ~ . - three)
#
# 結果
summary(one_two_mod)
#
# one + two + three
df_123 <- df_123 |>
mutate(all = one + two + three)
df_123 |>
ggplot(aes(x = all)) +
geom_histogram(color = "white") +
theme_minimal()
#
# allとコントロール変数
all_mod <- lm(crime ~ all + pop + I(pop^2) + l_pop + I(l_pop^2) +
mitsudo + year + pref, data = df_123)
#
# 結果
summary(all_mod)
#
(冒頭の画像は、Bing Image Creator で生成しました。プロンプトは Natural, beautiful landscape of global field, close up of white Lily and yellow Lily and red Lily flowers, photo です。)
