特集「がん統計データおよびその解析」［研究ノート］

統計数理（2011）
第 59 巻 第 2 号 205–215
c 2011 統計数理研究所


変化係数を用いたがん死亡危険度の年次変動要因の推測
1

1

2

3

冨田 哲治 ・佐藤 健一 ・中山 晃志 ・片野田 耕太 ・3

祖父江 友孝 ・大瀧 慈

1（受付 2010 年 12 月 24 日；改訂 2011 年 3 月 15 日；採択 4 月 28 日）

要旨

がんなどの疾患の罹患・死亡リスクは，様々な背景要因の変化に伴って経年的に変動している．このような背景要因と経年変動の関係を解析することは，我が国のがんリスクの動向を把握し，今後のがん対策へ繋げていくために基礎となる重要なテーマである．がん死亡数は地域ごとに毎年集計されており，このような繰り返し観測されたデータは特に経時データと呼ばれている．本論文では，観測値駆動型モデルを用いて経時データの特性を考慮しつつ変化係数により説明変数の経時的な変化を推測するための解析方法を提案する．実データへの適用例として，日本における男性の大腸がん死亡数データに提案手法を適用し解析を行った．
キーワード： 観測値駆動型モデル，一般化線形 ARMA，変化係数，一様性検定．
1.

はじめに

がんなどの疾患の罹患・死亡リスクは，様々な背景要因の変化に伴って経年的に変動している．このような背景要因と経年変動の関係を解析することは，我が国のがんリスクの動向を把握し，今後のがん対策へ繋げていくために基礎となる重要なテーマである．がん死亡数は地域ごとに毎年集計され，その一部はウェブ上で一般公開もされている．このような，繰り返し観測されるデータは特に経時データと呼ばれており，経時データを扱う場合には観測時点間の相関を考慮することが重要である．また，経時データに対する回帰モデルにおいて，変化係数は時間とともに変化する背景要因の影響を評価でき，がん死亡危険度の年次変動に関わる背景要因の影響を定量的に評価し，視覚的に理解するのに有用である．
がん死亡数といった離散値の経時データの解析には，一般化線形モデル（Dobson, 1990）にランダム変数を導入することで時点間相関をモデル化した一般化線形混合モデル（例えば, Breslowがよく使われている．しかし，尤度関数が陽に表せないため未知母数推定
and Clayton, 1993）に必要な計算量が膨大となり，観測時点や説明変数の数が比較的多い場合には適用が困難である（Chan amd Ledolter, 1995; Durbin and Koopman, 2000; Jung and Liesenfeld, 2001）．観測値
駆動型モデル（例えば, Zeger and Qaqish, 1988）は，ランダム変数の代わりに過去の観測値を用いたモデルであり，推移モデルとも呼ばれている（Diggle et al., 1994）．観測値駆動型モデルは，
123広島大学 原爆放射線医科学研究所：〒734–8551 広島県広島市南区霞 1–2–3
国際医療福祉大学 福岡看護学部：〒810–0072 福岡県福岡市中央区長浜 1–3–1
国立がんセンター がん対策情報センター：〒104–0045 東京都中央区築地 5–1–1

統計数理

206

第 59 巻

第 2 号 2011

尤度関数が陽に得られ扱いやすいという利点があり，特に，Davis et al.（2003）
は，自己回帰移動平均（AutoRegressive Moving Average, ARMA）に基づく観測値駆動型モデルを提案し，推定量の性質などについて報告している．
変化係数は，Hastie and Tibshirani（1993）
によって提案され，時間とともに変化する説明変数の効果を視覚的に評価できるため，経時測定データの解析に極めて有用である．変化係数の推定法は，固定した時間の近傍データを用いて回帰を繰り返すカーネル平滑化法（例えば, Hoover
et al., 1998; Tonda et al., 2011）
が一般的だが，この方法では変化係数の同時信頼区間の構築や検定が困難であった．そこで，Satoh and Yanagihara（2010）は変化係数の形状に線形構造を導入することで，連続値の経時データに対する変化係数の推測法を提案した．また，佐藤 他（2009）
は Liang and Zeger（1986）による一般化推定方程式を用いた離散値の経時データにおける変化
係数の推測法を提案し，冨田 他（2010）は空間データに対して位置とともに変化する回帰係数（変化係数曲面）の推測法を提案した．
本稿では，観測値駆動型モデルに変化係数を導入し，時点間相関を考慮しつつ，がん死亡数
危険度の年次変動要因を変化係数として視覚的に評価するための解析法の開発を試みる．変化係数を用いて，がんの危険・予防要因となる因子の効果の経時的な変化を捉えることは，がん対策の立案において有用であると考えられる．2 節で，観測値駆動型モデルと変化係数を組み合わせることでモデルを記述し， 3 節および 4 節で推測法を提案する． 5 節において，日本における男性の大腸がん死亡数データを用いた解析例を紹介する．
2.

