Python: 回帰モデルで市場反応分析

Python 統計学 マーケティング

これまで緑本などで学んできた統計モデルを、マーケティングに応用するための勉強を行っています。

今回は市場反応分析を線形回帰モデルとポアソン回帰モデルで行います。
市場反応分析に関する理論や使用するデータは、マーケティングの統計モデル (統計解析スタンダード)の第3章「消費者の市場反応のモデル化」を参考にさせていただいてます。

マーケティングの統計モデル (統計解析スタンダード)

実装はこちらです。

statistical_modeling_for_marketing/chapter3.ipynb at master · ohke/statistical_modeling_for_marketing · GitHub

市場反応分析

マーケティング施策は、顧客の行動(購入や予約など)の変化を促すために行われます。

  • 例えば、ある商品を値下げすると、その商品の販売個数が増える
  • 例えば、ある商品を山積み陳列すると、その商品の販売個数は増えるが、類似商品の販売個数は減る ...など

施策を効率的に実施するためには、上のように施策の効果を定性的に捉えるだけではなく、定量化する必要があります。

  • ある商品をX%値下げすると、その商品の販売個数がY%増える
  • ある商品を山積み陳列すると、その商品の販売個数はX%増えるが、類似商品の販売個数はY%減る

このように、施策を受けた顧客の行動変化を定量的にモデル化することを市場反応分析といいます。

価格と陳列によるPI値の変化を分析

それでは、擬似的なデータを使って、具体的に市場反応分析をしていきます。データはこちらからダウンロードできます。

最初の事例として、2つのブランドの醤油(醤油Aと醤油Bとする)を販売するケースを考えます。

両方の商品について「価格を変更する」「山積み陳列する」というそれぞれ2つの施策が組み合わせて行われています。
上の施策を実施した結果、醤油AのPI値の対数がどのように変化するのかをモデル化します。つまり醤油AのPI値の対数が目的変数となります。

データを概観する

まずはデータを概観します。

  • LogPI_Aが醤油AのPI値の対数
    • PI値とは、その商品の販売点数÷客数で計算される値で、その商品の人気度・支持度を表す指標で、
    • LogPI_Aは、 \log \frac{\mbox{販売点数}\times \mbox{客数}}{1000} で計算されています
  • LogPriceIndex_AとLogPriceIndex_Bはそれぞれ、 \log \frac{\mbox{商品Aの売価}}{\mbox{商品Aの期間最大売価}} 、 \log \frac{\mbox{商品Bの売価}}{\mbox{商品Bの期間最大売価}} で計算される値で、ここでは価格掛率と呼びます
    • 値が小さいほど安売りされていることを表す
  • Display_AとDisplay_Bは、それぞれの商品が山積みされていれば1となるダミー変数

f:id:ohke:20180324132509p:plain

醤油Aのと価格掛率の相関を見てみます。

  • 醤油Aが安売りされている(価格掛率LogPriceIndex_Aが小さい)と、PI値が大きくなる
  • 醤油Bが安売りされている(価格掛率LogPriceIndex_Bが小さい)と、PI値が大きくなる
    • 醤油Aと醤油Bが競合関係にあることがわかります

f:id:ohke:20180324134021p:plain

山積み陳列の有無による醤油AのPI値の分布を箱ひげ図で見てみます。

  • 醤油Aが山積み陳列される(Display_A=1)と、PI値が大きくなる
  • 醤油Bが山積み陳列される(Display_B=1)と、PI値が小さくなる

f:id:ohke:20180324134518p:plain

2標本t検定で醤油Bの山積み陳列の効果を確認

上の箱ひげ図では、醤油Aの山積み陳列によるPI値の増加ははっきり現れていますが、醤油Bの山積み陳列による影響はそれより弱く、誤差という見方ができるかもしれません。

モデル化の前に、醤油Bの山積み陳列によるPI値の変化が誤差ではない(起きにくい)ことを、念のため検定で確認しておきます。ここでは、価格掛率は一旦無視します。なおこの判断は誤りであることが後にわかります。線形回帰モデルを作る時にそれも見ていきます。

なお、検定では統計学入門 (基礎統計学?)の第11・12章を参考にしています。

平均値の差が有意かどうかをt検定で確認します。2つのグループで分散が異なっている(Display_B=1では0.765、Display_B=0では0.910)ため、ウェルチのt検定を行います。

まずはt統計量を計算します。 s が不偏分散、 n がデータ数です。下付き添字が1ならDisplay_B=1、0ならDisplay_B=0を表します。

$$ t=\frac{\overline{PI_1} - \overline{PI_0}}{\sqrt{\frac{s_1^2}{n_1}+\frac{s_0^2}{n_0}}} $$

t分布の自由度 \nu は以下の式で近似できます。 \nu から最も近い整数 \nu^{*} を選び、自由度 \nu^{*} のt分布の5%の下側信頼限界の t_{0.05}(\nu^{*}) 値を算出します。上で求めた値がこの値よりも小さければ(つまり t \mbox{＜} t_{0.05}(\nu^{*}) )、有意差があると判定します。これは左片側検定です。

$$ \nu=\frac{(\frac{s_1^2}{n_1}+\frac{s_0^2}{n_0})^2}{\frac{(s_1^2 / n_1)^2}{n_1-1}+\frac{(s_0^2 / n_0)^2}{n_0-1}} $$

