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

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

都道府県別の交通事故発生件数のデータの分析2 - R でデータをグラフにする。The 5NG (1-散布図, 2-線グラフ, 3-ヒストグラム, 4-箱ひげ図, 5-棒グラフ)

Bing Image Creator で生成: Close up of Kalanchoe flowers, background is blue sky, white clouds, long river in the forest, photo

 

www.crosshyou.info

の続きです。今回はデータを R でグラフにしてみます。

Chapter 2 Data Visualization | Statistical Inference via Data Science

 https://moderndive.com/v2/viz.html#FiveNG

このサイトにある5つのグラフを描こうと思います。

まずは scatterplots 散布図ですね。

交通事故と人口密度の散布図を描いてみます。

横軸が人口密度で、縦軸が交通事故です。色の濃淡は老年化指数を表しています。

明確な関連性は見られないですね。

The 5NG の2つ目は、lineargraphs です。線グラフですね。

年度と交通事故の線グラフを描きます。

ある年までは年々増加していますが、それから先は減少していることがわかります。

3つめの The 5NG は histograms ヒストグラムです。交通事故のヒストグラムを描いてみます。

500の当りが一番、度数が多いですね。右側の裾野が広い分布です。

The 5NG の4つ目は boxplots 箱ひげ図です。都道府県別の交通事故の箱ひげ図を描いてみます。

mutate(pref = reorder(pref, jiko, median)) のコードで pref を jiko の中央値の順に並び替えています。

静岡県が人口10万人当たり交通事故発生件数が多いんですね。意外でした。

さて、最後の The 5NG は barplots book 棒グラフです。交通事故の合計値を年度ごとに棒グラフにしてみます。

線グラフのところでもわかりましたが、交通事故は最近は減少傾向ですね。

今回は、

moderndive.com

Welcome to ModernDive (v2) | Statistical Inference via Data Science

というウェブサイト上の本の中の The 5NG を参考にして、5つの基本的なグラフ、散布図、線グラフ、ヒストグラム、箱ひげ図、棒グラフを描きました。

今回は以上です。

初めから読むには、

 

www.crosshyou.info

です。

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

#
# The 5NG 1 - scatter plots
# 交通事故と人口密度
df |> 
  ggplot(aes(x = mitsudo, y = jiko)) +
  geom_point(aes(color = elder))
#
# The 5NG 2 - linergraphs
# 交通事故と年度
df |> 
  group_by(year) |> 
  ggplot(aes(x = year, y = jiko)) +
  geom_line(aes(group = pref))
#
# The 5NG 3 - histograms
# 交通事故のヒストグラム
df |> 
  ggplot(aes(x = jiko)) +
  geom_histogram(color = "white")
#
# The 5NG 4 - boxplots
# 都道府県別の交通事故の箱ひげ図
df |> 
  mutate(pref = reorder(pref, jiko, median)) |> 
  ggplot(aes(x = jiko, y = pref)) +
  geom_boxplot(aes(group = pref))
#
# The 5NG 5 - barplots
# 年度ごとの交通事故の合計件数
df |> 
  group_by(year) |> 
  summarize(total_jiko = sum(jiko)) |> 
  ggplot(aes(x = year, y = total_jiko)) +
  geom_col()
#

 

都道府県別の交通事故発生件数のデータの分析1 - R にデータを取り込む

Bing Image Creator で生成: Close up Crassula ovata flowers, background is plain green grass field under blue sky, photo

今回からしばらくは、都道府県別の人口密度と交通事故発生件数のデータを分析してみようと思います。

政府統計の総合窓口(www.e-stat.go.jp)からデータを取得します。

可住地面積1平方キロメートル人口密度と人口10万人当たりの交通事故発生件数をデータとして選択しました。前回の分析でも使った老年化指数も選択しました。

CSVファイルにこのような形式でデータをダウンロードできました。

5行目に変数名を挿入しています。

mitsudo: 可住地面積1平方キロメートル当たり人口密度【人】

elder: 老年化指数(65歳以上人口/15歳未満人口を指数化)

jiko: 人口10万人当たり交通事故発生件数【件】

です。

