人口モデルに基づくウクライナの人口予測

神戸大学大学院保健学研究科国際保健学領域（国際協力研究科兼任）
中澤 港<minato-nakazawa@umin.net>
30 Nov. 2013 @ 統計数理研究所＊本研究は「ポスト人口転換期におけるオプティマルな対処方策の研究」（研究代表者：東京大学・梅﨑昌裕准教授）の一部として実施した研究を発展させたものである
1

背景とアイディア■ 人口予測の方法

たいていの場合，コーホート要因法かそのバリエーション➔ 女性の年齢別出生率，年齢別死亡率（と年齢別移動率）を予測➔ 年齢別女性人口に年齢別出生率を掛け，（両性の場合は出生性比により男女別に）翌年のゼロ歳人口を計算➔ 男女各歳人口がその年の年齢別死亡率に従って死亡した時に生残する人数（＋純移動）を翌年の１歳上の人口とする➔繰り返し。が基本形

将来の年齢別死亡率や年齢別出生率については，固定するか線形補外，あるいはシナリオとして与えるのが普通
固定するとか線形補外には根拠が薄いが，合理的なトレンド予測は困難。とくに出生率には政策の影響や突発的な変動の影響も大きいので，予測は難しい■ 本研究での予測方法

基本形は移動を無視したコーホート要因法

将来の年齢別死亡率と年齢別出生率の予測に工夫➔
現在まで複数年の年齢別死亡率と年齢別出生率に，比較的単純なパラメタライズドモデルを当てはめ➔ 個々のパラメータを arfima モデルなどで時系列解析➔
パラメータの推移トレンドを予測➔予測されたパラメータをモデルに代入し，年齢別死亡率と年齢別2出生率の推移を予測基本情報■ 元にするデータ

人口構造データ（コーホート要因法は，ある時点の人口構造を元にして，出生，死亡，移動による変化を足していくので，元にする性・年齢別人口が必要）：ウクライナの国勢調査は 2001 年にされたものが公開済み
http://www.ukrcensus.gov.ua/eng/ 年齢別死亡率： Human Mortality Database(HMD) に含まれている， 1961 年から 2005 年までの男女別生命表から
 女性の年齢別出生率： Human Fertility Database(HFD) に含まれている 1959 年から 2006 年までの女性の年齢別出生率
 HMD/HFD に登録されている国のデータなら同じ方法で分析可■ ウクライナ人口の歴史的背景
 1986 チェルノブイリ原発事故
 1991 旧ソ連からの独立
 2004 オレンジ革命を通した欧米への接近 2010 親ロシア政権に復帰■ 既存の人口予測
 US Census Bureau が年齢５歳階級で実施，発表
3
Human Mortality Database■■■

年齢別死亡率は生命表という形で発表されている
UCB の Dr. Wilmoth と Max-Planck の研究者がメンテしている "Human MortalityDatabase" (HMD) は，世界 37 の国と地域について，長年にわたる生命表を提供
オーストラリア，オーストリア，ベラルーシ，ベルギー，ブルガリア，カナダ，チリ，チェコ，
デンマーク，エストニア，フィンランド，フランス，ドイツ，ハンガリー，アイスランド，アイル
ランド，イスラエル，イタリア，日本，ラトビア，リトアニア，ルクセンブルグ，オランダ，
ニュージーランド，ノルウェー，ポーランド，ポルトガル，ロシア，スロバキア，スロベニア，スペイン，スウェーデン，スイス，台湾，英国，米国，ウクライナ

4Human Fertility Database■ 女性の年齢別出生率は， Human Ferility Database (HFD)

から生データとして提供されている■ HFD は Max-Planck 研究所とウィーン人口研究所によって維持管理されている
http://www.humanfertility.org/cgi-bin/main.php

5HMD/HFD を使う方法■ New User としてメールアドレスとパスワードを登録■ データを得る

