お断り

本内容は、2024年6月9日にTokyoRのLTで発表した内容の概要となります。内容について発表の際よりもブラッシュアップしています。特にメモリ割り当て量については、本稿ではmutate(cumsum())処理のメモリ割り当て量も追加しているため、発表の時よりインパクトの弱い（差のない）内容となっておりますが、こちらの方が正しい値ということでご理解いただければ幸いです。

Rで日本の年齢中央値を算出する

令和2年（2020年）の国勢調査のデータをもとに市町村別の年齢中央値を算出し、地図にプロットしたいと考えました。その際に、中央値算出のプログラムでちょっとした工夫をしたのでご紹介します。R以外のプログラミング言語でも参考になるようであれば幸いです。

データセット

まずデータセットですが、政府統計の総合窓口E-STATからダウンロードしました。
https://www.e-stat.go.jp/dbview?sid=0003445139
このデータを加工して、全国の年齢別人口を出したのが下記コードです。どこかからダウンロードするというのも面倒、と思われる方もいらっしゃると思いますので、実際の値を直接記述してあります。なお、国勢調査では100歳以上の人口はひとまとまりになっていますので、正確な平均値は算出できません。中央値が100歳以上にならなければ、中央値の方は正しく算出できます。

median.r

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()関数で実現できます。

median.r

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が奇数か偶数かで場合分けしようか？など、うだうだ考えるうちに下記の数行をコーディングするのに半日もかかってしまいました。しかし、やっていることは単純なので、コード（コーディング者の能力の低さ）を味わっていただければ幸いです。

median.r

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に貢献した気になれる。自己満足に浸れる。