統計数理（2004）
第 52 巻 第 1 号 117–134
c 2004 統計数理研究所


特集「極値理論」［原著論文］

年齢時代区分データによる
超高齢者寿命分布の推測
1

渋谷 政昭 ・華山 宣胤

2（受付 2003 年 8 月 4 日；改訂 2004 年 3 月 19 日）

要

旨

年齢時代区分別に観測された超高齢者の生存者統計および死亡者統計を極値理論を用いて解析し，寿命分布の限界を議論する．使用したデータは厚生労働省（厚生省）より発表されている
全国高齢者名簿および人口動態統計より得た．連続変数に対して議論している極値理論を，年
齢時代区分別に作表したデータへ適用することが必要である．そのため，一般パレート分布関数から計算した死亡確率を用いた多項分布モデルを当てはめる方法，乱数を加えて得られる連続量擬似確率標本を用いる方法を提案する．また，既存の手法としてレキシス図上の Piecewise
Constant Intensity Model も適用する．推定の結果，日本人男性寿命分布の有限な上限値が得られた．
キーワード：コホート分析，最尤法，人口動態統計，全国高齢者名簿，Piecewise ConstantIntensity Model，平均余命関数．
1.

まえがき

超高齢者の数が十分に多くなり，その年齢も信頼できるようになった．ここでは人口動態統計データ（死亡統計）および全国高齢者名簿（生存統計）を用いて，世代とともに高齢者の寿命分
布がどのように変化しているかを解析する．日本では平均寿命が着実に増加しているが，高齢者寿命分布ではどのような変化が起きているか調べたい．特に極値理論を利用して分布限界が有限，無限を含むモデルを考え，寿命分布の上限を推定する．
超高齢者（the oldest-old）の寿命については，種々の側面から興味がもたれている．生物（動物）学では加齢，老化現象を生理学的に，あるいは遺伝子の機能まで視野に入れて研究が行われている．このような研究ではヒトを対象とする必要はない．たとえば Carey et al.（1992）は，生活環境を変えることにより地中海ハエの寿命が延びた実験結果を示し，「生物の寿命が有限であることを立証しようとしたすべての研究は失敗している」と断言している．人間社会では，医療，医療制度への影響が直接的である．さらに高齢者の介護，健康，生活が問題であり，関連した労働，社会，ビジネスが考えられている．これら社会経済問題の場合には，高齢者の数，
年齢構成，その将来予測が主要な関心であるが，ここでは人口学的研究には立ち入らず，むしろ好奇心だけから寿命分布の上限を調べる．国勢調査などの作表では超高齢者を年齢区分の最高区分としているが，その限度を何歳にするかは関連した問題である．
1
2

高千穂大学 経営学部：〒168–8508 東京都杉並区大宮 2–19–1; sibuyam@takachiho.ac.jp
尚美学園大学 芸術情報学部：〒350–1153 埼玉県川越市下松原 655; nob-hanayama@jcom.home.ne.jp

118

統計数理 第 52 巻

第 1 号 2004

記録的な自然災害にたいする対策，製品の寿命または信頼性の保証のために極値理論が利用されている．確率標本の，ある閾値以上の条件付分布は，適当に標準化し閾値を大きくした極限で一般パレート分布に従うという，理論的結果がある．人口を世代（いわゆる出生時コホート）ごとに区分すれば，高齢まで生存した条件のもとで寿命が一般パレート分布に従うと仮定できるのではないか，というのが動機である．
死亡統計の粗い近似により，世代ごとの標本平均余命を計算すると図 1 のようにほぼ直線的に減少する．各線が 1 世代の標本平均余命であり，左端の点の横軸が 1980 年（もっとも古いデータ）における，その世代の年齢である．すべての世代について言えるのではないし，100 歳を越えると不安定になるが，直線性が特徴付ける一般パレート分布を検討する根拠となる．図

0

1

2

3

Mean residual life Ma 1880 - 1892

90

95

図 1.

100

105

110

死亡統計の世代ごと標本平均余命（年齢 年）
．

0

5000

10000

15000

#deaths older than x-axis Ma 1880 - 1892

90

95

図 2.

100

105

死亡統計の世代ごと生存者数（年齢 年）
．

110

119

年齢時代区分データによる超高齢者寿命分布の推測

2 は同じ死亡統計で，各世代の 1980 年以降の生存数である．この図のように超高齢者数は着実に増えているが，図 1 を見るとその平均余命はほとんど変わっていない．（これらの図については第 3 節の始めで説明する．
） 極値理論を用いた解析には Kaufmann（2001）などがあるが，結論が明確でなかった．人口学で利用しているモデル（たとえば Kannisto（1999））では寿命分
布に上限がない．
このように，高齢者の寿命を極値理論に基づいて議論することは有意義であるように思われるが，実際に我々が解析しようとしている高齢者寿命に関するデータは，年齢時代区分で作表されたデータ（いわゆる不完全データの一種）である．したがって，連続変量を前提として議論されている極値理論をこのデータの解析へ直接適用することはできない．そこで，本論文の要点は，連続変量に対して議論されている極値理論を年齢時代区分データにおけるコホート分析（Osmond and Gardener（1982））と結びつけるところにある．また，副次的に死亡統計と生存統計のコホートの違いを議論する．
2.