統計モデル

地域 i の t 年における観測死亡数を dit とする（i = 1, . . . , N ; t = 1, . . . , T ）．過去の死亡数の観測履歴 Fi,t−1 = {dis |s ≤ t − 1} が与えられたとき，dit の条件付き分布は平均 µit のポアソン分布に従うとする，(2.1)

dit |Fi,t−1 ∼ Poisson(µit ) .

また，平均母数 µit は，k 個の説明変数 ai = (ai1 , . . . , aik ) を用いて，(2.2)µit = eWit ,Wit = ai β(t) + Zit ,で記述されているとする．ここで，β(t) は時間とともに変化する説明変数の効果を表す回帰係数であり，変化係数（Hastie and Tibshirani, 1993）と呼ばれている．一方，Zit は観測値の観測時点間の相関をモデル化するための変数である．本研究では，尤度関数が陽に得られ扱いやすいことから，観測値駆動型モデルを用いて，過去の観測値の関数として Zit を記述する．Zit の構造の違いによって様々な観測値駆動型モデルが提案されており，Davis et al.（1999）
で詳しくレビューされている．ここでは，Zit は過去全ての観測値の加重和で表されるものとする．つ
まり，
(2.3)

Zit =

∞


π eit− ,

eit =

=1

dit − µit
.
√
µit

ここで，式
（2.3）
は無限次元の移動平均（Moving Average, MA）と呼ばれ，加重和の重みである π
は未知母数である．式（2.3）
における無限個の未知母数 π は自己回帰移動平均（AutoRegressiveを用いて有限個の未知母数で表現可能であることが知られている．こMoving Average, ARMA）
こでは，ARMA(p, q) フィルター：π(B) =∞

=1π B  =

θ(B)
− 1,
φ(B)

変化係数を用いたがん死亡危険度の年次変動要因の推測

207を用いて加重和の重み π を有限個の未知母数で表すことにする．ただし，B は Bdit = dit−1 で定義されるラグ作用素であり，pqφ(B) = 1 −φ B  , θ(B) = 1 +θ B  ,=1

=1

は ARMA(p, q) における自己回帰と移動平均を規定するラグ作用素の多項式である．これより，
Zit は p + q 個の未知母数 φ = (φ1 , . . . , φp ) ，θ = (θ1 , . . . , θq ) を用いて，
(2.4)

Zit =

p


φ (Zit− + eit− ) +

=1

q


θ eit− ,=1と表現される．
3.

変化係数の推定における変化係数 β(t) の推定法について紹介する．推定法の概要は，β(t)本節では，式（2.2）の形状に線形構造を導入することで，式（2.2）
を Davis et al.（2003）の提案した一般化線形 ARMA
モデルの枠組みに帰着させ，尤度関数を陽に求めることにある．
まず，変化係数 β(t) の形状として β(t) = Ξx(t) で表現されるクラスを考える．ここで，Ξ =(ξ1 , . . . , ξ k ) は k × r の未知母数行列であり，x(t) は r 個の既知の基底関数からなるベクトルである．例えば，変化係数の形状として r − 1 次多項式を用いるなら，x(t) = (1, t, . . . , tr−1 ) とな
る．より滑らかな形状を仮定したい場合には，B スプライン基底やガウス基底などの基底を用いることも可能である．基底については，例えば，Satoh et al.（2003）, Ruppert et al.（2003）の3 章，小西・北川（2004）の 4 章などを参照されたい．これにより，変化係数 β(t) を推定することは，未知母数行列 Ξ を推定することと同値となる．今，xit = ai ⊗ x(t)，ξ = vec(Ξ ) と定義
すれば，過去の観測値の履歴：Fi,t−1 が与えられたときの dit の条件付き期待値は，µit = eWit ,Wit = xit ξ + Zit ,と書け，Davis et al.（2003）の提案した一般化線形 ARMA モデルで表される．モデルに含まれる未知母数は，変化係数の形状を規定する ξ = vec(Ξ )，時点間相関を記述する p 次の自己回帰に関する母数 φ = (φ1 , . . . , φp ) と q 次の移動平均に関する母数 θ = (θ1 , . . . , θq ) である．
Davis et al.（2003）の提案した尤度関数の構成法を参考に，最尤法により未知母数 δ = (ξ , θ , φ )
の推定を行う．式（2.1）
より，過去の観測値の履歴 Fi,t−1 が与えられたときの dit の分布はポアソ
 Tン分布に従い，密度関数を f (dit |Fi,t−1 ) とすれば，対数尤度関数は N
