時間の長さを分析する

粕谷 英一（九大・理・生物）生態学会（2015年3月）自由集会『データ解析で出会う統計的問題: 時間の長さ・時系列データの統計モデリング』

時間
状態の持続時間
回帰的枠組みなら
目的変数

時間
状態の持続時間 duration
ある現象が起こるまでの経過時間
latency

瞬間的な現象について 時間間隔
interval

時間
状態の持続時間 duration
例.餌を食べ始めてから食べ終わるまで
例.交尾開始から交尾終了まで
例.静止してから動き始めるまで状態→生存

例.寿命

例.余命

時間
ある現象が起こるまでの経過時間
latency

例.病気の潜伏期間
例.刺激を与えてから反応が現れるまでの時間

時間
時間間隔

interval

ある現象とある現象の間の時間間隔
例.捕食者が攻撃してくる間隔
例.ある餌場を訪れる間隔
例.バス停に人が来る間隔

時間

生存時間

文字通りの生存時間でないものも機械やシステムの故障までの時間など

状態の持続時間

duration

ある現象が起こるまでの経過時間
latency

時間間隔interval

時間 生存時間文字通りの生存時間でないものも機械やシステムの故障までの時間など
信頼性（reliability）

時間

データの特徴負の値がない正（正かゼロ）

連続量
尺度水準

比率尺度（比尺度）

重さ、長さ、体積などと同じ尺度水準

時間

データと分析の特徴

注目する特徴

時間

生存時間

生存率

ハザード

打ち切り
説明変数の値の時間的変化
よく使われる”変わった”モデル Cox比例ハザードモデル

時間以降の内容の目次

時間依存的な説明変数
打ち切り

Kaplan-Meier推定量
logrank検定

生存関数とハザード
比例ハザードモデル
パラメトリックなモデル

時間
広義の回帰の枠組み

説明変数

目的(応答)変数

回帰式の関数型
目的変数のばらつき具合
誤差構造と確率分布

時間の経過説明変数→目的（応答）変数
Cox比例ハザードモデルでは共変量

説明変数の値が途中で変化しない

定数型（fixed）
観察開始時の値がとられることが多い

説明変数の値が途中で変化する

Cox比例ハザードモデルでは普通

時間依存的（time-dependent）
時間のラグがある（時間により説明変数の係数が変化）

時間時間依存的な説明変数
パッケージsurvivalなど

比例ハザードモデルの関数coxph
説明変数指定のオプション
データ整形の関数

tt

tmerge

データ読み取りの関数

Surv

交互作用が使えることもある

打ち切り

censor

時間の長さが正確にわからない
生存時間
ある時点では生きていた
死亡時点は不明
しばしば「観察をやめた」ことによるある時点では生きていた
死亡時点は不明

打ち切り

censor

値はわからない

right

右側

left

左側

ある値以上であることはわかる
値はわからないある値以下であることはわかる

値はわからない
interval
ある範囲に入っていることはわかる

打ち切り

censor

値はわからないright

ある値以上であることはわかる

右側パッケージsurvivalの関数Surv

例.
2000頭
1年目生存、2年目生存
784
1年目生存、2年目死亡
16
1年目死亡
440
1年目生存、2年目打ち切り 760

1年目までの生存率
２年目までの生存率は？

例.
2000頭
1年目生存、2年目生存
784
1年目生存、2年目死亡
161年目死亡440
1年目生存、2年目打ち切り 760
1年目の生存率＝(784+16+770)/2000 =0.78
2年目の生存率＝784/(784+16) =0.98
2年目までの生存率＝1年目の生存率×2年目の生存率=0.7644

KaplanーMeier推定量(product limit推定量)
ある時点までの生存率の推定値＝区間の生存率の積
区間
死亡が起こったらその時点で区間を区切る

死亡が起こった時点で区切った区間の生存率

KaplanーMeier推定量
死亡が起こった時点で区切った区間の生存率
その積
例．
死亡あるいは打ち切り時点

1 3 4 5 5 6 7 7 78時点

個体数

死亡数

0357109740
122

区間の生存率
1.00.89
0.71
0.5

KaplanーMeier推定量
死亡が起こった時点で区切った区間の生存率

例．
時点
0357

区間の生存率
1.00.890.71
0.5