一般パレート分布

定義と諸性質
一般極値分布の吸引領域に属する確率分布の，ある閾値以上の条件付分布は，適当に標準化し，閾値を大きくした極限で一般パレート分布



(2.1)

F (x) = F (x; γ, a) =に従う. その分布範囲は

1 − (1 + γx/a)−1/γ ,1 − exp (−x/a) ,0 < x < ∞,γ = 0 ,γ = 0,γ ≥ 0,

0 < x < a/(−γ) ,γ < 0,

である. (2.1) の分布を GPrt (γ, a) で表わす．γ = 0 の場合は指数分布である．一般パレート分
布については，たとえば Coles（2001）参照．
確率変数 X ∼ GPrt(γ, a) の期待値は E(X) = a/(1 − γ), γ < 1, である．γ ≥ 1 では期待値が存在しない．災厄関数は
h(x) =f (x)1=,
1 − F (x)a + γx∀γ ∈ R ,である．X ≥ x0 の条件の下で，X ∼ GPrt(γ, a) の生存関数は
1 − F (x0 + y)=
P {X − x0 ≥ y|X ≥ x0 } =1 − F (x0 )



γy
1+
a + γx0

−1/γ,

y ≥ 0;

γ = 0 ,

ただし，γ < 0 のとき，y < a/(−γ) − x0 ．つまり X − x0 | X ≥ x0 ∼ GPrt(γ, a + γx0 ) であり，左裾を打ち切ったとき同じ型の分布に従う．γ = 0 の場合も同様である．
閾値を変えたときに分布の型が変わらないことは，一般極値分布の標本最大値の分布型が変わらない最大値安定性と対応している．この性質から条件付期待値（平均余命 mean residual
life，平均超過関数 mean excess function）m(x0 ) および条件付生存関数の逆関数 y(p; x0 )（p は生存確率）が次のように閾値 x0 の 1 次関数となる．
a + γx0, γ < 1,m(x0 ) := E(X − x0 |X ≥ x0 ) =1−γ←
1 − F (x0 + y)y(p; x0 ) :== (p−γ − 1)(a + γx0 )/γ , γ = 0 , ∀γ,
1 − F (x0 )





120

統計数理 第 52 巻

第 1 号 2004

ただし，{}← は，{} 内の関数の逆関数を表わす．逆に条件付期待値 m(x0 ) が閾値 x0 の 1 次関数であれば，積分方程式を解いて，分布関数 (2.1) を導ける．

­ が負か零か
データ {xj ; j = 1, 2, . . . , n} に (2.1) を当てはめたとき，γ < 0 か γ = 0 かが興味の焦点の一つである．それを見る簡便な方法は次の標本平均余命 m̂ のグラフである．これは閾値を決定するときの判断にも利用される.
m̂(x0 ) :=n

(xj − x0 ) I[xj ≥ x0 ]

j=1


n

I[xj ≥ x0 ] ,( I[true] = 1, I[false] = 0) .

j=1

γ < 0 の場合には，(2.2)

F (x; η, ω) = 1 − (1 − x/ω)η ,0 ≤ x < ω;

η = −1/γ , ω = a/(−γ) ,となる．つまり X/ω がベータ分布 Be(1, η) に従う．期待値は E(X) = ω/(1 + η), 条件付期待値（平均余命）は



γ
1−γ
−γ
ω − x0=(ω − x0 ) =1−γ
1+η

m(x0 ) = E(X − x0 |X ≥ x0 ) =である．負の勾配
1
−γ
=
=E
α :=
1−γ
1+η

X − x0ω − x0x0 +aγ| X ≥ x0

は生続率（perseverance parameter）と呼ばれている．災厄関数は
h(x) =f (x)η
=
1 − F (x)
ω



1−

x
ω

,

である．
寿命分布に上限 ω < ∞ があれば，極限分布で γ < 0 となるとは限らない．有界の上限でハザード関数 h(x) が必ず無限大となるが，その速度が遅いときだけ γ < 0 となる．
人口学では “force of mortality”（ここでは災厄関数）が指数的に増加する，無限区間上に分布する Gompertz 曲線が用いられている．これは負のグンベル分布で，最小値極限分布のひとつで γ = 0 である．高齢者の寿命分布に限ると極値理論の観点から Gompertz 曲線は不自然である．後述の Aarssen and de Haan（1994）も，この点を議論している．
最尤推定
一般パレート分布 (2.1) からの大きさ n の確率標本の対数尤度は
(2.3)

l=−

n


1+γ
log(1 + τ xj ) − n log a ,γ j=1である．これから最尤方程式はn