i=1
t=1 log f (dit |Fi,t−1 )で与えられる．このとき，未知母数を含まない項を除いた対数尤度関数 (δ) は次式で与えられる．
(δ) =TN yit Wit (δ) − eWit (δ) ,Wit (δ) = xit ξ + Zit (δ) .

i=1 t=1

Wit と Zit は δ に依存する関数であるが，表記を簡単にするため以下省略して説明する．対数尤度関数 (δ) の最大化には Newton-Raphson 法といった数値最適化を用いて行う．数値最適化に必要な (δ) の 1 階および 2 階の導関数は，Wit の偏微分，Zit と eit に関する漸化式で得られ，詳細な記述は付録にまとめた．Newton-Raphson 法の初期値は，ξ = ξ̂GLM ，θ = 0，φ = 0，{eit = 0|t ≤ 0} とした．ξ̂GLM は，観測時点間の相関を独立とした通常の一般化線形モデルにおける推定値である．Davis et al.（2003, 2005）の結果から，δ の最尤推定量 δ̂ は漸近的に，

 2
−1
∂ (δ) (3.1),
δ̂ ∼ Nrk+p+q (δ, Ω) , Ω̂ = −∂δ∂δ δ=δ̂

統計数理

208

第 59 巻

第 2 号 2011が成り立つ．これより，変化係数 β(t) = (β1 (t), . . . , βk (t)) の推定値およびその漸近分布は，次式で与えられる．
β̂j (t) = ξ̂j x(t) ∼ N(βj (t), λj (t)) ,(3.2)

j = 1, . . . , k ,ただし，λj (t) = x(t) Ωjj x(t)，Ωjj は ξ̂j の分散共分散行列に対応する Ω の部分行列である．
4.

信頼区間の構築と検定が提案した連続値の経時データに対する変化係数のここでは，Satoh and Yanagihara（2010）
推測法を参考に，変化係数 β(t) = (β1 (t), . . . , βk (t)) の各 βj (t) に関する信頼区間の構築と検定問題について扱う． ξ̂ の分散共分散行列の推定量 Ω̂ = {Ω̂ij } が得られているとき，変化係数の推定量 β̂j (t) の分散は Ωjj の推定量 Ω̂jj を用いて λ̂j (t) = x(t) Ω̂jj x(t) で与えられる．このとき，Rao（1973）の 1 章 1f 節にある Cauchy-Schwarz の不等式を用いて，(4.1){(ξ̂ j − ξ j ) x}2{β̂j (t) − βj (t)}2≤ (ξ̂j − ξ j ) Ω̂−1≤ supjj (ξ̂ j − ξ j ) ,t∈Rx∈Rrx Ω̂jj x
λ̂j (t)supが成り立つことが，Satoh and Yanagihara（2010）により示された．この不等式により，変化係
数 βj (t) の信頼領域の構成問題は，変化係数 βj (t) の形状を規定する未知母数 ξ j の問題に置き換えられ，時間 t の関数である変化係数 βj (t) が時間 t に依存しない上限で評価することが可能となる．また，ξ̂ j の漸近正規性により，右辺は自由度 r のカイ二乗分布 χ2r に近似的に従う．
したがって，cr,α を χ2r の上側 100α％点，すなわち Pr(χ2r ≥ cr,α ) = α とすると，変化係数 βj (t)の 100α％信頼領域は次式で与えられる．
Ij,α (t) = β̂j (t) − λ̂j (t)cr,α , β̂j (t) + λ̂j (t)cr,α .
(4.2)
このとき，式（4.2）の信頼区間の被覆確率は漸近的に Pr(βj (t) ∈ Ij,α (t) : ∀ t ∈ R) ≥ 1 − α が成り立つ．
次に，変化係数に関する検定問題を考える．まず，変化係数 βj (t) の値がすべての時間 t に
おいて 0 であるという帰無仮説は，(4.3)