tidyverse パッケージの読み込みをします。

read_csv() 関数でCSVファイルのデータを読み込みます。

glimpse() 関数を使って読み込んだデータの内容を確認します。

na.omit() 関数を使って、NAのある行を削除します。

summary() 関数を使って各変数の基本統計値を確認します。

year, pref の Length が2256ということは、このデータフレームは2256行のデータフレームです。

人口密度は、1番低いところで、1平方キロメートル当たり226人、1番高いところで9873人です。

老年化指数は、一番若いところで、19.1、一番老年のところで417.4です。

人口10万人当たりの交通事故発生件数は、一番少ないところで109.8件、多いところで1328.3件です。

summary() 関数では計算されてない、標準偏差と変動係数を算出します。

変動係数は、標準偏差 / 平均値 で計算します。

人口密度が一番バラツキが大きく、交通事故が一番バラツキが小さいことがわかります。

今回は以上です。

次回は、

www.crosshyou.info

です。

 

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

#
# tidyverse パッケージの読み込み
library(tidyverse)
#
# CSVファイルの読み込み
df_raw <- read_csv("koutsuujiko.csv",
                   skip = 4)
#
# glimpse()
glimpse(df_raw)
#
# NA行の削除
df <- na.omit(df_raw)
#
# summary()
summary(df)
#
# 標準偏差と変動係数
df |> 
  summarize(
    std_mitsudo = sd(mitsudo),
    std_elder = sd(elder),
    std_jiko = sd(jiko),
    cv_mitsudo = sd(mitsudo) / mean(mitsudo),
    cv_elder = sd(elder) / mean(elder),
    cv_jiko = sd(jiko) / mean(jiko)
  )
#

 

読書記録 - 「日本政治学史 丸山眞男からジェンダー論、実験政治学まで」 酒井 大輔 著 (中公新書)

この本の巻末には参考文献のリストがありますが、これが新書の参考文献リストしては破格の量で圧倒されます。

著者は大学の先生かと思ったら、国家公務員ということなので、公務員の仕事をしながらこんなにも大量の参考文献を読み、本書を書きあげたということに感動しました。

内容は、私がまったく政治学という分野への土地勘が無いので、正しく理解できたどうかはわかりませんが、第二次世界大戦が終わってから、それまでの大学内でのサークルのように乱立していた政治学の研究会が、日本政治学会という大学の垣根を超えた全国規模の学会が設立されて、だんだんと発展してきたということはわかりました。

政治学者は実際の政治団体にどこまで関与するのか、とか難しい問題だと思いました。

 

都道府県別の老年化指数と県内総生産額対前年増加率のデータの分析 - R で Fixed Effects Regression

Bing Image Creator で生成: Close up Schlumbergera flowers, background is high mountains and white clouds in the sky, photo

 

www.crosshyou.info

の続きです。

前回は、gdp_growth と elder_index を単回帰分析しました。今回は都道府県別の効果を考慮して回帰分析をしてみます。

www.econometrics-with-r.org

https://www.econometrics-with-r.org/10.3-fixed-effects-regression.html

このサイトを参考にして分析します。

Fixed Effects Regression は個々の都道府県特有の効果、これを固定効果(Fixed Effects)と名付けて、モデルに組み込んだものです。

参考ウェブサイトの図を拝借しますと、

今回の分析では、Yが gdp_growth: 県内総生産額対前年増加率で、X1 が elder_index: 老年化指数です。X はこの一つだけです。D2, D3, D4... D45, D46, D47 が pref: 都道府県の各県、ということです。

早速、やってみます。

結果をみてみます。

elder_index の係数は -0.009818 です。都道府県特有の固定効果を考慮してもなお、elder_index が大きいと gdp_growth は小さくなります。

参考サイトには、plm パッケージの plm() 関数を使ったやり方も記載されていますので、そちらのほうも実践してみます。

まずは、plm パッケージを読み込みます。

次に、plm() 関数でモデルの係数を推定します。

plm() 関数では、index = で都道府県と年と指定して、model = "within" で Fixed Effects Regression と指定しています。

結果をみてみます。

係数が -0.0098179 とlm() 関数での推定と一致しました。