Pythonではこんな感じで計算できます。

m = len(treatment_A)
n = len(control_A)
s_1 = treatment_A.var()
s_2 = control_A.var()

t = (treatment_A.mean() - control_A.mean()) / math.sqrt(s_1**2 / m + s_2**2 / n)
print('t:', t)

# 自由度νを近似的に計算
nu = (s_1 ** 2 / m + s_2**2 / n)**2 / ((s_1**2 / m)**2 / (m-1) + (s_2**2 / n)**2 / (n-1))
nu_star = round(nu) # 整数ν*に丸める

print('-t_0.95(ν*):', -stats.t.ppf(0.95, nu_star))

t\mbox{＜}-t_{0.05}(\nu^{*}) \iff -8.47\mbox{＜}-1.65 となるため、グループ間の平均値に有意差があることも確認できました。

線形回帰でモデル化

PI値は連続変数なので、正規分布を仮定した線形回帰でモデル化します。

$$ LogPIA_t=\beta_0+LogPriceIndexA_t\beta_1+LogPriceIndexB_t\beta_2+DisplayA_t\beta_3+DisplayB_t\beta_4+\epsilon_t $$

import statsmodels.api as sm

model = sm.OLS(data['LogPI_A'],
               sm.add_constant(data[['LogPriceIndex_A', 'LogPriceIndex_B', 'Display_A', 'Display_B']]))
result = model.fit()

結果は以下のとおりです。いずれも各施策の効果が係数となっています。

  • LogPriceIndex_Aの係数が負数(-5.02)、Display_Aの係数が正(0.774)となっており、醤油Aの販売にプラスに働いていることが数値化されています
  • LogPriceIndex_Bの傾きは正数(1.93)で、競合商品の値引きの影響を受けることが確認できます

f:id:ohke:20180324152547p:plain

一方でDisplay_Bの傾きがわずかですが正数(0.0374)となっています。95%信頼区間では0を跨いでいるので負になる可能性もありますが、検定とは反した結果となっています。
LogPriceIndex_BとDisplay_Bの相関が高い(-0.641)にも関わらず、LogPriceIndex_Bを無視して検定したため、Display_Bの影響を過大評価してしまいました。安売りする商品を山積み陳列するというのは、相関係数的にも直感的にもよくあるパターンに見えます。

  • 他の条件が異なるグループ間での比較にt検定を使うべきではありません
    • LogPriceIndex_Bも揃えたグループ間で平均値を比較すべきです
  • こうした複数の変数がからむ場合、初めから回帰モデルに落とした方が良いかと思います

モデルにもLogPriceIndex_B×Display_Bの交互作用項を追加する必要がありそうです。
LogPriceIndex_B×Display_Bを交互作用項(LogPriceIndex_Display_B)としてパラメータに追加して線形回帰した結果が以下です。Display_Bの係数は-0.330(95%信頼区間は[-0.550, -0.110])で、Display_B単体では醤油Aの販売に負に影響していることがわかります。

f:id:ohke:20180326083747p:plain

価格と陳列による販売個数の変化を分析

同じく醤油Aと醤油Bを販売する事例について考えますが、今回は販売個数(つまりカウントデータ)のモデル化を行います。

データを概観する

同様にデータを概観します。

  • Sale_Unit_Aは醤油Aの販売個数
  • LogPriceIndex_AとLogPriceIndex_Bは、先の事例と同じ価格掛率ですが、今回は対数化されていません(0以上1以下の値域をとる)
  • Display_AとDisplay_Bも同じく山積み陳列の有無
  • Visitorsは来店者数

f:id:ohke:20180324155858p:plain

販売個数と価格掛率・来店者数の分布をプロットします。

  • 醤油Aが安くなれば販売個数は増え、醤油Bが安くなれば販売個数が減る傾向が見えます
  • 来店者数と販売個数も緩い正相関がありそう

f:id:ohke:20180324160455p:plain

山積み陳列の有無毎に、販売個数を箱ひげ図でプロットします。ここでもやはり、競合商品(醤油B)を山積み陳列した時は醤油Aの販売個数が減る傾向が伺えます。

f:id:ohke:20180324160506p:plain

ポアソン回帰でモデル化

販売個数はカウントデータなので、ポアソン回帰(+対数リンク関数)でモデル化します。ポアソン分布の平均値 \mu は以下で計算します。

  • Visitorsはオフセット(傾きを持たない定数項)として扱っています

$$ \log(\mu_t)=\log(Visitors_t)+\beta_0+PriceIndexA_t\beta_1+PriceIndexB_t\beta_2+DisplayA_t\beta_3+DisplayB_t\beta_4 $$

Pythonでの実装は以下となります。

  • offsetオプションで、オフセット値を指定できます
import statsmodels.formula.api as smf

model = smf.glm('Sale_Unit_A ~ PriceIndex_A + PriceIndex_B + Display_A + Display_B',
                data=data,
                offset=np.log(data['Visitors']),
                family=sm.families.Poisson(link=sm.families.links.log))
result = model.fit()

結果は以下です。今回はDisplay_Bの傾きが負数となっており、わずかながらですが競合商品の山積み陳列が悪影響を与えることがわかります。

f:id:ohke:20180324161921p:plain