H0 : βj (t) ≡ 0(∀ t ∈ R) ,と書ける．素朴な検定統計量として Tj (t) = β̂j (t)2 /λ̂j (t) を考えると，Pr(Tj (t) ≤ x : ∀ t ∈ R2 ) を評価する必要があるが，(4.4)

Pr(Tj (t) ≤ x : ∀ t ∈ R) = Pr max Tj (t) ≤ x ,t∈Rが成り立つことに注意すれば，信頼領域の構成と同様に maxt∈R Tj (t) ≤ Tj ，Tj = ξ̂j Ω̂−1jj ξ̂ j が成2
り立つことが分かる．Tj は式（4.3）の帰無仮説のもとで漸近的に χr に従うので，帰無仮説 H0は，Tj > cr,α ならば有意水準 100α％で棄却される，もしくは，Pr(χ2r > Tj ) から p 値を求めることが可能である．
また，変化係数 βj (t) の形状を構成する基底 x(t) の成分に 1 を含めれば，変化係数が時間 tに依らず定数であるという一様性の仮説(4.5)

H0 : βj (t) ≡ const.

(∀ t ∈ R)の検定も可能である．今，基底を x(t) = (1, x(g) (t) ) ，x(g) (t) = (x2 (t), . . . , xr (t)) とすれば，
(g)
(g)
(g)
(g)
β j (t) = ξj1 + β j (t) と書ける．ただし，β j (t) = ξj x(g) (t)，ξj = (ξj2 , . . . , ξjr ) ．これより，

変化係数を用いたがん死亡危険度の年次変動要因の推測

209

変化係数に関する一様性の帰無仮説（4.5）は，(g)

H0 : βj (t) ≡ 0(4.6)(∀ t ∈ R) ,
(g)

(g)

で表される．式（4.3）の帰無仮説の検定と同様に，ξ (g)の分散共分散行列のj の推定量を ξ̂ j ，ξ̂ j
(g)
(g)
(g)
(g)
(g)
(g)
推定量を Ω̂jj ，λ̂j (t) = x(g) (t) Ω̂(g)
x
(t)
とすれば，
素朴な検定統計量T(t)= β̂j (t)2 /λ̂j (t)
jj
j
(g)

−1

(g)

(g)
に対して，maxt∈R Tj(g) (t) ≤ Tj(g) ，Tj(g) = ξ̂j Ω̂(g)
（4.6）
の帰無仮説
ξ̂ j が成り立つ．Tj が式
jj(g)
2のもとで漸近的に χr−1 に従うので，帰無仮説 H0 は，Tj > cr−1,α ならば有意水準 100α％で棄却される，もしくは，Pr(χ2r−1 > Tj(g) ) から p 値を求めることが可能である．

5.

実データへの適用例

