本記事は、ドコモアドベントカレンダー17日目の記事になります。

こんにちは。NTTドコモサービスイノベーション部3年目社員の九島です。
主な業務は、マーケティング分野への機械学習・AI活用です。
今回は、これからのマーケティングではより重要となる指標のLife Time Value(以下、LTV)を機械学習・AIを用いて予測する方法についてご紹介したいと思います。

※本記事ではPythonが準備されていることを前提としています。また、Kaggleのデータセットを利用します。Kaggleについては本記事内でデータの取得方法について解説しています。

この記事を読んでわかること

LTVとは何か
LTVを機械学習・AIで予測する方法と検証方法(コード＆解説付き)
Kaggleデータの利用方法

LTVとは

LTVとは、マーケティングにおいて重要視すべき指標の一つでして、「企業に対してその顧客が生涯で生み出す利益」を表します。
和訳では顧客生涯価値と表現されます。

一般的にLTVは、商品購入における単価×頻度×期間で表されるため、
例えば、高額な商品を高頻度で、長い期間購入し続ける顧客のLTVは、高い数値になります。
このことから企業は、LTVを向上させるために、単価を上げるアップセルや、商品・サービスの利便性を上げてロイヤリティ(満足度)を高めることで、頻度や期間を増やしていく等、様々なアプローチを考えます。

LTVの種類

前章では、一般的なLTV算出式として、単価×頻度×期間をご紹介しましたが、その他にもいくつか種類がありまして、一般的な算出式を含めて以下で計３種類の式をご紹介します。
※基本的に各要素は顧客全体の平均値を用います。

購入単価 × 購入頻度 × 継続期間

購入単価が5,000円、購入頻度が月2回、継続期間が2年間(=24ヶ月)の場合は以下のようにLTVが算出されます。

LTV = 5,000(円) × 2 × 24 = 240,000(円)

最も一般的なLTVの算出式です。

顧客一人あたりの年間取引額 × 収益率 × 顧客一人あたりの継続年数

顧客一人あたりの年間取引額が150万円、収益率が70％、顧客一人あたりの継続年数が6年の場合は以下のようにLTVが算出されます。

LTV = 150万(円) × 0.7 × 6 = 630万(円)

こちらは収益率を考慮できます。

購入単価 ÷ 解約率

購入単価が15,000円、解約率(= 月あたりの解約した顧客数 ÷ 月の顧客全体数)が20％の場合は以下のようにLTVが算出されます。

LTV = 15,000(円) ÷ 0.2 = 75,000(円)

継続期間の割り出しが困難な場合や、サブスクリプション型の販売形態を持つSaaS系企業が活用できます。
このLTVの意味を補足すると、この解約率だと1/0.2=5ヶ月間しか顧客が持たないので、一人あたり15,000の5ヶ月分の75,000円しか収益が見込めない計算になります。

また、上記それぞれの式で算出した値から、顧客獲得コストを示すCost Per Acquisition(CPA)を差し引いた値をLTVと定義する場合もあり、企業ごと、商品・サービスの形態ごとにそれぞれ適したLTV算出式を使用してLTVを測ります。

LTV予測手法

本記事では、LTVを予測する機械学習の手法として、2019年にGoogleから発表された「A Deep Probabilistic Model for Customer Lifetime Value Prediction」に掲載されているモデルを紹介します。
前章までのLTVの説明では、顧客全体の平均値を使用していますが、本モデルでは顧客ごとにLTVを算出・予測します。

概要

本モデルは、Deep Neural Network(以下、DNN)を用いて顧客とLTV(初回購入から一年間の購入合計金額)のデータを学習します。
LTV予測モデリングの課題として、一度の購入で離反してしまう顧客や、LTVの分布がヘビーテールであることを考慮する必要があることから、本モデルではzero-inflated lognormal (以下、ZILN)を新たに導入し、損失関数として活用します。
ZILNは、離反の確率とヘビーテールの性質を両方考慮することができます。
また、ZILNの損失関数は、DNNモデルだけでなく、線形モデルでも使用することができます。
詳細な理論を知りたい方は元論文をご覧ください。