1log(1 + τ xj ) = γ ,
n j=1



n


1
xj
1+γ 1
= ,
γ n j=1 1 + τ xjτ

τ = γ/a ,121

年齢時代区分データによる超高齢者寿命分布の推測となる．第 1 式を第 2 式に代入すれば，τ だけが未知の方程式となる．γ < 0 のとき分布範囲がパラメータに依存して，最尤推定量の正則条件が満たされないが，このような場合を Smith（1985）が議論している．−1/2 < γ ならば Fisher 情報行列
I(γ, a) =n(2γ + 1)(γ + 1)



2
1/a

1/a
(γ + 1)/a2,

γ > −1/2 (or η > 2) ,が有限で，その逆行列，つまり最尤推定量の漸近分散共分散行列は
γ+1
Var((γ̂, â)) =
n



γ+1
−a

−a
2a2



,

γ > −1/2 (or η > 2) ,となり，正則条件と同じ議論が成り立つ．これから γ < 0 のときの分布上限 ω = a/(−γ) の（漸近）推定誤差はVar(ω̂) =a2 (γ + 1)2 (2γ + 1)
a2
≤
,
4
nγ
nγ 4

1
1
η
SD(ω̂) ∼ √
= √ ,
ω
n|γ|
n

γ → 0,となる（情報行列の計算は付録を参照）．γ = −1/2 の場合も最尤推定量の漸近有効性，漸近正規性が成り立つが，収束速度が速くなり (n log n)−1/2 の位となる. −1/2 < γ < −1 の場合，
最尤推定量が存在するが漸近正規性は成り立たない．漸近有効性も疑わしいが未解決である．
γ ≤ −1 の場合は一般に最尤推定量は存在しない．なお，Smith（1987）は，一般パレート分布の吸引領域に属する分布からの確率標本について，最尤推定量の精度を議論している．
閾値を決定するには，標本平均余命関数のグラフ，および閾値を変えた時の推定値・推定誤差のグラフを描いて判断する．
パラメータ γ の値に依存する一般パレート分布の特徴を表 1 にまとめておく．また，この節で述べたパラメータの意味と相互関係を付録でまとめておく．
3.

年齢時代区分データとそのモデル

3.1 年齢時代区分データ
一般に年齢時代区分によって与えられた統計表については，項目（セル）の番号を (i, j), i =1, . . . , I; j = 1, . . . , J で表わすが，この節では添字 i により年齢を，j により年度を表わす．つまり，例えば 1970–2000 年の 100–120 歳について与えられた統計表では，(i, j), i = 100, . . . , 120; j =1970, . . . , 2000 によってセルを表わす．利用可能なデータの詳細についても後で述べる．性別，
出生年度別コホートを一様な集団とみなすことを，再び強調しておく．

生存者数統計
j 年度の一定の調査日に i 歳で生存する人々を考え，その人数を Mij により表わすことにする．
（時間，年齢）平面の上で傾斜 45 度の線分で 1 人の生涯を表わすレキシス図（Lexis diagram,
Keiding（1999））を用いると，Mij は調査時点を表わす縦線と，斜線の交点の数である（図 3）．
表 1.

一般パレート分布の特性．

122

統計数理 第 52 巻

第 1 号 2004

図 3.

生存統計のレキシス図．

図 4.

死亡統計のレキシス図．

このグループの人が 1 年以内に死亡する事象を Vij とする．その人数は(3.1)

Mij − Mi+1,j+1 =: Dij ,である．k = j − i で出生コホート年度を表わす．
Mij に数えられる人々の年齢は [i, i + 1) に分布している，つまり，j 年度の調査時点までの生存時間が同一ではない．このため，厳密には「出生年度別コホートを一様な集団とみなす」という仮定では Dij に数えられる人々の死亡が同一分布に従うと考えられない．しかしここでは，Dij を，k − j − i 年度に生まれ満 i 歳以上 i + 1 歳未満で死亡した人数の近似値と考える．
高齢生存者数を毎年記録しているので Mij = 0 を満たす j を確認できれば，コホート全体の死
亡データが確定する．
死亡者数統計人口動態統計では，第 j 年度に i 歳以上 i + 1 歳未満で死亡した，という事象 Wij の人数
Nij が記録されている．レキシス図で見ると，Nij は個人の死亡を表わす斜線の端点が正方形の領域に落ちる数である（図 4）．
このときには出生コホートが不明になる．もっとも簡単な解決法は Nij を 2 等分して，あい続くコホートに振り分けることである．生存者数統計の，下記の確率 q(k, i) の事象の数を

123

年齢時代区分データによる超高齢者寿命分布の推測(Nij + Ni,j+1 )/2 により近似する．まえがきで述べた図 1，図 2 はこのような近似計算で作成した．他の方法については，推測と関連して述べる．コホート全体が死亡したかどうかは，データから判断するが，生存統計が判断に役立つ．
3.2 生存者数統計のモデル
j 年度に i 歳であるコホート（この集団を k = j − i で区別する）に関するデータ Mij にたい
して，P (Vij ) の近似である確率q(k, i) := Fk (i + 1) − Fk (i) = F k (i) − F k (i + 1) ,(3.2)