実データへの適用例として，日本における男性の大腸がん死亡数データに提案手法を適用し解析を行った．データは 1975 年から 2002 年までの期間に都道府県別に毎年集計された 28 時
点の経時データである．期待死亡数算出のための基準集団には，1975 年から 1994 年の期間における全国平均を用いた．説明変数には，人口密度 [persons/km2 ]，日照量 [KWh/m2 /day] お
よび自然空間放射線 [nGy/h] を用いた，ただし全ての説明変数に関して各平均値でシフトしてある．2 節の表記に従って書けば，ai = (ai1 , ai3 , ai2 , ai4 ) ，ai1 = 1，ai2 = 人口密度，ai3 = 日照量，ai4 = 自然空間放射線，i = 1, . . . , 47 となる．人口密度は都会度を表す指標として用いた．日照量に関して，Mizoue（2004）が横断的研究により消化器系がん死亡危険度との負の相関関係を報告しており，同様な日照量データ（1961 年から 1990 年の間の平年値）を用いた．自然空間
放射線は，Abe et al.（1981）の調査に基づく放射線医学総合研究所の公開データを用いた．
変化係数 β(t) の形状を規定する基底関数として r − 1 次の多項式を用い，AIC を用いて最適な次数 (r, p, q) の選択を行った．表 1 は，r = 4（つまり，3 次多項式を適合）の場合に，(p, q) をそれぞれ 0 から 4 まで動かしたときの AIC の値である．
AIC が最小となるのは (r, p, q) = (4, 1, 1) であった．他の r の値についても検討したが， AIC
が最小となるのは (r, p, q) = (4, 1, 1) であった．図 1 は，(r, p, q) = (4, 1, 1) の場合の適合曲線である．黒色の太線が適合曲線であり，灰色の細線は観測値をつないだ折れ線である．
図 2（a）から図 2（d）
は，(r, p, q) = (4, 1, 1) の場合の変化係数 β(t) の推定曲線である．図 2（a）
が全ての説明変数の値が各平均値の場合の死亡危険度の時間変化の推定曲線，図 2（b）が人口密度の効果，図 2（c）
が日照量の効果，図 2（d）
が自然空間放射線の効果を表す変化係数の推定曲線に対応する．上下の破線は，各推定曲線の 95％信頼区間である．
変化係数の有意性の検定および一様性の検定の結果を表 2 に示す．表 2 の各セルの数値は検定統計量の値であり，括弧の数値は対応する p 値である．
図 2（b）から，人口密度の効果は常に正値であることから，地方（人口密度の低い地域）に比
べて都会（人口密度の高い地域）の死亡リスクが経年的に高いことが分かる．また，地方と都会
表 1.

AIC による (p, q) の選択（r = 4 の場合）．

統計数理

210

第 59 巻

第 2 号 2011

図 1. (r, p, q) = (4, 1, 1) における適合曲線．

図 2.

(r, p, q) = (4, 1, 1) における変化係数 β(t) の推定曲線（実線）と 95％信頼区間（破線）．

変化係数を用いたがん死亡危険度の年次変動要因の推測
表 2.

211

(r, p, q) = (4, 1, 1) における変化係数に関する検定．

の差は経年的に見て減少傾向であるが，表 2 の一様性の検定結果から，5％有意とはならなかった．図 2（c）から，日照量の効果は常に負値であることから，日照量の低い地域は高い地域に比べて死亡リスクが経年的に高く，この結果は Mizoue（2004）の報告と一致する．また，表 2の一様性の検定結果から，日照量の効果は経年的に定数であることが示唆された．図 2（d）から，自然空間放射線の効果は常に 0 近くにあること，表 2 の検定結果から，自然空間放射線の死亡危険度への寄与はみられなかった．
以上の解析結果をまとめると，都会および日照量の少ない地域での死亡危険度が高い傾向があることが示唆された．また，このような傾向は経年的に変化していないと考えられ，その結果，死亡危険度の適合曲線はほぼ平行となった．日照量のがん死亡危険度に関する効果は，
日光浴により生成されるビタミン D が関与していると考えられている（Mizoue, 2004; Garland．
and Garland, 2006; Boscoe and Schymura, 2006）
6.

おわりに