最初を基準にした生存率
1.0
0.89
0.630.32

パッケージsurvivalの関数survfit

logrank検定
２つのサンプル（処理が異なる間）での生存曲線の比較）

３つ以上のサンプル
拡張
ブロック（層別）あり

logrank検定
死亡が起こった時点で区切り、両サンプルの生存個体数、死亡個体数を計算
サンプル1の⽣生存個体数 　サンプル1の死亡個体数
サンプル2の⽣生存個体数 　サンプル２の死亡個体数

複数の2×2分割表
Mantel-Haentzel検定

logrank検定
区間の数だけ

2×2分割表サンプル1の⽣生存個体数 　サンプル1の死亡個体数
サンプル2の⽣生存個体数 　サンプル２の死亡個体数

Mantel-Haenszel検定

ブロック内の相関

共通の傾向（共通オッズ比）
死亡の順序（生存時間の順序）
パッケージsurvivalの関数survdiff

Wilcoxon順位和検定の拡張版(Mann-WhitneyのU検定)

２つのサンプル

生存日数の比較

時間データは順序尺度も使える確率分布

distribution

正（正かゼロ） 連続量指数分布
ワイブル分布
ガンマ分布

アーラン分布

カイ二乗分布
対数正規分布
ロジスティック分布
対数ロジスティック分布
ゴンペルツ分布
正規分布

確率分布
確率密度関数

ヒストグラムある値以下のデータの割合

分布関数
確率密度関数の積分ある値より大きいデータの割合
生存関数=(１−分布関数)

そのデータが生存時間なら生存率になる
生存曲線

正規分布(平均=0)の場合

生存関数

確率密度関数
normal distribution − pdf

0.6
0.20.4

1 − pnorm(x)

0.6
0.4

0.0

0.20.0dnorm(x)

0.8

0.8

1.0

1.0

normal distribution − survival−3−2−10x123−3−2−10x123ハザード hazard

(危険度)
✖️危険率
✖️検定での有意水準

時点yで生存していた個体の瞬間的死亡率

f ( y)h ( y) =S ( y)
ハザード関数

確率密度関数生存関数

ハザード hazard(危険度)時点yで生存していた個体の死亡率
生存関数

線の傾き

S ( y)S ( y)

y y+dydy→0
時点yで生存していた個体の瞬間的死亡率

ハザード hazard(危険度)

h ( y) =

生存関数

接線の傾き

S ( y)S ( y)y
時点yで生存していた個体の瞬間的死亡率
f ( y)h ( y) =S ( y)

1.0

normal distribution − pdf

0.6
0.4
0.00.2dnorm(x)0.8

確率密度関数生存関数−3−2−10123

ハザード関数x

normal distribution − hazard

0.210050.40.6

dnorm(x)/(1 − pnorm(x))

0.8151.0normal distribution − survival

0.01 − pnorm(x)正規分布(平均=0)の場合−3−2−10x123−3−2−10x123指数分布
確率密度関数f ( y) = λe−λy

ポアソン過程生存関数

S ( y) = e−λy

ハザード関数h ( y) = λ起こる確率が一定

指数分布
状態の終了が起こる確率が一定のときの、状態持続時間瞬間的な事象の起こる時間間隔
指数分布

生存率一定の生存曲線

指数分布

生存関数

確率密度関数exponential distribution − survival

0.60.20.4

1 − pexp(x, rate = 1)

0.6
0.4
0.2

0.00.0468100x246810x0.60.8

1.0exponential distribution − hazard

0.4

ハザード関数0.2

20.00

dexp(x, rate = 1)/(1 − pexp(x, rate = 1))

dexp(x, rate = 1)

0.8

0.8

1.0

1.0

exponential distribution − pdf

0246x810024x68

0.2

0.4

0.6

0.8

1.0

0.4

0.6

0.8

1.2

dweibull(x, shape = 1, scale = 1)/(1 − pweibull(x, shape = 1,scale = 1))
0.0
0.5
1.0
1.5
2.0

1.0

Weibull distribution − pdf

0.00.2

確率密度関数1 − pweibull(x, shape = 1, scale = 1)

0.0

dweibull(x, shape = 1, scale = 1)

ワイブル分布Weibull distribution − hazard

0240624x

ハザード関数パラメーターにより、ハザードが時間的変化