k = j − i,の多項分布を当てはめる．Fk (t) は GPrt (γk , ak ) の左打切り分布，F k (t) はその生存関数である．簡単のために添え字 k を省略し，γ < 0 の場合に限って，(2.2) のパラメータを用いると
1,q(i) = F (i) − F (i + 1) ,

対数尤度は
l=

k


F (t) = F (t; ω, η) :=mi log q(i) ,t < 0,(1 − t/ω) ,

0 ,

η

0 ≤ t < ω,
ω ≤ t.

mi = Mij − Mi+1,j+1 ,i=1

尤度方程式は

k


∂l
mi ∂
=
(F (i) − F (i + 1)) ,
∂ω
q(i)
∂ω
i=1
ηt
∂
F (t) = 2
∂ω
ω



1−

t
ω

η−1

,

k


∂l
mi ∂
=
(F (i) − F (i + 1)) ,
∂η
q(i)
∂η
i=1



∂
t
F (t) = F (t) log 1 −
∂η
ω



,

0 < t < ω,となる．最尤法が有効であるためには，ω が 1 に比べて十分大きいことが必要である．
3.3 死亡者数統計のモデル
k 年度に誕生した人が，j 年度に死亡したとし，その誕生日（時点），死亡日（時点）を，それ
ぞれ
k+v,

j +w,

ただし v, w は, 独立な（0,1）一様確率変数と仮定する．連続化した死亡年齢 y は(3.3)

y = j + w − (k + v) = (j − k) + (w − v) = j − k + u ,となる．u = w − v は (−1, 1) 上の対称三角分布である．したがってこの人の死亡年齢は j − kまたは j − k − 1 である．つまり j 年度に i 歳で死んだという事象 Vij には，j − i コホートの人と j − i + 1 コホートの人が混じっている．このため、コホート毎に i.i.d. を仮定するだけでは，一般パレート分布のパラメータを推定することはできない．
以下では，以上の問題に対処する方法として，ポアソン点過程モデルを考える方法，擬似確率標本により推定する方法および表から直接推定する方法を提案する．
ポアソン点過程モデル
死亡点から出生点が自動的に定まる．死亡点をポアソン過程とみなし，その強度の変化が緩やかであると想定する．レキシス図の上で死亡強度が断片（piece）ごとに一定であることを仮定

124

統計数理 第 52 巻

第 1 号 2004

した Piecewise Constant Intensity Model（PCI モデル）を適用する（Keiding（1990））．ただし断片の一点での強度が一般パレート分布のハザード関数から定まるとする．これも近似法であるが，競合リスクを分離して推定できる点が重要である．
年齢時代区分データによる寿命分布解析における PCI モデルを構成するため，1 年度 1
歳単位の正方形領域 9 個からなる断片達 Vrs を次のように定義する．まず 2 重添字の集合
N 2 = {(i, j) : i, j ∈ N } の中の部分集合を



Mrs =

(r + 1 , s + 2 ) ,s = k + r; r, s = 0, 3, 6, . . .

1 ,2 ∈{0,1,2}で定義する．次にレキシス図上の上記断片を
Vrs =



Wij ,

j = k + i; r, s = 0, 3, 6, . . .

(i,j)∈Mrsにより定義する．この断片は k(= r − s = 0, 3, . . . ) 世代 k + 1 世代を完全に含み，k − 1, k + 2
世代の半分を含む．
断片 Vrs 上の強度が一様であるとし，λrs で表わす. 一方コホートごとの寿命分布が一定で，閾値以上の条件付分布が一般パレート分布に従うならば，その k 世代の災厄関数（強度）は前述の通り，h(x; k) = 1/(ak + γk x) である．その Vrs の中心での値を
λrs = h(r + 1, s − r) ,(3.4)

r, s = 0, 3, 6, . . .

とする．Vrs で僅かでも生存した全ての人々の生存時間の累和 xrs+ を，person-years と呼ばれる各年度各歳の生存人口数 Mij の累和で近似すると，点過程の尤度は以下のように表わされる．
L∝ d

λrsrs exp(−xrs+ λrs )

 d
r,s

≈

λrsrs exp(−λrsr,sMij ) ,

(i,j)∈Mrs

drs =



Nij

(i,j)∈Mrs

