お断り
本内容は、2024年6月9日にTokyoRのLTで発表した内容の概要となります。内容について発表の際よりもブラッシュアップしています。特にメモリ割り当て量については、本稿ではmutate(cumsum())処理のメモリ割り当て量も追加しているため、発表の時よりインパクトの弱い（差のない）内容となっておりますが、こちらの方が正しい値ということでご理解いただければ幸いです。
Rで日本の年齢中央値を算出する
令和2年（2020年）の国勢調査のデータをもとに市町村別の年齢中央値を算出し、地図にプロットしたいと考えました。その際に、中央値算出のプログラムでちょっとした工夫をしたのでご紹介します。R以外のプログラミング言語でも参考になるようであれば幸いです。
データセット
まずデータセットですが、政府統計の総合窓口E-STATからダウンロードしました。
https://www.e-stat.go.jp/dbview?sid=0003445139
このデータを加工して、全国の年齢別人口を出したのが下記コードです。どこかからダウンロードするというのも面倒、と思われる方もいらっしゃると思いますので、実際の値を直接記述してあります。なお、国勢調査では100歳以上の人口はひとまとまりになっていますので、正確な平均値は算出できません。中央値が100歳以上にならなければ、中央値の方は正しく算出できます。
library(tidyverse)
library(conflicted)
#日本の年齢別人口(2020年国勢調査)。下記URLのデータから加工した結果。
#https://www.e-stat.go.jp/dbview?sid=0003445139
df_all_japan <- dplyr::tibble(age = seq(from = 0, to = 100),
value = c(831824, 866525, 910005, 934063, 973665,
998664, 996576, 1020657, 1024325,
1048871, 1057736, 1063185, 1083113,
1077379, 1069104, 1070370, 1113159,
1123237, 1151389, 1159285, 1177049,
1174456, 1193935, 1193983, 1191883,
1210380, 1212167, 1192761, 1209983,
1206673, 1234225, 1258775, 1299569,
1335672, 1356353, 1407719, 1456774,
1476522, 1478920, 1491632, 1558309,
1596338, 1656969, 1699648, 1779813,
1852394, 1956602, 1991137, 1954205,
1895955, 1837431, 1807894, 1764215,
1758955, 1371356, 1689290, 1582910,
1542017, 1491947, 1461318, 1469672,
1494495, 1450282, 1407740, 1475001,
1516441, 1512118, 1598839, 1679855,
1768015, 1886160, 2052521, 2014143,
1893386, 1165585, 1234832, 1488292,
1416571, 1430462, 1360771, 1206058,
1017635, 1048747, 1035057, 989231,
893754, 794162, 739149, 662252, 580506,
493824, 425042, 358386, 277613, 224151,
169988, 123057, 90075, 64897, 44707,
79523))
median()を使ったメモリ割当量は約2GBに
Rで中央値を算出するにはmedian()関数を使えば良いです。もちろん私も実際にやってみました。すると、Rstudioのメモリ消費量を表す値が甚だしく大きくなる（全市町村分map()で繰り返し処理を行うと10GBくらい）ことに気付きました。当然のことながらnが大きくなると（今回の場合は1億2千万件以上）メモリ割り当て量が膨大になります。それではどのくらい使用しているのかモニタリングしてみました。profmemというライブラリのprofmem()関数で実現できます。
library(profmem)
p <- profmem({
vec <- rep(df_all_japan$age, df_all_japan$value)
med <- median(vec)
})
アウトプットは以下のようになります。
> p
Rprofmem memory profiling of:
{
vec <- rep(df_all_japan$age, df_all_japan$value)
med <- median(vec)
}
Memory allocations:
what bytes
1 alloc 544
2 alloc 1760
3 alloc 1760
4 alloc 1072
5 alloc 184
6 alloc 312
7 alloc 184
8 alloc 492857096
9 alloc 456
10 alloc 1400
11 alloc 6488
12 alloc 6488
13 alloc 1072
14 alloc 352
15 alloc 872
16 alloc 1688
17 alloc 552
18 alloc 232
19 alloc 872
20 alloc 492857096
21 alloc 492857096
22 alloc 492857096
23 alloc 528
24 alloc 1648
25 alloc 1648
26 alloc 1072
27 alloc 256
28 alloc 456
29 alloc 216
30 alloc 256
total 1971460752
calls
1 $() -> get()
2 $() -> get()
3 $() -> get()
4 $() -> get()
5 $() -> get()
6 $() -> get()
7 $() -> get()
8 <internal>
9 <internal>
10 median()
11 median()
12 median()
13 median()
14 median()
15 median()
16 median()
17 median()
18 median()
19 median()
20 median() -> median.default()
21 median() -> median.default() -> sort() -> sort.default() -> sort.int()
22 median() -> median.default() -> sort() -> sort.default() -> sort.int()
23 <internal>
24 <internal>
25 <internal>
26 <internal>
27 <internal>
28 <internal>
29 <internal>
30 <internal>
total
行が分離していて見づらいため、head()で5行目まで絞ってアウトプットの見方を説明します。
> head(p)
Rprofmem memory profiling of:
{
vec <- rep(df_all_japan$age, df_all_japan$value)
med <- median(vec)
}
Memory allocations:
what bytes calls
1 alloc 544 $() -> get()
2 alloc 1760 $() -> get()
3 alloc 1760 $() -> get()
4 alloc 1072 $() -> get()
5 alloc 184 $() -> get()
6 alloc 312 $() -> get()
total 5632
アウトプットは3列になっており、それぞれ"what","bytes","calls"という列名になります。Rのヘルプからprofmem()の説明を見ると、
what: (character) type of memory event; either "alloc" or "new page"
bytes: (numeric) number of bytes allocated or NA_real_ (when what is "new page")
となっています。whatはメモリが割り当てられたか、bytesはメモリに割り当てられたバイト数と考えて間違いないでしょう。
median()が内部で行っている処理
となると、今回問題になるのは8行目、20行目、21行目、22行目の約500MB(492,857,096バイト)の部分です。まず8行目はとなっているのでよくわかりませんが、1億2千万件以上のベクトルを作っているので、そりゃそうだと思います。20行目～22行目はmedian()の処理に費やしています。中でも21,22行目は"sort"処理に費やしているようです。つまりmedian()関数にはsort()関数の処理を含んでいるということがわかります。結果的にmedian()処理の部分で約500MB*3 = 約1.5GB、全体では約2GBのメモリを費やしていることがわかります。
人間的に中央値を算出する方法をプログラミングしてみる
では、全件（今回の場合は1億2千万件以上）を読み込まずに中央値算出できないか？というより算出できるよね？というのが、今回取り組んだことです。人間はどうやって中央値を出すか考えたとき、昇順（降順）に並べ、nの半分の位置がどこにあるかを判断しています。このとき、累計(cumsum)列があると、半分の位置がわかりやすくて便利ですよね。
ただしnが偶数のとき（例えばn=6の時）は、半分の位置の前後（n=6だったら、3列目と4列目）の値を足して2で割ります。この辺、どこまで厳密にやろうか？nが奇数か偶数かで場合分けしようか？など、うだうだ考えるうちに下記の数行をコーディングするのに半日もかかってしまいました。しかし、やっていることは単純なので、コード（コーディング者の能力の低さ）を味わっていただければ幸いです。
p_2 <- profmem({
df_human_med <- mutate(.data = df_all_japan, cumsum = cumsum(value))
total <- tail(df_human_med$cumsum, 1)
half_n <- (total + 1) / 2
lower_n <- floor(half_n)
upper_n <- ceiling(half_n)
lower_row <- min(which(lower_n <= df_human_med$cumsum))
upper_row <- min(which(upper_n <= df_human_med$cumsum))
med <- (df_human_med$age[lower_row] + df_human_med$age[upper_row]) / 2
})
このp_2のアウトプットがなんと676行になってしまいました！しかし、そのうち662行はmutate()処理、つまり累計の計算に費やされていました。そして最後の14行で、total～medまでの7行のコードの計算に費やしています。そんなわけで、676行のbytesの値のサマリーを以下に示します。
メモリ割り当て量が1000分の1以下に！
> library(summarytools)
> dfSummary(p_2$bytes)
p_2$bytes was converted to a data frame
Data Frame Summary
p_2
Dimensions: 676 x 1
Duplicates: 490
--------------------------------------------------------------------------------------------------
No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing
---- ----------- ------------------------------ --------------------- ------- ---------- ---------
1 bytes Mean (sd) : 2147.9 (13223.8) 186 distinct values : 676 0
[numeric] min < med < max: : (100.0%) (0.0%)
184 < 560 < 237216 :
IQR (CV) : 808 (6.2) :
:
--------------------------------------------------------------------------------------------------
> sum(p_2$bytes)
[1] 1451968
676行のbytesのうち最大値は約237KBで、平均2KBという結果でした。そして総バイト割り当ては1451986バイトですので、1.4～1.5MBのメモリ割り当てですみました。つまり、
> 1971460752/1451968
[1] 1357.785
メモリ割り当て量はmedian()関数を使ったときの1358分の1ですんだということになります。
結論という名の自己満足的まとめ
- 日本人の年齢の中央値のように、nが非常に大きいものの中央値を算出する場合は、(n+1)/2 と累計値とを比較して中央値を算出すると、PCへの負荷が大幅に低減する。
- Ecoである。あるいはEcoに貢献した気になれる。自己満足に浸れる。