Pythonでデータサイエンス(生存時間解析)

ビジネスの現場では、縦軸に発生日、横軸に経過日の表での残存をみる俗にいうコホート分析が使われることがよくあります。

これはビジネスのデータは日々更新され、始点が異なることが多いのと、終点がはっきりしないことが多いためです。

ロジスティック回帰などを用いる場合、終点がはっきりしないデータは扱いにくくなり、個人的には生存時間解析の方が使いやすいと思うことが多くなります。

使い方

シンプルにイベントの発生期間までをモデリングしたい場合、カプランマイアーを使います。

このとき、データの作り方をよく考える必要があります。(後述)

from lifelines import KaplanMeierFitter as kmf

mdl = kmf()
mdl.fit(経過期間,イベント発生有無)

#残存を可視化してくれる。
mdl.plot()

対して、回帰のように各変数の影響をみたい場合は、CoxPHを用います。

from lifelines import CoxPHFitter as cph

df
cph_mdl = cph()
cph_mdl.fit(データ,経過期間カラム名,イベントカラム名)
#一般的な回帰のようなsummary
cph_mdl.print_summary()

#変量の影響度の可視化
cph_mdl.plot()

使ってみる

データをつくる

今回は、二つのアプリケーションにログインするユーザが、解約するというシナリオを考えます。

この場合、観察できる期間=ログインし続けているか、と解約する(=イベント)か、というのは別の概念になっています。

これが逆にログインしなくなる場合は、イベント発生はログインが見えなくなった時ですし、そもそもログイン観察ができない場合は、前の行動からの期間を期間と設定する必要がでてきます。
ここのデータの作り方がモデリングに与える影響が強いのが生存時間解析の癖があるところです。

import numpy as np
import pandas as pd
import random

#経過期間はワイブル分布するとする
a = 1.5
s = np.random.weibull(a, 1000000)

#ユーザークラス
class User:
  def __init__(self):
    self.lifetime = random.choice(s) 
    self.application_type = random.choice(['a','b'])
    #アプリケーションによって解約発生率を変える
    self.churn = random.choice([0,0,0,1,1]) if self.application_type == 'a' else random.choice([0,0,1,1,1])

#データフレームにする

buf = []
for i in range(30):
  u = User()
  buf.append([i,u.lifetime,u.application_type,u.churn])

df = pd.DataFrame(buf,columns = ['id','lifetime','application_type','churn'])

フィッティングしてみる

from lifelines import KaplanMeierFitter as kmf
mdl = kmf()

mdl.fit(df['lifetime'],df['churn'])

mdl.plot()

f:id:esu-ko:20200719125517p:plain

横軸は経過期間、縦はそれに応じた、イベントが発生しなかった残存率です。今回は経過期間こそワイブル分布にしましたが、チャーンの発生は経過期間と連動していないため、なかなか下げ止まらなくなっています。

# coxphは数値のみ
d2 = df[['lifetime','churn']]
d2['is_a'] = df['application_type'].apply(lambda ) 

from lifelines import CoxPHFitter as cph

cph_mdl = cph()
cph_mdl.fit(d,'lifetime','churn')
cph_mdl.print_summary()

coxphもやってapplication_typeの影響をみてみます。coxphは数値のみで、かつワンホットエンコーディングするとエラーを返すので、1つに減らした状態でやります。

▼get_dummiesあたりの話はここ esu-ko.hatenablog.com

結果はprint_summaryでみることができます。
全体数では30ですが、イベント発生は13ユーザーのみ、といった基本情報と、回帰係数の情報をみることができます。

application_type_bはexp(coef)が2.61なので、bだと(aにくらべて)churnが2倍以上起きやすくなるようですね。ただしp値は0.15なので有意とはいえなさそうです。

変量間の違いの可視化

今回は説明変数が一つしかないため、あまり意味がありませんが、一応見ておきます。

cph_mdl.plot(hazard_ratios=True)

for app in df['application_type'].unique():
  mdl = kmf(label=app)
  tmp = df[df['application_type'] == app]
  mdl.fit(tmp['lifetime'],tmp['churn'])

  mdl.plot()

f:id:esu-ko:20200719125756p:plain

hazard_ratios=Trueをわたすことでexp(coef)の95%区間を表示します。係数間を一度に比較する時に便利です。

f:id:esu-ko:20200719125422p:plain

属性間の残存を見比べたいならkmfでplotの出し分けをします。確かにaの方がよく残存しているように見えます。