Smith, T.R., Wakefield, J. (2016) A Review and Comparison of Age-Period-Cohort Models for Cancer Incidence. Statistical Science, 31(4), 591-610.
いわゆるコウホート分析についてのレビュー。どうやら私の能力と根性を超える論文だが、都合により無理やり読んだ。主旨は以前に読んだ松本(2019)と似ている。
1. イントロダクション
人間の疾患についてのAge-Period-Cohortモデル(APCモデル)について考える。時代効果とコホート効果は外的要因への曝露の代理変数である。
[APCモデルの歴史。メモ省略。人口学・社会学からはじまったわけだが、ガン発症率への適用は80年代からなのだそうだ]
2. Motivating Data
Washington State Cancer Registryから得たガン発症数データを用いる。ガンのリスクは基本的に年齢と共に増大する。cohortは発がん性物質への曝露の代理変数となる。年代は診断手続きの異質性を捉えている。
乳がん・肺がん(男性)については、年代は11水準(30-34歳から80-84歳まで)、時代は4水準(1992-1996年から1997-2001年まで)。子宮頚がんについては、年代は13水準(20-24歳からはじめる)、時代は4水準。
3. APCモデルの基礎
3.1 表記
ふつうAPCモデルは、発症率なり死亡率なりの対数について加法的なモデルである。
年齢集団\(i=1,\ldots,A\), 時点\(j=1,\ldots,T\)の発症数を\(y_{ij}\)、その行列を\(Y\)とする。コホートのインデクス\(k\)は年齢と時代の関数となる。たとえば、年齢と時代が同じ時間尺度で測られていたら(5年刻みとか)、\(k =A-i+j\)と定義すれば\(k \in \{1, \ldots, A+T-1\}\)となる。年齢x時代の行列において、同じコホートは左上から右下に斜めに並ぶ。
リスク集合のサイズを\(N_{ij}\)としその行列を\(N\)とする。なお、よく考えてみるとリスク集合のサイズの決め方は一筋縄ではいかない。たとえば、2007-2011年の30-34歳のサイズとは？ 単純には、2007年の7月1日人口、2008年の7月1日人口, … の合計であろう。しかしそこには2007年の前半で死んだ人とか、2007年の後半で30歳になった人とかが含まれていないし、年の後半で35歳になったり死んだりした人が完全に含まれてしまっている。出生率とか死亡率とか移民とかが年齢-時代セルのなかで一定であれば、行って来いになるから気にしないでよいんですけどね。本論文では深入りしません。
3.2 APCモデルの識別
基本的なAPCモデルは次の通り。$$ \log E\left[ \frac{y_{ij}}{N_{ij}} \right] = \mu_{ij} = \delta + \alpha_i + \beta_j + \gamma_k$$ \(\delta\)は全体的な生起率の対数、\(\alpha_i, \beta_j, \gamma_k\)は相対リスクの対数として表現された年齢効果・時代効果・コホート効果だ… と解釈したくなるのが人情だが、そういう直接的な解釈は難しい。なぜならモデルがover parametrizedだからである。理由はふたつある。
- 因子を持っているモデルはいつもそうなのだが、切片のあるモデルでは推定可能な水準数よりひとつ多い水準があることになるので制約が必要である[言い方がよくわからんのだけど、前後の文脈からいえば、たとえば1要因A水準の分散分析モデル\(y_{ij} = \mu + a_i + e_{ij}\)で、\(\mu, a_1, \ldots, a_A\)をすべて自由推定するのは無理だという話であろう]。典型的な解決策は合計ゼロ制約である。
- ３つの因子の間に線形従属性がある。もし上の式のパラメータを直接に解釈したいなら、なんらかの仮定を置かないといけない。それらの仮定はローデータだけでは検証できない。
パラメータ\(\delta, \alpha_1, \ldots, \alpha_A, \beta_1, \ldots, \beta_T, \gamma_1, \ldots, \gamma_k\)を縦に積んで \(\theta\)としよう。適切な計画行列\(X\)を定義すれば\(\mu = \mu(\theta) = X’\theta\)と書ける。ところが\(X\)はランク落ちしている。ひとつには一般的な因子問題のせい、もうひとつにはコホート効果は年齢効果と時代効果に線形従属しているせいである。
というわけで、年齢効果・時代効果・コホート効果をすべて識別することはできない。
[どういう線形従属性があるのかという説明。メモ省略]
3.3 識別問題に対してこれまでに提案されてきた解決策
最初に試みられたのは、年齢・時代・コホートの効果をなんらか制約して識別に持っていこうというアプローチである。
\(\sum_i \alpha_i = \sum_j \beta_j = \sum_k \gamma_k = 0\)とすれば\(\delta\)は識別できる。しかし\(X\)はまだ識別できない。\(X\)のランクが4つ落ちているのが1つ落ちているのに変わっただけだからだ。もうひとつ制約しないといけない。
そこで、時代効果の一階差分 \(\beta_2 – \beta_1, \ldots, \beta_T – \beta_{T-1}\)の合計を0とするとか(結局\(\beta_T = \beta_1\)となる)、どこかの隣り合うパラメータを等しくする(\(\gamma_1 = \gamma_2\)とか)といった制約が提案された。どこを制約するかで結果が全然変わってきちゃうのが弱点だが[実例を示している]、実質的知識で制約を決められる場合もある。子宮頚がんの発症率は45歳以上では変わらないはず、とか。
ICAR事前分布に基づくベイジアンアプローチ(6.2で後述)では、合計ゼロ制約と共に、３つの効果における線形トレンドに対する暗黙的な罰則を含んでいる。層別死亡率のモデリングへの拡張においては、一部の効果が層を通じて同じだと仮定する限りにおいて、一部の時間効果の対比を識別できるということが示されている(4.3で後述)。
年齢・時代・コホートの関数を因子とみるのではなくて連続関数とみることもできる。一般的な因子問題による識別不能性は回避できるが、線形従属性の問題は残る。[具体的にどういう整形従属性があるかという例。メモ省略]
たとえば、コホート効果がないと仮定して年齢と時代だけの線形モデルを組めば識別はできる。こういうモデルを「ドリフト」モデルという。年齢・時代のモデルは年齢・コホートのモデルと同じ適合度となる。
連続関数の選択肢のひとつに罰則付きスプライン関数がある。RのEpiパッケージでは、リファレンスとする特定のコホートについて効果ゼロとし、時代の関数を平均ゼロ・傾きゼロと仮定している。
このように、個々の時間尺度上で連続関数をモデル化するというのは魅力的ではあるが、なにが推定可能なのかについての透明性が失われる点、そして仮定がアドホックに思われるのが難点である。
識別問題を「解決する」アプローチは他にもいろいろあって…
- Robertson & Boyle (1986): 個人レベルの記録にアクセスできる場合の方法。
- Lee & Clater (1992): 年齢別の時代トレンドと年齢別のコホートトレンドを許す。しかし一意な解を得るためには線形制約が必要。
- Yang, Fu, Land(2004): 線形制約を決めるために\(X\)のnull spaceを使うintrinsic 推定アプローチ[←さっぱりわからん]。科学的に正当化するのは依然として難しいという批判もある。
- O’Brien(2000): 2因子モデル。ないし、予測子つきの2因子モデル。たとえば肺がんの発症率で、年齢、時代、喫煙率を予測子にする、とか。
識別不能性という問題を踏まれると、３つの効果についての推定値を図表で示す際には十分な注意が必要である。いかなる推定値もなんらかの制約に条件づけられており、その制約はたいてい実質科学的基盤を持たず、データだけから検証することもできない。ただの探索ツールと割り切るか、推定可能な時間効果の関数である部分だけに焦点をあてるとよい。たとえば、効果の二階差分は変換に対して不変であり識別可能である。二階差分は効果のlocal curvatureを示し、指数尺度上で二時点間の相対リスク比になっている。
APCモデルの識別可能性問題を「解決」するアプローチのひとつとして、３つの効果のパラメータについての推定可能な関数という観点のみによってモデルを表現するというアプローチがある。たとえばHolford(1983)は、それぞれの効果の集合を、線形効果と曲線効果に分割している。曲線効果の線形結合は推定可能である。また、３つの効果の傾きからなるある種の関数も推定可能である。３つの効果の線形傾きを\(s_\alpha, s_\beta, s_\gamma\)として、\(u s_\alpha + v s_\beta + (v-u) s_\gamma\)という形式の関数は推定可能である。
Rosenberg & Anderson (2011)は、縦断研究やクロスセクション研究における年齢トレンドのような多くの疫学的要約が、Holfordのモデルでいうところの、パラメータの推定可能関数として表現できることを示している。
[よくわからんが、パラメータを識別できなくても研究の役には立つよという話かしらん]
4. APCモデルの拡張
4.1 予測
[パス]
4.2 年齢の時間間隔と時代の時間間隔が異なる場合
たとえば年齢が5年単位で5水準、時代が年単位で10水準としよう。1年目の年代1が6年目の年代2ということになる。こういう場合、時代効果に対してM個おきに定数を足し、コホート効果に対してM個おきに同じ定数を引いても、全体の生起率対数は変わらないことになる。このように、年齢の時間間隔と時代の時間間隔が異なるとさらなる識別不能性が生じる。
Holford(2006)はこの問題に扱うために、ミクロ時間尺度(年齢と時代との間で時間間隔が等しい)とマクロ時間尺度という考え方を導入している。
[よくわからんけど、なんかうまく解けそうな問題であるような気もしますね。状態空間モデルにして、状態変数はすべて年単位とし、観察方程式のほうで年齢の5年集計を表現すればよいのではなかろうか…]
4.3 層別APC
[性別のような時間独立な層別変数を取り込むという話。必要になったら読もう、ということでパス]
5. 近年のAPCモデルの適用についての批評
[論文の実例を挙げて、誤った解釈を批判している… 怖い…]
6. APCモデルのパラメータ化
[さあ、ここからが本番だ。気を取り直してがんばろう]
6.1 Nielsenのパラメータ化
Martinez Miranda, Nielsen and Nielsen (2015 JRSS)の定式化を紹介する。
\(\Delta\)を差分オペレータとする(\(\Delta \alpha_i = \alpha_i – \alpha_{i-1}\))。[原文の表記だと目が回るので、以下しばらくの間 \( s_\alpha \equiv \mu_{A1} – \mu_{A-1,1}, s_\beta \equiv \mu_{A2} – \mu_{A1} \)と略記する] $$ \mu_{ij} = \mu_{A1} + (i-A) s_\alpha + (j-1)s_\beta + a_{ij} $$ $$ a_{ij} = \sum_{t=i}^{A-2} \sum_{s=t}^{A-2} \Delta^2 \alpha_{s+2} + \sum_{t=3}^{j} \sum_{s=3}^{t} \Delta^2 \beta_{s} + \sum_{t=3}^{A-i+j} \sum_{s=3}^t \Delta^2 \gamma_s$$ [落ち着け、そんなに難しいことはいってないはずだ。
第1式では、まだAPCのモデル化は始まっていない。第1-3項は切片で、年齢と時代を床にした空間に張った平面になっている。
第2式はやたらにややこしそうなので、心の準備をしておくと、$$ \sum_{i=p}^{q} \Delta x_i = (x_p – x_{p-1}) + (x_{p+1} – x_{p}) + \cdots + (x_q – x_{q-1}) = x_q – x_{p-1}$$ だから、差分系列の足し上げは始点の直前から終点までの増分を表す。オーケー、一番簡単そうな第2項について考えてみよう。\(j=\{1,2\}\)ではゼロである。\(j=3\)以降では$$ \sum_{t=3}^{j} \sum_{s=3}^{t} \Delta^2 \beta_{s} = \sum_{t=3}^{j} (\Delta \beta_t – \Delta \beta_2) = (\beta_j – \beta_2) – (j-2) (\beta_2 – \beta_1)$$ だから、\(\beta_2\)からみた\(\beta_j\)の高さと「もし\(\beta_1\)と\(\beta_2\)を通る線が時点2から時点\(j\)までまっすぐ伸びていったら\(\beta_j\)はどこまで高くなっていたか」との差、つまり最初の式に書いた切片平面からみた\(\beta_j\)の偏差になっているわけだ。
えーと、まあ要するに、もともとの発想である \( \mu_{ij} = \alpha_i + \beta_j + \gamma_k \) を、切片平面の高さ\(\mu_{A1}\), 傾き\(s_\alpha, s_\beta\), ３つの効果の二階差分系列\(\Delta^2 \alpha_i, \Delta^2 \beta_j, \Delta^2 \gamma_k\)を使って書き直したってことなんでしょうね]
これを整理すると、結局\(\mu = M’ \xi\)と書ける。\(\xi\)は切片を表す3つのパラメータ\(\mu_{A1}, s_\alpha, s_\beta\)と、二階差分を表す\(A-2\)個と\(T-2\)個と\(A+T-3\)個のパラメータを縦に積んだもの。\(M\)は… [式は省略するけど] 計画行列である。
これは識別できる。Rでいうとapcパッケージがこの定式化をとっている。
[なるほど… apcパッケージのドキュメントを読んでもぴんとこなかったし、これを読んでもやっぱり狐につままれたような気分だが、まあとにかく二階差分をパラメータとすると識別できるってわけね。そういってしまえば凄そうだが、要は識別できなさそうな部分を\(s_\alpha, s_\beta\)にギュッと押し込めただけであって、この傾きを定義する3点 \(\mu_{A1}, \mu_{A-1, 1}, \mu_{A2}\) の選び方に依存しちゃうんじゃない？ これって、床の上のものをすべて押し入れにぶち込んで部屋を片づけたと称するようなもんなんじゃないでしょうか]
[…と不思議に思っていた話が、読み飛ばしていた次の段落で説明されていたことに気が付いた。難しくなってくるのでほぼ逐語訳する]
仮に \(\mu_{A1} – \mu_{A-1,1}\)を\(\Delta \alpha_A – \Delta \gamma_2\)に、\(\mu_{A2} – \mu_{A1}\)を\(\Delta \beta_2 + \Delta \gamma_2\)に書き換えれば、このパラメータ化は、Holford(1983)が提案した線形効果と曲線効果というモデルに類似したモデルとなる。Holfordのモデルでは、すべての曲線効果と、\(u s_\alpha + v s_\beta + (v-u) s_\gamma\)という形式をとる傾きの線形結合を推定できる。たとえば\( (u,v) = (1,0) \)とすれば\(s_\alpha – s_\gamma\)が識別可能であり、\( (u,v) = (0,1) \)とすれば\( s_\beta + s_\gamma \)が識別可能である。Holfordのモデルにおいては線形トレンドの推定可能な関数は無限個あることになる。
同様に、Nielsenのパラメータ化では、３つの初期点をどう選ぶか(つまり最初の差分の関数をどう選ぶか)は一意に決まらない。3点のインデクスが線形従属でない限り、\(\mu_{A1}, \mu_{A-1, 1}, \mu_{A2}\)の代わりにどの3点を選んでもよい。Nielsen & Nielsen (2014) は初期点を、年齢とコホートの中央値、そして外れ値でないこと、という基準で選んでいる。このガイドラインを用いれば、ベースライン・レートは端ではなくて中央の率を反映することになる。
これをベイジアンの枠組みに持っていくこともできる。初期点としてなんらか3点選ぶとして、$$ \pi(\mu_{i_1, j_1}, \mu_{i_2, j_2}, \mu_{i_3, j_3}) \propto 1$$ $$ \Delta^2 \alpha_3, \ldots, \Delta^2 \alpha_A | \tau^2_\alpha \sim N(0, \tau_\alpha^{-2} I) $$ $$ \Delta^2 \beta_3, \ldots, \Delta^2 \beta_T | \tau^2_\beta \sim N(0, \tau_\beta^{-2} I) $$ $$ \Delta^2 \gamma_3, \ldots, \Delta^2 \gamma_{A+T-1} | \tau^2_\gamma \sim N(0, \tau_\gamma^{-2} I) $$ という風に。
6.2 ベイジアンIntrinsic自己回帰モデル
3つの効果のそれぞれについて、一次ないし二次のランダムウォーク事前分布を考える。二次RWの場合、二階差分がiidな正規確率変数だと捉えているわけで、これは、効果の平均が前の2点と後の2点における効果で条件づけられているようなガウシアン確率場(GMRF)を意味している[なるほどおおお？？？]
つまりだ。$$ \alpha | \tau^2_\alpha \sim RW2(\tau^2_\alpha) $$ ということは $$ \Delta^2 \alpha_i | \tau^2_\alpha \sim N(0, \tau_\alpha^{-2} I) $$ だから、\(i=3, \ldots, A-2\)について$$ \alpha_i | \alpha_{-i}, \tau^2_\alpha \sim N \left( \frac{2}{3} (\alpha_{i-1} + \alpha_{i+1}) – \frac{1}{6} (\alpha_{i-2} – \alpha_{i+2}), \frac{1}{6 \tau^2_\alpha} \right) $$ となる。[ひゃあああ… ほんとにGMRFだ… まじか…]
[ここから理解できなくなるのでほぼ逐語訳]
RW2事前分布は正則なMVNにならない。精度行列のランクは\(n-2\)である。こういうクラスの事前分布をintrinsic GMRFという。RW2モデルに切片と傾きを含める場合、非正則性はふつう、示唆された精度行列のnull spaceに張った線形制約の集合によって条件付けることで、その空間での非正則な密度関数がより低次元な空間での正則な密度関数と等価になるようにすることによって扱われる[原文: the impropriety is usually addressed by conditioning on a set of linear constraints spanning the null space of the implied precision matrix so that the improper density is equivalent to a proper density on a lower dimensional space]。特に、RW2モデルでは合計ゼロ制約と傾きゼロ制約がある。たとえば年齢の効果では、\(i\)番目の列が\(\{1,i\}, i=1,\ldots,A\)である行列\(L\)について\(L \alpha = 0\)であることと等価である。APCの文脈では、これらの制約はover-parameterizedでないモデルを与えるが、切片や傾きは解釈できない。データからは切片と傾きがわからないからである。
[予測という観点からはRW1よりRW2のほうが良いという話。略]
RW2事前分布を使ったベイジアンAPCモデルは、ふつう以下のように指定される。$$ y_{ij} | \mu_{ij} \sim Poi(N_{ij} \exp(\mu_{ij})) $$ $$ \mu_{ij} = \delta + \alpha_i + \beta_j + \gamma_k + z_{ij}$$ $$ \delta \sim M(m_\delta, s^2_\delta)$$ $$ \alpha | \delta^2_\alpha \sim RW2(\tau^2_\alpha)$$ $$ \beta | \tau^2_\beta \sim RW2(\tau^2_\beta)$$ $$ \gamma|\tau^2_\gamma \sim RW2(\tau^2_\gamma)$$ $$ z|\tau^2_z \sim N(0, \gamma^{-2}z I)$$ $$ \tau_\theta \sim \pi_\phi, \, \phi = \alpha, \beta, \gamma, z$$ \(z_{ij}\)は構造化されていないランダム効果で、これを付け加えることで制約された時間的効果の周りの異質性が許容されるし、MCMCによる計算も簡単になる(Gibbsサンプリングが使えるようになるから)。３つの効果に合計ゼロ制約はかけるが線形トレンドゼロ制約はかけないことが多い(線形の項はないんだから筋がとおっている)。以下ではこれを古典的RW2モデルと呼ぶ。
いっぽう、線形の傾きも追加するやり方もある。年齢と時代に傾きを追加すれば$$ \mu_{ij} = \delta + v_{1i} + v_{2j} + \alpha_i + \beta_j + \gamma_k + z_{ij}$$ これを新型RWモデルと呼ぶことにする。
RW2のAPCモデルは、計算が楽だという点、そして年齢と時代の間隔が等しくなくてもよいという点が魅力的である。Nielsenのパラメータ化の場合、間隔が等しくないときにどうするんだという課題が残されている。
6.3 分散成分による事前分布の指定
[RW2のモデルで\(\tau\)の事前分布をどうするかという話だと思う。力尽きたのでスキップ]
6.4 実装
RW1ないしRW2のベイジアンAPCモデルはBAMPパッケージ, WinBUGS, RのBAPCパッケージ(これはINLAを使っている)で推定できる。Nielsenのパラメータ化のベイズ推定はINLAかMCMCでできるが、INLAではちょっとややこしい。
[…後略]
7. モデルの比較
[データ例について、Nielsenのパラメータ化(MLE, INLA, MCMC)とRW2のモデル(古典的、新型。どちらもINLA)の推定結果を比較している。補足資料にコードがあるのかな。疲れたのでスキップ]
8. 考察
識別可能なパラメータ化に基づいたモデルを使え、ないし、変換に対して不変なパラメータ化から得られた要約のみに注目せよ。
本論文では、RW2に基づくモデルも、Nielsenのパラメータ化によるベイズ推定も、結果はだいたい一致することを示した。
本論文は識別可能性の問題に焦点をあてたが、他にも研究はたくさんある。個人レベルのデータが利用可能な場合とか。また、対数リンクじゃなくてpowerリンクを使うという話題もある。
—–
途中で力尽きてしまったが、apiパッケージがなにをやっておるのかがわかったので、これで良しとしよう。(自分に甘い)
自分の問題に当てはめると、自力でMCMCするとかR-INLAを使うとか、正直なところ勘弁してほしい。でもBAMPパッケージってRW1でさえめったに収束しないと思うんですよね。いっぽうapcパッケージは使いやすいんだけど、調査は年次で年代は5歳刻みという場合に困る。うううううむ。