この -0.0098179 が有意にゼロと違うかどうかを検定します。

AER パッケージを読み込みます。

coeftest() 関数で有意かどうか検定できます。

p値が0.1042 と 10%よりも少し大きいです。5%の有意水準では、0と有意な違いがあるとは言えないです。

老年化指数が大きいほうが県内総生産額対前年増加率は低いようだけど、自信をもって断言はできないな~、という感じでしょうか?

前回、前々回は infer パッケージを使って、シミュレーション・ベースでの検定をしました。今回は、infer パッケージは使えないので、for loop でシミュレーションしてみましょう。

1000回シミュレーションをするとして、1000回の結果を格納する箱を作ります。

for loop で1000回シミュレーションします。

結果をデータフレームに格納します。

keisu の平均値を確認しておきます。

実際の係数は、-0.0098179 なので、若干の違いはあります。シミュレーションの回数をもっと多くすると、もっと実際の係数の値に近づくでしょう。

95%信頼区間を求めます。

信頼区間がゼロを含んでいます。

グラフにして視覚化します。

赤い2つの垂線が95%信頼区間の水準です。緑の垂線がゼロの位置です。オレンジ色の垂線が推定された係数の位置です。緑の垂線が信頼区間の中にありますので、95%の信頼度では、推定値が0と違うとは言えないですね。

今回は以上です。

初めから読むには、

www.crosshyou.info

です。

今回のコードは

# https://www.econometrics-with-r.org/10.3-fixed-effects-regression.html
# Fixed Effect Model
#
# gdp_growth(it) = Beta1 * elder_index(it) + prefFixedEffect + u(it)
growth_fe_lm_mod <- lm(gdp_growth ~ elder_index + pref - 1,
                       data = df)
#
# 結果
growth_fe_lm_mod
#
# plm パッケージを読み込む
library(plm)
#
# plm() で Fixed Effects Regression
growth_fe_mod <- plm(gdp_growth ~ elder_index,
                     data = df,
                     index = c("pref", "year"),
                     model = "within")
#
# 結果
growth_fe_mod
#
# AERを読み込む
library(AER)
#
# growth_fe_modの係数
coeftest(growth_fe_mod,
         vcov. = vcovHC,
         type = "HC1")
#
# Bootstrap 法
# elder_index の係数を格納するベクトルオブジェクトを作成
keisu <- numeric(1000)
#
# for loopで1000回実行
set.seed(333)
for (i in 1:1000) {
    index <- sample(1:1000, replace = TRUE)
    model <- lm(gdp_growth ~ elder_index + pref - 1,
                data = df[index, ])
    keisu[i] <- model$coefficients[1]
}
#
# keisuからデータフレームを作成
boot_dist <- tibble(
  replicate = 1:1000,
  stat = keisu
)
boot_dist
#
# keisu の平均値
mean(keisu)
#
# 95%信頼区間
percentile_ci95 <- boot_dist |> 
  get_confidence_interval(level = 0.95,
                          type = "percentile")
percentile_ci95
#
#
# 視覚化
boot_dist |> ggplot(aes(x = stat)) +
  geom_histogram(color = "white") +
  geom_vline(xintercept = percentile_ci95$lower_ci, color = "red",
             linewidth = 1.5) +
  geom_vline(xintercept = percentile_ci95$upper_ci, color = "red",
             linewidth = 1.5) +
  geom_vline(xintercept = growth_fe_lm_mod$coefficients[1],
             color = "orange", linewidth = 1.5) +
  geom_vline(xintercept = 0, color = "green",
             linewidth = 1.5)
#

です。

 

都道府県別の老年化指数と県内総生産額対前年増加率のデータの分析6 - R で単回帰分析 (infer パッケージのワークフロー, Two numerical vars - SLR)

Bing Image Creator で生成: Close up of yellow Cymbidium flowers, background is green grass field and blue sky, photo

www.crosshyou.info

の続きです。今回は R で単回帰分析をしましょう。infer パッケージのワークフローに沿ってやってみます。

https://infer.netlify.app/articles/observed_stat_examples#two-numerical-vars---slr

このサイトを参考にしています。

