2016/05/22

生存時間モデルの推定
with R
まず，生存時間分析用ライブラリを下記コマンドでインストールして下さい

install.packages("survival")福田研・朝倉研合同ゼミ

交通行動分析ゼミ 生存時間モデル
山下さん解説資料

基礎概念

ある基準の時刻からある事象が生起、• 生存時間モデルとは あるいは終了するまでの時間の解析を対象とするモデル［期間モデル］• 打ち切り• 生存関数とハザード関数
特徴
解析対象となる事象が生起するまでの時間は必ず正の値
時間の分布の裾が右に長くなることが多い
分布の正規性を仮定することが適切でない場合が多い
来訪者
C
A

B
D

観測者

休憩

2

福田研・朝倉研合同ゼミ

交通行動分析ゼミ 生存時間モデル
山下さん解説資料

基礎概念• 生存時間モデルとは• 打ち切り• 生存関数とハザード関数

対象とする事象の生起するまでの時間 生存関数 𝑆 𝑡
Tの分布を以下のように表すことが多い 𝑆 𝑡 = 𝑃, 𝑇 ≥ 𝑡 = 1 − 𝑃, 𝑇 ≤ 𝑡 = 1 − 𝐹 t
ハザード関数 ℎ 𝑡生存関数 S（t）
－対象とする事象がある時点 ｔ においてまだ生起していない確率
ハザード関数 h（t）