8x

Weibull distribution − survival

6生存関数

8ガンマ分布ハザード関数

確率密度関数

生存関数

00.2

0.4

0.6

0.8

1.00.00.00214

対数正規分布

2x638

lognormal distribution − pdf

確率密度関数

450.4

0.6

0.8

0.2

0.4

0.60.8dchisq(x, df = 1)/(1 − pchisq(x, df = 1))

0.2

dchisq(x, df = 1)

確率密度関数dlnorm(x, meanlog = 0, sdlog = 1)/(1 − plnorm(x, meanlog = 0,sdlog = 1))
0123
450.0

dlnorm(x, meanlog = 0, sdlog = 1)

1.0

1.0

カイ二乗分布chi−square distribution − pdf

chi−square distribution − hazard

10ハザード関数

00214x623x8410x

lognormal distribution − hazard

ハザード関数5

ロジスティック分布

logistic distribution − hazard
2.00.6logistic distribution − pdf

ハザード関数

1.5
1.0
0.5

0.3

0.0

0.0

0.1

0.2

dlogis(x)0.4

dlogis(x)/(1 − plogis(x))

0.5

確率密度関数−3−2−1012

3−3−2−10x123x対数ロジスティック分布
log−logistic distribution − pdf

0.2

0.4

0.60.8ハザード関数0.0dllog(x, shape = 1, scale = 1)

0.2

0.4

0.60.8

確率密度関数0.0dllog(x, shape = 1, scale = 1)

1.0

1.0

log−logistic distribution − pdf

0123x450

123
x45

0.5

1.0

1.50011223344Gompertz distribution − pdf

確率密度関数5ertz(x, shape = 1, rate = 1)/(1 − pgompertz(x, shape = 1,rate = 1))051015200.0dgompertz(x, shape = 1, rate = 1)

0.2

0.4

0.6

0.8

1.0

1.2

inverse Gaussian distribution − pdf

確率密度関数

5ゴンペルツ分布
x0.0

0.5

1.0

1.5