はじめに、slope_hat を求めます。

散布図を描いてみます。

geom_smooth() 関数で回帰直線を追加しました。この傾きが有意に0と違うかどうかが問題です。

null distribution を生成します。

null distribution と slope_hat(-0.00719) をグラフにします。

-0.00719 ってほとんど 0 かと思いましたが、けっこう 0 とは離れていることがグラフにするとわかります。

p値を計算します。

p値は0.05よりも小さい0.014です。なので、elder_index と gdp_growth は5%の有意水準で関係性があることがわかります。

95%信頼区間も求めてみます。

95%信頼区間を求めます。

-0.0128 ~ -0.00148 となりました。0を含んでいません。

グラフ表示します。

95%信頼区間が0を含んでいないことがわかります。

今回の分析で、老年化指数が大きいほど、県内総生産額対前年増加率は低い、という関係があることがわかりました。

今回は以上です。

次回は、

www.crosshyou.info

です。

 

初めから読むには、

www.crosshyou.info

です。

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

#
# 傾きの算出
slope_hat <- df |> 
  specify(gdp_growth ~ elder_index) |> 
  calculate(stat = "slope")
slope_hat
#
#  散布図
df |> 
  ggplot(aes(x = elder_index, y = gdp_growth)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE)
#
# null distribution の生成
set.seed(222)
null_dist <- df |> 
  specify(gdp_growth ~ elder_index) |> 
  hypothesize(null = "independence") |> 
  generate(reps = 1000, type = "permute") |> 
  calculate(stat = "slope")
#
# null distribution と slope_hat
visualize(null_dist) +
  shade_p_value(obs_stat = slope_hat,
                direction = "two-sided")
#
# p値の計算
null_dist |> 
  get_p_value(obs_stat = slope_hat,
              direction = "two-sided")
#
# ブートストラップ法
set.seed(222)
boot_dist <- df |> 
  specify(gdp_growth ~ elder_index) |> 
  generate(reps = 1000, type = "bootstrap") |> 
  calculate(stat = "slope")
#
# 95%信頼区間
percentile_ci <- boot_dist |> 
  get_confidence_interval(level = 0.95,
                          type = "percentile")
percentile_ci
#
# ブートストラップ法と信頼区間
boot_dist |> 
  visualize() +
  shade_confidence_interval(endpoints = percentile_ci) +
  geom_vline(xintercept = 0.0, color = "red", linewidth = 1.5)
#

 

都道府県別の老年化指数と県内総生産額対前年増加率のデータの分析5- R で比率の差の検定 (infer パッケージのワークフローで Two categorical (2 level) variables)

Bing Image Creator で生成: Close up of chrysanthemum flowers, background is Orion in the night sky, photo

 

www.crosshyou.info

の続きです。今回は、2つのカテゴリカル変数があって、変数の比率の差の検定を infer パッケージのワークフローでやってみます。

https://infer.netlify.app/articles/observed_stat_examples#two-categorical-2-level-variables

が参考ページです。

まずは、2つのカテゴリカル変数を作成します。

一つめは、gdp_growth >= 0 なら positive, gdp_growth < 0 なら negative という変数です。

もう一つは、elder_index が平均よりも小さければ young, 平均よりも上ならば elder となるカテゴリカル変数です。

これで、2つのカテゴリカル変数が揃いました。table() 関数でクロス集計表をみてみます。

negative: gdp_growth < 0 のとき、elder は58、young は47とe lder のほうが多いです。

positive: gdp_growth >= 0 のとき、elder は125、young は146と elder のほうが多いです。

prop.table() 関数で比率をみてみます。

negative のとき、elder の比率は55%、positive のとき、elder の比率は46%です。

比率の差は、

9%ぐらいです。この9%の差は有意な差なのかを調べましょう、ということです。

infer パッケージのワークフローをする前に、せっかくですのでグラフを描いてみます。

positive の件数のほうがだんぜん多いことがわかります。

それでは、infer パッケージのワークフローに入ります。まずは、d_hat: 比率の差を求めます。

0.0911 さきほど求めた値と(あたりまえですが)同じです。

null distribution を生成します。