本論文では，都道府県別がん死亡数の経時測定データに対して，観測値駆動型モデルを導入することで観測時点間相関を考慮しつつ，変化係数を用いてがん死亡危険度の変動要因の推測を行った．変化係数の推定は，カーネル平滑化法（Hoover et al., 1998 など）を用いるのが一般的であるが，本論文では，Satoh and Yanagihara（2010）の推測法を参考に，変化係数の形状に線形構造を仮定した推測法を適用した．線形構造を仮定することで，変化係数に関する同時信
頼区間の構築や時間変化の有無に関する仮説検定が可能となった．経時データにおける観測値の時点間相関を考慮したモデルとして，ランダム変量を用いた一般化線形混合モデル（GLMM）に代表されるランダム効果モデル（例えば, Breslow and Clayton, 1993）と，一般化推定方程式（GEE）に代表される周辺モデル（例えば, Liang and Zeger, 1986）がよく使われている．GLMM
では，導入するランダム変量の数だけ陽に解けない多重積分を数値計算で評価する必要があり，観測時点や説明変数の数が比較的多いデータに対する適用は困難となる．GEE は，計算量が少なく適用は容易であるものの，観測時点が多い場合には分散共分散に関する未知パラメータ数が増大し推定が不安になり易く，また，尤度が得られないことから，AIC といったモデル選択
規準も利用できない．本論文では，ランダム変量の代わりに過去の観測値を用いた観測値駆動型モデルを適用した．観測値駆動型モデルを用いる利点は，少数の未知パラメータで時点間相関のモデル化が可能であり，尤度関数が陽に表現可能であることから計算量が少なく，実デー
タへの適用が容易であることにある（Davis et al., 1999, 2003, 2005）．観測値駆動型モデルの導入により，より効率よく変化係数の推測が可能となった．また，尤度関数が陽に得られ容易に計算可能であることから，AIC といった従来の変数選択の手法を用いることができた（表 1）．
次に，提案手法を日本における男性の大腸がん死亡数データに適用した結果について考察を行う．本論文では，がん死亡危険度の変動要因として人口密度，日照量，自然空間放射線について検討を行った．変化係数に関する一様性の検定結果は 5％有意ではないものの（表 2），人口密度の効果に経年的な減少傾向が見られた（図 2（b））．これは，近年になって都会と田舎の

統計数理

212

第 59 巻

第 2 号 2011

地域差が小さくなっていると考えられる．日照量のがん死亡危険度への影響に関する報告は，
Mizoue（2004）, Garland and Garland（2006）
および Boscoe and Schymura（2006）らにより報告されているが，これらの報告は横断データの解析報告である．本論文では，経時データに対して変化係数を推定し，効果の時間変化の有無についても検討を行った．経時データ解析により，日照量の効果に時間変化が認められず経年的にほぼ一定であることが分かった（表 2，図 2（c））．
最後に，今回扱ったがん統計データは，都道府県という比較的大きな構成単位の集計データであったが，より細かな地域差に関心がある場合には，市区町村や第 2 次医療圏などを構成単位としたがん統計データの解析が必要となる．このような解析を行う際には，時点間相関に加えて空間相関を適切にモデルに導入することが重要となり，空間相関にも対応した解析法の開発は今後の課題である．

付録：(δ) の導関数に関する漸化式
対数尤度関数 (δ) の導関数は次式で与えられる：
T
N 


 ∂Wit∂
yit − eWit
=
,
∂δ
∂δ
i=1 t=1


T
N 

 2

∂2
Wit ∂ Wit
Wit ∂Wit ∂Wit
=
−
e
−
e
y
,
it
∂δ∂δ 
∂δ∂δ ∂δ ∂δ
i=1 t=1

ただし，



xit
∂Wit
∂Zit
=
,
+
∂δ
∂δ
0p+q

∂ 2 Wit∂Zit=.
∂δ∂δ 
∂δ∂δ 

Zit の 1 階微分は以下のような漸化式から得られる：
 

p
q
∂Zit−
∂Zit  ∂φ∂eit−
=(Zit− + eit− ) + φ
+
+
∂δ
∂δ
∂δ
∂δ
=1

=1

∂θ
∂eit−
eit− + θ∂δ∂δ

各成分は，ある添字 a に対して，

 
p
q

∂eit−
∂Zit∂eit−
∂Zit−
+
=
φ
+
θ
,
∂ξa
∂ξa
∂ξa
∂ξa
=1
=1

 
p
q

∂Zit−
∂Zit∂eit−
∂eit−
+= (Zit−a + eit−a ) +
φ
+
θ
,
∂φa
∂φa
∂φa
∂φa
=1
=1


p
q


∂Zit−
∂Zit∂eit−
∂eit−=φ+θ.
+ eit−a +
∂θa
∂θa
∂θa
∂θa
=1
=1
Zit の 2 階微分も以下のような漸化式から得られる：
 


p

∂eit− ∂φ∂eit−
∂ 2 Zit∂φ ∂Zit−
∂Zit−
+
=
+
+
∂δ∂δ
∂δ
∂δ 
∂δ 
∂δ
∂δ
∂δ 
=1
 2

∂ Zit−∂ 2 eit−+φ+∂δ∂δ 
∂δ∂δ q∂θ ∂eit−
∂eit− ∂θ∂ 2 eit−
+
+ θ
+
.
∂δ ∂δ∂δ ∂δ
∂δ∂δ 
=1

