の続きです。今回はcrime: 刑法犯認知件数と人口の関係をみてみます。
まず、人口だけのデータフレームを作ります。
このデータフレームに前回作ったcrimeだけのデータフレームを合体させます。inner_join()関数を使います。yearとcodeをキーにして一致させます。
札幌市は188万人で、人口千人当たりの刑法犯認知件数は17.0件とわかりますね。
人口のようなデータは対数変換しておくといいと、良く言われますので対数変換したデータも追加しておきます。
popとcrimeの散布図、l_popとcrimeの散布図を描いてみましょう。
L字型の分布でしょうか？
crimeをpop, l_pop, pop^2, l_pop^2で回帰分析してみます。
lm()関数でモデルを作ったら、summary()関数でモデルをみます。
モデル全体のp-valueは < 2.2e-16となりました。全ての係数でp値は0.05でした。
決定係数は0.6392で調整済み決定係数は0.6366でした。
モデルの予測値と実際のcrimeの値を比べてみます。
crimeが50以下の件数の少ないところはあんまりうまく予測できていないような感じですね。
この部分を拡大してみます。corrd_cartesian()関数を使います。
やっぱり、予測値と実際の値には違いがありますね。
yearとprefも加えた線形回帰モデルを作ってみましょう。
summary()関数で結果をみてみます。
決定係数が0.812に上がっています。
このモデルの予測値と実際のcrimeの値の散布図を描いてみます。
左下のエリアもそこそこフィットしている感じですね。この部分を拡大してみましょう。
だいぶ改善されていることがわかります。
あらためて、2番目の回帰モデルの係数をみてみます。
popの係数は -1.264e-04, I(pop^2)の係数は 1.149e-11 ということなので、
popは 1.264e-04 / (2 * 1.149e-11) = 5500435人になるまではpopが増えるとcrime は減少し、そこから先はpopが増えるとcrimeは減少する、ということです。
l_popの係数は -8.414e+02, I(l_pop^2)の係数は 3.444e+01 ということなので、
l_popは 8.414e+02 / (2 * 3.333e+01) = 16.62226 になるまでは、l_popが増えるとcrimeは減少して、そこから先はl_popが増えるとcrimeは減少するということです。
ただし、popが変化すれば、l_popも変化するから解釈は難しいですね。
curve()関数で、popの動きとcrimeの動きの曲線を描いてみましょう。
基本的には、popが増えるとcrimeも増える、ということですね。crimeは人口千人当たりの件数ですから、大都会ほど犯罪件数が多いということですね。
今回は以上です。
次回は、
です。
はじめから読むには、
です。
今回のコードは以下になります。
#
# year, code, popのデータフレームを作る
df_pop <- df_raw |>
select(year, code, pop) |>
na.omit()
df_pop
#
# df_popにdf_crimeを合体させる
df_pop <- df_pop |>
inner_join(df_crime, by = join_by(year, code))
df_pop
#
# popの対数変換値を作る
df_pop <- df_pop |>
mutate(l_pop = log(pop)) |>
relocate(year, code, pop, l_pop)
df_pop
#
# popとcrime
df_pop |>
ggplot(aes(x = pop, y = crime)) +
geom_point() +
theme_minimal()
#
# l_popとcrime
df_pop |>
ggplot(aes(x = l_pop, y = crime)) +
geom_point(color = "red") +
theme_minimal()
#
# crimeをpop, l_pop, crime^2, l_pop^2で回帰分析
pop_mod <- lm(crime ~ pop + I(pop^2) + l_pop + I(l_pop^2),
data = df_pop)
#
# 回帰分析結果を見る
summary(pop_mod)
#
# 予測値と実際のcrimeの散布図
df_pop |>
mutate(estimate = predict(pop_mod)) |>
ggplot(aes(x = crime, y = estimate)) +
geom_point() +
geom_abline(color = "red") +
theme_minimal()
#
# 予測値と実際のcrimeの散布図
# 左下の部分を拡大
df_pop |>
mutate(estimate = predict(pop_mod)) |>
ggplot(aes(x = crime, y = estimate)) +
geom_point() +
geom_abline(color = "red") +
coord_cartesian(xlim = c(0, 50), ylim = c(0, 40)) +
theme_minimal()
#
# year, prefも加えたモデル
pop_mod2 <- lm(crime ~ pop + I(pop^2) + l_pop + I(l_pop^2) +
year + pref, data = df_pop)
#
# 結果を見る
summary(pop_mod2)
#
# 予測値と実際の値の散布図
df_pop |>
mutate(estimate = predict(pop_mod2)) |>
ggplot(aes(x = crime, y = estimate)) +
geom_p5oint() +
geom_abline(color = "red") +
theme_minimal()
#
# 左下の部分を拡大
df_pop |>
mutate(estimate = predict(pop_mod2)) |>
ggplot(aes(x = crime, y = estimate)) +
geom_point() +
geom_abline(color = "red") +
coord_cartesian(xlim = c(0, 50), ylim = c(0, 50)) +
theme_minimal()
#
# pop_mod2の係数の確認
pop_mod2
#
# popの変曲点
1.264e-04 / (2 * 1.149e-11)
#
# l_popの変曲点
8.414e+02 / (2 * 3.333e+01)
#
# popとcrimeの関係
curve(5.185e+03 - 1.264e-04 * x + 1.1492e-11 + x^2 - 8.414e+02 * log(x) +
3.444e+01 * (log(x))^2,
from = min(df_pop$pop),
to = max(df_pop$pop))
#
(冒頭の画像は、Bing Image Creator で生成しました。プロンプトは Landscape of red Clivia flowers, full of flowers on the grass field, some yellow Dandelion flowers, under the blue sky, Photo です。)