null distribution と d_hat のグラフを描いてみます。

赤い垂線のところが d_hat の 0.911 です。

p値を計算します。

p値は0.142と0.05よりも大きな値ですね。0.0911という比率の差は、0とは有意に違いがあるとは言えません。

95%の信頼区間も求めてみます。

ブートストラップ法で、1000個の d_hat を作ります。

get_confidence_interval() 関数で95%信頼区間を求めます。

グラフにします。

赤い線が0.0の位置の垂線です。95%の信頼区間がこの赤い垂線を含んでいるので、比率の差が0とは有意に異なるとは言えません。

今回は以上です。

次回は、

www.crosshyou.info

です。

 

初めから読むには、

www.crosshyou.info

です。

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

#
# gdp_growth がプラス成長ならpositive, マイナス成長ならnegative
# の変数を作る
df <- df |> 
  mutate(growth = if_else(gdp_growth >= 0, "positive", "negative"))
df
#
# elder_index が平均より下ならyoung, 上ならelderの変数を作る
df <- df |> 
  mutate(age = if_else(elder_index < mean(elder_index), "young", "elder"))
df
#
# growth と age のクロス集計
table(df$growth, df$age)
#
# negative, positive 別の elder, young の比率
prop.table(table(df$growth, df$age), margin = 1)
#
# 比率の差
0.5523810 - 0.4612546
#
# growth と age のグラフ
df |> 
  ggplot(aes(x = growth, fill = age)) + 
  geom_bar()
#
# infer パッケージのワークフローで比率の差(0.0911264)が有意に0と違うかどうか
# を調べる
# d_hat(比率の差)
d_hat <- df |> 
  specify(age ~ growth,
          success = "elder") |> 
  calculate(stat = "diff in props",
            order = c("negative", "positive"))
d_hat
#
# null distributionの生成
set.seed(321)
null_dist <- df |> 
  specify(age ~ growth, success = "elder") |> 
  hypothesize(null = "independence") |> 
  generate(reps = 1000, type = "permute") |> 
  calculate(stat = "diff in props", order = c("negative", "positive"))
null_dist
#
# null distribution のグラフ
visualize(null_dist) +
  shade_p_value(obs_stat = d_hat, direction = "two-sided")
#
# p値の計算
null_dist |> 
  get_p_value(obs_stat = d_hat, direction = "two-sided")
#
# ブートストラップ法
set.seed(333)
boot_dist <- df |> 
  specify(age ~ growth, success = "elder") |> 
  generate(reps = 1000, type = "bootstrap") |> 
  calculate(stat = "diff in props",
            order = c("negative", "positive"))
boot_dist
#
# 信頼区間の算出
percentile_ci <- boot_dist |> 
  get_confidence_interval(level = 0.95, 
                          type = "percentile")
percentile_ci
#
# 信頼区間とブートストラップ法の d_hat
visualize(boot_dist) +
  shade_confidence_interval(endpoints = percentile_ci) +
  geom_vline(xintercept = 0, color = "red", linewidth = 2)
#

 

都道府県別の老年化指数と県内総生産額対前年増加率のデータの分析4 - R の infer パッケージを利用して、平均値の検定をする。

Bing Image Creator で生成: Close up of yellow Chimonanthus praecox flowers, background is white clouds in the sky, photo 

www.crosshyou.info

の続きです。

これまでの結果、沖縄県が一番、老齢化指数が低く、秋田県が一番、老齢化指数が高いことがわかりました。今回は、R の infer パッケージを利用して、二つの県の老年化指数が統計的に有意に違うのかを確認します。

まず infer パッケージを library() 関数を使って読み込みます。

秋田県と沖縄県の elder_index: 老齢化指数の平均値を調べるためのデータフレームを作成します。

2つの県の平均値を確認します。

秋田県の平均値は、328で、沖縄県の平均値は、116なので差は、212です。

infer パッケージのワークフローに従って差の平均値を求めます。

https://infer.netlify.app/articles/observed_stat_examples#one-numerical-variable-one-categorical-2-levels-diff-in-means

を参考にしています。

212となっています。