dinvgauss(x, mu = 1, lambda = 1)/(1 − pinvgauss(x, mu = 1, lambda = 1

0.0

dinvgauss(x, mu = 1, lambda = 1)

逆ガウス分布inverse Gaussian distribution − hazard

ハザード関数012
345x

Gompertz distribution − hazard

ハザード関数ある値より大きいデータの割合
生存関数

S ( y ) = 1− F ( y )目的変数＝生存関数
説明変数＝時間などという回帰的な分析 要注意
生存曲線など二重に評価
生存関数

データ点の独立性
時間時間比例ハザードモデル
proportional hazard model

Cox回帰

D.R.Cox(1924−)

比例ハザードモデル
http://www.stats.ox.ac.ukより(1972) J.R.Statist.Soc.B34:187-220.

パラメトリック・ブートストラップ(1961) Proc. 4th Berkeley Symposium, 105-123.
(1962) J. R. Statist. Soc. B 24, 406-424.
(2013) J. R. Statist. Soc. B75: 207-215.

D.R.Cox(1924−)

比例ハザードモデル
http://www.stats.ox.ac.ukより(1972) J.R.Statist.Soc.B34:187-220.

家畜などのウシ型結核菌感染とアナグマ
D. R. Cox et al. (2005) PNAS 102:17588–17593
J.R.Krebs

パラメトリック・ブートストラップ(1961) Proc. 4th Berkeley Symposium, 105-123.
(1962) J. R. Statist. Soc. B 24, 406-424.
(2013) J. R. Statist. Soc. B75: 207-215.

T.Clutton-Brock
R. Andersonなど

比例ハザードモデル
proportional hazard model

Cox回帰

中村『(1980年のこと)初めて聞く名前であった。…どうせ重回帰の変法だろうと想像して……だが届いた論文と原論文(Cox,1972)を読んでみても、まったく理解できず困ってしまった』比例ハザード性
説明変数

共変量(covariate)とも呼ぶ

説明変数が xの個体の時刻tにおけるハザードλ (t x ) = λ (t 0 ) ⋅ r ( x )相対危険度
ベースラインハザード
説明変数がすべて0の個体の、時刻tにおけるハザード

比例ハザード性
説明変数共変量(covariate)とも呼ぶ

説明変数が xの個体の時刻tにおけるハザードλ (t x ) = λ (t 0 ) ⋅ r ( x )対数線形性()r ( x ) = exp β x = exp ( β1 x1 + β 2 x2 +!)T比例ハザードモデル
説明変数が ( x1 , x2 ,!) の個体の時刻tにおけるハザード

λ ( t x1 , x2 ,!) = λ0 ( t ) ⋅ exp ( β1 x1 + β 2 x2 +!)λ0 ( t )ベースラインハザード

比例ハザードモデル
説明変数が ( x1 , x2 ,!) の個体の時刻tにおけるハザード

時間の要素

説明変数の効果λ ( t x1 , x2 ,!) = λ0 ( t ) ⋅ exp ( β1 x1 + β 2 x2 +!)λ0 ( t )ベースラインハザード

比例ハザードモデルの仮定
説明変数が異なる個体のハザード
比例ハザード性ハザードの比

時刻tにかかわらずいつでも同じ比例ハザード性＋対数線形生
ハザードの比

時刻tにかかわらずいつでも同じで、説明変数の一次式の指数関数

比例ハザードモデルの仮定
説明変数が異なる個体（個体たち）の生存関数
S ( t x1 , x2 ,!) = S0 ( t )exp( β1x1 + β 2 x2 +!)生存関数の対数の比は一定
log S ( t x1 , x2 ,!) = ( β1 x1 + β 2 x2 +!) ⋅ log S0 ( t )生存曲線は交わらない比例ハザードモデルの仮定が成り立っているとき

0.6
0.4
0.2
0.0

生
存率

0.8

1.0

proportional hazerd model− survival

0123x

時間

45

比例ハザードモデルの仮定が成り立っているとき

0.8

拡大図

0.0

0.2

0.4

0.6

0.8
0.6
0.4
0.2
0.0

生
存率

1.0

proportional hazerd model− survival

1.0

proportional hazerd model− survival

0123x時間

450.0

0.5

1.0x時間

1.5

比例ハザードモデルの仮定
説明変数が異なる個体（個体たち）の生存関数{}

log − log S ( t x1 , x2 ,!) = ( β1 x1 + β 2 x2 +!) + log {− log S0 ( t )}生存関数の対数の符号をかえたものの対数を時間に対してプロット形が同じで平行移動したもの

二重対数プロット34

proportional hazerd model− log−log

012log{-log(生存率)}

0123x時間45

比例ハザードモデル
説明変数が ( x1 , x2 ,!) の個体の時刻tにおけるハザード

λ ( t x1 , x2 ,!) = λ0 ( t ) ⋅ exp ( β1 x1 + β 2 x2 +!)こちらには手をつけずこちらだけ最尤法で推定（部分尤度、partial likelihood）

比例ハザードモデル
生存関数や生存曲線の計算

λ ( t x1 , x2 ,!) = λ0 ( t ) ⋅ exp ( β1 x1 + β 2 x2 +!)こちらを真の値だとして計算

適用例

比例ハザードモデル
再犯までの時間

432人

Rossi et al.(1980) Money, Work and Crime. Academic Press.
Allison(1995) Survival Analysis Using the SAS System.

目的変数
釈放後の逮捕までの時間（週単位）
説明変数
age 釈放時の年齢
prio それ以前の収監回数
fin

釈放後に経済的援助があったか（あり=1 　なし=0）データの一部
逮捕か打ち切り
経済的援助
時間

年齢

収監回数
就労(前) 仮釈放結婚

学歴week arrest fin age wexp mar paro prio educ個人
423 15
11 220 0
13 4424 5200 180 0
13 3425 19
10 18
0 00
2 3426 52
00 24
1 0
12 4427 12
11 22
1 112 4428 52
0
1 31
1 0
1
3 3429 52
00 20
0 0
1
1 4430 52
0
1 20
1 11
1 3431 52
00 29
1 0
1
3 4432 52
0
1 24
1 0
1
1 4

データの一部
目的変数

説明変数

逮捕か打ち切り

経済的援助

年齢

収監回数

時間

個人

week arrest fin age wexp mar paro prio educ423 151
1 220 01
3 4424 52
00 18
0 0
1
3 3425 19
10 18
0 0
0
2 3426 52
00 24
1 012 4427 1211 22
1 112 4428 5201 31
1 013 3429 52
00 200 011 4430 520
1 20
1 111 3431 52
00 29
1 013 4432 52
0
1 24
1 01
1 4

パッケージsurvival

関数coxphcoxph(Surv(week,arrest)~fin+age+prio, data=Rossi)
目的変数部

説明変数

データフレーム名

res.recidiv01<-coxph(Surv(week,arrest)~fin+age+prio, data=Rossi)

関数coxph基本的な回帰式の表現が使える

関数coxphの結果> res.recidiv01<-coxph(Surv(week,arrest)~fin+age+prio, data=Rossi)> res.recidiv01
Call:coxph(formula = Surv(week, arrest) ~ fin + age + prio, data = Rossi)
coef exp(coef) se(coef) zp
fin -0.3470 0.707 0.1902 -1.82 0.06800
age -0.0671 0.935 0.0209 -3.22 0.00130
prio 0.0969 1.102 0.0273 3.56 0.00038

切片の項はない

Likelihood ratio test=29.1 on 3 df, p=2.19e-06 n= 432, number ofevents= 114

関数coxphの結果にsummary関数を使う> summary(res.recidiv01)
Call:coxph(formula = Surv(week, arrest) ~ fin + age + prio, data = Rossi)
n= 432, number of events= 114
coef exp(coef) se(coef)z Pr(>|z|)fin -0.34695 0.70684 0.19025 -1.824 0.068197 .
age -0.06711 0.93510 0.02085 -3.218 0.001289 **prio 0.09689 1.10174 0.02725 3.555 0.000378 ***--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1exp(coef) exp(-coef) lower .95 upper .95
fin 0.7068 1.4148 0.4868 1.0263
age 0.9351 1.0694 0.8977 0.9741
prio 1.1017 0.9077 1.0444 1.1622
Concordance= 0.63 (se = 0.027 )Rsquare= 0.065 (max possible= 0.956 )Likelihood ratio test= 29.05 on 3 df, p=2.189e-06
Wald test= 27.94 on 3 df, p=3.741e-06
Score (logrank) test = 29.03 on 3 df, p=2.203e-06

1.0

関数survfitによる予測値の算出とプロット

0.6
0.4
0.2
0.0（予測値）

Prop. not arrested

0.8

逮捕されない率

0102030
weeks

時間4050関数cox,zphによる比例ハザード性の検定
残差の時間との相関> cox.zph(res.recidiv01)
rho
chisq
p
fin -0.00657 0.00507 0.9433
age -0.20976 6.54147 0.0105
prio -0.08004 0.77288 0.3793
GLOBAL NA 7.13046 0.0679

パラメトリックな回帰モデル
加速モデル
accelerated failure time model
略してAFT model

ポアソン回帰モデル

加速モデル
i番⽬目の個体の⽣生存時間

yi

i番⽬目の個体の説明変数の値

x1,i , x2,i ,!

log e ( yi ) = β 0 + β1 x1,i + β 2 x2,i +!+ ( error term )生存時間での対数線形性
説明変数と係数の部分の効果
時間の進み方が速くなる／遅くなる

加速モデル
i番⽬目の個体の⽣生存時間

yi

i番⽬目の個体の説明変数の値

x1,i , x2,i ,!

log e ( yi ) = β 0 + β1 x1,i + β 2 x2,i +!+ ( error term )生存時間での対数線形性
分布はいろいろワイブル分布、指数分布、ロジスティック分布、対数ロジスティック分布など加速モデル
i番⽬目の個体の⽣生存時間

yi

i番⽬目の個体の説明変数の値

x1,i , x2,i ,!

log e ( yi ) = β 0 + β1 x1,i + β 2 x2,i +!+ ( error term )パッケージsurvivalの関数survreg
基本的な回帰式の表現が使える
ワイブル分布、指数分布、ロジスティック分布、対数ロジスティック分布などの他、ユーザー定義も可能ポアソン回帰モデル
目的変数

死亡というイベントの数
説明変数の値が同じ個体のグループ
観察期間が同じ
説明変数の値が同じ個体のグループ
観察期間が同じ
説明変数の値が同じ個体のグループ
観察期間が同じ

死亡数
目的変数

説明変数
説明変数

個体数×観察期間
オフセット