[Rによるデータ分析入門]離散選択モデルの様々（５）：生存分析
本コラムは「離散選択モデルの様々」では、
（１）多項ロジット・モデル
（２）順序ロジット・モデル（Ordered Logit Model）
（３）ヘックマンの二段階推定モデル
（４）カウント・データ
を紹介してきましたが、（５）では生存分析で使われるKaplan-Mayer生存曲線の推計とCox比例ハザードモデルを紹介します。（１）～（４）を読みたい人は以下のリンクを参照してください。
生存分析, Survival Analysisとは
生存分析とは、元々、疫学や生物学分野で開発された手法で、あるイベントが発生するまでの時間（あるいは期間の長さ）を分析対象とする手法です。たとえば、ある病気の患者が完治するまでの時間や実験用マウスが死亡するまでの時間などを分析する際に使われます。最近では、社会科学分野でも用いられており、失業者が次の仕事を見つけて就職までの時間、結婚から妊娠するまでの年数、企業が倒産したり事業所が併産されたりするまでの時間、取引開始から取引終了までの時間などの分析に用いられています。
こうした期間の長さの分析には、次の二つの理由で、通常の回帰分析では問題が生じることが知られています。第一の理由は、あるイベントが発生する発生するまで期間の長さは必ず正の値をとりますが、通常の回帰分析を用いると予測値がマイナスになりうるという問題が生じます。ただし、この点はTobitモデルを適用することで対応が可能です。第二の理由は、この分析では、通常、観察期間の中にイベントが起こったサンプルとイベントが起こらなかったサンプルが混在していることが多いという点です。以下の図は、A、B、Cの3人の失業者がいつ就職したかを示すデータだと考えてください。AさんとBさんは分析期間中に就職しているのに対して、Cさんは分析対象期間終了時点でまだ失業状況にあります。AさんとBさんのサンプルのみで回帰分析を行うと分析結果に歪が生じることが知られており、こうしたデータを分析するための手法が開発されています。
なお、本コラムでは理論的な紹介は最小限にとどめRの使用方法を中心に紹介していますので、生存分析の背景等は以下の文献を参照してください。
山本勲「実証分析のための計量経済学」
ウィラワンードニ－ダハナ・勝又壮太朗「Rによるマーケティング・データ分析」新世社 第９章 ハザードモデル
サンプルデータ、Rのコードも付随しています。
では、具体的な事例を見ながらRによる分析方法を紹介します。ここで紹介するのは Cameron and Trivedi (2005)の第17章11節(17.11)で紹介されたMcCall (1996)の米国の労働者の失業期間の分析例です。データは米国のCurrent Polulation Surveyに付随する失職労働者のデータ（Displaced Worker Supplements）で、1986, 1988, 1990と1992年のデータです。
※ Cameron adn Trivedi (2005)では、多くの変数をコントロール変数として加えられていますが、本コラムでは簡略化のため変数の数を制限しています。
データは、著者のColin CameronのWEBサイトからダウンロードできます。
なお、本コラムで上記サイトのemma1996.dtaをCSV形式に変換したデータを使用します。本コラムの結果を再現する使用するRのスクリプトもおいておきます。
このデータに含まれる変数のリストが本コラムの一番下の補論に記載してあるので必要に応じて参照してください。
まずデータをread_csv()で読み込んで基本構造を確認しましょう。ここで分析対象にするのは失業期間spell（2週間単位）です。そして、調査期間中に（フルタイム雇用者として）再就職したかどうかを示すダミー変数がcensor1です。以下、spellとcensorの最初の10件を表示しています。
たとえば1行目の労働者のデータは、失業期間が５（単位が2週間なので10週間）でcensor1が1なので調査期間中にフルタイム雇用者として再就職していることが分かります。5行目の労働者は失業期間が９（18週間）ですが、censor1が0なので調査終了までに再就職できていないことがわかります。
spellの記述統計量も見ておきましょう。失業期間の平均が５、最大値は28になっています。
フルタイム再就職ダミーcensor1をみると、調査期間中に再就職できた人(censor1=1)が1073人、調査期間中に再就職できていない人（censor1=0)が2270人であることがわかります。
spellのヒストグラムも見ておきましょう。大半の人が失業期間10（20週）までに含まれていることがわかります。ただし、このヒストグラムでは調査終了時点で再就職できていない人が混在していることに注意をしてください。
記述分析：生存曲線
まず記述分析としてKaplan-Mayerの生存曲線を描いてみましょう。生存曲線とは、今回の例では（次の図１の実線）、失業期間に応じてどの程度の人が失業状態を継続しているかを示すグラフです。
このグラフの読み方ですが、横軸の５は失業期間が５、単位が2週間なので失業から10週間経過すると3割ぐらいの人が再就職する一方で、7割ぐらいの人が失業状態を継続していることを示します。そして、失業期間が25（50週間）でもまだ4割ぐらいの人が再就職できないでいることが分かります。なお、実線の上下の点線は95％信頼区間を示します。
この生存曲線を失業保険の受給の有無で比較してみましょう。このデータには失業保険の受給を申請したかどうかのダミー変数uiが含まれています。この変数を使って失業率の需給の有無（w/o UI, without UI、w UI, with UI）で失業期間が異なるかをグラフで示します。
このグラフから、たとえば失業期間が10（＝10週間）でみると、失業保険の給付を受けていない人で失業状態にある人は約50％、失業状態にない人は約73％ぐらいで、失業保険を受給している人の方が失業期間が長いことがわかります。
Rによる生存曲線グラフの作成
生存曲線を描くにはパッケージ{survival}のsurvfitを使用します。予めインストールし、library(survival)を呼び出しておきます。survfit()は以下のように使います。
result_haz1 <-survfit(Surv(期間を示す変数, イベントダミー)~X,data=データフレーム)
plot(result_haz1)
今回の場合、「期間を示す変数」は失業期間のspell、イベントダミーはフルタイム再就職ダミーであるcensor1を使います。Xはグループ分けするためのカテゴリー変数で、グループ分けしない図1のようなグラフを作成する場合は1と書きます。つまり、図1の場合は以下のスクリプトで作成できます。plot()に、xlab=とylab=というオプションがついていますが、これにより縦軸と横軸のラベルを付けることができます。
# Kaplan-Mayer survival curve
result_haz1 <-survfit(Surv(spell,censor1)~1,data=dataf)
plot(result_haz1,xlab="Duration",ylab="Survival ")図2は以下のスクリプトで作成できます。
# 失業保険受給の有無での比較
survdiff(Surv(spell,censor1)~ui,data=dataf)
result_haz2 <-survfit(Surv(spell,censor1)~ui,data=dataf)
plot(result_haz2,lty=c(1,2),xlab="Duration",ylab="Survival Probability")
legend("topright",c("w/o UI","w UI"),lty=c(1,2))なお、plot()関数のオプションltyはline typeの略で線の種類を番号で指定します。lty=c(1,2)で線種を最初の線は1番実線、次の線は2番点線で示せ、というオプションです。この辺りを詳しく知りたい人は以下のコラムを参照してください。
legend()関数はplot()に凡例を付ける関数で、"topright"で右上部に、"w/o UI"、”w UI”という凡例をつけltyで線種を指定しています。
なお、survfitの結果オブジェクトをsummary(result_haz2)で表示させると、以下のように生存確率の推計値が表示されます。
Cox propotional hazardモデル
ここで「期間」の決定要因を回帰分析のような方法で分析するCox比例ハザードモデル（Cox proportioanl hazard model）モデルの推計方法を紹介します。使用する関数はcoxph()関数で、使い方は以下の通りです。
coxph(Serv(期間を示す変数, イベントダミー)~x1+x2+x3, data=データフレーム)
例として、ui（失業保険給付受給ダミー）、logwage(前職の賃金)、tenure(前職の勤続年数)、houshead（家計の主たる稼ぎ手ダミー）、marrid（結婚ダミー）、female（女性ダミー）、child（子供の有無ダミー）を説明変数とするCox比例ハザードモデルを推定してみましょう。
結果の見方には注意が必要で、係数がプラスであれば当該変数が大きければ「イベントが発生しやすく（再就職しやすい）、期間が短い」ことを示します。UIの係数は-1.076とマイナスなので「再就職しにくい、つまり再就職までの期間が長い」ことを示します。前職の賃金が高い、主たる稼ぎ主である、結婚している、女性であれば再就職までの期間が短いと解釈できます。
各変数のインパクトですが、各変数が変化したときに翌期にイベントが発生する確率、すなわちハザード率がどの程度変化したか、に注目します。そして、係数を指数化したexp(coef)に注目します。Cox比例ハザードモデルでは変数がダミー変数の場合は、exp(coef)-1が「ダミー変数が０と１のときの生存率の差」を示します。uiの場合、exp(coef)は0.3409なのでexp(coef)-1=0.66となり、ハザード率（翌期にイベントが発生する確率）が66％低い（再就職しにくい）と解釈できます。一方、marriedのexp(coef)は1.4033なのでexp(coef)-1=0.40で40%ほど翌期にイベントがは発生する確率が高い（再就職しやすい）ことを示します。
本コラムは「Rによるデータ分析入門」のWEBサポートページとして作成されました。WEBサポートの一覧は以下を参照してください。
WEBサポートの一覧は以下を参照してください。
補論：変数名の一覧
spell periods jobless: two-week intervals
censor1 1 if re-employed at full-time job
censor2 1 if re-employed at part-time job
censor3 1 if re-employed but left job: pt-ft status unknown
censor4 1 if still jobless
ui 1 if filed UI claim
reprate eligible replacement rate
logwage log weekly earnings
tenure years tenure in lost job
disrate eligible disregard rate
slack 1 if lost job due to slack work
abolpos 1 if lost job due to position abolished
explose 1 if expected to lose job
stateur state unemployment rate in year job lost
houshead 1 if head of household
married 1 if married
female 1 if female
child 1 if have children
ychild 1 if have children under 6
nonwhite 1 if nonwhite
age age at time of survey
schlt12 1 if less than 12 years completed schooling
schgt12 1 if more than 12 years completed schooling
smsa 1 if resides in smsa
bluecoll 1 if lost bluecollar job
mining 1 if lost job in mining
constr 1 if lost job in construction
transp 1 if lost job in transportation
trade 1 if lost job in trade
fire 1 if lost job in finance,insurance, real estate
services 1 if lost job in services
pubadmin 1 if lost job in public administration
year85 1 if lost job in 1985
year87 1 if lost job in 1987
year89 1 if lost job in 1989
midatl 1 if reside in mid-atlantic region
encen 1 if reside in east north central region
wncen 1 if reside in west north central region
southatl 1 if reside in south atlantic region
escen 1 if reside in east south central region
wscen 1 if reside in west south central region
mountain 1 if reside in mountain region
pacific 1 if reside in pacific region