ログインして欲しいデータをダウンロードし，テキスト保存して Perl などで欲しい形式に加工するか， R で直接読む
 HMD については， Dr. Rob J. Hyndman がメンテナをしている demography パッケージからアクセス可能。
library(demography)➔ukraine <- hmd.mx("UKR", "username", "passwd")➔とすると， ukraine というオブジェクトにウクライナの生命表の mx が得られる。パスワードを平文でスクリプトに書いてしまうのが気になるが，便利ではある。 str(ukraine) でデータ構造をみると， year, age, pop, rate からなるリスト。 popと rate はそれぞれ female, male, total からなるリスト
 female, male, total は，それぞれ 111x51 の matrix になっている
 HFD はまだ対応したパッケージが無い
6

2001 年センサスによるウクライナの人口構造

x <- read.delim("https://sites.google.com/site/postdemotrans/progress/ukrainepop2001.txt)
library(pyramid) # CRAN からインストール可能
pyramid(x, Cstep=10, Laxis=0:5*100000, AxisFM="d")7

ウクライナ生命表への死亡モデルの fitting■ 生命表へのモデルの当てはめの検討課題

どの死亡モデルを使うか？
➔ Siler モデルがメカニズムを考えている中ではパラメータが 5 つと比較的少ない➔ Denny モデルが 3 パラメータで適合度がよい
 当てはめは dx, lx, qx のどれでやるか？

最小化対象とする関数は，年齢別死亡について，モデルとデータの差の二乗和をとるのでよいか？
 optim （）の method は何がいいか？
8

方法（年齢別死亡率データの準備）■ HMD に登録し，ウクライナの生命表データをダウンロード■ 元データを Perl のスクリプトにより変換する■ Denny のモデルを当てはめ■ （発表では Siler モデルを当てはめて，異常値が出た話を含めたが省略）

Perl の23 行
script9

モデルの fitting 方法■ Denny モデルそのものを関数 Denny() として定義。年齢別の lx を返す■ データ (rlx) とモデル (mlx) の年齢ごとの差の二乗和を返す
関数 fDenny() を定義■ optim() で fDenny() を最小化するパラメータ a,b,c を得る関数 fitDenny() を定義■ データをファイルから読み込み， 1 年分ずつ fitDenny() を実行してパラメータを行列に保存していく

10

Denny の年齢別死亡モデル■ fmsb パッケージに定義済み■ より単純な Denny のモデルを使ってみる> library(fmsb)
 > Denny function (a, b, c, t) {
lx <- ifelse(t < 105, 1/(1 + a * (t/(105 - t))^3 + b *sqrt(exp(t/(105 t)) - 1) + c * (1 - exp(-2 * t))), 0)
lx <- ifelse(lx > 1, 1, lx)
return(lx * 1e+05)
 } <environment: namespace:fmsb>
11

Denny モデルの当てはめlibrary(fmsb)■ ukm <- read.delim("./data/ukmdx.txt") # ウクライナ男性の dx の読込■ ryears <- 1961:2005■ fittedpar <- matrix(rep(0, 10*length(ryears)), length(ryears), 10)■ colnames(fittedpar) <- c("fa", "fb", "fc", "fRMSE", "fConv", "ma","mb", "mc", "mRMSE", "mConv")■ fittedlxs <- matrix(rep(0, 111*length(ryears)), 111, length(ryears))■ colnames(fittedlxs) <- paste("Y", ryears, sep="")■ # calc for females■ for (yy in ryears) {■ lx <- dxtolx(subset(ukf, Year==yy)$dx)■ params <- rep(0,5)■ params <- fitDenny(, lx, Method="BFGS", mode=2)■ fittedpar[yy-1960, 1:5] <- params■ fittedlx <- Denny(params[1], params[2], params[3], 0:(length(lx)-1))■ fittedlxs[, paste("Y", yy, sep="")] <- fittedlx■ }■ write.table(fittedlxs,"./ukflx-modeldenny.txt",quote=FALSE,sep="\t",row.names=FALSE) # ファイル保存■ # 女性も同様に推定■12

Denny のモデルを lx に当てはめる■ 1961-2005 の間，男女別々に→だいたい良い近似■ パラメータは forecast() で使うためにファイルに保存13時系列解析により係数を求め，信頼区間付きで予測
forecast パッケージを用い， Arfima モデルを当てはめ，将来のパラメータをテキスト保存■ どのモデルを使うか，当てはめ方はどうすべきかは，まだ検討が不十分dennyparams <- read.delim("./ukraine-params-denny.txt")
library(forecast)# females parameter a
ffa <- ts(dennyparams$fa, freq=1, start=1961)
fitffa <- arfima(ffa)print(projffa <- forecast(fitffa, level=c(0.95), h=50))
# females parameter b
ffb <- ts(dennyparams$fb, freq=1, start=1961)fitffb <- arfima(ffb)print(projffb <- forecast(fitffb, level=c(0.95), h=50))
# females parameter cffc <- ts(dennyparams$fc, freq=1, start=1961)fitffc <- arfima(ffc)print(projffc <- forecast(fitffc, level=c(0.95), h=50))
write.table(data.frame(YEAR=2006:2055,PFA=projffa$mean,
PFB=projffb$mean, PFC=projffc$mean),"projected-dennyfemales.txt",sep="\t",row.names=FALSE)
# 男性も同様に実行。
■14

時系列解析により係数を求め，信頼区間付きで予測■ Denny モデルをウクライナの死亡データに当てはめて得られた係数 a, b, c について男女別に Arfima モデルを当てはめ，予測した結果（左が男性，右が女性）

15

女性の年齢別出生率の単純なモデル化■ Hadwiger のモデル（下式）■ >library(fmsb)■ >Hadwiger■ function (a = 3.4, b = 2.5, c = 22.2)■ {ages <- 15:54■
rval <- a * b/c * (c/ages)^1.5 * exp(-b^2 * (c/ages +ages/c - 2))■
names(rval) <- ages■
return(rval)■ }■ <environment: namespace:fmsb>■

16

Hadwiger モデルの当てはめ
■
■
■
■
■
■
■
■
■
■
■
■
■
■
■
■
■
■
■

library(fmsb)
asfr <- read.delim("./ASFR.txt")
ryears <- 1961:2005
ages <- 15:54
fittedpar <- matrix(rep(0, 5*length(ryears)), length(ryears), 5)
colnames(fittedpar) <- c("a", "b", "c", "RMSE", "Conv")
fittedasfr <- matrix(rep(0, length(ages)*length(ryears)), length(ages),length(ryears))colnames(fittedasfr) <- paste("Y", ryears, sep="")# calc for femalesfor (yy in ryears) {pasfr <- subset(asfr, ((Year==yy)&(Age %in% ages)))$ASFRparams <- c(3.4, 2.5, 22.2, 0, 0)
params <- fitHad(params[1:3], pasfr, Method="Nelder-Mead")fittedpar[yy-1960, ] <- params
fittedpasfr <- Hadwiger(params[1], params[2], params[3])fittedasfr[, paste("Y", yy, sep="")] <- fittedpasfr}
write.table(fittedasfr,"./ukfittedasfrhadwiger.txt",quote=FALSE,sep="\t",row.names=FALSE)
write.table(fittedpar,"./ukasfr-paramshadwiger.txt",quote=FALSE,sep="\t",row.names=FALSE)

