はじめに
Sales cycle（見込み客を見つけて結果が出るまでの期間）を時変 Weibull 分布の生存時間分析で推定する記事です。
もっと Stan ゴリゴリの内容によせて Stan アドベントカレンダーに載せたほうがよかったかも
想定する読者とニーズを満たせない方は以下を想定しています。
-
想定する読者
- ビジネスサイドの人でいろいろなDXのネタを集めている方
- ビジネスサイドの人でデータサイエンスについて知りたい方
- Sales の人で Sales cycle に課題を感じている方
- データサイエンティストを志している方
- Stan のコードを見たい方
-
ニーズを満たせない方
- ゴリゴリの数式が知りたい方
- 収束のテクニックなど Stan の細かなテクニックを知りたい方
（少し長くなってしまったので、ビジネスサイドの方は前半部分をメインに、データサイエンティストを志している方は後半部分をメインに読んでいただくのが良さそうです）
背景: Sales cycle 改善の難しさ
前に Web マーケをやっていたことがあるのですが、ABテストを高速・大量に実施することで数ヶ月でCVRが2.5倍になるなど、PDCA を高速で回していくことで成果に結びつくことを肌で感じていました。
世の中全体でクロックスピードが上がっていくなかで、その流れに取り残されないためには Web マーケの領域だけでなく会社のあらゆるところでクロックスピードを上げていくことが重要です。
ただ、Web マーケでは高速に PDCA を回して改善していくという文化が比較的根付いている印象ですが、営業で同じようなことを実施できている企業は多くはないのではないでしょうか。
PDCA のスピードを上げていくためにもキャッシュフローを良くするためにも、意図的に Sales cycle にフォーカスして施策を回すできると理想です。
ただ、それをやろうと思っても、Sales cycle は計算することが難しいため、数字をもとに改善していくのは難しいです。
特に既存業務を設計しなおす必要があるソリューション（SFAなど）や比較的金額の大きい商品（家）などは、意思決定までの期間がながく、かつデータ（=商談数）も少ないので、計算が難しいことが多いです。
具体的に Sales cycle の計算の難しさを見てみましょう。
正確な値が出るまでの時間が長すぎる問題
例えば、Sales cycle が長いもので1年かかると仮定します。
すると、正確な値が出るのは1年後です。
1年前の流入や施策、営業人員がどのように影響していたかを振り返るのは難しいですし、たとえ正確に思い出して振り返ったとしても改善は1年後にしか実施できず、クロックスピードとしては致命的です。
もう少し具体的なデータ・計算をもとにこの問題を深ぼっていきます。
ここからはイメージしやすいように5件のデータを例に計算していきます。
| No. | 商談発生時期 | 商談結果時期 | Sales cycle |
|---|---|---|---|
| 1 | 2020/12 | 2021/1 | 1ヶ月 |
| 2 | 2020/12 | 2021/3 | 3ヶ月 |
| 3 | 2020/12 | 2021/6 | 6ヶ月 |
| 4 | 2020/12 | 2021/9 | 9ヶ月 |
| 5 | 2020/12 | 2021/12 | 12ヶ月 |
（平均すると6.2ヶ月）
結果が出ているものに絞った計算した場合
「暫定でいいから、早めに計算しよう」と、一定期間で結果が出たものだけで計算したとします。
ここでは仮に3ヶ月で計算したとします。つまり、2020年1月に発生した商談の結果を2020年3月で判断するとします。
すると計算に使うデータは......
| No. | 商談発生時期 | 商談結果時期 | Sales cycle |
|---|---|---|---|
| 1 | 2020/12 | 2021/1 | 1ヶ月 |
| 2 | 2020/12 | 2021/3 | 3ヶ月 |
$$
\frac{1+3}{2}=2
$$
となり、平均 Sales cycle は2ヶ月と算出され、本当の6.2ヶ月からは程遠い結果となってしまいます。
結果出ていないものを含めた計算した場合
結果が出ているものだけで計算すると、セレクションバイアスが発生してしまい正確な値がでないので、対策として結果が出ていないものを含めて計算してみましょう。
具体的には分母の件数に結果が出ていないものを含めます。
| No. | 商談発生時期 | 商談結果時期 | Sales cycle | 備考 |
|---|---|---|---|---|
| 1 | 2020/12 | 2021/1 | 1ヶ月 | |
| 2 | 2020/12 | 2021/3 | 3ヶ月 | |
| 3 | 2020/12 | --- | -ヶ月 | 計算の段階では不明 |
| 4 | 2020/12 | --- | -ヶ月 | 計算の段階では不明 |
| 5 | 2020/12 | --- | -ヶ月 | 計算の段階では不明 |
平均すると
$$
\frac{1+3}{5}=0.8
$$
と、分母が増えてしまったので、前の計算より短く出てしまいました笑
（ここまであからさまな計算ミスをする人はいないとは思いますが）
結果が出ていないものを今結果が出たと仮定して計算した場合
そこで、計算段階で結果が出ていものについては、その瞬間に結果が出たと仮定して計算してみます。
| No. | 商談発生時期 | 商談結果時期 | Sales cycle | 備考 |
|---|---|---|---|---|
| 1 | 2020/12 | 2021/1 | 1ヶ月 | |
| 2 | 2020/12 | 2021/3 | 3ヶ月 | |
| 3 | 2020/12 | *2021/3 | *3ヶ月 | 結果が出たと仮定 |
| 4 | 2020/12 | *2021/3 | *3ヶ月 | 結果が出たと仮定 |
| 5 | 2020/12 | *2021/3 | *3ヶ月 | 結果が出たと仮定 |
$$
\frac{1+3+3+3+3}{5}=2.6
$$
となり、最初の計算結果（2ヶ月）よりは良くなりましたが、やはり本当の6.2ヶ月からは程遠い結果となってしまいます。
一定期間で結果が出る商談の割合を計算した場合
そもそも Sales cycle の平均値を KPI としているので、どうしても正確な値にならないです。
なので、平均値を KPI にすることを捨てて、「一定期間で結果が出る商談の割合」を KPI にするとかなりましな意思決定になります。
具体的には
| No. | 商談発生時期 | 商談結果時期 | Sales cycle | 結果の有無 |
|---|---|---|---|---|
| 1 | 2020/12 | 2021/1 | 1ヶ月 | 有 |
| 2 | 2020/12 | 2021/3 | 3ヶ月 | 有 |
| 3 | 2020/12 | --- | -ヶ月 | 無 |
| 4 | 2020/12 | --- | -ヶ月 | 無 |
| 5 | 2020/12 | --- | -ヶ月 | 無 |
$$
3ヶ月で結果が出る商談の割合 = 0.2
$$
となります。
ただ、多くの人・組織は平均値で KPI を追っていることが多いので、上記方法だと理解を得られるかが問題となります。
また、どの値が適正値なのかも判断が難しくなるので、過去データとの比較をして感覚を養う必要があります。
解決方法
このように単純な計算だけで Sales cycle を算出しようとすると正確な値を算出するのが難しいです。
今回は時変 Weibull 分布の生存時間分析を用いることで（完璧ではないですが）より正確に推定していきます。
理論的な話
今回使う知識は
- 生存時間分析
- ワイブル分布
- 打ち切り
- 時変モデル
ですので、順を追って説明していきます。
（細かく説明するとそれだけで1つの記事になるので、簡単に説明します）
また、小難しい内容になるので面倒な人は「推定」のパートまで飛ばしてください。
生存時間
生存時間分析は、あるイベントが発生するまでの時間を分析する手法です。
例えば、治療後の患者の生存や機械の故障の発生といったもので扱われます。
簡単に言うと、時間とともに結果が出る確率が変化するので、それを明確にする分析です。
今回は商談の結果が出る（イベント）までの時間を分析します。
生存時間分析を理解するためには
- 確率密度関数
- 生存関数
- ハザード関数
という3つの関数を理解することが重要です。
1. 確率密度関数: f(t)
イベント（商談結果が出る）までの時間の分布を表す関数です。
先の例で言う Sales cycle を発生させる確率関数です。
| No. | 商談発生時期 | 商談結果時期 | Sales cycle |
|---|---|---|---|
| 1 | 2020/12 | 2021/1 | 1ヶ月 |
| 2 | 2020/12 | 2021/3 | 3ヶ月 |
| 3 | 2020/12 | 2021/6 | 6ヶ月 |
| 4 | 2020/12 | 2021/9 | 9ヶ月 |
| 5 | 2020/12 | 2021/12 | 12ヶ月 |
2. 生存関数: S(t)
確率密度関数を積分すると、ある時間tまでにイベントが発生する確率が算出されます。
「3ヶ月で商談結果が出る確率」などです。
生存関数S(t)はこの逆で、ある時間tまでにイベントが発生しない確率です。
なので、1 - (確率密度関数を積分した値)で計算できます。
$$
S(t) = 1 - \int_0^t f(x)dx = \int_t^\infty f(x)dx
$$
確率密度関数と生存関数の関係をグラフで見てみましょう。
上の図では上部が確率密度関数で、下部が生存関数です。
t=50~60ぐらいで確率密度関数が高くなり（イベントが発生する確率が高くなり）、一方で生存関数は急激に下がっています。
3. ハザード関数: h(t)
ある時間tまでイベントが発生しなかった状態で、その瞬間にイベントが発生する確率を表す関数をハザード関数h(t)と言います。
確率密度関数との関係がややこしいですね笑
整理すると、確率密度関数がどの時間でイベントが発生するかを表す関数なのに対して、ハザード関数はイベントが発生しなかった前提でその瞬間にイベントが発生を表す関数です。
Weibull 分布
確率密度関数の中身としてはいろいろな関数が使われますが、今回は Weibull 分布を使います（他にも指数分布などが使われます）。
この分布は時間に対する劣化現象や寿命に関する分布で、鎖を引っ張る場合において最も弱い輪が破壊することにより鎖全体が破壊したとするモデルと言われています。
具体的には以下の確率分布をもつ分布です。
$$
f(t) = \frac{m}{\eta}\left(\frac{t}{\eta}\right)^{m-1} \exp\left[-\left(\frac{t}{\eta}\right)^m\right]
$$
mは形状パラメーター、$\eta$は尺度パラメーターと呼ばれます。
mは形状とついているとおり、これにより分布の形状が変わります。
- m < 1 のときは、時間とともにイベントが発生する確率が小さくなる
- m = 1 のときは、時間に対して一定確率でイベントが発生する（指数関数）
- m > 1 のときは、時間とともにイベントが発生する確率が高くなる
よりイメージしやすく、故障の例で説明すると
- m < 1: 初期故障（最初に故障が多い）
- m = 1: 偶発的故障（故障する確率は時間に関係ない）
- m > 1: 摩耗的な故障（時間が経つにつれ摩耗し故障が発生しやすくなる）
具体的に m を変化させた場合の生存関数と確率密度関数を見てみましょう。
- m=0.2 は初期故障型なので、最初はイベントが多く発生しますが、その後は急激にイベントが発生しなくなっている
- m=1は0.2の場合と比べて、確率密度関数が緩やかに減っていているのが分かる
- m=3はt=50あたりにピークが来ており、一定時間が経った際にイベントが発生するのが分かる
また、もう一つのパラメーターである$\eta$はイベント発生するまでの全体的な時間の長さを表しています。
$\eta$が大きくなると全体的にイベントが発生するまでの時間が長くなります（グラフで右側にズレていく）。
このように Weibull 分布はいろいろな形を取ることができるので、柔軟に状況を推測することができます。
打ち切り
このまま確率密度関数をもとに計算しようとすると、すべての結果が出るまで待つ必要があります。
これだと、前半部分で説明した問題（全部の結果が出るまで待っていられない）が発生してしまいます。
そこで使うテクニックが「打ち切り」です。
「まだ途中だけど、ここまではイベント（商談結果が出る）が発生しなかった」ということを計算に組み入れる手法です。
ここを解説していると長くなるので、結果が出ていなくても、「ここまで結果が出なかった」という情報を推定に使っているんだな、程度の理解をしていただければ大丈夫です。
時変モデル
今回、営業活動によって Sales cycle が時間とともに変化していく、と仮定しています。
例えば、2019年1月での Sales cycle は長かったけど努力して 2020年12月時点では短くなっている、という状況です。
なので、モデルでもそれを反映させる必要があります。
具体的には Weibull 分布のパラメーター（mと$\eta$）が時間とともに変化していく（時変）というモデルを使います。
ここで計算の工夫として、t時点でのパラメーターはt-1時点でのパラメーターから少しズレたものであるという仮定を用います。
先月まで Sales cycle が6ヶ月だったのが、今月になって劇的に改善して1ヶ月になったということがない、という仮定です（なだらかに変化するという仮定）。
具体的にはt時点のパラメーターとt-1時点のパラメーターの差が平均0の正規分布に従うというモデルにします。
$$
m(t)-m(t-1) \sim \text{Normal}(0, \sigma_m) \\
\eta(t)-\eta(t-1) \sim \text{Normal}(0, \sigma_{\eta})
$$
要は急激に変動することなく、1個前の値に近くなるというモデルです。
推定
ここからは具体的なデータをもとに推定を行っていきます。
本当は Python で書くつもりだったのですが、時間がなかったので前に書いたコードを流用できる R にしました。
（暇ができたら Python ver. を追記するかもしれないです）
使うパッケージは tiyverse, lubridate, rstan, ggmcmc あたりです。
# tidyverse
library(tidyverse)
library(lubridate)
# stan
library(rstan)
library(ggmcmc)
# others
library(here)
データの生成
1ヶ月あたり50個の商談が発生する状態が24ヶ月続いたと仮定してデータを作成します。
Sales cycle のデータを生成する前に、まずは Weibull 分布パラメーターを作成します。
para <-
tibble(t = 1:24, m = seq(1.2, 0.8, by=-0.4/23)) %>%
mutate(eta = (t - 12) ^ 2 / 4 + 60) %>%
return()
para
## # A tibble: 24 x 3
## t m eta
## <int> <dbl> <dbl>
## 1 1 1.2 90.2
## 2 2 1.18 85
## 3 3 1.17 80.2
## 4 4 1.15 76
## 5 5 1.13 72.2
## 6 6 1.11 69
## 7 7 1.10 66.2
## 8 8 1.08 64
## 9 9 1.06 62.2
## 10 10 1.04 61
## # … with 14 more rows
月ごとに Weibull 分布のパラメーターが変化する状態を作っています。
具体的にはこんな形にしました。
- $\eta$: 2次関数的に真ん中で凹んでいる
- 最初 Sales Cycle が長かった
- 次第に Sales Cycle が短くなった
- しかし、最近また Sales Cycle が長かった
- という状態です。
- m: 1.2から0.8に変化
- 最初は摩耗的な故障タイプだった
- 最近は初期故障タイプになってきている
次に実際のデータを作ります。
data <- NULL
for (t in 1:24) {
data_t <- tibble(sales_cycle = rweibull(50,
shape = para[t, ]$m,
scale = para[t, ]$eta)) %>%
mutate(t = t)
data <- data %>% bind_rows(data_t)
}
data <-
data %>%
mutate(sales_cycle = sales_cycle %>% round() + 1) %>%
bind_cols(tibble(days = runif(nrow(data), min = 0, max = 30) %>% round())) %>%
mutate(start_date = ymd("2018/12/01") + months(t - 1) + days(days),
end_date = start_date + days(sales_cycle),
censored_sales_cycle = ifelse(
end_date > ymd("2020/11/30"),
as.integer(ymd("2020/11/30") - start_date), sales_cycle
),
censored_flg = ifelse(end_date > ymd("2020/11/30"), 1, 0)) %>%
filter(start_date < ymd("2020/11/30")) %>%
select(t, start_date, end_date, sales_cycle, censored_sales_cycle, censored_flg) %>%
return()
data
## # A tibble: 1,196 x 6
## t start_date end_date sales_cycle censored_sales_cycle censored_flg
## <int> <date> <date> <dbl> <dbl> <dbl>
## 1 1 2018-12-20 2019-06-11 173 173 0
## 2 1 2018-12-11 2019-01-29 49 49 0
## 3 1 2018-12-23 2019-02-12 51 51 0
## 4 1 2018-12-08 2019-01-26 49 49 0
## 5 1 2018-12-02 2018-12-22 20 20 0
## 6 1 2018-12-30 2019-02-15 47 47 0
## 7 1 2018-12-25 2019-11-16 326 326 0
## 8 1 2018-12-27 2019-05-01 125 125 0
## 9 1 2018-12-16 2019-01-29 44 44 0
## 10 1 2018-12-23 2019-02-26 65 65 0
## # … with 1,186 more rows
Stan の実装
今回はモデルの推定のために Stan という言語を使います。
時変性を組み込もうとすると階層ベイズモデルが便利なのですが、Stan はその階層ベイズを実装するための言語です。
まずは Stan に渡すためにデータを加工しておきます。
stan_data <-
list(
N = data %>% nrow(),
t = data %>% pull(t),
sales_cycle = data %>% pull(censored_sales_cycle),
censored_flg = data %>% pull(censored_flg),
t_max = data %>% pull(t) %>% max(),
sigma_m_max = 5,
sigma_eta_max = 150
)
map(stan_data, head)
## $N
## [1] 1196
##
## $t
## [1] 1 1 1 1 1 1
##
## $sales_cycle
## [1] 173 49 51 49 20 47
##
## $censored_flg
## [1] 0 0 0 0 0 0
##
## $t_max
## [1] 24
##
## $sigma_m_max
## [1] 5
##
## $sigma_eta_max
## [1] 150
つぎに Stan でモデルを記述します。
data {
int N; // レコード数
int t[N];
vector<lower = 0>[N] sales_cycle;
vector<lower = 0, upper = 1>[N] censored_flg; // 打ち切りしたか否か
int t_max;
// パラメーターの変動性の最大値を設定
real<lower = 0.0001> sigma_m_max;
real<lower = 0.0001> sigma_eta_max;
}
parameters {
// weibull のパラメーターの宣言
vector<lower = 0.001>[t_max] m;
vector<lower = 0.001>[t_max] eta;
// 時変のパラメーターを宣言
real<lower = 0> sigma_m;
real<lower = 0> sigma_eta;
}
model {
// Weibull 分布
for (n in 1:N) {
if (censored_flg[n] == 1) {
// 打ち切りの場合
target += weibull_lccdf(sales_cycle[n] | m[t[n]], eta[t[n]]);
} else {
// 観測された場合
target += weibull_lpdf(sales_cycle[n] | m[t[n]], eta[t[n]]);
}
}
// パラメーターの時変性
for (month in 2:t_max) {
m[month] ~ normal(m[month - 1], sigma_m);
eta[month] ~ normal(eta[month - 1], sigma_eta);
}
// 時変の変動性
sigma_m ~ normal(0, sigma_m_max);
sigma_eta ~ normal(0, sigma_eta_max);
}
このコードをcensored-weibull.stanという名前で保存しておきます。
このコードをもとに学習していきます。
# モデルの設定
weibull_model <- stan_model(here::here("scr/stan/censored-weibull.stan"))
# 学習
weibull_fit <- sampling(weibull_model, data = stan_data)
Stan 結果の可視化・考察
結果を確認・確認するために、学習結果を取得します。
weibull_para <-
weibull_fit %>%
ggs() %>%
group_by(Parameter) %>%
summarise(
q05 = quantile(value, probs = 0.05),
q25 = quantile(value, probs = 0.25),
q50 = quantile(value, probs = 0.50),
q75 = quantile(value, probs = 0.75),
q95 = quantile(value, probs = 0.95)
) %>%
mutate(
parameters = str_remove(Parameter, "\\[[0-9]+\\]"),
parameter_num = str_extract(Parameter, "[0-9]+") %>% as.integer()
) %>%
select(-Parameter) %>%
select(parameters, parameter_num, everything())
この結果を可視化すると以下のようになります。
推定値を黒で表示、その推定範囲をグレーの帯で表示、真の値を赤い点線で表示しています。
考察
結果を見ながら少し考察してみます。
- $\eta$: イベント発生までの時間に影響するパラメーター
- 凹状態になっているところは再現できていそうです。
- 一方で、t=17あたりからズレが大きくなり、実際の値より低く見積もられています。
- これは打ち切りの影響が多いと考えられます。
- 打ち切りは「これまではイベントが発生しなかった、これ以降にイベントが発生する」を推定に活用しています。
- ただ、打ち切りタイミングが早いと情報量が少なくなるので、真の値よりも低めに推定されてしまいます。
- m: 分布の形状に影響するパラメーター
- 右肩下がりという形状は再現できていそうです。
- 一方でt=20あたりから乖離が大きくなっています。
- おそらくこれは、$\eta$の推定値がズレたしわ寄せがmに来ているのだと思われます。
Weibull 分布の平均値は
$$
\mu = \eta \Gamma(1 + \frac{1}{m})
$$
で計算されます。
これで最後の月の Sales cycle の平均値を計算すると、以下のようになります。
| 真の値 | 時変Weibull分布での推定値 | 結果が出ていないものは月末に結果が出たと仮定して平均 |
|---|---|---|
| 108日 | 79日 | 11日 |
となり、真の値から少し外れているものの、最初のほうで紹介した簡単な計算方法よりは精度が高くなっています。