次に null distribution を求めます。これは、pref と elder_index をバラバラに組み合わせて、秋田県と沖縄県の平均値の差を計算することを何回も繰り返すことです。

1回目のバラバラの組み合わせでは、差は-3でした。2回目は-55.6、3回目は6.58となっています。

実際の差の212と、この null distibution を視覚化して比較します。

ヒストグラムが null distribution です。秋田県と沖縄県の elder_index の値をバラバラにシャッフルして差を計算していますので、0を中心とした分布になっています。赤い垂線が実際の差の212の線です。elder_index の値は秋田県、沖縄県と全く関係がなければ212という大きな差は起こりえない感じです。

p値を確認します。

p値は0です。

以上が infer パッケージのワークフローを使ったシミュレーションベースの統計的推定です。

昔からある理論ベースの t検定もやってみます。t.test() 関数を使います。

akita, okinawa とうオブジェクト名で、それぞれの elder_index のベクトルを作成して、t.test() 関数で処理しています。

p値は 1.966e-07 とほぼ0です。差の 95%信頼区間は、181 ~ 243 です。

infer パッケージのワークフローでも信頼区間は求めることができます。

やってみましょう。

https://infer.netlify.app/articles/observed_stat_examples#one-numerical-variable-one-categorical-2-levels-diff-in-means-1

を参考にしています。

d_hat は既に求めていますので、ブートストラップ法による 差の distribution を求めます。

get_ci() 関数で信頼区間を求めます。

グラフにしましょう。

緑の垂線が、infer パッケージのワークフローで求めたシミュレーションベースの95%信頼区間で、赤い垂線が、t.test() 関数で求めた理論ベースの95%信頼区間です。

シミュレーションベースの信頼区間のほうが幅が狭いですね。

数値を確認しておきます。

186 ~ 238 です。

今回は以上です。

次回は、

www.crosshyou.info

です。

初めから読むには、

www.crosshyou.info

です。

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

#
# inferパッケージを読み込む
library(infer)
#
# 沖縄県と秋田県のelder_index: 老齢化指数の比較
df2 <- df |> 
  filter(pref == "秋田県" | pref == "沖縄県") |> 
  mutate(pref = factor(pref, levels = c("秋田県", "沖縄県")))
df2
#
# 秋田県と沖縄県の平均値
df2 |> 
  group_by(pref) |> 
  summarize(avg = mean(elder_index))
#
# 平均値の差
d_hat <- df2 |> 
  specify(elder_index ~ pref) |> 
  calculate(stat = "diff in means",
            order = c("秋田県", "沖縄県"))
d_hat
#
# null distribution
set.seed(123)
null_dist <- df2 |> 
  specify(elder_index ~ pref) |> 
  hypothesize(null = "independence") |> 
  generate(reps = 1000, type = "permute") |> 
  calculate(stat = "diff in means",
            order = c("秋田県", "沖縄県"))
null_dist
#
# null distribution と d_hat
null_dist |> 
  visualize() + 
  shade_p_value(obs_stat = d_hat, direction = "two-sided")
#
# p値
null_dist |> 
  get_p_value(obs_stat = d_hat, direction = "two-sided")
#
# 理論ベースの t検定
akita <- df2 |> 
  filter(pref == "秋田県") |> 
  pull(elder_index)
akita
#
okinawa <- df2 |> 
  filter(pref == "沖縄県") |> 
  pull(elder_index)
okinawa
#
t.test(akita, okinawa)
#
# ブートストラップ法
set.seed(123)
boot_dist <- df2 |> 
  specify(elder_index ~ pref) |> 
  generate(reps = 1000, type = "bootstrap") |> 
  calculate(stat = "diff in means",
            order = c("秋田県", "沖縄県"))
boot_dist
#
# 信頼区間
percentile_ci <- boot_dist |> 
  get_ci(type = "percentile")
#
# boot_dist と信頼区間のグラフ
visualize(boot_dist) +
  shade_confidence_interval(endpoints = percentile_ci) +
  geom_vline(xintercept = 181, color = "red", linewidth = 2) +
  geom_vline(xintercept = 243, color = "red", linewidth = 2)
#
# 信頼区間の数値
percentile_ci
#