17女性の年齢別出生率の単純なモデル化■ 1961-2005 年の女性の年齢別出生率に当てはめた結果■ どの年もまあ当てはまる

18

Hadwiger モデルパラメータからの将来予測
■
■
■
■
■
■
■
■
■
■
■
■
■
■
■
■
■
■
■
■
■
■
■

params <- read.delim("./Ukraine-fertility/ukasfr-params-hadwiger.txt")
library(forecast) # Rob Hyndman が開発した予測用パッケージ
library(fmsb)# parameter a
ffa <- ts(params$a, freq=1, start=1961)
fitffa <- arfima(ffa)
print(projffa <- forecast(fitffa, level=c(0.95), h=50))# parameter b
ffb <- ts(params$b, freq=1, start=1961)
fitffb <- arfima(ffb)
print(projffb <- forecast(fitffb, level=c(0.95), h=50))# parameter c
ffc <- ts(params$c, freq=1, start=1961)
fitffc <- arfima(ffc)
print(projffc <- forecast(fitffc, level=c(0.95), h=50))
proja <- as.vector(projffa$mean)
projb <- as.vector(projffb$mean)
projc <- as.vector(projffc$mean)
TFRs <- rep(0, length(proja))for (i in 1:length(proja)) {TFRs[i] <- sum(Hadwiger(proja[i], projb[i], projc[i]))}write.table(data.frame(YEAR=2006:2055, TFRs, proja, projb, projc),"projectedhadwiger-asfrs.txt",sep="\t",row.names=FALSE)

