学習のまとめと例題のデータを使ってRで回帰診断図を出してみました。

学習のまとめ

回帰診断を使う理由

線形回帰モデルは、

外れ値がある場合
誤差項の
- 独立性が満たされない場合
- 等分散性が満たされない場合
- 正規性の仮定が満たされない場合

について、最小二乗法による予測の結果には誤解を生じさせる可能性がある。これらを判断するために回帰診断(regression diagnostics)を行う。

どんな手法を用いて診断するのか

調べたいこと	手法
外れ値	残差プロット、てこ比、Cookの距離
独立性	残差プロット、自己相関、DW比
等分散性	残差プロット、標準化残差の絶対値の平方根プロット*1
正規性	残差プロット、正規Q-Qプロット

回帰診断図

例1の回帰診断図

47都道府県の乗用車所有率に関する回帰診断図

解説・説明についてはワークブックを参照。図に現れる赤い線は図中の点に対するLOESS(LOWESS?)、局所回帰らしい。点たちの大まかな傾向を読み取ることができる。

問17.1の回帰診断図

44都道府県の乗用車所有率に関する回帰診断図

データの取得

さて図の出典を見ると、

都道府県・市区町村のすがた(社会・人口統計体系) 2017年度
都道府県別・車種別自動車保有台数(軽自動車含む) 令和元年8月末現在

と書かれている。ググればわかるのだが、データを探すのも加工するのも大変面倒である。細かい注意点を備忘録がてら記載しておく。

「車種別自動車保有台数(軽自動車含む)」については、最新の月末データもしくは各年度の年度末のデータにしかオンラインではアクセスできないっぽいので、ワークブックと完全同一のデータを用いることはおそらく不可能と思われる。先駆者の記事『統計検定準1級対応統計学実践ワークブック』をR, Pythonで解く～第17章回帰診断法～ - Qiitaでは令和二年２月現在のデータを用いたと記載あり。
各変数の単位について言及すると、
- 乗用車台数の単位は[台]を用いる
- 人口は[人]単位であって、[万人]単位ではない
  - データ元では[万人]単位だったりするので、変換を忘れずに
- 可住地面積は[ヘクタール]単位であって、[平方キロメートル]単位ではない
  - データ元では[平方キロメートル]単位だったりするので、変換を忘れずに
先駆者の記事では『フォーマットが複雑なので手作業で準備した』と書いてある通り、フォーマットが複雑(といっても、集計値をExcelにベタ貼りした人間にはやさしく機械が読み取るにはちょっと面倒なレベル)。Excel的な表計算ソフトがあるならそちらを使って手作業で加工するのが断然良い。
- ただし、データラングリング・前処理の練習としては全てRで加工する・全てPythonで加工するのも良いと思う
車種別自動車保有台数(軽自動車含む)は都道府県単位での集計値ではなく、地方運輸局x支局 + 沖縄という単位での集計値となっている
- 北海道と沖縄の取り扱いに注意。いわゆるcase式やdecode式を用いて都道府県の列を派生させる必要がある

今回利用したデータを自分が二次加工して作ったデータ、をgistに載っけておきます。

gist.github.com

ソースコード

前提：各データはURLへアクセスの上、Excelファイルをローカルに落とすこと

library(tidyverse)
library(readxl)

# 都道府県・市区町村のすがた(社会・人口統計体系)
# データセット情報
# 社会・人口統計体系 / 統計でみる都道府県のすがた2019 / 社会生活統計指標
# https://www.e-stat.go.jp/stat-search/files?page=1&layout=datalist&toukei=00200502&tstat=000001124677&cycle=0&tclass1=000001124955&stat_infid=000031788971&tclass2val=0
# データを見ると2017年度のデータが入っているのでこれを利用する
path1 <- "/Users/axjack/Downloads/a201.xls"
dat1 <- read_xls(path1)
dat1 |> 
  slice(13:59) |> 
  select(`...8`, `...9`, `1`, `13`) |> 
  rename(
    PrefNumber = `...8`
    , PrefName = `...9`
    , Pop = `1`
    , p.density = `13`
  ) |> 
  mutate(Pop = as.numeric(Pop) * 10000 ) |> # 万人を人に変換
  mutate(p.density = as.numeric(p.density) * 0.01 ) |>  #平方キロメートルをヘクタールに変換
  mutate(PrefName = str_remove_all(PrefName,"[都府県]$")) -> dat11


# 都道府県別・車種別自動車保有台数
# 平成31年(2019年)
# https://www.airia.or.jp/publish/statistics/number.html
path2 <- "/Users/axjack/Downloads/r5c6pv000000mdpm.xlsx"
dat2 <- read_xlsx(path2)
dat2 |> select(1:3) |> slice(-(1:3)) -> dat21
names(dat21) <- paste0('col',1:3)
dat21 |> 
  mutate(col1 = str_remove_all(col1,'\\s')) |> 
  mutate(col2 = str_remove_all(col2,'\\s')) |> 
  replace_na(list(col1 = '',col2 = '')) |> 
  filter(col2 != '計') |> filter(col1 != '合計') |> 
  mutate(PrefName = case_when(
    col2 %in% c('札幌','函館','旭川','室蘭','釧路','帯広','北見') ~ '北海道',
    col1 == '沖縄' ~ '沖縄',
    TRUE ~ col2
  )) |> 
  mutate(col3 = as.numeric(col3)) |> 
  group_by(PrefName) |> summarise(NumCar = sum(col3) ) -> dat22


# 結合
dat22 |> left_join(dat11) |> mutate(c.ratio = NumCar/Pop ) |> 
  select(PrefNumber, PrefName, Pop, NumCar, p.density, c.ratio) |>
  arrange(PrefNumber) -> dat

#write_csv(dat,"tjw_ch17_dat.csv")
#dat <- read_csv("https://gist.githubusercontent.com/axjack/84505a444f76436229f456fe915e44ff/raw/a8f41565a2986f7eedb9c1c9ead70ace388b7492/tjw_ch17_dat.csv")


# 重回帰分析
## 47都道府県の場合
dat.lm <- lm(c.ratio ~ sqrt(p.density), dat)
summary(dat.lm)

par(mfrow = c(2,2))
plot(dat.lm)


# 44都道府県の場合
## 東京・大阪・神奈川を抜く
dat |> filter( !(PrefName %in% c('東京','大阪','神奈川')) ) -> dat.44
dat.44.lm <- lm(c.ratio ~ sqrt(p.density), dat.44)
summary(dat.44.lm)

par(mfrow = c(2,2))
plot(dat.44.lm)