元論文では、二つの公開データセット(Kaggle、KDDCUP)で実験を行い、その有効性を示しています。
本記事ではその内の一つのKaggleのデータセットを使用してLTVを予測する手順を次章以降記載します。

環境構築

実際の実行手順に入る前に、開発環境を準備します。

Pythonの準備

自分はPython3.8を利用しました。

python3 -V

出力結果:
Python 3.8.10

必要なパッケージのインストール

以下のpipコマンドでパッケージをインストールします。

pip install -q git+https://github.com/google/lifetime_value

詳しくは、公式Githubをご覧ください。
https://github.com/google/lifetime_value

Kaggle API

Kaggleのデータセットを使用するためには、KaggleのHPでユーザ登録をし、API Tokenの情報が含まれたjsonファイルを取得する必要があります。
以下の手順で進みます。

ユーザ登録
KaggleHP(https://www.kaggle.com/)の右上にあるRegisterからユーザ登録します。
API Tokenを取得
KaggleHP右上のアイコン->Your Profile->Accountタブ->API項目の順番で進み、Create New API Tokenでkaggle.jsonをダウンロードします。
kaggleAPIを利用するためのパッケージをインストール
pip install kaggle
2番で取得したkaggle.jsonを.kaggleディレクトリへ移動
mv {kaggle.jsonがあるディレクトリ}/kaggle.json ~/.kaggle/
実行権限を付与
chmod 600 ~/.kaggle/kaggle.json
ダウンロードしたいデータの元コンペのRulesをAccept
KaggleHPのCompetitions->コンペを選択->Rulesタブ->I Understand and Acceptをクリックします。
この作業により、コンペデータをダウンロードできるようになります。

実行手順

以降で、LTV予測の実行手順をコードベースで記載し、それぞれの手順について説明します。
※作業するディレクトリをworkとしています。

Kaggleデータの準備

kaggleコマンドでコンペデータをダウンロードします。
前述の通り、事前にkaggleのHPでRulesをAcceptしておく必要がありますので、注意してください。

kaggle competitions download -c acquire-valued-shoppers-challenge

また、ダウンロードしたコンペデータは、そのままではzipとgzipが二重でかかっていますので、下記のコードを実行して解凍します。

import shutil
import gzip

shutil.unpack_archive("/work/acquire-valued-shoppers-challenge/acquire-valued-shoppers-challenge.zip", "/work/acquire-valued-shoppers-challenge")

source_file = "/work/acquire-valued-shoppers-challenge/transactions.csv.gz"
target_file = "/work/acquire-valued-shoppers-challenge/transactions.csv"

with gzip.open(source_file, mode="rb") as gzip_file:
    with open(target_file, mode="wb") as decompressed_file:
        shutil.copyfileobj(gzip_file, decompressed_file)

パッケージ・モジュールのインポートと初期設定

プログラム実行に必要なパッケージ・モジュールをインポートします。
こちらでエラーが出る場合は、各種パッケージをpip installでインストールしてください。
また、各種初期設定も実施します。

import os

import numpy as np
import pandas as pd
from scipy import stats
import seaborn as sns
from sklearn import model_selection
from sklearn import preprocessing
import tensorflow as tf
from tensorflow import keras
from tensorflow.keras import backend as K
import tensorflow_probability as tfp
import tqdm
from typing import Sequence

import lifetime_value as ltv

tfd = tfp.distributions
%config InlineBackend.figure_format='retina'
sns.set_style('whitegrid')
pd.options.mode.chained_assignment = None  # default='warn'

グローバル変数の定義

後ほど記載する関数内で使用するグローバル変数を定義します。
Githubにあるコードでは、COMPANYがstring型で定義されていますが、後述のCSVデータ読み込みでこの値がint型でないとqueryで引っかからなかったので、変更して記載しています(これ以降に記載する関数内のCOMPANYもint型であることを考慮して一部変更しています)。

COMPANY = 103600030  # @param { isTemplate: true, type: 'string'}
LOSS = 'ziln'  # @param { isTemplate: true, type: 'string'} ['mse', 'ziln']
MODEL = 'dnn'  # @param { isTemplate: true, type: 'string'} ['linear', 'dnn']
LEARNING_RATE = 0.0002  # @param { isTemplate: true}
EPOCHS = 400  # @param { isTemplate: true, type: 'integer'}
OUTPUT_CSV_FOLDER = '/work/result'  # @param { isTemplate: true, type: 'string'}

CATEGORICAL_FEATURES = ['chain', 'dept', 'category', 'brand', 'productmeasure']
NUMERIC_FEATURES = ['log_calibration_value']

ALL_FEATURES = CATEGORICAL_FEATURES + NUMERIC_FEATURES

CSVデータの読み込み

上記で解凍したCSVデータを読み込む関数です。
データから特定のcompanyを指定して抽出します。

def load_transaction_data(company):
  all_data_filename = '/work/acquire-valued-shoppers-challenge/transactions.csv'
  one_company_data_filename = (
      '/work/acquire-valued-shoppers-challenge/transactions_company_"{}".csv'
      .format(COMPANY))
  if os.path.isfile(one_company_data_filename):
    df = pd.read_csv(one_company_data_filename)
  else:
    data_list = []
    chunksize = 10**6
    # 350 iterations
    for chunk in tqdm.tqdm(pd.read_csv(all_data_filename, chunksize=chunksize)):
      data_list.append(chunk.query("company=={}".format(company)))
    df = pd.concat(data_list, axis=0)
    df.to_csv(one_company_data_filename, index=None)
  return df

前処理

学習する前段階の処理を実行する関数です。
学習＆検証の際に使用するデータを作成し、元データにマージしています。
holdout_valueがLTVです。

def preprocess(df):
  df = df.query('purchaseamount>0')
  df['date'] = pd.to_datetime(df['date'], format='%Y-%m-%d')
  df['start_date'] = df.groupby('id')['date'].transform('min')

  # Compute calibration values
  calibration_value = (
      df.query('date==start_date').groupby('id')
      ['purchaseamount'].sum().reset_index())
  calibration_value.columns = ['id', 'calibration_value']

  # Compute holdout values
  one_year_holdout_window_mask = (
      (df['date'] > df['start_date']) &
      (df['date'] <= df['start_date'] + np.timedelta64(365, 'D')))
  holdout_value = (
      df[one_year_holdout_window_mask].groupby('id')
      ['purchaseamount'].sum().reset_index())
  holdout_value.columns = ['id', 'holdout_value']

  # Compute calibration attributes
  calibration_attributes = (
      df.query('date==start_date').sort_values(
          'purchaseamount', ascending=False).groupby('id')[[
              'chain', 'dept', 'category', 'brand', 'productmeasure'
          ]].first().reset_index())

  # Merge dataframes
  customer_level_data = (
      calibration_value.merge(calibration_attributes, how='left',
                              on='id').merge(
                                  holdout_value, how='left', on='id'))
  customer_level_data['holdout_value'] = (
      customer_level_data['holdout_value'].fillna(0.))
  customer_level_data[CATEGORICAL_FEATURES] = (
      customer_level_data[CATEGORICAL_FEATURES].fillna('UNKNOWN'))

  # Specify data types
  customer_level_data['log_calibration_value'] = (
      np.log(customer_level_data['calibration_value']).astype('float32'))
  customer_level_data['chain'] = (
      customer_level_data['chain'].astype('category'))
  customer_level_data['dept'] = (customer_level_data['dept'].astype('category'))
  customer_level_data['brand'] = (
      customer_level_data['brand'].astype('category'))
  customer_level_data['category'] = (
      customer_level_data['category'].astype('category'))
  customer_level_data['label'] = (
      customer_level_data['holdout_value'].astype('float32'))
  return customer_level_data

CSVデータの読み込みと前処理の適用を実行し、学習・検証用のデータを作成する関数です。

def load_customer_level_csv(company):
  customer_level_data_file = (
      '/work/acquire-valued-shoppers-challenge/customer_level_data_company_"{}".csv'
      .format(company))
  if os.path.isfile(customer_level_data_file):
    customer_level_data = pd.read_csv(customer_level_data_file)
  else:
    customer_level_data = preprocess(load_transaction_data(company))
  for cat_col in CATEGORICAL_FEATURES:
    customer_level_data[cat_col] = (
        customer_level_data[cat_col].astype('category'))
  for num_col in [
      'log_calibration_value', 'calibration_value', 'holdout_value'
  ]:
    customer_level_data[num_col] = (
        customer_level_data[num_col].astype('float32'))

  return customer_level_data

COMPANYを引数として学習・検証用のデータを作成します。

# Processes data. 350 iteration in total. May take 10min.
customer_level_data = load_customer_level_csv(COMPANY)

以下で、顧客ごとのLTVの分布を確認できます。
横軸がLTV(対数)、縦軸が顧客のカウントのヒストグラムです。

customer_level_data.label.apply(np.log1p).hist(bins=50)

学習用＆検証用データの準備

上記で作成したデータを学習用と検証用に分ける関数です。
線形モデル、DNNモデルそれぞれで定義しています。

def linear_split(df):
  # get_dummies preserves numeric features.
  x = pd.get_dummies(df[ALL_FEATURES], drop_first=True).astype('float32').values
  y = df['label'].values
  y0 = df['calibration_value'].values

  x_train, x_eval, y_train, y_eval, y0_train, y0_eval = (
      model_selection.train_test_split(
          x, y, y0, test_size=0.2, random_state=123))

  return x_train, x_eval, y_train, y_eval, y0_eval


def dnn_split(df):
  for key in CATEGORICAL_FEATURES:
    encoder = preprocessing.LabelEncoder()
    df[key] = encoder.fit_transform(df[key])

  y0 = df['calibration_value'].values
  df_train, df_eval, y0_train, y0_eval = model_selection.train_test_split(
      df, y0, test_size=0.2, random_state=123)

  def feature_dict(df):
    features = {k: v.values for k, v in dict(df[CATEGORICAL_FEATURES]).items()}
    features['numeric'] = df[NUMERIC_FEATURES].values
    return features

  x_train, y_train = feature_dict(df_train), df_train['label'].values
  x_eval, y_eval = feature_dict(df_eval), df_eval['label'].values

  return x_train, x_eval, y_train, y_eval, y0_eval

学習

線形モデル、DNNモデルを作成する関数です。

def linear_model(output_units):
  return tf.keras.experimental.LinearModel(output_units)


def embedding_dim(x):
  return int(x**.25) + 1


def embedding_layer(vocab_size):
  return tf.keras.Sequential([
      tf.keras.layers.Embedding(
          input_dim=vocab_size,
          output_dim=embedding_dim(vocab_size),
          input_length=1),
      tf.keras.layers.Flatten(),
  ])


def dnn_model(output_units, df):
  numeric_input = tf.keras.layers.Input(
      shape=(len(NUMERIC_FEATURES),), name='numeric')

  embedding_inputs = [
      tf.keras.layers.Input(shape=(1,), name=key, dtype=np.int64)
      for key in CATEGORICAL_FEATURES
  ]

  embedding_outputs = [
      embedding_layer(vocab_size=df[key].nunique())(input)
      for key, input in zip(CATEGORICAL_FEATURES, embedding_inputs)
  ]

  deep_input = tf.keras.layers.concatenate([numeric_input] + embedding_outputs)
  deep_model = tf.keras.Sequential([
      tf.keras.layers.Dense(64, activation='relu'),
      tf.keras.layers.Dense(32, activation='relu'),
      tf.keras.layers.Dense(output_units),
  ])
  return tf.keras.Model(
      inputs=[numeric_input] + embedding_inputs, outputs=deep_model(deep_input))

損失関数の定義とモデルの作成を実行します。

if LOSS == 'mse':
  loss = keras.losses.MeanSquaredError()
  output_units = 1

if LOSS == 'ziln':
  loss = ltv.zero_inflated_lognormal_loss
  output_units = 3

if MODEL == 'linear':
  x_train, x_eval, y_train, y_eval, y0_eval = linear_split(customer_level_data)
  model = linear_model(output_units)

if MODEL == 'dnn':
  x_train, x_eval, y_train, y_eval, y0_eval = dnn_split(customer_level_data)
  model = dnn_model(output_units, customer_level_data)

作成したモデルを学習用データに適用し、学習します。

model.compile(loss=loss, optimizer=keras.optimizers.Adam(lr=LEARNING_RATE))

callbacks = [
    tf.keras.callbacks.ReduceLROnPlateau(monitor='val_loss', min_lr=1e-6),
    tf.keras.callbacks.EarlyStopping(monitor='val_loss', patience=10),
]

history = model.fit(
    x=x_train,
    y=y_train,
    batch_size=1024,
    epochs=EPOCHS,
    verbose=2,
    callbacks=callbacks,
    validation_data=(x_eval, y_eval)).history

以下で、損失関数のログを確認できます。

pd.DataFrame(history)[['loss', 'val_loss']][2:].plot()

検証(ゲインチャート、正規化ジニ係数、デシルチャート、スピアマンの順位相関係数)

以降では、複数の指標でモデルの有効性を検証します。
下記で各指標で確認できることと有効性の見方を一覧で示します。

指標	確認できること	有効性の見方
ゲインチャート	LTVの高い顧客と低い顧客の判別性	対角線45度に近いほど判別性が低い
正規化ジニ係数	同上	[0,1]の範囲で、1に近いほど判別性が高く、 0に近いほど判別性が低い
デシルチャート	正解LTVと予測LTVが上手くキャリブレーションされているか	同じ分位で各平均値が密接であるほど、キャリブレーションが上手くいっている
スピアマンの順位相関係数	正解LTVと予測LTVの相関	[-1,1]の範囲で、1、-1に近いほどそれぞれ正の相関、負の相関があり、0に近いほど無相関である

まず、検証用データに対して学習したモデルを適用し、推論します。

if LOSS == 'mse':
  y_pred = model.predict(x=x_eval, batch_size=1024).flatten()

if LOSS == 'ziln':
  logits = model.predict(x=x_eval, batch_size=1024)
  y_pred = ltv.zero_inflated_lognormal_pred(logits).numpy().flatten()

ゲインチャート

モデルによって、LTVの高い顧客とそうでない顧客をどの程度判別できているかをゲインチャートで確認します。
対角線45度に近いほど、判別性が低くなります。
以下で、ゲインチャートを作成します。

df_pred = pd.DataFrame({
    'y_true': y_eval,
    'y_pred': y_pred,
})

gain = pd.DataFrame({
    'lorenz': ltv.cumulative_true(y_eval, y_eval),
    'baseline': ltv.cumulative_true(y_eval, y0_eval),
    'model': ltv.cumulative_true(y_eval, y_pred),
})

num_customers = np.float32(gain.shape[0])
gain['cumulative_customer'] = (np.arange(num_customers) + 1.) / num_customers

ax = gain[[
    'cumulative_customer',
    'lorenz',
    'baseline',
    'model',
]].plot(
    x='cumulative_customer', figsize=(8, 5), legend=True)

ax.legend(['Groundtruth', 'Baseline', 'Model'], loc='upper left')

ax.set_xlabel('Cumulative Fraction of Customers')
ax.set_xticks(np.arange(0, 1.1, 0.1))
ax.set_xlim((0, 1.))

ax.set_ylabel('Cumulative Fraction of Total Lifetime Value')
ax.set_yticks(np.arange(0, 1.1, 0.1))
ax.set_ylim((0, 1.05))
ax.set_title('Gain Chart')

正規化ジニ係数

上記で算出したゲインから正規化ジニ係数を求めます。
この数値が1に近いほど、顧客の判別性が高いと判断できます。
以下で、正規化ジニ係数を算出します(normalizedのカラムが正規化ジニ係数です)。

gini = ltv.gini_from_gain(gain[['lorenz', 'baseline', 'model']])
gini

デシルチャート

正解LTVに対して、予測LTVが上手くキャリブレーションされているかを測るため、デシルチャートを用います。
デシルチャートは、正解値を10分位で分け、分位ごとに平均値を正解と予測で算出します。
同じ分位で各平均値が密接であるほど、キャリブレーションが上手くいっていると言えます。
以下で、デシルチャートを作成します。

df_decile = ltv.decile_stats(y_eval, y_pred)

ax = df_decile[['label_mean', 'pred_mean']].plot.bar(rot=0)

ax.set_title('Decile Chart')
ax.set_xlabel('Prediction bucket')
ax.set_ylabel('Average bucket value')
ax.legend(['Label', 'Prediction'], loc='upper left')

スピアマンの順位相関係数

正解LTVと予測LTVの相関を確認します。
以下で、スピアマンの順位相関係数を算出します。

def spearmanr(x1: Sequence[float], x2: Sequence[float]) -> float:
  """Calculates spearmanr rank correlation coefficient.

  See https://docs.scipy.org/doc/scipy/reference/stats.html.

  Args:
    x1: 1D array_like.
    x2: 1D array_like.

  Returns:
    correlation: float.
  """
  return stats.spearmanr(x1, x2, nan_policy='raise')[0]


spearman_corr = spearmanr(y_eval, y_pred)
spearman_corr

出力結果：0.3142630651229492

出力結果の保存

最後に、出力結果を保存します。

df_metrics = pd.DataFrame(
    {
        'company': COMPANY,
        'model': MODEL,
        'loss': LOSS,
        'label_mean': y_eval.mean(),
        'pred_mean': y_pred.mean(),
        'label_positive': np.mean(y_eval > 0),
        'decile_mape': df_decile['decile_mape'].mean(),
        'baseline_gini': gini['normalized'][1],
        'gini': gini['normalized'][2],
        'spearman_corr': spearman_corr,
    },
    index=[0])

output_path = os.path.join(OUTPUT_CSV_FOLDER, str(COMPANY))

if not os.path.isdir(output_path):
  os.makedirs(output_path)

output_file = os.path.join(output_path,
                           '{}_regression_{}.csv'.format(MODEL, LOSS))

df_metrics.to_csv(output_file, index=False)

まとめ

本記事では、マーケティングにおいて重要な指標であるLTVを機械学習・AIを用いて予測する方法についてご紹介しました。
異なるデータセットや問題設定に対して適用するとどうなるか、自分も色々試してみたいと思いますので、みなさまもぜひお試しください。

本記事で少しでもみなさまのお役に立てれば幸いです。

参考文献

LTV関連

Wang, Xiaojing, Liu, Tianqi, and Miao, Jingang. (2019). A Deep Probabilistic Model for Customer Lifetime Value Prediction.
https://arxiv.org/abs/1912.07753
元論文の公式Github
https://github.com/google/lifetime_value
LTVの計算方法【モデル別3大計算式】…サブスクEC/SaaS事業者向け
https://www.tsuhannews.jp/shopblogs/detail/66619
LTV(顧客生涯価値)とは？3種類の計算方法と活用方法をご紹介
https://collabos-service.jp/blog/management/4231/
LTVとは？5つの向上の利点や計算方法、分析の仕方を解説
https://botchan.chat/base/ltv2