.

変化係数を用いたがん死亡危険度の年次変動要因の推測

213

各成分は，ある添字 a，b に対して，
 2
 
p
q

∂ Zit−∂ 2 eit−∂ 2 Zit∂ 2 eit−
+
=
φ
+
θ
,
∂ξa ∂ξb∂ξa ∂ξb∂ξa ∂ξb∂ξa ∂ξb
=1
=1
 2
 
p
q
∂ Zit−
∂ 2 Zit∂Zit−b∂eit−b ∂ 2 eit−∂ 2 eit−
=
+
+
φ
+
θ
,
+
∂ξa ∂φb∂ξa∂ξa∂ξa ∂φb∂ξa ∂φb∂ξa ∂φb
=1
=1
 2

p
q

∂ 2 Zit∂ 2 eit−∂ Zit−∂ 2 eit−b  ∂ 2 eit−
=
φ
+
+
θ
,
+
∂ξa ∂θb∂ξa ∂θb∂ξa ∂θb∂ξa∂ξa ∂θb
=1
=1∂Zit−a∂eit−a∂Zit−b∂eit−b
∂ 2 Zit=+++∂φa ∂φb
∂φb
∂φb
∂φa
∂φa
 2
 
p
q

∂ Zit−∂ 2 eit−∂ 2 eit−
+
+
φ
+
θ
,
∂φa ∂φb∂φa ∂φb∂φa ∂φb
=1
=1
 2

p
∂ 2 Zit∂Zit−a∂eit−a ∂ 2 eit−∂ Zit−=++φ+∂φa ∂θb∂θb∂θb∂φa ∂θb∂φa ∂θb=1q∂eit−b  ∂ 2 eit−
+
θ
,
∂φa
∂φa ∂θb
=1
 2

p
q

∂ 2 eit−
∂eit−b  ∂ 2 eit−∂ Zit−∂eit−a∂ 2 Zit
=
φ
+
+
+
θ
.
+
∂θa ∂θb∂θa ∂θb∂θa ∂θb∂θb∂θa∂θa ∂θb

+

=1

=1

eit の 1 階微分と 2 階微分は次式で与えられる：


1
∂Wit
1
∂eit
= − e 2 Wit + eit
,
∂δ
2
∂δ


1
1
1 ∂Wit ∂Wit
∂ 2 Wit∂ 2 eit= − e 2 Wit + eit
+ eit
.

∂δ∂δ
2
∂δ∂δ 
4∂δ ∂δ 謝辞本論文作成にあたって，ニューサウスウェールズ大学の William Dunsmuir 教授には貴重なコメントを頂いた．記して深く感謝申し上げたい．また，本誌査読者の方からは論文改訂にあたって有益なコメントを頂いた．ここに感謝いたします．本研究は，厚生労働省科研費「がん罹患・死亡動向の実態把握の研究」からの援助を受けた．

参 考 文 献Abe, S., Fujitaka, K., Abe, M. and Fujimoto, K.（1981）. Extensive ﬁeld survey of natural radiationin Japan, Journal of Nuclear Science and Technology, 18, 21–45.
. Solar ultraviolet-B exposure and cancer incidence and
Boscoe, F. P. and Schymura, M. J.（2006）
mortality in the United States, 1993–2002, BMC Cancer, Art. No. 264.
Breslow, N. E. and Clayton, D. G.（1993）. Approximate inference in generalized linear mixed models,
Journal of the American Statistical Association, 88, 9–25.
. Monte Carlo EM estimation for time series models involving
Chan, K. S. and Ledolter, J.（1995）
counts, Journal of the American Statistical Association, 90, 242–252.
. Modelling time series of count data, AsympDavis, R. A., Dunsmuir, W. T. M. and Wang, Y.（1999）

統計数理

214

第 59 巻

第 2 号 2011