ただし Mij は Vij で僅かでも生存した人の数を表わし，(3.1) と少し違う．これをパラメータに関して最大化する（Keiding（1975））．
擬似確率標本による推定
一般パレート分布からの確率標本に基づく推定法（たとえば最尤法）を適用する簡便法として，
区分計数データから ‘復元シミュレーション’ により，擬似確率標本を生成する．(3.3) 式より，
事象 Wij の出現数 Nij の個数の三角分布乱数 u を独立に生成し，寿命 y = j − i + u を u > 0
ならば k = j − i コホート，u < 0 ならば k − 1 コホートに属するデータとする．死亡日が既知であっても季節変動があるならば，それを調整するために擬似化したデータが必要である．
表からの直接推定
事象の確率をより詳しく書き，その組合せで推定方程式を構成できる．ただしその解法は容易ではなく検討中である．式 (3.3) に戻って，出生年度 k，死亡年齢（日計算）y の人（死亡年齢（年計算）は，整数部分 [y] = i）の死亡年度がいつになるかを考えると死亡日が k + v + y であるから，死亡年度が j となる条件は
j ≤k+v+y <j+1⇔

i ≤ y + v < i + 1.

125

年齢時代区分データによる超高齢者寿命分布の推測

これと [y] = i を合わせると，y についての条件

1

死亡年度 = j ⇔ i ≤ y < i + 1 − v ,(3.5)

ρ0 (k, i) :=0

Fk (i + 1 − v)dv − Fk (i) ,が得られる．Fk (t) は (3.2) の下で定義した．ρ0 (k, i) は，条件式の確率である．同様に(3.6)

死亡年度 = j + 1 ⇔ i + 1 − v ≤ y < i + 1 ,ρ1 (k, i) := Fk (i + 1) −

10

Fk (i + 1 − v)dv ,となる．ρ1 (k, i) は，条件式の確率である．出生年度 k − 1 の人が i 歳で j 年度に死亡する確率は，この議論で k を k − 1 に変えるだけでよく ρ1 (k − 1, i) である．
注意.

この記号を用いると P (Vij ) は (3.2) の q(k, i) ではなくP (Vij ) = Rk ρ1 (k, i) + Rk−1 ρ0 (k − 1, i + 1) ,である．Rk−1 , Rk は k − 1 世代 k 世代の閾値における生存比である．
(3.5), (3.6) の積分は陽に表わせる: 指数分布 (γ = 0) の場合，

10F (t − v)dv = 1 −

1
0



ベータ分布 (γ < 0) の場合，

1ωF (t − v)dv = 1 −
1+η
0

(3.7)

1

exp(−(t − v)/a)dv = 1 − a exp

a



 t

− 1 exp −a,
1+η  ω − t 1+η 
ω+1−t= F (t) + O(ω

−

ω
−2

),

ω

ω → ∞.

尤度の計算には，閾値における各世代の生存比が必要であるが，観測比率をそのまま既知パラ
メータとする．
簡単のため k − 1, k の 2 世代を含む 3 時代区分データ (Nij , Ni,j+1 , Ni,j+2 ), k = j − i, i =0, 1, . . . , から 3 組のパラメータ (γk−1 , ak−1 ), (γk , ak ), (γk+1 , ak+1 ) を推定する．この 2 世代のi 歳での死亡確率を 4 分割し，対応する観測数を用いた対数尤度

