ある意味，インパクトのあるデータだったのか
> 大阪府西成区（P10の最小値）のヤバさが数字ではっきり出ちゃってますね
> 昨日のセンター試験数学「男の市区町村別平均寿命の箱ひげ図」で最小値が極端に小さいP10は大阪府（大阪市西成区73.5）か
> それはそうと、大阪の異常に平均寿命の短い市町村区はなんなんでしょう...
のようないくつかの反応があるようですが，二番目に衝撃的（？）だった神奈川県については誰も何も言っていないのかな？
神奈川県での平均寿命ワーストは「神奈川県川崎市川崎区 男 78.2」，「神奈川県横浜市中区 男 78.5」ですね。第三位の「神奈川県横浜市南区 男 80」よりも格段に低いです。
ある意味，大阪市西成区と同じような特性を持っているところかもしれません（曖昧に述べておきます）。
かなり前に，簡易宿泊所の火災で多くの人命が犠牲になったというニュースが記憶に残っています。
xlrd は，セルごとに文字列か数値かを判断して読み込んでくれたり，R では空セルも列のタイプによって NA だったり だったりだが，この関数では空セルは '' になる。
というように，R での読み込み関数より扱いやすい。
# https://www.e-stat.go.jp/stat-search/file-download?statInfId=000031693271&fileKind=0
# 図表データのダウンロード
# ファイル名: ckts-lifetable2015.xlsx
# 生命表1〜生命表7 が該当
# 7 列目が「平均余命」の 1 行前の 2 列目が「都道府県名+市区町村名」
# 3，25 行目先の 7 列目が目的値（男，女）
import pandas as pd
import xlrd
seireisi = ['札幌市', '仙台市', 'さいたま市', '千葉市', '東京都区部',
'横浜市', '川崎市', '相模原市', '新潟市', '静岡市', '浜松市',
'名古屋市', '京都市', '大阪市', '堺市', '神戸市', '岡山市',
'広島市', '北九州市', '福岡市', '熊本市']
input_file = 'ckts-lifetable2015.xlsx'
output_file = 'all-python.csv'
wb = xlrd.open_workbook(input_file)
sheet_names = wb.sheet_names()[1:]
field0 = []
field1 = []
field2 = []
field3 = []
for sheet_name in sheet_names:
sheet = wb.sheet_by_name(sheet_name)
n = sheet.nrows
print(f'処理中:{sheet_name}, {n} 行')
data = [sheet.row_values(row) for row in range(sheet.nrows)]
for i in range(2, n, 49):
if data[i][6] != '平均余命':
raise Exception('フィールド 6 が "平均余命" ではない')
name = data[i - 1][1]
name2 = name.split(' ')
name2_length = len(name2)
if name2_length == 1 or\
(name2_length == 2 and name2[1] in seireisi) or\
data[i + 3][6] == '…':
print(f'{name2} は除去')
else:
male = data[i + 3][6]
female = data[i + 25][6]
field0.append(name2[0])
field1.append(''.join(name2[1:]))
field2.append(male)
field3.append(female)
df = pd.DataFrame({'都道府県': field0, '市区町村': field1, '男': field2, '女': field3})
df.to_csv(output_file, index=False)
print(len(field0)) # 1888 件（福島県の 8 町村を除いた件数）
「素の R で，データスクレイピングとプロット（その1）」で使用した Excel ファイルは *.xls 形式で，かつ 3 段組みの表になっており処理がちょっと面倒であった。
これとは別の，市区町村別の生命表のエクセルファイルが同じく政府統計ポータルサイトからダウンロードできる。
https://www.e-stat.go.jp/stat-search/file-download?statInfId=000031693271&fileKind=0
ダウンロードしたファイルは，ckts-lifetable2015.xlsxである。
拡張子が .xlsx なので，openxlsx パッケージで入力できるので，readlx::read_excel よりストレスフリーである。
この Excel ファイルはしっかりとしたひな形にしたがって作成されているので，プログラムで例外処理をさせることなく，簡単に必要なデータを取り出せる。
「素の R で，データスクレイピングとプロット（その1）」で示したプログラムより遙かに短くて済む。
# https://www.e-stat.go.jp/stat-search/file-download?statInfId=000031693271&fileKind=0
# 図表データのダウンロード
# ファイル名: ckts-lifetable2015.xlsx
# 生命表1〜生命表7 が該当
# 7 列目が「平均余命」の 1 行前の 2 列目が「都道府県名+市区町村名」
# 3，25 行目先の 7 列目が目的値（男，女）
# install.packages("openxlsx")
library(openxlsx)
# 文字列を数値化。数値以外の場合はメッセージを出力し NA を返す。
get_value = function(txt, str) {
if (grepl("[^0-9.]", str)) {
cat("非数値データ：", txt, "\n", str, "を NA に変換\n")
return(NA)
}
return(as.numeric(str))
}
seireisi = c("札幌市", "仙台市", "さいたま市", "千葉市", "東京都区部",
"横浜市", "川崎市", "相模原市", "新潟市", "静岡市", "浜松市",
"名古屋市", "京都市", "大阪市", "堺市", "神戸市", "岡山市",
"広島市", "北九州市", "福岡市", "熊本市")
con = file("all3.csv", "w") # 結果の出力ファイル
count = 0
for (page in 1:7) { # 生命表1〜生命表7
data = read.xlsx("ckts-lifetable2015.xlsx", sheet = page + 1)
for (i in 1:nrow(data)) {
if (!is.na(data[i, 7]) && data[i, 7] == "平均余命") {
str1 = data[i - 1, 2] # 自治体名
str2 = data[i + 3, 7] # 男
str3 = data[i + 25, 7] # 女
name = unlist(strsplit(str1, " "))
if (length(name) == 1 || length(name) == 2 && name[2] %in% seireisi) {
cat("除外：", str1, "\n") # 除外した自治体名を表示
} else {
count = count + 1
cat(c(str1, get_value(str1, str2), get_value(str1, str3), "\n"),
file = con, append = TRUE)
}
}
}
}
close(con)
cat(count, "件") # 1896 件
センター試験の男女の都道府県別平均寿命の散布図は，小数点以下 2 桁までのデータが記載されている
https://www.mhlw.go.jp/toukei/saikin/hw/life/tdfk15/dl/tdfk15-09.xls
によるものと判明したので，やり直してみる。
表5-1 が平均寿命の推移（男），表5-2 が平均寿命の推移（女）で，それぞれが 7 段組みになっているので，最終的に以下のようなデータフレームにまとめ上げる。
pref.name male1965 female1965 male1975 female1975 male1985 female1985 male1995 … female2015
1 北海道 67.46 72.82 71.46 76.74 74.50 80.42 76.56 …
2 青森 65.32 71.77 69.69 76.50 73.05 79.90 74.71 …
3 岩手 65.87 71.58 70.27 76.20 74.27 80.69 76.35 …
4 宮城 67.29 73.19 71.50 77.00 75.11 80.69 77.00 …
5 秋田 65.39 71.24 70.17 75.86 74.12 80.29 75.92 …
6 山形 66.49 71.94 70.96 76.35 74.99 80.86 76.99 …
：
46 鹿児島 67.36 72.71 70.54 76.53 74.09 80.34 76.13 …
47 沖縄 NA NA 72.15 78.96 76.34 83.70 77.22 …
# https://www.mhlw.go.jp/toukei/saikin/hw/life/tdfk15/dl/tdfk15-09.xls
# 図表データのダウンロード
# 表5-1 が平均寿命の推移（男），表5-2 が平均寿命の推移（女）
# install.packages("readxl")
library(readxl)
# データが 7 段組み（S40, S50, S60, H7, H17, H22, H27）になっているので分解と結合
pref.name = c("北海道", "青森", "岩手", "宮城", "秋田", "山形", "福島",
"茨城", "栃木", "群馬", "埼玉", "千葉", "東京", "神奈川", "新潟",
"富山", "石川", "福井", "山梨", "長野", "岐阜", "静岡", "愛知",
"三重", "滋賀", "京都", "大阪", "兵庫", "奈良", "和歌山", "鳥取",
"島根", "岡山", "広島", "山口", "徳島", "香川", "愛媛", "高知",
"福岡", "佐賀", "長崎", "熊本", "大分", "宮崎", "鹿児島", "沖縄")
year = c(1965+0:4*10, 2010, 2015)
get.data = function(name, value) {
name = gsub(" ", "", name)
return(sapply(pref.name, function(pn) value[which(pn == name)]))
}
data3 = matrix(0, 47, 14)
for (page in 1:2) {
data = as.data.frame(read_excel("tdfk15-09.xls", sheet = page+5,
col_types = rep(c("skip", "text", "numeric"), 7)))
data2 = data[c(6:53),]
data2 = data2[-ifelse(page == 1, 52, 53),]
data2[47, 1] = "沖縄" # 1965 年には沖縄は含まれない
for (i in 1:7) {
data3[, i*2-(2-page)] = get.data(data2[, i*2-1], data2[, i*2])
}
}
data4 = data.frame(data3)
colnames(data4) = c(outer(c('male', 'female'), year, paste0))
data5 = cbind(data.frame(pref.name), data4)
all2.data = data5 # 最終的なデータフレーム
write.csv(all2.data, "all2.csv", row.names=FALSE)
all2.csv を用いて，目的の散布図を描く。
par(las=1, bty="l", mar=c(4, 4, 2, 0.5), mgp=c(2, 0.4, 0), tck=-0.005)
df2 = read.csv("all2.csv")
plot(female2015 ~ male2015, data=df2, asp = TRUE, xlim=c(80, 80.5), ylim=c(85.5, 88.0),
xlab="", ylab="", main = "平成27年度の都道府県別平均寿命")
for (intercept in 1:5) {
abline(intercept*0.5+5, 1, col="gray")
}
abline(v=seq(78.0, 82, by=0.5), h=seq(85.5, 88.0, by=0.5), lty=3, col="gray")
mtext("男", side=1, line=1.5)
mtext("女", side=2, line=2.3, las=1)
mtext("https://www.mhlw.go.jp/toukei/saikin/hw/life/tdfk15/dl/tdfk15-09.xls に基づき作図", side=1, line=2.5, adj=1, cex=0.8)
正しく描けたようだ。
次は，男女差のヒストグラム。
選択肢の中にあるヒストグラムが描けた。
par(las=1, bty="l", mar=c(4, 4, 2, 0.5), mgp=c(1.5, 0.4, 0), tck=-0.01)
hist(df2$female2015 - df2$male2015, breaks=seq(5.5, 7.5, by=0.25), right=FALSE,
col="aliceblue", xlab="歳", ylab="度数",
main="平成27年の男女の都道府県別平均寿命の差")
さて，次は，都道府県別の男女の平均寿命の散布図が描かれているのだ。
これも，前に示したデータ all.csv を使えばすぐに描ける。重なるデータ点は黒丸で示し，何個重なっているかを数字で表す。
と思いきや。ん？？様子が変だ。センター試験の問題用紙にある図と違いがある。
しばらく悩んだが，何の説明もしていないが，どうやらセンター試験の散布図はデータに jitter() を掛けているいるのだ。と思ったが，小数点以下 2 桁まで表示されている別の Excel ファイルにもとづいたもののようだ。
なので，以下は削除します。
言いがかりを付けてしまい，ご迷惑をおかけした関係各位にお詫び申し上げます。申し訳ありませんでした。
で，元のプログラムに付け加えた。JITTER = TRUE にしてから以下のプログラムを実行するとそれらしきものになるが，jitter は乱数で決めているので，毎回結果は少しずつ変わるし，jitter の factor の大きさ指定によっても違いの大きさも変わる。以下の図ではデータの重なりは全くないようになっている。
「厚生労働省の Web ページにより作成」と書いてあるのに，しかも入学試験で何の断りもなくjitter を加えた図を描くのはいかがなものか？
中澤さんが言及している図であるが，
この図から「フィンランドの出生率の2015年以降の急低下に実質的な意味があるのか，今後持続するのかどうかはわからない。」というのを導くって，図の読み取りが難しい？
横軸と年代を入れ替えて以下のような図にすると，わかりやすいような気がするのだけど。
20-24 はずっと前から下がって，1/4 にもなっている。25−29 は 1970 年頃が谷でその後持ち直したが 1990-1994 をピークに下がっている。30−34 は 1970 年頃が谷でその後持ち直していたが 2005-2009 から下がっている。35-39 は2010-2014 から下がっている。
素人だからよく分からないのではあるが。もとの図よりはわかりやすいような気がする。
if (!require(wpp2019)) { install.packages("wpp2019", dep=TRUE); library(wpp2019) }
data(tfr)
data(percentASFR)
SY <- 1950+0:13*5
EY <- 1955+0:13*5
TFRFinland <- subset(tfr, name=="Finland")[, sprintf("%4d-%4d", SY, EY)]
pASFRFinland <- subset(percentASFR, name=="Finland")[, c("age", sprintf("%4d-%4d", SY, EY))]
ASFRFinland <- pASFRFinland
for (i in 0:13) {
ASFRFinland[, i+2] <- pASFRFinland[, i+2]/100*TFRFinland[1, i+1]
}
color = c("black", "red", "blue", "brown", "aquamarine4 ", "magenta", "purple")
old = par(las=1, tck=-0.01, mar=c(3, 3, 1.5, 2), mgp=c(1.2, 0.3, 0), bty="l",
cex.axis=0.6, cex.lab=0.8)
data = ASFRFinland[, 2:15]
colnames(data) = sprintf("%4d-%4d", SY, EY-1)
matplot(t(data), type="l", lwd=2,
ylim=c(0, 1), xaxt="n",
col= color, lty=1,
xlab="year", ylab="ASFR")
title(main=list("ASFR by age in Finland from 1950 to 2020", cex=0.9))
legend = sprintf("%d-%d", 3:9*5, 3:9*5+4)
legend[4] = paste0("age\n", legend[4])
text(14, data[, 14], legend, cex=0.6, col=color, xpd=TRUE, pos=4)
axis(1, at=1:14, label=colnames(data))
par(old)
であるが，ggplot2 というのが気に入らない。
まず，データの取得であるが，Excel で 12 のシートがあるからと怖じ気づいたようだが，以下のようにすればよい。簡単とはいえないかもしれないが，手作業するより気は楽だ。なお，tidyverse なども使わない。
2020/01/25:追記 誤って，一箇所だけバグのあるソースコード（下記に赤字で示した箇所）を掲載してしまいました。
より簡単な方法（プログラム）を投稿しました。
「素の R で，データスクレイピングとプロット（その4）」
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
# https://www.mhlw.go.jp/toukei/saikin/hw/life/ckts15/dl/ckts15-06.xls
# 図表データのダウンロード
# 統計表1-1〜1-12 が「市区町村別平均寿命」
# install.packages("readxl")
library(readxl)
# データが 3 段組み（最終ページは 2 段組み）になっているので分解と結合
for (page in 1:12) {
data = as.data.frame(read_excel("ckts15-06.xls", sheet = page+4))
# stopifnot(data[3,1] == "市区町村") # データ開始行の確認
data = data[-(1:3),]
colnames(data) = rep(letters[1:5], 3)[1:ncol(data)]
# stopifnot((page != 12 && ncol(data) == 14) || # 段組み数の確認
# (page == 12 && ncol(data) == 9))
if (page == 1) {
data2 = rbind(data[1:4], data[6:9], data[11:14])
} else if (page != 12) {
data2 = rbind(data2, data[1:4], data[6:9], data[11:14])
} else {
data2 = rbind(data2, data[1:4], data[6:9])
}
}
# 列名は，都道府県・市郡 name1, 区町村 name2，男 male，女 female
colnames(data2) = c("name1", "name2", "male", "female")
# 都道府県名が の行は除いて良いか確認
# data2[is.na(data2$name1),]
data2 = data2[!is.na(data2$name1),]
data2 = data2[!is.na(data2[,3]),] # 米原市の後の不正な1行を除く
n = nrow(data2) # 1965 行
pref.code = integer(n) # 都道府県コード
pref.count = 0
for (i in 2:n) {
cap =substr(data2$name1[i], 1, 1)
if (cap != " ") { # 先頭1文字が全角空白でなければ都道府県名
pref.count = pref.count+1
# cat(pref.count, data2$name1[i], "\n") # 47都道府県コードが表示されることを確認
} else {
pref.code[i] = pref.count # 該当行に都道府県コードを付与
}
}
data3 = data2
# 各行に都道府県コード pref.code を追加
data3$pref.code = pref.code
# 文字データを数値に変換。但し，欠損値 "…" があることに注意。
data3$male = as.numeric(gsub("…", NA, data3$male))
data3$female = as.numeric(gsub("…", NA, data3$female))
# 都道府県・市郡 name1, 区町村 name2 の正規化
data3$name1 = gsub(" ", "", data3$name1)
data3$name2 = gsub("0", "", data3$name2)
data3$name2 = gsub(" ", "", data3$name2)
data3$name2[is.na(data3$name2)] = "" # name1 == 宮城県 で，name2 = を空に
data4 = data3
# 特別区を持つ市をまとめたデータ行の区分として，pref.code = -1 とする
rownames(data4) = seq_len(nrow(data4))
special = NULL # 特別区を持つ市の名前を記録
for (i in seq_len(nrow(data4))) {
if (data4$name1[i] == "東京都区部") { # これは特別中の特別
data4$pref.code[i] = -1
# print(data4[i, ])
} else if (grepl("市", data4$name1[i]) && grepl("区", data4$name2[i])) {
special = c(special, data4$name1[i])
# print(data4[i, ])
}
}
special = unique(special) # ユニークな特別区のリスト
for (i in seq_len(nrow(data4))) {
# name1 が特別区のある市の名前で，name2 が空の場合が特別区をまとめたデータ行
if (data4$name1[i] %in% special && data4$name2[i] == "") {
data4$pref.code[i] = -1
# print(data4[i,])
}
}
all.data = data4 # 最終的なデータフレーム
write.csv(all.data, "all.csv", row.names=FALSE)
必要のないデータも含んでいるので，以下のように必要なデータ（男の市区町村別データ）だけを取り出して boxplot で描画する。
par(las=1, bty="l", mar=c(4, 4, 2, 0.5), mgp=c(2, 0.4, 0), tck=-0.005)
df = read.csv("all.csv")
df2 = df[df$pref.code == 0,]
pref.name = unique(factor(df2$name1))[-1]
df2 = df[df$pref.code > 0,]
means = by(df2$male, df2$pref.code, mean, na.rm=TRUE)
boxplot(df2$male ~ df2$pref.code, horizontal=TRUE,
xlab= "", ylab="",
main = "平成27年度の都道府県別（市区町村別）平均寿命",
at=rank(-means), names=pref.name,
cex.axis=0.6, col="gray", range=0, whisklty=1, outline=TRUE)
points(means, rank(-means), pch = 19, cex=0.3, col="red")
mtext("平均寿命（歳）", side=1, line=1.5)
mtext("都道府県", side=2, line=2.3, las=0)
mtext("https://www.mhlw.go.jp/toukei/saikin/hw/life/ckts15/dl/ckts15-06.xls に基づき作図", side=1, line=2.5, adj=1, cex=0.8)
これで綺麗な図が描ける。図中の赤丸はおまけで示した平均値。
今すぐにでも1千万円を儲ける方法は，宝くじを 2千141万円分ほど買うことです。
還元率が46.7パーセントなので，欲しい金額のほぼ倍額の宝くじを買えば良いだけです。
ただし，かならずしも投資額の 46.7% が得られるのではありません。
それ以上の金額を回収できることもあるでしょうし，回収額がほとんど0ということもあるというだけの話です。
要するに，ギャンブルですね。
「Web ページからデータスクレイピング」の続き。
Web ページに，作表されたデータが表示されているとき，これを簡単に取り出すことができる。
1. R, Excel を立ち上げておく。
2. 目的のページを表示する
たとえば，
成人喫煙率（JT全国喫煙者率調査）
http://www.health-net.or.jp/tobacco/product/pd090000.html
を表示し，表の左上「年度」から表の右下「8.7」までのセルを，マウスドラッグで選択し，バッファにコピーする（command + c）。
3. Excel へペーストする。
4. 表の本体，左上の 「20歳代」 から 右下の 「8.7」 までのセルを，マウスドラッグで選択し，バッファにコピーする（command + c）。
5. R のコンソールで，
data = read.table(pipe("pbpaste", "r"))
と入力する。
> head(data)
V1 V2 V3 V4 V5 V6
1 80.5 84.7 86.7 81.4 74.6 82.3
2 6.6 13.5 19.0 23.0 23.0 15.7
3 83.5 84.8 87.3 83.4 78.0 83.7
4 10.6 14.3 22.0 24.1 24.1 18.0
5 83.2 84.1 85.8 82.3 73.3 82.3
6 11.0 16.4 20.9 23.1 20.3 17.7
:
のようになる。
6. データは男女が交互になっているので，まず奇数行 data[1:54*2-1, ]，次に偶数行 data[1:54*2, ] を取りだしたものを cbind() する。
data = cbind(data[1:54*2-1, ], data[1:54*2, ])
> head(data)
V1 V2 V3 V4 V5 V6 V1 V2 V3 V4 V5 V6
1 80.5 84.7 86.7 81.4 74.6 82.3 6.6 13.5 19.0 23.0 23.0 15.7
3 83.5 84.8 87.3 83.4 78.0 83.7 10.6 14.3 22.0 24.1 24.1 18.0
5 83.2 84.1 85.8 82.3 73.3 82.3 11.0 16.4 20.9 23.1 20.3 17.7
7 78.0 79.3 82.5 81.3 70.8 78.5 8.1 13.6 17.8 21.1 20.4 15.4
9 78.5 80.6 83.7 80.3 71.1 79.1 9.9 13.1 16.8 20.7 19.8 15.4
11 79.9 78.4 81.0 78.3 67.8 77.5 9.8 13.0 16.1 23.3 20.0 15.6
:
のようになる。
7. あとは，西暦年を列に付加するとか，変数名を付けるとか，どうしてもやりたければ整然データにするとか，どうにでも。
もう一つ例を挙げておく。
Wikipedia にある，「日本の自殺」の中の「自殺者数および人口10万人中の自殺率の推移」
これを上に書いたように範囲を指定してコピーし Excel などにペーストすると，数値のはいっているセルがだんだん右の方にズレている。
Excel ではなく，普通のエディタにペーストすればよい。数値がカンマで 3 桁区切りになっているので，全文置換でカンマを削除する。また，戦中のデータがないときの "" を NA に変換し，ヘッダー行を整えた後に，データ部分をバッファにコピーし，R で data = read.table(pipe("pbpaste", "r")) する。
ここでの教訓は，なんでもかんでもプログラム（特に tidyverse ?? なんか）に頼らないこと。得手不得手があるので，単純なエディタで十分な作業もあるんだということ。
あるところで，代金 200 円也でデータスクレイピングしたデータファイルが売られている。
成人喫煙率（JT全国喫煙者率調査）
http://www.health-net.or.jp/tobacco/product/pd090000.html
つまり，以下は少なくとも 200 円の価値のある作業だということか（笑）。
上に示した URL を開き，ページを「別名で保存...」する。
次に以下の R プログラムで処理する。
# 対象ファイル
file.name = "最新たばこ情報|統計情報|成人喫煙率（JT全国喫煙者率調査）.html"
# 保存された html ファイルを読み込む（エンコーディングに注意）
con = file(file.name, encoding = "cp932")
text = readLines(con) # 入力
data = character(858 - 47) # データ保管用のベクトル
count = 0 # 数値データの存在行のカウント
for (txt in text) { # 各行について
txt = gsub("<.*?>", "", txt) # html タグを除去する
if (nchar(txt) > 0) { # 本文もしくは数値（文字）など
count = count + 1
if (count >= 47 && count <= 857) { # この範囲なら目的のデータ
data[count - 46] = txt # 保存
}
}
}
data = data[data != "平成元年"] # 不要な要素を除去する
data = as.data.frame(matrix(data,
ncol = 15, byrow = TRUE), stringsAsFactors = FALSE) # データフレームに変換する
data = data[-c(2, 9)] # 2 列目と 9 列目は性別データなので消去
colnames(data) = c("year", # 列名（変数名）を付ける
"male20-29", "male30-39", "male40-49", "male50-59", "male60-69", "male all",
"female20-29", "female30-39", "female40-49", "female50-59", "female60-69", "female all")
data$year = 1965:2018 # 1 列目を西暦年にする
data = sapply(data, as.numeric) # 文字列を数値に変換する
これで，以下のようなデータフレームができあがる。
> head(data)
year male20-29 male30-39 male40-49 male50-59 male60-69 male all
[1,] 1965 80.5 84.7 86.7 81.4 74.6 82.3
[2,] 1966 83.5 84.8 87.3 83.4 78.0 83.7
[3,] 1967 83.2 84.1 85.8 82.3 73.3 82.3
[4,] 1968 78.0 79.3 82.5 81.3 70.8 78.5
[5,] 1969 78.5 80.6 83.7 80.3 71.1 79.1
[6,] 1970 79.9 78.4 81.0 78.3 67.8 77.5
female20-29 female30-39 female40-49 female50-59 female60-69 female all
[1,] 6.6 13.5 19.0 23.0 23.0 15.7
[2,] 10.6 14.3 22.0 24.1 24.1 18.0
[3,] 11.0 16.4 20.9 23.1 20.3 17.7
[4,] 8.1 13.6 17.8 21.1 20.4 15.4
[5,] 9.9 13.1 16.8 20.7 19.8 15.4
[6,] 9.8 13.0 16.1 23.3 20.0 15.6
長いなあ。面倒だなぁ。とお思いの方，次の記事で超簡単なデータスクレイピングをお見せします。
それはさておき，早速 graphics::matplot() で図を描いてみる。
old = par(mar = c(3.5, 3, 1.5, 5), mgp = c(1.5, 0.4, 0), tck = -0.01,
bty = "n", las = 1)
colors = c("black", "red", "blue", "brown", "mediumseagreen", "magenta")
matplot(data[, -1], type = "l", lwd = 2, xaxt = "n", ylim = c(0, 100),
lty = 1, col = colors,
xlab = "年", ylab = "喫煙率", main = "性別・年代別喫煙率の推移")
mtext("JT全国喫煙者率調査", xpd=TRUE, side = 3, line = -0.8, cex=0.8)
axis(1, at = 1:nrow(data), label = data[, 1], pos = 0)
delta = c(0, -1, 0, 0, 0, 0, -0.1, -0.1, 0.8, 0.5, -0.7, 0)*strheight("H")
text(55, data[54, 2:13]+delta, colnames(data)[2:13], col = colors, pos = 4,
cex = 0.7, xpd = TRUE)
par(old)
西原氏の作業手順は以下の通り。
1. 整然データを作るという前提
2. すべての年のファイルに以下の処理を行う
注釈的要素の除去
空行や空列の除去
整然データを作るために転置する
転置すると行列になり扱いにくいのでデータフレームに戻す
変数の名前を付ける
「公 立」のような余計な空白の削除
結合が解除されてできる NA を、そのセルの上にある内容で穴埋めする
合計を示すデータの削除
データは文字列として扱われているので整数に変換
元のデータに含まれない調査年を示す列の不可
3. 結果を一つのデータフレームにする
ずいぶんと処理内容が多い。西原氏自身が「この記事は非常に長いものになっている。この長さは、データクリーニングの繁雑さに比例したものである。つまり、データクリーニングが容易ではなく、うんざりするほどのものであることを反映している。自らの手でデータを扱わない人は、この分量を見てデータクリーニングの大変さを感じていただければと思う。」というほどだ。
しかし，手順を逆にすれば，話は比較的簡単になる。
西原氏は，余計なものを除くこと，余計なものでなくすることを目標にしている。しかし，官製 Excel データファイルは，余計なものばかりなのだ。
これに対して，必要なもの（データ）は全体から見ればほんのわずかだ。
やるべき事は各ファイルを読み込み，必要なデータだけを取りだし，一つのデータフレームを作る。
これだけで，データ分析（例えば図を描く）ができるようになる。
整然データにしたければ，そのデータフレームを変換すればよい（多くの場合はそれさえも不要）。
まあ，西原氏の書いたとおり，一定の書式で作られたデータではないが，それでもいくつかの書式に準拠しているのでその規則を利用して，必要なデータを取り出そう。
1986 年のデータ（csv ファイル）は以下のようになっている。
"60 男女別学校数","...2","...3","...4","...5","...6","...7","...8","...9","...10","...11"
NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA
NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA
"区分","計",NA,NA,"国立","公立",NA,NA,"私立",NA,NA
NA,"計","本校","分校","本校","計","本校","分校","計","本校","分校"
"計","5491","5295","196","17","4178","3990","188","1296","1288","8"
"男女ともにいる学校","4362","4188","174","15","3865","3695","170","482","478","4"
"男のみの学校","400","393","7","1","124","120","4","275","272","3"
"女のみの学校","708","693","15","1","189","175","14","518","517","1"
"生徒のいない学校","21","21","0","0","0","0","0","21","21","0"
"この表は，男子校あるいは女子校という分類ではなく，現実に在学している生徒の状況により分類して集計した。",NA,NA,NA,NA,NA,NA,NA,NA,NA,NA
西原氏がデータ分析の一例として挙げた図を描くには，上の赤字で示した2つの数値（400, 708）があれば十分だ（それ以外の場所の数値も同時に取り出すことは容易だ）。
ファイルを読み取り，"男のみの学校"で始まる行とその次の行の，カンマで区切られた2番目の文字列を整数に変換する。これだけでできるかどうかやってみる。
get.value = function(txt) { # 必要なデータを取り出す関数
a = unlist(strsplit(txt, ",")) # カンマで区切られたデータをバラす
return(as.integer(gsub('"', '', a[2]))) # 2番目のデータを整数にして返す
}
dat = matrix(0, 31, 2)
for (year in 1986:2016) {
file.name = paste0(year, ".csv")
text = readLines(file.name)
for (i in 1:length(text)) {
txt = text[i]
if (grepl("男のみの学校", txt)) { # "男のみの学校"を含む行なら
a = get.value(text[i]) # その行からデータを取り出す
b = get.value(text[i+1]) # 次の行からデータを取り出す
cat(file.name, a, b, "\n")
dat[year-1985, ] = c(a, b)
break # 次のファイルの処理へ
}
}
}
2007 年までのデータはこれで取り出せた。2008 年からは必要なデータの位置が少し違うので，
データを取り出す関数を，以下のようにする。
get.value = function(txt) { # 必要なデータを取り出す関数
a = unlist(strsplit(txt, ",")) # カンマで区切られたデータをバラす
if (a[1] != "NA") { # 最初のデータが NA でなければ，2番目のデータを整数にする
return(as.integer(gsub('"', '', a[2])))
} else if (a[1] == "NA") { # 2008 年以降の場合は最初のデータが NA で，必要なデータは4番目にある
return(as.integer(gsub('"', '', a[4])))
}
}
私もビックリするくらい，実に簡単な規則にそってデータを取り出すことができた。
後はデータフレームにして，列に名前を付ける。必要ならこれを整然データにすれば良いが，そんなことはどうでもよい。
df = data.frame(dat)
colnames(df) = c("male", "female")
rownames(df) = 1986:2016
print(df)
male female
1986 400 708
1987 384 696
1988 366 690
1989 352 689
：
2014 125 320
2015 117 314
2016 112 311
西原氏と同じ図を描いてみよう。ただし，ggplot ではなく，graphics::matplot() で。
old = par(mar=c(4, 4, 1.5, 1), mgp=c(2.5, 0.8, 0), las=1, tck=-0.01)
matplot(df, type="l", ylim=c(0, 700),
xaxt = "n",
xlab = "年", ylab = "学校数",
main="日本における男子のみ・女子のみの高校（通信制除く）の数")
text(25, c(450, 200), c("男子のみ", "女子のみ"))
axis(1, at=1:31, labels=1986:2016)
par(old)
で，「まずは、e-Stat から男女別学校数をまとめた統計表の Excel ファイルをダウンロードしよう。先に述べたように、1年度につき1ファイルなので、必要となる年度の Excel ファイルを1つずつ手動でダウンロードする。この作業もできれば自動化したいところであったが、自動化のスクリプトを書く方がよほど時間がかかりそうだったので、手動で処理することにした。」
と，あった。ダウンロードすべきファイルは実に 1986 年から 2016 年までの 31 個である。私は，手動ではやりたくない。結果として，プログラムもそんなに複雑ではなかった。
1. 「学校基本調査，男女別学校数」のページを開く。以下の 2 ページ。
https://www.e-stat.go.jp/stat-search/files?page=1&query=%E7%94%B7%E5%A5%B3%E5%88%A5%E5%AD%A6%E6%A0%A1%E6%95%B0&layout=dataset&toukei=00400001&tstat=000001011528&metadata=1&data=1
https://www.e-stat.go.jp/stat-search/files?page=1&query=%E7%94%B7%E5%A5%B3%E5%88%A5%E5%AD%A6%E6%A0%A1%E6%95%B0&layout=dataset&toukei=00400001&tstat=000001011528&metadata=1&data=2
このページを表示させて，1 ファイルずつ指定して，保存したのだろうか。
しかも，「ファイル名を西暦年としたことには理由がある。おそろしいことに、e-Stat から手に入る Excel ファイルには、それが何年の調査結果なのかを示す情報がまったく入っていない。このため、ファイル名に西暦年の情報を入れておくことで、あとでデータクリーニングをするときにファイル名から西暦年の情報が得られるようにしておくのだ。」というように，名前を変えながら（間違えないように）保存しないといけなかったのか。大変だ。
2. このページの html ファイルには，ダウンロードすべきファイルの情報が書かれている。
ページを別ファイルで保存する。2 ページあるので，page1.html, page2.htmlと して保存する。
3. 年度とそ年度の結果の url を取り出す。
調査年度は '<span class="stat-sp">調査年月 </span>2019年</li>' のような行にある。
ここから 2019 を取り出すのは簡単だ。
ダウンロードすべきファイル情報は，キーワード"EXCEL" を含む行にある。
'<a href="https://www.e-stat.go.jp/stat-search/file-download?statInfId=000031893802&fileKind=0" class="stat-dl_icon stat-icon_0 stat-icon_format js-dl stat-download_icon_left" data-file_id="000008560862" data-release_count="1" data-file_type="EXCEL" tabindex="22">'
ここから 'https://www.e-stat.go.jp/stat-search/file-download?statInfId=000031893802&fileKind=0' のような部分を取り出す。ただし，実際にアクセするときには，& は 単に一文字の & にしないといけないので注意。
4. utils::download.file() でファイル(*.xls, *.xlsx) をダウンロードする。
実は，2015 年以降は *.xlsx なので openxlsx::read.xlsx で読めるのだが，2014 年以前は *.xls なので読めない。
xls, xlsx の両方が読めるのは eadxl::read_excel() だが，これは URL を指定してネットから読むことができない。
そこでまず，utils::download.file() を使って，上記 3. で取り出した url の内容を保存する。保存時に年度により適切な拡張子を付けてやる。
5. *.xls, *.xlsx を readxl::read_excel() で読み込み，write.csv() で csv ファイルに書き出す。
プログラムは以下のようになる。使用するしないにかかわらず，1963 年〜2019 年までのデータをダウンロードした。
library(readxl) # read_excel() を使うため
url.info.pages = c("page1.html", "page2.html") # ダウンロードするファイル情報を含むページの html ファイル
last.file = FALSE # 年次が同じ2種類のファイルが掲載されているので，必要な最後のファイルが来たらおしまいにする
for (page in url.info.pages) { # 2 ページについて
text = readLines(page) # テキストとして読み込む
for (i in 1:length(text)) { # 各行を見ていく
txt = text[i] # 当該行
if (grepl("調査年月", txt)) { # 「調査年月」を含む行から 4 桁の西暦年 year を残す
year = sub("^.+</span>", "", txt)
year = sub("年.+", "", year)
}
if (grepl("data-file_type=\"EXCEL\"", txt)) { # キーワード EXCEL を含む行から，URL を残す
url = sub("<a href=\"", "", txt)
url = sub("&fileKind=0.*", "&fileKind=0", url)
if (year < 2015) { # 2015 年以前は *.xls
file.name = paste0(year, ".xls")
} else { # 2015 年以降は *.xlsx
file.name = paste0(year, ".xlsx")
}
download.file(url, file.name) # url からファイルを（そのまま）ダウンロードし，保存
df = read_excel(file.name) # Excel ファイルを読み込む
write.csv(df, file= paste0(year, ".csv")) # csv ファイルとして書き出す。ファイル名は 2019.csv のように
cat('output to ', paste0(year, ".csv"), "\n\n")
if (year == 1970) { # 処理したフィあるが 1970 年のものならこれでおしまい（調査年順に掲載されていない）
last.file = TRUE # 最後のファイルだとして内側の for ループを抜ける
break
}
}
if (last.file) { # 最後のファイルなら外側の for ループも抜ける（終了）
break
}
}
}
「e-Stat から入手したデータの利用にあたり，やるべき事は単純だ」 へと続く