19

Hadwiger モデルパラメータの時系列解析将来予測■ Arfima モデル（要検討）■ forecast パッ
ケージを使った当てはめ結果が右の図■ c がほぼ初婚または初産年齢に対応すると言われている

20保存しておいたパラメータをまとめて人口予測
■
■
■
■
■
■
■
■
■
■
■
■
■
■
■
■
■
■
■
■
■

# 予測された人口は人口ピラミッドとして animation パッケージを使って動画にして表示
pop <- read.delim("./ukrainepop2001.txt")
dfp <- read.delim("./projected-denny-females.txt")
dmp <- read.delim("./projected-denny-males.txt")
hp <- read.delim("./projected-hadwiger-asfrs.txt")
library(pyramid); library(fmsb);# Settingwx <- 101; mp <- fp <- matrix(rep(0, wx*50), wx, 50) # For 50 years
mp[,1] <- pop$Males; fp[,1] <- pop$Females# Start simulation calculationfor (j in 2:50) {asfr <- c(rep(0, 14), Hadwiger(hp$proja[j-1], hp$projb[j-1], hp$projc[j-1]), rep(0, wx54))
mmx <- qxtomx(lxtoqx(Denny(dmp$PMA[j-1], dmp$PMB[j-1], dmp$PMC[j-1], 1:(wx1))))
fmx <- qxtomx(lxtoqx(Denny(dmp$PMA[j-1], dmp$PMB[j-1], dmp$PMC[j-1], 1:(wx1))))baby <- sum(fp[, j-1]*asfr)
babym <- as.integer(baby*1.06/2.06+0.5); mp[1, j] <- babym
babyf <- as.integer(baby*1/2.06+0.5); fp[1, j] <- babyf
mp[2:wx, j] <- as.integer(mp[1:(wx-1), j-1]*(1-mmx))
fp[2:wx, j] <- as.integer(fp[1:(wx-1), j-1]*(1-fmx))}# これは期待値だけなので，今後， 95% 信頼区間の下限と上限の組合せを試す予定
21

保存しておいたパラメータをまとめて人口予測■ # Draw graphs■ library(animation) # animation パッケージを使う■ ani.options(interval = 0.5, outdir=getwd())■ saveVideo({ par(cex=0.8, las=1)■ pyramid(data.frame(M=pop$Males/1000,

F=pop$Females/1000, A=pop$Age), Laxis=0:5*100,
AxisFM="d", Cstep=10)■ for (i in c(2:50)) { pyramid(data.frame(M=mp[, i]/1000,F=fp[, i]/1000, A=pop$Age),■ Llab="Males\n(x1000)", Rlab="Females\n(x1000)",Clab="",■ Laxis=0:5*100, AxisFM="d", Cstep=10,main=sprintf("Year %d projected", i+2005))■ ani.pause()■ }}, video.name="changing-pyramid2.mp4", other.opts="-b300k")■ # mp4 で保存。 RStudio から実行すると，ファイル生成後，22自動的に再生される

シミュレーション結果（これはアニメ GIF ）23