(Ni,k+i log(Rk ρ0 (k, i)) + Ni,k+i+1 (log(Rk ρ1 (k + 1, i)) + log(Rk+1 ρ0 (k + 1, i))i+Ni,k+i+2 log(Rk+2 ρ1 (k + 2, i))) ,から推定する．ただし Rk , Rk+1 , Rk+2 は k, k + 1, k + 2 世代の閾値における生存比であり，

ρ0 (k, i) =
ρ1 (k, i) =1
1 + ηk(i + 1 + ηk − ωk ) 1 −


ω −i

η
i

1 + ηk

ωk

k

1−

k

−

i
ωk

η


ω −i+η
kk

1 + ηkk− (i + 1 − ωk ) 1 −

1−

η
i+1

i+1
ωk

η 
k

,

k

ωk

,

である．
世代数を増やすと，同時に解かねばならないパラメータ数が増え数値解法が困難になる．これはコホートモデルで絶えず生じる困難である．

126

4.

統計数理 第 52 巻

第 1 号 2004

実データの解析

この論文ではデータの理解と解析法の検討に重点を置き，男女の比較には至っていない．女性のほうが寿命が長いために観測数が多く，そのために死滅を確認できる世代が限定されるという矛盾がある．以下の解析は男性データに限る．
ギネス認定の，最近の世界最高齢者は泉重千代さん（1865–1986, 120 歳），ジャンヌ・カルマ
ンさん（1875–1997, 122 歳），本郷かまとさん（1887–, 115 歳女性最高），中願寺雄吉さん（1889–,
114 歳男性最高）などである．以下の解析で泉さんは，同世代の記録が少ないために含まれていない．（本稿査読中に中願寺さん，本郷さんが亡くなられた．
）
4.1 生存者数統計の解析
厚生労働省（厚生省）は 1964 年以来毎年 9 月に，100 歳以上の全国高齢者名簿，いわゆる長寿番付，を発表している．その付表で 100 歳 以上 1 歳刻みで男女別高齢者数を数えている．表2 は 1964–2002 年の男性生存者数である．ただし，これは名簿（氏名，住所，誕生日）で，厚生労働省の要請にたいする市町村の自由な回答の収集である．各市町村がどのように調べている
か，統一されておらず不明である．東京都の市では，住民基本台帳から自動的に選出して報告しているようであり，死亡統計と比較すると，最近ではこのような処理が支配的であると予想される．世代ごとに分けると，約 50 人以上のデータが得られる世代の誕生年度は 1873–1902年である．全死亡者を死亡統計と比較できる世代は 1880–1890 年に限られる．
死亡統計は歴年単位の集計である．生存統計は 7–8 月頃の調査である．そのために (3.2) 式の Dij ; j = k + i, は Nij , Ni+1,j , Ni,j+1 , Ni+1,j+1 のそれぞれの一部である．これらの間の違いを見るために，図 5 では 1880 年誕生世代男性について Dij ; j = k + i, i = 100, 101, . . . を横軸に,Nij , Ni+1,j , Ni,j+1 , Ni+1,j+1 を縦軸にプロットした．それぞれの点を ‘+ - * ％’ で表わす．調査時点が年後半であるため，Dij が Nij より Ni+1,j+1 に近いのは当然かもしれない．少なくともこの時期では地方自治体の報告が，完全ではなかったとも判断できる．

compare death numbers, 1890

100

150

*
+

N[i,j]

-

*
+

50
0

%

*+

%

*

+
**
*+%
-+%
*-+%
*-+%
%
0

*

+ %
-

%
+
-

%

50

100

150

dM[i,j]

図 5.

生存統計と死亡統計の比較．+ : N (ij), − : N (i + 1, j),  : N (i, j + 1), ％ : N (i + 1, j +1)．

年齢時代区分データによる超高齢者寿命分布の推測
表 2.

年齢時代別 100 歳以上男性生存者数．全国高齢者名簿より．

127

128

統計数理 第 52 巻
表 3.

第 1 号 2004

生存者統計を用いたパラメータ推定結果．

1

2

3

4

Mean residual life, Male , 1880 - 1890

100

102

図 6.

104

106

108

110

112

生存統計の標本平均余命（年齢 年）．

表 3 に，100 歳以上の各人の寿命が同一のパレート分布に従うと仮定し，生存者統計にモデル (3.2) を当てはめた結果を示す．ただし，110 歳以上の死亡者数，死亡確率を統合し推定した．解析に用いた 1885–1887 年の世代は，死亡統計に PCI モデルを当てはめた場合と対応している．この世代について，γ の推定値およびその標準偏差（フィッシャー情報行列から算出）
から γ = 0 の仮説は棄却される．また，−a/γ の推定値から寿命分布の限界は 134 歳と推定された．先に述べた実際の長寿記録よりかなり大きな値であるが，利用できる最も新しい世代である 1890–1892 年の世代を用いた推定もほとんど同様の結果となる．一方，最も古い 3 つの世代である 1864–1866 年の世代を用いた推定では，γ の推定値が正となり，標準偏差から γ = 0の仮説は棄却されている．図 6 の標本平均余命でも，減少傾向があるが水平に近く，これらの解析結果と一致している．
4.2 死亡者数統計の解析
厚生労働省（厚生省）が毎年発表している人口動態統計では，85 歳以上（年度によって 90 歳
以上）の原因別死亡者数が 1 つの年齢区分にまとめられている．しかし，同省収蔵の保管第 3
表を閲覧すれば原因別各歳死亡者数が得られる．表 4 は 1980–2001 年の 84–117 歳の死亡者数である．図 1 を見ると，83 歳以下についても，平均余命関数が直線を示すことが予想されるため，83 歳以下のデータも加えるべきだろう．

年齢時代区分データによる超高齢者寿命分布の推測
表 4.

年齢時代別男性死亡者数．

129

130

統計数理 第 52 巻

第 1 号 2004

表 5. Piecewise Constant Intensity Model によるパラメータ推定結果．

表 6. 擬似確率標本による GPrt の推定．

世代に分けると，約 50 人以上のデータが得られる，死亡確定できる世代の誕生年度は，男性が 1880–1892 年，女性が 1878–1896 年である．
表 5 に，レキシス図上の断片 Vrs ; r = 95, 98, . . . , 113, s = 1980, 1983, . . . , 1998 に PCI モデル(3.4) を当てはめた結果を示す. この断片領域は全員が死亡している 3 年連続した世代のなかで最も新しいもので，1885–1887 年度世代である．γ の推定値およびその標準偏差から γ = 0 の仮説は棄却され，分布の上限 ω = −a/γ の推定値から寿命分布の限界は 125–127 歳と推定された．この結果は生存者数統計から推定されたものより小さいものである．
先に述べたように，PCI モデルの適用は，競合リスクを分離して推定できる利点がある．そこで，癌死亡について同様のモデルの適用を試みた結果，γ > 0 と推定され，その標準偏差か

年齢時代区分データによる超高齢者寿命分布の推測

131

ら γ = 0 の仮説は棄却された（表 5 を参照）．つまり，超高齢者に関しては，癌による死亡危険率は年齢の増加とともに減少すると推定された．これは，ある高年齢を超えると，癌での死亡率は減るという一般の認識と一致する．
擬似確率標本による解析結果
前述のように Nij 個の 3 角乱数により j −i 年度，j −i−1 年度のデータを生成する．1880–1892
年度世代より 1880–1891 年度世代データを得る．このデータの標本平均余命も，図 1 と同様に
100 歳頃までは常に直線的であるが，それ以上になると不安定であって，水平に近くなる世代もある．
表 6 に，Coles（2001）の関数 gpd.fit（）による推定の結果をまとめる．閾値は，利用可能なす
べてのデータを用いる場合と，高齢者名簿とも対応する 100 歳を選んだ場合を表に含めた．同じ世代の 2 列の，上段と下段である．推定値の標準偏差は漸近分散に推定値を入れて求めた．
p 値の計算も同様である．
p 値がかなり大きい場合もあるが，γ の推定値はすべて負である．ω の推定値がかなり大き
いもの小さいものがあるが，閾値を 100 歳に制限しない場合はだいたい 118–135 に納まっている．古い世代はデータが少ないこともあり，推定区間の幅が広い．
なお 1889 年度世代は，中願寺さんが生存中であるため対象外であるのを，あえて存命を除外して解析した．また推定誤差には，擬似確率標本作成の影響を考慮していない．
5.

終わりに

Aarssen and de Haan（1994）は 1877–1881 年に生まれ，1971 年 1 月 1 日に生存しており，オ
ランダ住民として死亡した合計約 1 万人の性別，日単位享年のデータを解析した．コホートを作るため，1881（, 80, . . . , 77）年生まれで 1971 年 1 月 1 日に 90（, 91, . . . , 94）歳に達していない人を除いた．また 1 年内の高齢者死亡率に季節性があり，冬に高く夏に低いので，一様分布となるように調整した．順序統計量の対数に基づくモーメント推定量を用いて，γ = 0 の仮説を棄却し，ω の推定量 113–124, α の推定量 0.07–0.14 を導いている．
Kaufmann（2001）は 1993 年西ドイツ女性 90 歳以上の年齢別死亡数を解析して，γ = 0 の仮説を否定できないと主張している．
このように極値理論を適用した研究結果でも γ < 0, γ = 0 の両論がある．ただし Kaufmann（2001）は年度のデータを用い，世代の寿命を調べていない．Aarssen and de Haan（1994）の結論は本稿と同じである．しかし著者たちは彼らの結論をさらに裏付けたと確信しているわけでもない．
現在，高齢者の死因では癌と循環器病が大きな割合いを占めており，その割合が死亡年齢により異なる．全死亡者中の癌死亡者の割合を見ると年齢とともに減少の傾向がある．しかも前
述のように循環器病死亡者数の割合が多い 100 歳以上では平均余命が減少しない．
このように対抗する死因（competing risks）が混合している場合，極値理論に立つならば，さらに閾値を大きくして，平均余命関数が直線となるように選ぶことになるが，現実にはデータが不足で確認できないことになる．医療技術の発達により癌死亡（の年齢分布）が変化するだ
ろうし，このような意味で，一般パレート分布（ベータ分布）が当てはまるという結論も一時期の判断となる可能性がある．
この論文で使用した元のデータおよび世代別に分け直したデータを，要求があれば配布することができる．

132

統計数理 第 52 巻

第 1 号 2004

付 録
一般パレート分布の情報量
次のモーメントを計算しておく．
E((1 + τ X)r ) = (τ /γ)∞0(1 + τ x)r−1−1/γ I [1 + τ x ≥ 0]dx

= 1/(1 − γr) ,r < 1/γ,

τ = γ/a .

特にE(1/(1 + τ X)) = 1 − E(log(1 + τ X))/(γ + 1) = 1/(γ + 1) ,
2

2

2

2

γ > −1 ,2

E((X/(1 + τ X)) ) = E((1 − 1/(1 + τ X)) /τ ) = 2γ /(τ (2γ + 1)(γ + 1)) .

対数尤度 (2.3) より，n = 1 のとき，
γ
1+γ
∂l
1
x
= 2 log 1 + x −,∂γγaaγ 1 + γx/a

x
∂l
1+γ
1
=
− ,
∂a
a2 1 + γx/a
a

x
2
γ
2
1+γ
∂2l
= − 3 log 1 + x +
+ 2
∂γ 2
γ
a
aγ 2 1 + γx/a
a γ
1
1+γ
x
∂2l
= 2
−
∂γ∂a
a 1 + γx/aa3

2(1 + γ)
x
∂2l
=−
+
∂a2
a3
1 + γx/ax
1 + γx/a2γ(1 + γ)a4,x
1 + γx/ax
1 + γx/a

2
+

2

,

1
,
a2

これらの x を X に変え，上記のモーメントを用いて期待値を計算する．
一般パレート分布のパラメータ
一般パレート分布 (2.1) のパラメータ (γ, a) は種々の別表現（parametrization）があり，第 2
節でも 6 つのパラメータについて述べている．ここでそれらの意味と相互関係をまとめておく．
まず γ は極値理論で裾の指数（tail index），極値指数（extreme index）と呼ばれる重要なパラメータである．分布関数を表わすには η の方が自然であるが，指数分布を γ = 0 の場合として含むことと，γ の大きい方が裾が重いというわかりやすさで γ が標準的となった．a はこの表現での尺度パラメータである．γ < 0 の場合に限ると分布範囲が有限区間となるので，上限
表 7.

3 パラメータ間の関係．

表 8.

4 パラメータ間の関係．

年齢時代区分データによる超高齢者寿命分布の推測

図 7.

133

4 パラメータ間の関係．

ω と，確率変数を ω で割って得られるベータ分布のべき乗パラメータ η が分かりやすい．平
均超過関数（平均余命）の直線性が一般パレート分布の特徴なので，その勾配（の負数）α は応用
上分かりやすいパラメータであり，推定方法によっては役立つであろう．最尤方程式の計算では，和の中のパラメータを一つにする τ を導入すると計算しやすくなる．
表 7 は 3 つのパラメータ γ, η, α の間の関係で，その 1 つで他の 2 が定まる．これらの代表を γ とする．表，図では γ によって，代表であることを強調している．(γ, a) のうちの一方を表 8 にしたがって ω ，または τ に変えることができる．言い換えると，4 つのパラメータγ, a, ω, τ から 2 つを択ぶ 6 通りの組合せのうち，図 7 の示すように (ω, τ ) を除く 5 通りの組
合せにより分布は定まる．

参 考 文 献Aarssen, K. and de Haan, L.（1994）. On the maximal life span of humans, Mathematical Population
Studies, 4, 259–281.
Carey, R. J., Liedo, P. Orozco, D. and Vaupel, W. J.（1992）. Slowing of mortality rates at older ages
in large medfly cohorts, Science, 258, 457–461.
Coles, S.（2001）. An Introduction to Statistical Modeling of Extreme Values, Springer-Verlag, London.
Kannisto, V.（1999）. Trends in the mortality of the oldest-old, Statistics, Registries and Science:
Experiences from Finland （ed. J. Alho）, 177–194, Statistics Finland.
Kaufmann, E.（2001）. About the longevity of humans, Section 14.2, Statistical Analysis of Extreme
Values（written by R. D. Reiss and M. Thomas）, 2nd ed., Birkhaueser, Basel.
Keiding, N.（1975）. Maximum likelihood estimation in the birth-and-death process, Annals of Statistics, 3, 363–372.
Keiding, N.（1990）. Statistical inference in the Lexis diagram, Philosophical Transactions of theRoyal Society of London, Series A , 332, 487–509.
Keiding, N.（1999）. Lexis diagram, Encyclopedia of Statistical Sciences（eds. S. Kotz, C. B. Readand D. L. Banks）, Update Volume 3, 411–413, Wiley, New York.
Osmond, C. and Gardener, M. J.（1982）. Age, period and cohort models applied to cancer mortality,
Statistics in Medicine, 1, 245–259.
Smith, R. L.（1985）. Maximum likelihood estimation in a class of nonregular cases, Biometrika, 72,
67–90.
Smith, R. L.（1987）. Estimating tails of probability distributions, Annals of Statistics, 15, 1174–1207.

134

Proceedings of the Institute of Statistical Mathematics Vol. 52, No. 1, 117–134 (2004)

Estimation of Human Longevity Distribution Based on
Tabulated Statistics
Masaaki Sibuya(Faculty of Business Management, Takachiho University)

Nobutane Hanayama(Faculty of Art and Information, Shobi-gakuen University)(Age, period)-specific data for the oldest-old survivors and deaths are analyzed using
the extreme value theory and the limit of longevity dstribution is discussed. The data
were obtained from National Oldest-old Survivors List and Population Movement Statistics by the Ministry of Health and Labor in Japan. In applying the theory of extreme value
statistics of continuous variables to the analysis of (age, period)-specific tabulated data,
we propose a procedure for applying the multinomial distribution model based on probabilities calculated from the generalized Pareto distribution, and another procedure using
a continuous pseudo random sample generated by adding random numbers. Furthermore,
the piecewise constant intensity model on the Lexis diagram, widely used for such data,is also applied. ML estimates show the finite upper limit of longevity distribution.

Key words: Cohort analysis, maximum likelihood estimation, mean residual life function, national
oldest-old survivors list, Piecewise Constant Intensity Model, population movement statistics.