【和訳の一部】Buy ’Til You Die - A Walkthrough
UPDATE: 2022-10-28 20:22:31
はじめに
LTVの話題でよく出てくるDaniel McCarthy, Edward Wadsworth (2014)Buy ’Til You Die - A Walkthroughの一部分(1章~2章)の和訳。理論は下記参照。
- Counting Your Customers: Who Are They and What Will They Do Next?
- 将来発生するトランザクション数を予測する方法
- ID-POSデータから、顧客の将来価値を測定する
- Rを用いたLTV(Life Time Value)の推定
- Introduction how to calculate LTV
- pLTV: 顧客生涯価値を予測する
1 イントロダクション
BTYDパッケージには、顧客の非契約購買行動(non-contractual purchasing behavior)を捉えるモデル、より簡単に言えば、「顧客が死ぬまで購買する(=顧客として活動しなくなるまで購買する)」という顧客のストーリーを捉えるモデルが含まれています。パッケージで提示されている主なモデルは、Pareto/NBD、BG/NBD、BG/BBモデルであり、これらは、顧客がドロップアウトするまでの正確な時間を観察できない顧客の購買シナリオを説明するモデルを提示します。順番にそれぞれをカバーします。これらのモデルに慣れていない場合、Fader et al. (2004) はBG/NBDモデルの説明を提供し、Fader et al. (2005) はPareto/NBDモデルの説明を提供し、Fader et al. (2010)はBG/BBモデルの説明を提供しているのでそちらを参照してください。
- Fader, Peter S., and Bruce G.S. Hardie. “A Note on Deriving the Pareto/NBD Model and Related Expressions.” November. 2005. Web. http://www.brucehardie.com/notes/008/
- Fader, Peter S., Bruce G.S. Hardie, and Jen Shang. “Customer-Base Analysis in a Discrete-Time Noncontractual Setting.” Marketing Science, 29(6), pp. 1086-1108. 2010. INFORMS. http://www.brucehardie.com/papers/020/
- Fader, Peter S., Hardie, Bruce G.S., and Lee, Ka Lok. ““Counting Your Cus- tomers” the Easy Way: An Alternative to the Pareto/NBD Model.” Marketing Science, 24(2), pp. 275-284. 2005. INFORMS. http://brucehardie.com/papers/018/
2 Pareto/NBD
Pareto/NBDモデルは、顧客がいつでも購入できる非契約の状況(non-contractual situations)で使用されます。 4つのパラメタを使用して、顧客の購入確率と、顧客のドロップアウト確率を説明します。CDNOWデータセットを使用して、BTYDパッケージが提供するPareto/NBDの機能を見ていきましょう。図1に示すように、Pareto/NBDモデルはこのデータセットを非常によく記述しています。
2.1 データの準備
Pareto/NBDモデルのパラメタを推定するのに必要なデータは驚くほど少ないです。モデルの顧客ごとのアプローチは維持されますが、必要な情報はすべての顧客について3つだけです。キャリブレーション期間に行った取引の回数(frequency)、最後の取引の時間(recency)、観察された合計時間です。BTYDパッケージで使用される顧客別十分統計量行列(customer-by-sufficient-statistic matrix)は、すべての顧客(行数)と上述の統計量(列数)のを持つ行列です。
イベントログとして利用可能なデータがあります。dc.ReadLines()は、カンマ区切りのファイル内のイベントログをRのデータフレ
ームに変換するための関数です。以下の例では、2列目に顧客ID(cust)、3列目に日付(date)、
5
列目に売上(sales)を持つファイルcdnowElog.csvからイベントログを作成しています。
library(BTYD)
cdnowElog <- system.file("data/cdnowElog.csv", package = "BTYD")
elog <- dc.ReadLines(
cdnowElog,
cust.idx = 2,
date.idx = 3,
sales.idx = 5
)
elog[1:3, ]
cust date sales
1 1 19970101 29.33
2 1 19970118 29.73
3 1 19970802 14.96dc.Readlines()は、カンマ区切りの元のファイルに表示されている通りに、
日付を文字として保存します。しかし、BTYDのデータ変換関数(data-conversion
functions)の多くでは、日付は互いに比較する必要があります。そのため、イベントログの日付をDate型のオブジェクトに変換します。
elog$date <- as.Date(elog$date, "%Y%m%d")
elog[1:3,]
cust date sales
1 1 1997-01-01 29.33
2 1 1997-01-18 29.73
3 1 1997-08-02 14.96イベントログの日付は正しいフォーマットになりましたが、もう少し整理する必要があります。Pareto/NBDのようなトランザクション・フロー・モデルは、
購入間の時間(interpurchase
time)に関係しています。取引タイミングの日までしか正確ではないので、同じ日に発生したすべてのトランザクションをマージする必要があります。このために、dc.MergeTransactionsOnSameDate()を使用します。この関数は、1日に1人の顧客につき1つのトランザクションのみのイベントロ
グを返し、その日の消費額の合計が売上になるように変換します。
elog <- dc.MergeTransactionsOnSameDate(elog)
Started merging same-date transactions...
... Finished merging same-date transactions.モデルが機能することを検証するために、データをキャリブレーション期間とホールドアウト期間に分割する必要があります。これは、イベントログまたは顧客別時間マトリックスのいずれかで比較的簡単にでき、すぐに作成する予定です。この時点(39週間)でデータセットが半分に分割されるので、1997年9月30日をカットオフ日とします。今、このように分割する理由は、顧客別時間行列(customer-by-time matrix)から顧客別十分統計量行列(customer-by-sufficient-statistic matrix)を作成するときに明らかになるでしょう。
end.of.cal.period <- as.Date("1997-09-30")
elog.cal <- elog[which(elog$date <= end.of.cal.period), ]最後のクリーンアップステップは非常に重要なステップです。キャリブレーション期間中、Pareto/NBDモデルは一般的にリピート取引、つまり最初の取引は無視されます。これは、少なくとも1回の取引を行ったすべての顧客を追跡することが容易であるため、実際にこのモデルを使用している企業にとって便利です(全く取引を行っていない顧客を考慮しようとするのとは対照的)。顧客の最初の取引を単純に取り除くことの1つの問題は以下の通りです。直近性(recency)および観測された総時間(total
time observed)の参照点として時間ゼロ(time
zero)を追跡しなければならない。このため、各顧客に関する重要な情報($cust.data)保存するだけでなく、フィルタリングされたイベントログ($repeat.trans.elog)を返すdc.SplitUpElogForRepeatTrans()を使用しています。
split.data <- dc.SplitUpElogForRepeatTrans(elog.cal)
Started Creating Repeat Purchases
Finished Creating Repeat Purchases
clean.elog <- split.data$repeat.trans.elog(補足START)
例えば、顧客ID=1は、本来①1997/01/01、②1997/01/18、③1997/08/02、④1997/12/12の4回の取引があるが、①は1回目なので除外され、④はカットオフで対象外なので、$repeat.trans.elogには②と③しかない。顧客ID=37は、本来①1997/01/02の1回のみの取引であるため、①は1回目なので除外されて、$repeat.trans.elogに顧客ID=37の情報はない。$cust.dataには顧客ごとのID(cuid)、初回取引日(birth.per)、初回取引合計金額(first.sales)、最終取引日(last.date)、最終取引合計金額(last.sales)が記録される。
split.data$repeat.trans.elog
cust date sales
416 1 1997-01-18 29.73
4364 1 1997-08-02 14.96
281 2 1997-01-13 11.77
251 6 1997-01-11 32.99
2648 6 1997-03-15 77.96
3374 6 1997-04-16 59.30
3455 6 1997-04-24 134.98
4005 6 1997-06-23 91.92
4242 6 1997-07-22 47.08
4285 6 1997-07-26 71.96
1039 7 1997-02-05 11.77
3529 9 1997-05-01 28.13
4644 9 1997-09-08 22.97
【略】
2098 21 1997-03-02 34.60
3345 21 1997-04-14 29.99
282 26 1997-01-13 227.14
4297 28 1997-07-27 80.46
1354 30 1997-02-12 40.47
4174 35 1997-07-11 31.14
407 38 1997-01-17 13.97
4000 45 1997-06-23 14.96
338 46 1997-01-14 14.96
split.data$cust.data
cust birth.per first.sales last.date last.sales
1 1 1997-01-01 29.33 1997-08-02 14.96
11 2 1997-01-01 63.34 1997-01-13 11.77
12 3 1997-01-01 6.79 1997-01-01 6.79
13 4 1997-01-01 13.97 1997-01-01 13.97
【略】
36 36 1997-01-02 14.96 1997-01-02 14.96
37 37 1997-01-02 34.97 1997-01-02 34.97
38 38 1997-01-02 47.00 1997-01-17 13.97(補足END)
次のステップは、顧客別時間行列(customer-by-time matrix)を作成することです。これは単純に、各顧客ごとに各日付の列を持つ行列です。これらの行列を作成するには、いくつかの異なるオプションがあります。
- Frequency:
行列は、その日にその顧客によって行われた取引の数が含まれます。
dc.CreateFreqCBT()を使用します。すでにdc.MergeTransactionsOnSameDate()を使用している場合、これは単に顧客ごとの時間リーチ行列になります。 - Reach:行列は、顧客がその日に取引を行った場合は1を含み、そうでない場合は0を含みます。
dc.CreateReachCBT()を使用します。 - Spend:行列は、その日にその顧客が支出した金額が含まれます。
dc.CreateSpendCBT()を使用します。is.avg.spendパラメタを変更することで、各日の消費額を使用するか、各日の平均消費額を使用するかを設定することができます。ほとんどの場合、is.avg.spendをFALSEのままにしておくことが適切です。
freq.cbt <- dc.CreateFreqCBT(clean.elog)
freq.cbt[1:10,1:10]
date
cust 1997-01-08 1997-01-09 1997-01-10 1997-01-11 1997-01-12 1997-01-13 1997-01-14 1997-01-15 1997-01-16 1997-01-17
1 0 0 0 0 0 0 0 0 0 0
2 0 0 0 0 0 1 0 0 0 0
6 0 0 0 1 0 0 0 0 0 0
7 0 0 0 0 0 0 0 0 0 0
9 0 0 0 0 0 0 0 0 0 0
11 0 0 0 0 0 0 0 0 0 0
17 0 0 0 0 0 0 0 0 0 0
18 0 0 0 0 0 0 0 0 0 0
19 0 0 0 0 0 0 0 0 0 0
21 0 0 0 0 0 0 0 0 0 0上のアウトプットから注意すべき点が2つあります。
- 顧客ID3、4および5がいないように見えますが、実際、いません。いないというのは、これらの顧客はリピート取引を行っておらず、
dc.SplitUpElogForRepeatTrans()を使用した際に削除されました。dc.MergeCustomersを使用した際にデータに戻されますので、心配しないでください。
- 顧客ID3、4および5がいないように見えますが、実際、いません。いないというのは、これらの顧客はリピート取引を行っておらず、
- コラムは1月8日から始まっていますが、これは最初の取引をすべて削除したためです(最初の週の間に2つの取引をした人はいません)。
少し問題があります。最初のトランザクションをすべて削除したので、顧客別時間行列(customer-by-time)には、
リピートがない顧客が含まれていません。実際、ほとんどのデータセットでは、他のどの数字よりも多くの顧客がリピート取引をゼロにしています。この問題の解決は非常に簡単です。すべてのトランザクションを使用して顧客別時間行列(customer-by-time)を作成し、フィルタリングされたCBTをこの合計CBTにマージします(フィルタリングされたCBTからのデータと合計CBTからの顧客IDを使用します)。
tot.cbt <- dc.CreateFreqCBT(elog)
cal.cbt <- dc.MergeCustomers(tot.cbt, freq.cbt)(補足START)
要するにリピートしていない顧客をfreq.cbtのデータを基準(取引日)に戻しているようです。
dim(freq.cbt)
[1] 946 266
dim(cal.cbt)
[1] 2357 266
tot.cbt[1:10,1:10]
date
cust 1997-01-01 1997-01-02 1997-01-03 1997-01-04 1997-01-05 1997-01-06 1997-01-07 1997-01-08 1997-01-09 1997-01-10
1 1 0 0 0 0 0 0 0 0 0
2 1 0 0 0 0 0 0 0 0 0
3 1 0 0 0 0 0 0 0 0 0
4 1 0 0 0 0 0 0 0 0 0
5 1 0 0 0 0 0 0 0 0 0
6 1 0 0 0 0 0 0 0 0 0
7 1 0 0 0 0 0 0 0 0 0
8 1 0 0 0 0 0 0 0 0 0
9 1 0 0 0 0 0 0 0 0 0
10 1 0 0 0 0 0 0 0 0 0
cal.cbt[1:10,1:10]
1997-01-08 1997-01-09 1997-01-10 1997-01-11 1997-01-12 1997-01-13 1997-01-14 1997-01-15 1997-01-16 1997-01-17
1 0 0 0 0 0 0 0 0 0 0
2 0 0 0 0 0 1 0 0 0 0
3 0 0 0 0 0 0 0 0 0 0
4 0 0 0 0 0 0 0 0 0 0
5 0 0 0 0 0 0 0 0 0 0
6 0 0 0 1 0 0 0 0 0 0
7 0 0 0 0 0 0 0 0 0 0
8 0 0 0 0 0 0 0 0 0 0
9 0 0 0 0 0 0 0 0 0 0
10 0 0 0 0 0 0 0 0 0 0(補足END)
キャリブレーション期間の顧客別時間行列(および先ほど保存した追加情報)から、最終的に先ほど説明した顧客別十分統計量行列(customer-by-sufficient-statistic
matrix)を作成することができます。使用する関数はdc.BuildCBSFromCBTAndDates()であり、これには、顧客別時間行列(customer-by-time
matrix)、
各顧客の開始日と終了日、およびキャリブレーション期間の終了時刻が必要です。また、終了期間も必要となります。例えば、日数(par = "day")を選択した場合、顧客別十分統計量行列(customer-by-sufficient-statistic
matrix)の直近性(recency)と観測された合計時間(total time
observed)の列は、39ではなく273になりますが、最終的には同じデータになります。この関数は、cbt.is.intern.cal.periodをFALSEに設定することで、ホールドアウト期間にも使用することができます。この関数をホールドアウト期間に使用した場合には若干の違いがありますが、異なる日付(単純にホールドアウト期間の開始日と終了日)を入力する必要があり、直近性(recency,
ホールドアウト期間ではほとんど価値がありません)を返しません。
# 顧客ごとの初回取引日ベクトル
birth.periods <- split.data$cust.data$birth.per
# 顧客ごとの最終取引日ベクトル
last.dates <- split.data$cust.data$last.date
# 顧客ごとの初回、最終取引データ
cal.cbs.dates <- data.frame(birth.periods, last.dates, end.of.cal.period)
cal.cbs <- dc.BuildCBSFromCBTAndDates(cal.cbt, cal.cbs.dates, per="week")
Started making calibration period CBS...
Finished building CBS.(補足START)
観測された合計時間の273というのは、データの開始日からカットオフまでの日数差のことで、起算日を含む計算なのでプラス1のことを言っているのかよくわからんが、実際の計算結果は38.85で272になる。私の英語能力の低さが原因か・・・？
as.Date("1997-09-30") - as.Date("1997-01-01")
Time difference of 272 days
head(cal.cbs.dates, 10)
birth.periods last.dates end.of.cal.period
1 1997-01-01 1997-08-02 1997-09-30
2 1997-01-01 1997-01-13 1997-09-30
3 1997-01-01 1997-01-01 1997-09-30
4 1997-01-01 1997-01-01 1997-09-30
5 1997-01-01 1997-01-01 1997-09-30
6 1997-01-01 1997-07-26 1997-09-30
7 1997-01-01 1997-02-05 1997-09-30
8 1997-01-01 1997-01-01 1997-09-30
9 1997-01-01 1997-09-08 1997-09-30
10 1997-01-01 1997-01-01 1997-09-30cal.cbsには、顧客別十分統計量行列(customer-by-sufficient-statistic
matrix)が入っており、
# week計算
head(cal.cbs, 10)
x t.x T.cal
1 2 30.428571 38.85714
2 1 1.714286 38.85714
3 0 0.000000 38.85714
4 0 0.000000 38.85714
5 0 0.000000 38.85714
6 7 29.428571 38.85714
7 1 5.000000 38.85714
8 0 0.000000 38.85714
9 2 35.714286 38.85714
10 0 0.000000 38.85714
# day計算
tmp <- dc.BuildCBSFromCBTAndDates(cal.cbt, cal.cbs.dates, per="day")
head(tmp, 10)
x t.x T.cal
1 2 213 272
2 1 12 272
3 0 0 272
4 0 0 272
5 0 0 272
6 7 206 272
7 1 35 272
8 0 0 272
9 2 250 272
10 0 0 272顧客別十分統計量行列(customer-by-sufficient-statistic matrix)の解釈として、1番目の顧客を例にすると、
xはFrequencyなので2回の取引したt.xはRecencyなので1997年01月01日基準で213日に購入した(=つまり、直近59日前に取引した)。これが大きいと直近に取引したことになる。小さいと直近は取引してないことになる。T.calは、観測された合計なので、1997年01月01日から1997年09月30日までの272日間が合計時間
a <- as.Date("1997-09-30") - as.Date("1997-01-01")
a
Time difference of 272 days
b <- as.Date("1997-09-30") - as.Date("1997-08-02")
b
Time difference of 59 days
T.cal <- a - b
T.cal
Time difference of 213 days(補足END)
上記の処理について、パッケージにはdc.ElogToCbsCbt()という、すべてを行う単一の関数が含まれていることを聞いて喜ぶことでしょ
う。しかし、状況によっては、顧客の初回トランザクションを削除したくない場合や、データがイベントログではなく顧客ごとのタイムマトリクスとして利用可能な場合など、小さな変更を加えたい場合があるかもしれないので、上記のプロセスを意識しておくと良いでしょう。しかし、ほとんどの標準的な状況では、dc.ElogToCbsCbt()で十分です。この関数のパラメタと出力については、パッケージのドキュメントを読むことで、ほとんどの目的でこの関数を使用するのに十分な理解を得ることができます。
2.2 パラメタ推定
これでデータが必要な形式になったので、モデルのパラメタを推定することができます。パラメタを推定するには、pnbd.EstimateParameters()を使用します。
これには、キャリブレーション期間の顧客別十分統計量行列(customer-by-sufficient-statistic
matrix)と、(オプションで)
初期値パラメタが必要です。(1,1,1,1,1)がデフォルトの初期値パラメタとして使用されます。対数尤度を最大化する関数はpnbd.cbs.LLであり、これは顧客別十分統計量行列(customer-by-sufficient-statistic
matrix)に対して与えられたパラメタのセットに対する対数尤度を返します。
params <- pnbd.EstimateParameters(cal.cbs)
params
[1] 0.5533971 10.5801985 0.6060625 11.6562237
LL <- pnbd.cbs.LL(params, cal.cbs)
LL
[1] -9594.976どんな最適化でもそうですが、最初に得られる結果に満足してはいけません。 最初の結果を出発点にして、さらに数回実行してみて、収束するかどうかを確認してみましょう。
p.matrix <- c(params, LL)
for (i in 1:2) {
params <- pnbd.EstimateParameters(cal.cbs, params)
LL <- pnbd.cbs.LL(params, cal.cbs)
p.matrix.row <- c(params, LL)
p.matrix <- rbind(p.matrix, p.matrix.row)
}
colnames(p.matrix) <- c("r" , "alpha" , "s" , "beta" , "LL")
rownames(p.matrix) <- 1:3.
p.matrix
r alpha s beta LL
1 0.5533971 10.58020 0.6060625 11.65622 -9594.976
2 0.5534354 10.57952 0.6060276 11.65666 -9594.976
3 0.5533710 10.57952 0.6059074 11.65812 -9594.976この最適化の推定は収束します。ここでは実証するつもりはありませんが、いくつか初期値を変更してから推定関数をテストするのは常に良いアイデアです。パラメタがわかったので、BTYDパッケージはそれらを解釈する関数を提供しています。ご存知のように、GammaとAlphaは、NBD取引過程のガンマ混合分布を表します。このガンマ分布は図2で見ることができ、pnbd.PlotTransactionRateHeterogeneity(params)を使ってプロットしています。また、sとBetaはPareto(またはガンマ指数関数)ドロップアウト過程のガンマ混合分布を説明するものです。このガンマ分布は図3で見ることができ、pnbd.PlotDropoutHeterogeneity(params)を使ってプロットしています。これらから、我々の顧客についてわかります。彼らは個々のポアソン取引過程のパラメタの値が低い傾向にあり、個々の指数関数的ドロップアウトプロセスのパラメタの値が低い傾向にあります。
2.3 個人レベルの推定
母集団のパラメタがわかったので、個人レベルで顧客の推定を行うことができます。まず、新規に獲得した顧客がある期間に行うと予想される取引の回数を推定することができます。例えば、新規に獲得した顧客が1年の期間に行うリピート取引の数に興味があるとします。12ヶ月、365日、1年ではなく、1年を表すために52 週を使用していることに注意してください。これは、パラメタが週単位のデータを使用して推定されたからです。
pnbd.Expectation(params, t=52)
[1] 1.473434(補足START) 新規顧客が52週間で取引する期待購入回数は1.47回くらいというニュアンス。上記の取引期間で、普通に計算すると2357人中946人がリピートしているので、リピート率は40%となる。
946/2357
[1] 0.4013577(補足END)
また、キャリブレーション期間中の顧客の購買行動に基づいて、特定の顧客の期待値を計算することもできます。pnbd.ConditionalExpectedTransactionsで、ホールドアウト期間中に顧客が行うと期待される取引数を計算し、2つ目はpnbd.PAlive()で、キャリブレーション期間の終了時に顧客がまだ生きている確率を表します。上記のように、使用される期間は、どの期間がパラメタの推定に使用されたかによって異なります。
# `cal.cbs`には、顧客別十分統計量行列(customer-by-sufficient-statistic matrix)が入っている。例、顧客ID=1516
cal.cbs["1516",]
x t.x T.cal
26.00000 30.85714 31.00000
# x(Frequency)、t.x(Recency)、T.cal(total observed time)
x <- cal.cbs["1516", "x"]
t.x <- cal.cbs["1516", "t.x"]
T.cal <- cal.cbs["1516", "T.cal"]
# 以後52週間の期待取引回数
pnbd.ConditionalExpectedTransactions(params, T.star = 52,x, t.x, T.cal)
[1] 25.45647
pnbd.PAlive(params, x, t.x, T.cal)
[1] 0.997874(補足START)
実際の顧客ID=1516のデータでは、キャリブレーション期間(9ヶ月)で27回の取引、ホールドアウト期間(9ヶ月だけど予測は12ヶ月)で15回の取引なので、条件を揃えて35週にすると、18回の取引が期待される結果が返る。ちょっと多めの取引回数が返されている。
pnbd.ConditionalExpectedTransactions(params, T.star = 35,x, t.x, T.cal)
[1] 18.3693
elog %>%
dplyr::filter(cust == 1516) %>%
dplyr::group_by(cust, date) %>%
dplyr::summarise(sales = sum(sales)) %>%
dplyr::ungroup() %>%
dplyr::mutate(point = if_else(date >= "1997-09-30", "holdout", "calibration")) %>%
dplyr::group_by(point) %>%
dplyr::summarise(sum_salse = sum(sales),
cnt_date = n()) %>%
dplyr::ungroup()
# A tibble: 2 x 3
point sum_salse cnt_date
<chr> <dbl> <int>
1 calibration 1051. 27
2 holdout 429. 15現在のアクティブ確率は99%。キャリブレーション期間の終了時に顧客がまだ生きている確率なので、1997-09-30がキャリブレーション終了で、1997-09-29に取引があるので、アクティブ率も高い。
# A tibble: 42 x 4
cust date sales point
<chr> <date> <dbl> <chr>
1 1516 1997-02-25 11.8 calibration
2 1516 1997-02-28 65.8 calibration
3 1516 1997-03-29 58.2 calibration
4 1516 1997-06-09 28.3 calibration
5 1516 1997-06-11 25.1 calibration
6 1516 1997-07-02 11.8 calibration
7 1516 1997-07-05 10.8 calibration
8 1516 1997-07-11 24.1 calibration
9 1516 1997-07-21 35.3 calibration
10 1516 1997-07-22 57.2 calibration
11 1516 1997-07-24 193. calibration
12 1516 1997-07-26 23.5 calibration
13 1516 1997-07-30 88.0 calibration
14 1516 1997-08-02 43.3 calibration
15 1516 1997-08-08 75.4 calibration
16 1516 1997-08-10 25.7 calibration
17 1516 1997-08-12 26.1 calibration
18 1516 1997-08-18 28.3 calibration
19 1516 1997-08-20 29.7 calibration
20 1516 1997-08-22 9.97 calibration
21 1516 1997-09-06 30.3 calibration
22 1516 1997-09-10 11.5 calibration
23 1516 1997-09-14 26.5 calibration
24 1516 1997-09-16 21.5 calibration
25 1516 1997-09-17 32.5 calibration
26 1516 1997-09-21 28.5 calibration
27 1516 1997-09-29 29.0 calibration
28 1516 1997-10-15 11.5 holdout
29 1516 1997-10-21 15.0 holdout
30 1516 1997-10-24 38.5 holdout
31 1516 1997-10-30 18.0 holdout
32 1516 1997-11-07 14.5 holdout
33 1516 1997-11-15 9.49 holdout
34 1516 1997-12-13 51.0 holdout
35 1516 1997-12-21 13.0 holdout
36 1516 1998-01-21 23.0 holdout
37 1516 1998-01-30 36.0 holdout
38 1516 1998-03-02 38.0 holdout
39 1516 1998-04-26 14.4 holdout
40 1516 1998-05-10 27.0 holdout
41 1516 1998-05-26 92.9 holdout
42 1516 1998-06-10 27.0 holdout
(補足END)
もう一つ注意すべき点があります。条件付き期待値関数を使用すると、「増加する頻度のパラドックス(increasing frequency paradox)」の作用を見ることができます。
for (i in seq(10, 25, 5)) {
cond.expectation <- pnbd.ConditionalExpectedTransactions(
params,
T.star = 52,
x = i,
t.x = 20,
T.cal = 39
)
cat ("x:", i, "\t Expectation:", cond.expectation, fill = TRUE)
}
x: 10 Expectation: 0.7062289
x: 15 Expectation: 0.1442396
x: 20 Expectation: 0.02250658
x: 25 Expectation: 0.003092672.4 当てはまりの良さ
私たちは、個々の顧客について推論を行う以上のことができるようになりた いと考えています。BTYDパッケージは、キャリブレーション期間とホールド アウト期間の両方で、期待される顧客の行動を実際の顧客の行動に対してプロッ トする機能を提供します。キャリブレーション期間内の実際の購買回数と期待購買回数の比較です。これは、以下のコードを使用して生成された図です。
pnbd.PlotFrequencyInCalibration(params, cal.cbs, 7)
この関数は明らかに(推定パラメタから)期待されるデータを生成できる必要があり、実際のデータ(キャリブレーション期間の顧客別統計量)を必要とします。また、しきい値(censor
number)と呼ばれる別の数字も必要です。プロットされるヒストグラムは右打ち切りになっており、ある数を超えると、すべての値が結合されます。しきい値として与えられた数字は、データがどこで切断されるかを決定します。
残念ながら、キャリブレーション期間の値を比較することでわかることは、我々のモデルとデータの間の適合は悪くないということだけです。モデルの適合性がホー
ルドアウト期間まで維持されていることを検証する必要があります。dc.ElogToCbsCbt ()は、キャリブレーション期間の顧客別十分統計量行列とホールドアウト期間の顧客別十分統計量行列の両方を生成し、これらを組み合わせて、ホールドアウト期間に各顧客が行った取引の数を求めることができます。しかし、dc.ElogToCbsCbt()を使用していないので、イベントログから直接情報を取得することにします。なお、総取引数からキャ
リブレーション期間中のリピートトランザクション数を差し引いています。初期トラ
ンザクションは関係ないので最初に削除しています。
elog <- dc.SplitUpElogForRepeatTrans(elog)$repeat.trans.elog
x.star <- rep(0, nrow(cal.cbs))
cal.cbs <- cbind(cal.cbs, x.star)
elog.custs <- elog$cust
for (i in 1:nrow(cal.cbs)) {
current.cust <- rownames(cal.cbs)[i]
tot.cust.trans <- length(which(elog.custs == current.cust))
cal.trans <-
cal.cbs[i, "x"]
cal.cbs[i, "x.star"] <- tot.cust.trans - cal.trans
}
cal.cbs[1:10, ]
x t.x T.cal x.star
1 2 30.428571 38.85714 1
2 1 1.714286 38.85714 0
3 0 0.000000 38.85714 0
4 0 0.000000 38.85714 0
5 0 0.000000 38.85714 0
6 7 29.428571 38.85714 8
7 1 5.000000 38.85714 0
8 0 0.000000 38.85714 2
9 2 35.714286 38.85714 2
10 0 0.000000 38.85714 0ここで、モデルがホールドアウト期間にどれだけうまく機能するかを見ることができます。図4は、以下のコードによって生成された出力を示 しています。これは、顧客をキャリブレーション期間の頻度に従ってビンに分割し、 これらのビンについて、実際のホールドアウト期間と条件付きで予想されるホ ールドアウト期間の頻度をプロットしています。
T.star <- 39 # length of the holdout period
censor <- 7 # This censor serves the same purpose described above
x.star <- cal.cbs[,"x.star"]
comp <- pnbd.PlotFreqVsConditionalExpectedFrequency(params, T.star,cal.cbs, x.star, censor)
rownames(comp) <- c("act", "exp", "bin")
comp
freq.0 freq.1 freq.2 freq.3 freq.4 freq.5 freq.6 freq.7+
transaction.actual 0.2367116 0.6970387 1.392523 1.560000 2.532258 2.947368 3.862069 6.359375
transaction.expected 0.1384724 0.5995607 1.195989 1.714041 2.398545 2.907467 3.818906 6.403484
bin.size 1411.0000000 439.0000000 214.000000 100.000000 62.000000 38.000000 29.000000 64.000000
上で見たように、グラフと行列出力も行います。BTYDパッケージのほとんどのプロット関数は、このような出力を生成します。これらの関数には、グラフには表示されていない追加の情報、つまり、グラフ内の各ビンのサイズが含まれているので、見る価値があります。例えば、このグラフでは、ビンのサイズが、ゼロでのギャップが6または7トランザクションでの精度よりもはるかに大 きな意味を持つことを示しているため、この情報は重要です。このグラフは、モデルがホールドアウト期間のデータに非常によく適合していることを示しています。キャリブレーション期間の頻度による集計は、それを行うための1つの方法にすぎません。 BTYD は、他のいくつかの尺度で集計するプロット関数も提供しています。ここで 実証するもう一つの方法は、時間による集計です。最初のステップは、もう一度、モデルを比較するために必要なデータを収集することです。顧客別時間マトリックスはすでに時間帯別のデータを収集しています。そして、毎日のトラッキングデータを週単位のデータに変換します。
tot.cbt <- dc.CreateFreqCBT(elog)
d.track.data <- rep(0, 7 * 78)
origin <- as.Date("1997-01-01")
for (i in colnames(tot.cbt)) {
date.index <- difftime(as.Date(i), origin) + 1
d.track.data[date.index] <- sum(tot.cbt[, i])
}
w.track.data <- rep(0, 78)
for (j in 1:78) {
w.track.data[j] <- sum(d.track.data[(j * 7 - 6):(j * 7)])
}
さて、実際の取引件数と週次の予想取引件数を比較したプロットを作成すると、図5のようになります。
n.periods.finalを78に設定していることに注意してください。これは、
週単位のデータを使用していることを示します。もしトラッキングデータが日次であれば、ここでは546を使うことになります。この概念は少し難しいかもしれませんが、pnbd.PlotTrackingInc()のドキュメントで説明されています。
合計期間に2つの数字(T.totとn.periods.final)があるのは、顧客統計量行列のマトリックスとトラッキングデータは異なる時間のポインがあるかもしれません。
T.cal <- cal.cbs[,"T.cal"]
T.tot <- 78
n.periods.final <- 78
inc.tracking <- pnbd.PlotTrackingInc(params, T.cal,T.tot, w.track.data, n.periods.final)
inc.tracking[,20:25]
[,1] [,2] [,3] [,4] [,5] [,6]
actual 73.00000 55.0000 70.00000 33.00000 56.00000 99.00000
expected 78.30848 76.4191 74.64776 72.98278 71.41403 69.93268
図5は、このモデルが時間の経過とともに顧客の購買傾向を確実に捉えている ことを示していますが、非常に乱雑な懐疑論者を納得させることはできないかもしれません。さらに、サンプルが示している行列は、購入はある週から次の週までに変化する可能性があるため、実際には多くの情報を伝えることができません。これらの理由から、図6に示されているように、時間をかけて累積するこ とでデータを平滑化する必要があるかもしれません。
cum.tracking.data <- cumsum(w.track.data)
cum.tracking <- pnbd.PlotTrackingCum(params, T.cal,T.tot, cum.tracking.data, n.periods.final)
cum.tracking
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
actual 0.000000 19.0000 42.00000 81.00000 119.0000 192.0000 261.0000 351.0000 428.0000 504.0000 610.0000 733.0000 828.000 914.0000 1005.000
expected 4.215215 16.5693 37.21379 66.71984 105.3896 153.7293 210.3993 275.4081 349.0571 431.0633 520.3901 616.0048 712.358 805.4379 895.532
[,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30]
actual 1078.0000 1149.00 1222.000 1286.000 1359.000 1414.000 1484.000 1517.000 1573.00 1672.000 1746.000 1800.000 1840.000 1885.000 1973.000
expected 982.8877 1067.72 1150.218 1230.548 1308.856 1385.275 1459.923 1532.906 1604.32 1674.253 1742.784 1809.986 1875.926 1940.666 2004.262
[,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45]
actual 2032.000 2069.00 2132.000 2192.000 2239.0 2290.000 2342.000 2395.000 2457.000 2507.000 2563.000 2612.000 2671.000 2718.000 2803.000
expected 2066.767 2128.23 2188.695 2248.206 2306.8 2364.515 2421.384 2477.442 2532.716 2587.237 2641.031 2694.122 2746.536 2798.294 2849.417
[,46] [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60]
actual 2864.000 2923.000 2968.000 3038.000 3104.0 3153.000 3183.000 3234.000 3276.000 3329.00 3357.000 3404.000 3444.000 3490.000 3533.000
expected 2899.928 2949.843 2999.182 3047.963 3096.2 3143.912 3191.111 3237.814 3284.032 3329.78 3375.071 3419.915 3464.324 3508.311 3551.884
[,61] [,62] [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74] [,75]
actual 3598.000 3677.000 3727.000 3786.000 3844.000 3879.000 3915.000 3962.000 3991.000 4035.000 4077.000 4112.000 4141.000 4183.000 4231.000
expected 3595.054 3637.832 3680.225 3722.244 3763.897 3805.193 3846.139 3886.743 3927.012 3966.954 4006.577 4045.885 4084.887 4123.588 4161.995
[,76] [,77] [,78]
actual 4279.000 4311.000 4339.000
expected 4200.112 4237.947 4275.504
##
コード
# ---- library ----
library(tidyverse)
library(lubridate)
library(BTYD)
# ---- Parate/NBD, Gamma-Gamma ----
# 読み込み
elog <- dc.ReadLines(system.file("data/cdnowElog.csv", package="BTYD"), cust.idx = 2, date.idx = 3, sales.idx = 5) %>%
dplyr::mutate(date = ymd(date)) %>%
dplyr::group_by(cast, date) %>%
dplyr::summarise(sales = sum(sales)) # bdc.MergeTransactionsOnSameDate(elog)
# キャリブレーション期間の設定
end_of_cal_period <- as.Date("1997-09-30")
elog_cal <- elog %>%
dplyr::filter(end_of_cal_period >= date)
# リピート取引のみに限定
split_data <- dc.SplitUpElogForRepeatTrans(elog.cal)
clean_elog <- split_data$repeat.trans.elog
# 顧客別時間行列(customer-by-time matrix)。取引があるとフラグ1が立つ。
# 対象の顧客の取引日データをベースに1回しか取引がない顧客のデータをマージ
freq_cbt <- dc.CreateFreqCBT(clean_elog)
tot_cbt <- dc.CreateFreqCBT(elog)
cal_cbt <- dc.MergeCustomers(tot_cbt, freq_cbt)
# 顧客の情報をまとめる
birth_periods <- split_data$cust.data$birth.per
last_dates <- split_data$cust.data$last.date
cal_cbs_dates <- data.frame(birth_periods, last_dates, end_of_cal_period)
# キャリブレーション期間の顧客別の平均売上を計算
m.x <- elog_cal %>%
dplyr::group_by(cust) %>%
dplyr::summarise(m.x=mean(sales)) %>%
dplyr::pull(m.x)
# 顧客別時間行列(customer-by-time matrix)と顧客別初回/最終取引行列を使って、
# 顧客別十分統計量行列(customer-by-sufficient-statistic matrix)を作成。それに顧客別の平均売上を結合。
# pnbd.EstimateParameters()がtibbleに対応していないのでas.data.frame()
cal_cbs <- dc.BuildCBSFromCBTAndDates(cal_cbt, cal_cbs_dates, per = "week") %>%
as.data.frame() %>%
dplyr::bind_cols(tibble(m.x = m.x))
cal_cbs %>%
head(10)
# # A tibble: 2,357 x 4
# x t.x T.cal m.x
# <dbl> <dbl> <dbl> <dbl>
# 1 2 30.4 38.9 24.7
# 2 1 1.71 38.9 30.3
# 3 0 0 38.9 21.6
# 4 0 0 38.9 15.0
# 5 0 0 38.9 38.8
# 6 7 29.4 38.9 11.8
# 7 1 5 38.9 63.8
# 8 0 0 38.9 27.9
# 9 2 35.7 38.9 21.6
# 10 0 0 38.9 14.4
# # … with 2,347 more rows
## param est
# palate / nbd
pnbd_params <- pnbd.EstimateParameters(cal_cbs[,1:3])
pnbd_params
# [1] 0.5533971 10.5801985 0.6060625 11.6562237
# gamma - gamma
spend_params <- spend.EstimateParameters(cal_cbs$m.x, cal_cbs$x)
spend_params
# [1] 9.958284e+03 2.809933e+00 6.069101e-03
## 個人レベルの推定
# 上から、キャリブレーション期間の期待取引数、キャリブレーション期間終了時でのアクティブ確率、期待平均金額、期待LTV
df_ltv <-
tibble::tibble(
pnbd_expected_trans = pnbd.ConditionalExpectedTransactions(pnbd_params, T.star = 52, cal_cbs$x, cal_cbs$t.x, cal_cbs$T.cal),
pnbd_prob_alive = pnbd.PAlive(pnbd_params, cal_cbs$x, cal_cbs$t.x, cal_cbs$T.cal),
expected_trans_val = spend.expected.value(spend_params, cal_cbs$m.x, cal_cbs$x)
) %>%
dplyr::mutate(pnbd_forecast_ltv = pnbd_expected_trans * pnbd_prob_alive * expected_trans_val)
df_ltv %>% slice(1516)
# # A tibble: 2,357 x 4
# pnbd_expected_trans pnbd_prob_alive expected_trans_val pnbd_forecast_ltv
# <dbl> <dbl> <dbl> <dbl>
# 1 1.85 0.869 24.7 39.7
# 2 0.218 0.168 30.3 1.11
# 3 0.136 0.295 33.4 1.34
# 4 0.136 0.295 33.4 1.34
# 5 0.136 0.295 33.4 1.34
# 6 4.72 0.749 11.8 41.7
# 7 0.331 0.255 63.8 5.38
# 8 0.136 0.295 33.4 1.34
# 9 2.04 0.959 21.6 42.4
# 10 0.136 0.295 33.4 1.34
# # … with 2,347 more rows