[mathjax]
購買の行動や意思決定にまつわるビッグデータの出現や消費者の多様化が進む中、企業のマーケティング活動においてはミクロな市場に対して理解を深め、的確なインサイトを突いた訴求を実現することが求められています。SEEDATAでは定性的なリサーチに加え、統計モデルを駆使した消費者の理解に取り組んでおり、今回はマーケティングにおける線形回帰モデルの例について見ていきます。
なお、分析対象のデータおよび回帰モデルは、『ビッグデータ時代のマーケティング』（佐藤忠彦, 樋口知之, 講談社 理工学専門書）の2.4節の手順に倣い、考察していきます。分析はPythonで行います。
線形回帰の数理的な知識については、次の記事を参照してください。
シミュレーションデータの生成
分析対象となるデータを人工的に生成します。
時刻(t)における、販売個数（目的変数）、価格掛け率および山積み陳列実施の有無（説明変数）をそれぞれ、(y_t, Z_t^1, Z_t^2)とし、データの生成過程として次を設定します。なお、価格掛け率のここでの定義は、定価を(k)円とした時、割引販売価格を(kZ_t^1)と表すようなものとします。
$$
log (y_t) = 0.5 – 2.5log(Z_t^1) + 0.8 Z_t^2 + epsilon_t big(epsilon_t sim mathbb{N}(0, 0.3) big) tag{1}
$$
シミュレーションデータは参考書籍に倣い、次のように生成します。データ数は100個です。
- (epsilon_t)を平均0, 分散0.3の正規分布から生成する。
- 価格掛け率(Z_t^1)を区間[0.5, 1.0]の一様乱数で生成する。これは、最大値引率を50%までとした設定に対応する。
- 区間[0, 1]の一様乱数(u_t)を発生させ、(u_t leq 0.3)かつ(Z_t^1 leq 0.8)のとき、(Z_t^2=1)とし、それ以外のときは(Z_t^2=0)とする。(Z_t^1 leq 0.8)は、20%以上値引きされた際に、山積み陳列がなされることが多いことを反映した設定で、また、(u_t leq 0.3)はその時の山積み陳列の実施頻度に対応する。
- (y_t)を式(1)の通り、生成する。
データを生成するPythonコードは次の通りです。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
np.random.seed(seed=100) # 乱数の種を設定def generate_data(n_samples=100):
epsilon_t = np.random.normal(loc=0.0, scale=0.3, size=n_samples)
z1_t = 0.5 + np.random.rand(n_samples) * 0.5
z2_t = np.zeros(n_samples, dtype=np.int)
u_t = np.random.rand(n_samples)
z2_t[(u_t <= 0.3) & (z1_t <= 0.8)] = 1
y_t = np.exp(0.5 - 2.5 * np.log(z1_t) + 0.8 * z2_t + epsilon_t)
return y_t, z1_t, z2_ty_t, z1_t, z2_t = generate_data()さて、生成したデータを次のようなコードでプロットしてみましょう。通常、2次元上に全てをプロットできるのは説明変数が1つしかない単回帰分析のみですが、今回は2つ目の変数がクラスであるため、凡例を用いて2次元上に可視化することができます。
fig, ax = plt.subplots(figsize=(8, 6))
for i in [0, 1]:
ax.scatter(z1_t[z2_t==i], y_t[z2_t==i], s=40, label=r'$Z_t^2={}$'.format(i))
ax.set_title('n_samples = 100')
ax.set_xlabel(r'$Z_t^1$', fontsize=15)
ax.set_ylabel(r'$y_t$', fontsize=15)
ax.grid(True)
ax.legend(fontsize=15)
値引率が大きい方（グラフの左側）が、あるいは山積み陳列を行った方（オレンジ）が、販売個数が増加する傾向が見てとれます。
補足
なお、販売個数が「非負整数」であるのに、なぜポアソン回帰ではなく線形回帰を適用するのかという疑問があるかもしれません。答えとしては、今回は「非負整数」という性質にこだわるメリットが大きくありません。例えば地震の発生回数や植物の種子数などに回帰を行う場合は「整数であること」および「現象発生頻度の少数性」が重要となりますが、今回のケースで最終的に知りたいのは売り上げに関する情報であり、販売個数が整数でなく連続値をとることによる不便はそれほど大きくありません。よって、よりシンプルなモデルである線形回帰を採用することになります。
データの分析
さて、シミュレーションデータの生成が完了したため、「データの生成過程を知らない」前提でデータに適した線形回帰モデルを構築するケースを考えましょう。
上記のグラフを見ていて第一に気がつくことは、グラフの右側にいくにしたがって目的変数(y_t)の散らばりが小さくなっていくことです。このようなパターンは、素朴に線形回帰を行うだけでは表現できません。そこで、線形回帰を適用できるように(y_t)に対する変数変換を考えてやる必要があります。
(x)が増えるにしたがって、(f(x))の変化が緩やかになっていくような(f)の代表的なものは対数関数*です。そこで、今回は(y)に対して対数変換を施すと、次のようなグラフになります。
*補足
もちろん、このような性質を持つ関数は対数関数に限りません。しかし、今回行いたいのは「真の関数」を知ることではなく、現象のモデル化です。よってまずはシンプルな関数によりモデリングを試みることが肝要だと言えるでしょう。
先ほどまでのパターンが緩和されたことが分かります。さて、次に考慮すべきは以下の2点です。
- 説明変数として採用可能な変数はどの変数か。
- 説明変数に変数変換を施すかどうか。
今回は、AIC（赤池情報量基準）と呼ばれるものを用いてモデル選択を行います。AICの詳細については今回は割愛しますが、追って議論します。上記の2点を考慮して、次の5つのモデル（モデル1~モデル5とします）のAICを比較します。
- (log(y_t) = beta_0 + epsilon_t)
- (log(y_t) = beta_0 + beta_1 Z_t^1 + epsilon_t)
- (log(y_t) = beta_0 + beta_1 Z_t^1 + beta_2 Z_t^2 + epsilon_t)
- (log(y_t) = beta_0 + beta_1 log(Z_t^1) + epsilon_t)
- (log(y_t) = beta_0 + beta_1 log(Z_t^1) + beta_2 Z_t^2 + epsilon_t)
今回は説明変数が2つしかないため、上記のようにモデルを全て列挙することができますが、説明変数が例えば10個ある場合は全てを列挙すると計算コストが高くなってしまいます。その場合はステップワイズ法やL1正則化と呼ばれる手法を取ることがあります。
Pythonでこのような分析を行うには、statsmodelsが便利です。
import statsmodels.formula.api as smfdata = pd.DataFrame(np.stack([y_t, z1_t, z2_t], axis=1), columns=['y_t', 'z1_t', 'z2_t'])例えば、モデル5を採用すると次のように回帰を実行することができます。
model = smf.ols('np.log(y_t) ~ np.log(z1_t) + z2_t', data=data).fit()今回は、5つのモデルに対してAICを比較したいので、次のように実行してみます。
predictors = ['1', 'z1_t', 'z1_t + z2_t', 'np.log(z1_t)', 'np.log(z1_t) + z2_t']
models = [smf.ols('np.log(y_t) ~ {}'.format(predictor), data=data).fit() for predictor in predictors]
aics = [model.aic for model in models]
print(aics)[233.13424268382613, 125.04124090867109, 40.710628861601634, 125.664754058125, 41.285362746380116]悩ましいのは、モデル3（(Z_t^1)に対数変換を施さないモデル）とモデル5（(Z_t^1)に対数変換を施すモデル）のAICが非常に近い値であるということです。今回はデータ数が100個しかないため、例えばシミュレーションデータ生成の過程ではじめに設定した乱数の種が異なれば、容易にこれらのAICは上下するでしょう。実際のデータ分析では必ずしもデータ数を増やすことができるとは限りません。100個という限られたデータ数から最適なモデルを選択するにはどうしたら良いでしょうか。
モデル3とモデル5は、今回のケースではモデルの性能自体に大きな優劣はありません（決定係数や尤度など他の指標を参考にしたとしても）。よって、この先は、モデルとして解釈がしやすく、応用がきくかどうかを根拠に最適なモデルを考えましょう。
マーケティングにおける重要な概念の一つに、「需要の価格弾力性」と呼ばれるものがあります。これは、（需要の変化率/価格の変化率）で表されます。結論としては、価格弾力性を定式化しやすいモデルはモデル5です。モデル比較の肝である説明変数(Z_t^1)は価格掛け率ですが、定価を(k)円とすれば売価は(k Z_t^1)で表されます。簡単な考察のため、(k=1)として価格弾力性がどのように表されるか、モデル5について確かめてみましょう。
$$
begin{eqnarray}
log(y_t) = beta_0 + beta_1 log(Z_t^1) + beta_2 Z_t^2 + epsilon_t tag{2}\[1.0em]
frac{1}{y_t} cdot frac{d y_t}{d Z_t^1} = beta_1 frac{1}{Z_t^1} tag{3} \[1.0em]
beta_1 = big(frac{d y_t}{y_t}big) big/ big(frac{d Z_t^1}{Z_t^1}big) tag{4}
end{eqnarray}
$$
式(3)は式(2)の両辺を売価である(Z_t^1)で微分したものです。式(4)の右辺が、価格弾力性の定義となっていることが確認できます。式(4)より、モデル5においては価格弾力性が(beta_1)と表されます（任意の定価(k)については(beta_1/k)）。よって、今回はモデルの応用性の観点から、モデル5を採用します。
results = smf.ols('np.log(y_t) ~ np.log(z1_t) + z2_t', data=data).fit()
print(results.summary()) OLS Regression Results
==============================================================================
Dep. Variable: np.log(y_t) R-squared: 0.859
Model: OLS Adj. R-squared: 0.856
Method: Least Squares F-statistic: 295.3
Date: Fri, 19 Jun 2020 Prob (F-statistic): 5.59e-42
Time: 01:30:50 Log-Likelihood: -17.643
No. Observations: 100 AIC: 41.29
Df Residuals: 97 BIC: 49.10
Df Model: 2
Covariance Type: nonrobust
================================================================================
coef std err t P>|t| [0.025 0.975]
--------------------------------------------------------------------------------
Intercept 0.4710 0.054 8.777 0.000 0.365 0.578
np.log(z1_t) -2.4162 0.172 -14.078 0.000 -2.757 -2.076
z2_t 0.8959 0.078 11.537 0.000 0.742 1.050
==============================================================================
Omnibus: 3.079 Durbin-Watson: 2.028
Prob(Omnibus): 0.215 Jarque-Bera (JB): 1.894
Skew: -0.065 Prob(JB): 0.388
Kurtosis: 2.338 Cond. No. 6.62
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.線形モデルの係数を参照すると、次のようにモデルがフィットされたことが分かります。
$$
log (y_t) = 0.4710 – 2.4162log(Z_t^1) + 0.8959 Z_t^2 + epsilon_t tag{2}
$$
これは、シミュレーションデータを生成した次のモデルと大まかに一致することが分かります。
$$
log (y_t) = 0.5000 – 2.5000log(Z_t^1) + 0.8000 Z_t^2 + epsilon_t
$$
最後に、現況再現性（モデルによる真値の再現性）をグラフで確認します。
fig, ax = plt.subplots(figsize=(12, 4))
x = np.arange(0, 100)
pred = results.predict(data)
ax.plot(x, np.log(y_t), label='True')
ax.plot(x, pred, label='Predicted')
ax.plot(x, np.log(y_t) - pred, label='Residual')
ax.set_title('Reproduction of Sales Numbers')
ax.set_xlabel(r'$t$', fontsize=15)
ax.set_ylabel(r'$log(y_t)$', fontsize=15)
ax.grid(True)
ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0, fontsize=15)
Trueは真値、Predictedはモデルによる再現値、Residualは残差を示します。販売個数がモデルにより十分に記述できていることを確認できます。
まとめ
今回は、マーケティングにおける線形回帰モデルの適用例について、Pythonを用いて簡単なケースを見ていきました。このようなモデルの発展形については、今後の記事でも触れていきます。
参考文献
『ビッグデータ時代のマーケティング』（佐藤忠彦, 樋口知之, 講談社 理工学専門書）