totics, Nonparametrics and Time Series（ed. S. Ghosh），63–114, Marcel Dekker, New York.
Davis, R. A., Dunsmuir, W. T. M. and Streett, S. B.（2003）. Observation driven models for Poissoncounts, Biometrika, 90, 777–790.
. Maximum Likelihood estimation for an
Davis, R. A., Dunsmuir, W. T. M. and Streett, S. B.（2005）
observation driven model for poisson counts, Methodology and Computing in Applied Probability, 7, 149–159.
Diggle, P. J., Liang, K-Y. and Zeger, S. L.（1994）. Analysis of Longitudinal Data, Oxford University
Press, New York.
. Introduction to Generalized Linear Models, Chapman and Hall, London.
Dobson, A.（1990）
Durbin, J. and Koopman, S. J.（2000）. Time series analysis of non-Gaussian observation based on，Journalstate space models from both classical and Bayesian perspectives（with discussion）
of the Royal Statistical Society: Series B, 62, 3–59.
. Do sunlight and vitamin D reduce the likelihood of colon
Garland, C. F. and Garland, F. C.（2006）
cancer?, International Journal of Epidemiology, 35, 217–220.
Hastie, T. and Tibshirani, R.（1993）. Varying-coeﬃcient models, Journal of the Royal StatisticalSociety: Series B, 55, 757–796.
. Nonparametric smoothing estimates
Hoover, D. R., Rice, J. A., Wu, C. O. and Yang, L. P.（1998）of time-varying coeﬃcient models with longitudinal data, Biometrika, 85, 809–822.
. Estimating time series models for count data using eﬃcient
Jung, R. C. and Liesenfeld, R.（2001）
importance sampling, Allgemeines Statistisches Archiv, 85, 387–407.
小西貞則，北川源四郎（2004）. 『情報量規準』，朝倉書店，東京．
Liang, K. Y. and Zeger, S. L.（1986）. Longitudinal data analysis using generalized linear models,
Biometrika, 73, 13–22.
. Ecological study of solar radiation and cancer mortality in Japan, Health Physics,
Mizoue, T.（2004）
87, 532–538.
Rao, C. R.（1973）. Linear Statistical Inference and Its Applications, 2nd ed., John Wiley, New York.
. Semiparametric Regression, Cambridge University
Ruppert, D., Wand, M. and Carroll, R.（2003）
Press, New York.
. Estimation of varying coeﬃcients for a growth curve model,
Satoh, K. and Yanagihara, H.（2010）
American Journal of Mathematical and Management Sciences（in press）.
Satoh, K., Yanagihara, H. and Ohtaki, M.（2003）. Bridging the gap between B-spline and polynomial
regression model, Communications in Statistics Theory and Methods, 32, 179–190.
佐藤健一，柳原宏和，加茂憲一（2009）. 離散分布の経時測定データにおける線形な変化係数の推測について，応用統計学, 38, 1–11.
冨田哲治，佐藤健一，柳原宏和（2010）. 空間データに対する交互作用モデルを用いた変化係数曲面の推測について，応用統計学, 39, 59–70.
. A nonparametTonda, T., Satoh, K., Nakayama, T., Katanoda, K., Sobue, T. and Ohtaki, M.（2011）
ric mixed-eﬀects model for cancer mortality, Australian and New Zealand Journal of Statistics（in press）.
. Marcov regression models for time series: A quasi-likelihood
Zeger, S. L. and Qaqish, B.（1988）
approach, Biometrics, 44, 1019–1031.

Proceedings of the Institute of Statistical Mathematics Vol. 59, No. 2, 205–215 (2011)

215

Statistical Analysis of Time Trend in Cancer Mortality
Using Time Varying Coeﬃcient
Tetsuji Tonda1 , Kenichi Satoh1 , Teruyuki Nakayama2, Kota Katanoda3,
Tomotaka Sobue3 and Megu Ohtaki1
1 Research Institute for Radiation Biology and Medicine, Hiroshima University
2 Faculty of Nursing, International University of Health and Welfare
3 Research Center for Cancer Prevention and Screening, National Cancer Center

Disease risk, which is usually quantiﬁed by incidence or mortality, may have patternsof change by region or time due to varying factors such as level of exposure to environmental hazards or lifestyle habits of residents. Our aim is to assess how regional variationsof cancer mortality change over time and what characteristics inﬂuence such changes using longitudinal cancer mortality data. In this paper, we develop a statistical method for
longitudinal cancer mortality data using an observation driven model with time-varying
coeﬃcient. As an illustration based on real data, we apply the proposed method to a setof prefecture-speciﬁc longitudinal data on mortality from males’ large bowel cancer inJapan.

Key words: Observation driven model, generalized linear ARMA, varying coeﬃcient, test of uniformity.