𝑃, 𝑡 ≤ 𝑇 < 𝑡 + ∆𝑡|𝑇 ≥ 𝑡∆(→*∆𝑡ℎ 𝑡 = lim

3 ( 43((6∆(); <=> ? @=−∆(∗3 (;(∆(→*= lim=A (3 (𝑓 𝑡 はTの確率密度関数－対象とする事象がある時点 ｔ においてまだ生起していないという条件の下 生存関数とハザード関数の関係(で、次の瞬間に事象が生起するというS 𝑡 = exp − K ℎ 𝑢 d𝑢
条件付き確率*3

福田研・朝倉研合同ゼミ

交通行動分析ゼミ 生存時間モデル
山下さん解説資料

生存時間の解析方法• 生存時間のモデル化

共変量がパラメー
ターとしてモデルに導入されている

生存時間の分布に特定の確率分布が仮定されている

交通行動分
析の分野での使用• Non-parametric model• Semi-parametric model• Parametric model

××△〇×〇〇〇〇

Kaplan-Meier法
𝑑O ＝時点𝑡O で事象が生起したケース数、事象が生起した時点を𝑡P , 𝑡R ,…
𝑛P ＝時点𝑡O の直前のリスク集合の大きさ(＝その直前に観測されていたケース数)
生存時間が𝑡O4P より長いという条件の下で生存時間が𝑡O 以上である条件付き確率の推定量
𝑆T 𝑡𝑛O − 𝑑O=𝑛O
𝑆T 𝑡O4P
V 4X
生存関数S（ｔ）のKaplan-Meier推定量 𝑆T 𝑡 = ∏( Z( W WW

VW4

福田研・朝倉研合同ゼミ

交通行動分析ゼミ 生存時間モデル
山下さん解説資料

生存時間の解析方法• 生存時間のモデル化

共変量がパラメー
ターとしてモデルに導入されている

生存時間の分布に特定の確率分布が仮定されている

交通行動分
析の分野での使用• Non-parametric model• Semi-parametric model• Parametric model

××△〇×〇〇〇〇

比例ハザードモデルでのハザード関数ℎ 𝑡|𝑥O = ℎ* 𝑡 exp 𝛽𝑥Oℎ 𝑡|𝑥O :共変量ベクトル𝑥O を持つケースiのハザード関数、ℎ* 𝑡 ：基準ハザード、
𝛽：未知パラメーターベクトル、𝑥O ：ケースiの共変量ベクトル(ケース間でのハザード関数の比は時点によらず一定と仮定)
パラメーターの推定：部分尤度法、Han and Hausman法5福田研・朝倉研合同ゼミ

交通行動分析ゼミ 生存時間モデル
6
山下さん解説資料

生存時間の解析方法• 生存時間のモデル化

共変量がパラメーターとしてモデルに導入されている

生存時間の分布に特定の確率分布が仮定されている

交通行動分
析の分野での使用• Non-parametric model• Semi-parametric model• Parametric model

××△〇×〇〇〇〇①比例ハザードモデル②加速故障モデルℎ 𝑡|𝑥O = ℎ* 𝑡 exp −𝛽𝑥O

exp −𝛽𝑥O

生存時間への共変量の影響の仕方が異なる→２つのモデルのうち対象とする問題に適したモデルを選択する

確率分布

ハザード関数ℎ* 𝑡生存関数𝑆* 𝑡

指数分布λexp(−λt)ワイブル分布𝛾λ (_4Pexp(−λ ( )対数ロジスティック分布

𝛾λ (_4P

1__1+λ(log 𝑡 − 𝜇1𝜑𝜎
対数正規分布 𝑡elog 𝑡 − 𝜇1−𝜑𝜎
ゴンベルツ分布_𝛾exp(−λt)

1+λ(1−𝜑log 𝑡 − 𝜇𝜎exp((γ/λ)(exp(−λt)−1)

生存時間モデルの推定 with R

適用事例• 活動時間・滞在時間の分析• 個々人の日常的な自由活動の活動時間 / 観光地における滞在時間• 詳しく知りたい人は下記二編の論文を参照：

Bhat, C.R. and Pinjari, A.R. (2008). Duration Modeling. In David A. Hensher , Kenneth J.
Button (ed.) Handbook of Transport Modelling (Handbooks in Transport, Volume 1) ,pp.105-131.
http://www.caee.utexas.edu/prof/bhat/ABSTRACTS/Duration_Modeling_Nov1.pdf
Hensher, D.A. and Mannering, F.L. (1994). Hazard-based duration models and theirapplication to transport analysis, Transport Reviews, Vol. 14, No. 1, pp. 63-82.
http://courses.washington.edu/cee500/Transport%20Reviews.pdf• 災害発生時（e.g. ハリケーン・豪雨）の避難行動開始時間の分析等にも応用• 自動車保有行動の分析• 自動車保有期間・自動車取り替え行動の分析• 詳しく知りたい人は山本先生（名大）の博士論文のレビュー章を見て下さい山本俊行 (2000) 連続時間軸上における世帯の自動車保有更新行動及び世帯内での配分・利用行動に関する研究．京都大学博士論文．
http://repository.kulib.kyotou.ac.jp/dspace/bitstream/2433/151521/3/D_Yamamoto_Toshiyuki.pdf
7Activity-Based Travel Choice Modeling•   •  1     • +'0
% &1
% !%$ " • '),(.,(-*0/1"# ( " 1
 →→ →→ ) • 個⼈の⼀⽇の⽣活⾏動すべてを包括した分析• トリップだけでなく活動（Activity）とそのつながりを明⽰的に考慮• 個⼈⼀⼈⼀⼈の⼀⽇の活動のマイクロな動的シミュレーション8erogeneity in the evacuation behavior. In contrast to the above finding related to the household with one child, the indicatorvariable representing households with more than one child is found to be fixed with a positive parameter. This suggests thatH., thanRodrigo,
M.R.would
and Satish,U.the(2013).
A random-parameter
hazard-basedmodel tohouseholdsSamiul,with more
one childperceiverisk uniformly
resulting into an
earlier evacuation.
understandhouseholdevacuation
behavior, Transportation
research
parttheC:trends
EmergingFig. 2 presents
the actualand predictedhazard timing
for the households.
The predicted hazard
capturesof the actual
hazard values.
The unusualVol.
drop27,in theis due to the hurricane landfall event when no households actually evacuate.
technologies,
pp. hazard
108-116.

Fig. 2. Actual and predicted hazard at different departure times.

Hurrycane Ivan (2004年) 襲来時のアラバマ州の避難行動調査を使った分析
9ඪͷಛ௃ʹ͍ͭͯड़΂Δɽ

2.1

ௐࠪ֓ཁ

本日の演習で用いるデータશࠃ౎ࢢަ௨ಛੑௐࠪ (‫چ‬শ:શࠃύʔιϯτϦοϓௐࠪ) ͸ɼௐࠪର৅ͷ̍೔ͷߦಈΛ௥͍͔͚Δ͜ͱʹΑΓਓʑͷަ௨ͷಛ௃Λଊ͑Δ͜ͱΛओ໨తͱ͢ΔௐࠪͰ͋Γɼত࿨ 62 ೥ͷୈ 1 ճௐ͔ࠪΒ਺͑ͯࠓճͷฏ੒ 22 ೥ͷௐࠪͰ 5 ճ໨Ͱ͋ΔɽҎલ͸ௐࠪһ͕‫ݸ‬ผʹ๚໰ͯ͠ௐࠪ΁ͷ‫ྗڠ‬ґཔ͓Αͼճऩ‫ۀ‬຿Λߦ͍ͬͯͨͷʹର͠ɼࠓճॳΊͯ༣ૹʹΑΔௐ͕ࠪߦΘΕͨ͜ͱ΋ಛ௃ͱ• 平成22年全国都市交通特性調査（全国PT調査）ͯ͠‫͛ڍ‬ΒΕΔɽௐࠪର৅೔͸ฏ੒ 22 ೥ͷ 10 ݄·ͨ͸ 11 ݄ͷฏ‫ٳ‬೔֤Ұ೔Ͱ͋Γɼฏ೔͸Ր༵͔Β໦༵ͷ͍ͣ
• 各都市圏で行われる通常Person
Trip調査と異なり，Ε͔ɼ‫ٳ‬೔͸࿈‫͍ͳ͍ͯͬͳʹٳ‬िͷ೔༵ͱ‫ݶ‬ఆ͞Ε͍ͯΔɽௐࠪର৅ͱͳΔॅຽ͸ҎԼͷද 2.1
規模や構造の異なる各都市で満遍なく調査を行い，都市ͷ֘౰౎ࢢͷॅຽ‫ج‬ຊ୆ாʹ‫͖ͮج‬ແ࡞ҝநग़͞Εͨੈଳͷ 5 ࡀҎ্ͷॅຽͰ͋Δɽͳ͓ຊ‫ڀݚ‬
構造が交通行動に及ぼす影響等を横断的に分析するのに，Ͱ͸ொଜ෦͸औΓѻΘͣɼ౎ࢢ෦ (ࢢ‫ͼٴ‬౦‫ژ‬ಛผ۠) ͷΈऔΓѻ͏ɽ
主に用いられる．
ද 2.1: ௐࠪର৅౎ࢢ10https://trafficnews.jp/post/62155/11

https://trafficnews.jp/post/62155/2/12౎ࢢަ௨ʹؔ͢Δҙࣝௐࠪɾɾɾަ௨΁ͷҙࣝ΍ར༻܏޲Λฉ͘΋ͷ

坂東徹 (2013). 都市構造が個人の交通行動に及ぼす影響の分析.
·ͨɼ۩ମతͳ࣭໰߲໨ΛҎԼͷද 2.2 ʹ‫͠ه‬ɼͦͷ͏ͪຊ‫Ͱڀݚ‬આ໌ม਺͓Αͼඃઆ໌ม਺ͱ͠
ͯ༻͍ΔσʔλΛද 2.3 ʹࣔ͢ɽ
東京工業大学工学部土木・環境工学科卒業論文.
ද 2.2: ࣭໰߲໨ද 2.3: આ໌ม਺͓Αͼඃઆ໌ม਺ͱͯ͠༻͍Δσʔλ

生存時間分析モデル
13

モデル構築の目的• 都市構造，個人のトリップ属性や社会経済属性が私事行動における活動範囲に及ぼす影響を分析する．
• より具体的には，「どれぐらいの距離まで，人々は買い物トリップを行うのか？」について，生存時間分析を適用し，その規定要因について実証的に検討する．
• 「都市交通に関する意識調査票」の回答結果を利用する

14

参考：交通実態調査票

意識調査票（表面）15参考：交通実態調査票

意識調査票（裏面）16公共交通利用者の日常買い物目的
トリップ長分布（n=36816）
Min.
0.100

1st Qu.
1.000

Median Mean
2.000
2.431

3rd Qu. Max.
3.000
12.000NA's
12141

17後期高齢者と非後期高齢者のトリップ長分布（経験分布）

18

Kaplan-Meier法により後期高齢者と非後期高齢者のトリップ長分布を求める

sa.km <- survfit(Surv(日常買い物距離) ~ 高齢D,data=traveldata)sa.km
summary(sa.km)

19トリップ長を決める要因は多数ある→Semi-parametric/Parametric models（Cox）比例ハザードモデル［セミパラメトリックモデル］

h(t|x) = h0 (t) exp( t x) = h0 (t) exp

MX

i xi

i=1

“ベースラインハザード”
ここはノンパラメトリック!

“共変量”の値に依存してハザードが変化する
ここはパラメトリック

ハザードの比率が時間 t に依存せず常に一定となる，強い仮定を置いたモデル．
どのような目的関数を最大化して最適パラメータ b を求めているかについては各自教科書等で確認すること．
20

21コンパクトシティ指標の作成1.

人口二次半径： RP (km) [都市単位]2.

相対的遠距離性： RF [大字(個人)単位]
𝑅𝐹 =都市Xijk

人口重心dkRP

大字 k22

2013/12/05

コンパクトシティ指標の作成
3.

人口集中度： CP [都市単位]

DID(mnm o,po)
𝐶𝑃 =qjkr

RP

4. 他地域人口密度集積性：CI (people/km2)[大字(個人)単位]∑|W}iwxW(yWzi ){u𝐶𝐼t = |~∑W}i(y{ HighWzi )23

人口密度高CILow

2013/12/05

低Cox比例ハザードモデルの当てはめ
sa.fit.cox<-survfit(sa.cox)
sa.fit.cox
summary(sa.fit.cox)
plot(sa.fit.cox)
残差分析，ハザードの比例性の妥当性分析 等が必要

24いろいろなハザード関数 h(t)(a)(b)

Weibullexponential model

tOt

O(c)

(d)

Weibullp >1 のとき

p <1 のとき

lognormal

Commonly Used Distributions and ParameterstOOt

Distributionf (t)
S(t)
h(t)
図 15-1: いろいろなハザード関数
Exponentialexp( t)exp( t)pt p 1 exp(t pσ)2 ) に従いexp( t,p )T の密度関数が次のようになる
pt p 1(iii) 対数正規分布 Weibull
log T が正規分布
N(µ,.
pt p 1pt p 1125
Log-logistic
1+ t p1+ t p(1+ t p )2

√1−(log t−µ)22σ 2Weibull ModelAssuming T ⇠ Weibull( , p) with probability density function
f (t) = pt p 1 exp(t p ), where p > 0 andWeibull型
Hazard-Function
the hazard function is given by> 0,h(t) = pt p 1 .
Hazard Function h(t)p<1p=1p>11.5

p is called shape parameter:

1.0

I If p = 1 the hazard is constanth(t)I If p > 1 the hazard increases

I If p < 1 the hazard decreases

0.5(exponential model)

02

46t

81026

パラメトリックモデルの推定
ワイブル分布のハザード関数：

p 1h(t) = pt
ワイブル加速時間故障モデル：

1⌘exp(+01/p1 x1n + · · · +M xM n )指数分布ー＞”exponential”
対数正規分布ー＞”log-normal”
対数ロジスティック分布ー＞”log-logistic”
27

Cox比例ハザードモデルの場合と多くの共変量の係数の符号が逆になっていることに注意
注意# survreg’s scale = 1/(rweibull shape)# survreg’s intercept = log(rweibull scale)であるので，𝛼€ =1/0.772=1.295
𝜆T = exp 0.669 = 1.952となる

どの確率分布を用いるのがベストかに関しては，
AIC = -2*最大対数尤度 + 2* パラメタ数を用いて評価することができる

28

表 5−3自動車保有
期間モデル（山本先生D論第五章）

ワイブル分布を想定
Rではないので符号は直接解釈してよい

定数項
自動車属性
保有形態

保有期間モデルの推定結果
重み付きモデルA
推定値t 値4.348-0.462
リース車-0.625
社用車-0.531
購入時
中古車
Mini,
Subcompact
0.169
車種
Large, Luxury
0.132
Sport car0.288
Pick-up truck
0.057
Van
0.008
Utility vehicle-0.070-0.387
年間走行距離/10,000 mile
世帯属性
0.089
世帯構成
夫婦と子供（0-15才）
0.152
夫婦と子供（16-20才）-0.012夫婦と子供とその他の大人
0.043
独身者
0.099
夫婦とその他の大人-0.026
子供数-0.088
大人数-0.113賃貸
0.167
自動車保有台数
運転者属性
1.697
年齢/100才-0.100性別
女性-0.528職業営業-0.189自営業-0.092
退職者
0.048
個人年収
$20,000未満-0.517$100,000以上
0.685
形状パラメータ(1/γ)
8143
サンプル数
L(C)
-6377
-5883
L(β)
2988(26)
χ (df)

重みなし
モデルB
推定値
t 値
4.821

t 値
-2.27
**

-4.63
-5.46
-15.77
4.22
2.85
5.75
1.22
0.13
-1.16
-3.52

-0.498-0.402-0.619
0.131
0.119
0.281
0.046
-0.004-0.098-0.546-3.19-1.78-12.55
1.73
1.32
2.94
0.51
-0.03
-0.85
-6.27

0.19
-0.88
1.48
0.44
0.13
0.06
0.11
0.09
0.22
1.13

1.56
2.23
-0.13
0.80
1.75
-1.13
-2.84
-3.12
8.81

0.119
0.213
-0.0240.265
0.128-0.051-0.110-0.1860.322

1.36
1.86-0.16
3.05
1.40
-1.46
-2.14
-3.40
9.16

-0.29
-0.46
0.07
-2.18
-0.27
0.60
0.38
1.11-3.87

11.78
-3.22
-4.00
-1.43
-1.51
1.30
-6.24
10.55*

1.702-0.080-0.698-0.252-0.0280.098-0.546
0.664
8143
-2973
-2706
536(26)

7.28
-1.60
-3.99
-1.27
-0.29
1.65
-4.27
20.45*

-0.02
-0.34
0.78
0.27
-0.54
-0.70
0.190.63*H0 : γ=1 に対する t 値，**H0 : βa =βb に対する t 値(βa はモデル A による推定値，βb はモデル B による推定値)

29