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は、 で計算されています
• LogPriceIndex_AとLogPriceIndex_Bはそれぞれ、 、 で計算される値で、ここでは価格掛率と呼びます
  • 値が小さいほど安売りされていることを表す
• Display_AとDisplay_Bは、それぞれの商品が山積みされていれば1となるダミー変数

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

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

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

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

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

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

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

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

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

まずはt統計量を計算します。 が不偏分散、 がデータ数です。下付き添字が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分布の自由度 は以下の式で近似できます。 から最も近い整数 を選び、自由度 のt分布の5%の下側信頼限界の 値を算出します。上で求めた値がこの値よりも小さければ(つまり )、有意差があると判定します。これは左片側検定です。

$$ \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))

となるため、グループ間の平均値に有意差があることも確認できました。

線形回帰でモデル化

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)で、競合商品の値引きの影響を受けることが確認できます

一方で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の販売に負に影響していることがわかります。

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

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

データを概観する

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

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

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

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

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

ポアソン回帰でモデル化

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

• 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の傾きが負数となっており、わずかながらですが競合商品の山積み陳列が悪影響を与えることがわかります。

Python: statsmodelsでガンマ回帰

これまで読み進めてきた緑本を改めて復習するために、statsmodelsで一般化線形モデル(GLM)を作り、糖尿病に関するデータの回帰分析を行います。

実装はこちらです。

blog/statsmodels_gamma.ipynb at master · ohke/blog · GitHub

データセット

今回はscikit-learn付属の糖尿病データセットを使います。

import pandas as pd
from sklearn.datasets import load_diabetes

diabetes = load_diabetes()

data = pd.DataFrame(diabetes.data, columns=("age", "sex", "bmi", "map", "tc", "ldl", "hdl", "tch", "ltg", "glu"))
data['y'] = diabetes.target

データは442レコード有り、各レコードは11個の変数を持ちます。
1年後の糖尿病の進行状況が目的変数 (値が大きいほど進行していることを表す)となっており、それ以外は説明変数です。説明変数は全て平均0で正規化されています。この10個の説明変数から1年後の進行状況 を予測する、回帰分析のタスクとなります。

各変数の相関係数を見ます。
との相関係数では、 (BMI)が最も高く、 (おそらく中性脂肪)が続いています。
また、 (年齢)は各変数と緩やかに相関があります。加齢とともに (総コレステロール)や が増加している傾向は見て取れます。

の分布も確認しておきます。平均152、標準偏差77となっています。

GLMの設計

それでは を予測するモデルをGLMで設計します。

GLMでは、確率分布・リンク関数・線形予測子を決める必要がありました (1/25と2/12の投稿も参考にしてください)。

• 0以上で上限の無い連続値のようですので、確率分布にはガンマ分布を使います
• ガンマ分布ですので、リンク関数には対数リンク関数を使います

ガンマ分布は、以下の確率密度関数で表される0以上の連続確率分布です。
はガンマ関数で、 が正の整数なら となります。ガンマ分布の平均は 、分散は です。

$$ p(y \mid s,\mu)=\frac{\exp(-\frac{y}{\mu})}{\Gamma (s)\mu^s}y^{s-1} $$

線形予測子で指定するパラメータは以下の5パターンで比較します。比較基準にはAICを使います(参考)。

1. : パラメータなし(切片のみ)モデル
2. : bmiモデル
3. : bmi+sexモデル
4. : ldl+ageモデル
5. : ldl+age+ldl×ageモデル(交互作用項有り)

まずはパラメータなしモデルで、ベースラインとなるAICを確認します。

• AICは5041
• このモデルの最大尤度は-2519.6 (残差逸脱度は126.8)

# pip install statsmodels
import statsmodels.api as sm
import statsmodels.formula.api as smf

# y~1とすることで、切片のみモデルを作る
model0 = smf.glm('y ~ 1', data=data, family=sm.families.Gamma(link=sm.families.links.log))
result0 = model0.fit(disp=0)

print('AIC:', result0.aic)
# AIC: 5041.186062839457

result0.summary()

次に、bmiモデルで最尤推定させたところ、AICは4884となりました。パラメータ無しのモデルと比較すると、-157改善されています。

model1 = smf.glm('y ~ bmi', data=data, family=sm.families.Gamma(link=sm.families.links.log))
result1 = model1.fit(disp=0)

print('AIC:', result1.aic)
# AIC: 4883.66124171387

もう1つ、bmi+sexモデルでも最尤推定します。AICは4885となり、bmiモデルよりも大きいため、複雑度が増すだけのモデルと言えます。したがって、bmiモデルの方が良い予測を行えるモデルです。

model2 = smf.glm('y ~ bmi+sex', data=data, family=sm.families.Gamma(link=sm.families.links.log))
result2 = model2.fit(disp=0)

print('AIC:', result2.aic)
# AIC: 4885.365291972465

ここまでの3つのモデルの予測値を比較してみます。

• パラメータ無し(model0)は、単一の平均値(152)で推定されています
• bmiモデル(model1)とbmi+sexモデル(model2)はほぼ同じ予測値となっています

次に、ldl+ageモデルとldl+age+ldl×ageモデルから、交互作用項の有無による差を見ていきます。
一般的に、ldl(LDLコレステロール)は加齢に伴って上がる検査値です(相関係数からも緩い正相関が見られます)。つまり、2つは完全に独立したパラメータではありません。そこで交互作用項ldl×ageを導入します。加齢によるLDLコレステロール値の上昇を打ち消す(つまり負の係数となる)ことを期待しています。

2つモデルのAICの差は5.5で、交互作用項の導入による予測精度の向上が期待できます。

model3 = smf.glm('y ~ ldl+age', data=data, family=sm.families.Gamma(link=sm.families.links.log))
result3 = model3.fit(disp=0)
print('AIC:', result3.aic)
# AIC: 5019.939140159227

model4 = smf.glm('y ~ ldl+age+ldl:age', data=data, family=sm.families.Gamma(link=sm.families.links.log))
result4 = model4.fit(disp=0)
print('AIC:', result4.aic)
# AIC: 5014.413561464431

ldl+age+ldl×ageの各パラメータを見ると、 ldlとageは正の係数(それぞれ1.49、1.64)となっていますが、 ldl×ageは負の係数(-28.92)となっており、期待通りの結果となっていることがわかります。

Pythonで実装しながら緑本を学ぶ (第10章 階層ベイズモデル)

データ解析のための統計モデリング入門(通称、緑本)を読み進めています。
述べられている理論を整理しつつ、Rでの実装をPythonに置き換えた際のポイントなども深掘りしていきます。

今回は第10章です。PyMC3を使って階層ベイズモデルを表現します。実装は以下で公開しています。

introduction_to_machine_learning_with_python/chapter10.ipynb at master · ohke/introduction_to_machine_learning_with_python · GitHub

10 GLMのベイズモデル化と事後分布の推定 -GLMMのベイズモデル化-

階層事前分布を使って一般化線形混合モデル(GLMM)を階層ベイズモデルとして扱うことが、本章のテーマです。

10.1 例題：個体差と生存種子数(個体差あり)

1個体あたりの種子数8個を持つ架空植物100個体について、生存種子数を計測した。個体 における生存種子数を とする。

8個の内 個というデータですので、二項分布が使えそうです。
しかし個体に由来するランダム効果によって過分散が生じています。平均4.03に対して分散9.93となっており、二項分布が期待する分散2.00(=8×(4.03/8)×(1-(4.03/8)))よりも、かなり大きいです。

10.2 GLMMの階層ベイズモデル化

第7章(こちらの投稿)では、GLMにこうした原因不明の個体差を組み込んだGLMMを導入し、パラメータを最尤推定で求めた。
パラメータが増えると最尤推定では困難になるため、MCMCサンプリングによる事後分布推定を行いたい。そこで、GLMMのベイズモデル化を行う。

まずリンク関数と線形予測子を決めます。二項分布ですので、リンク関数にはロジットを使います(第6章も参照)。 が個体差です。

$$ logit(q_i)=\beta+r_i $$

8章・9章で追ってきたとおり、事後分布は に比例します。

データ が得られる条件付き確率(尤度)は全個体の二項分布の積です。

$$ p({\bf Y} \mid \beta, {r_i})=\prod_i \binom{8}{y_i} q_i^{y_i}(1-q_i)^{8-y_i} $$

の事前分布は、特に制約がないので平均0・標準偏差100の正規分布に従うこととする(無情報事前分布)。

$$ p(\beta)=\frac{1}{\sqrt{2\pi \times 100^2}} \exp(\frac{-\beta^2}{2 \times 100^2}) $$

個体差 は、個体差のばらつきを表す変数 を導入し、平均0・標準偏差 の正規分布に従うことにする。

$$ p(r_i \mid s)=\frac{1}{\sqrt{2\pi s^2}} \exp(\frac{-r_i^2}{2s^2}) $$

は正であれば良いので、0から10000の連続一様分布に従う無情報事前分布とする。

$$ p(s)=U(0, 10^4) $$

このように、別の事前分布 によって定まる事前分布 を階層事前分布と呼ばれます。 これに対して、 を超パラメータ、 を超事前分布)と言います。

階層事前分布を使うベイズモデルが、階層ベイズです。

10.3 階層ベイズモデルの推定・予測

PyMC3を使い、10.2のモデルを構築します。超事前分布 、階層事前分布 の順に作っていることを確認してください。

• Normalで正規分布、Uniformで連続一様分布、Binomialで二項分布を設定できます
  • rは個体数分必要ですので、shape=len(data.y)としています
  • invlogitでロジット関数の逆関数(つまりロジスティック関数)を二項分布のpに設定しています
• NUTSで1600点(最初の100点は捨てる)をサンプリングしたサンプリング列を3つ作る
  • 合計4500点をサンプリング

import pymc3 as pm

# モデルの作成
with pm.Model() as model:
    # βの事前分布をN(0, 100)の正規分布で設定(無情報事前分布)
    beta = pm.Normal('beta', mu=0, sd=100)
    
    # 超パラメータsの(超)事前分布をU(0, 10000)の連続一様分布で設定(無情報事前分布)
    s = pm.Uniform('s', lower=0, upper=10000)
    
    # パラメータrの事前分布をN(0, s)の正規分布で設定(階層事前分布)
    r = pm.Normal('r', mu=0, sd=s, shape=len(data.y))
    
    # ロジットリンク関数を設定し、二項分布で推定する
    y = pm.Binomial('y', n=8, p=pm.invlogit(beta + r), observed=data.y.values)

# サンプリング
with model:
    # NUTSで101個目からサンプルを取得するチェインを3つ作る
    trace = pm.sample(1500, step=pm.NUTS(), tune=100, njobs=3, random_seed=0)
    
_ = pm.traceplot(trace)
pm.summary(trace).loc[['beta', 's']]

各パラメータの事後分布(95%信用区間)を見ると は-0.600〜0.723、 は0.010〜3.802となっており、いずれも は十分1に近い(収束している)。

の事後分布(下図中段左)を見ると、 毎に異なる分布となっており、個体差が表現されていることがわかる。

パラメータの事後分布から、生存確率を予測してみます。

1. 3つのサンプル列の平均値を求める( と のそれぞれ1500点ずつ)
2. と を1つずつ取り出し、 から100個体分の を正規分布で生成する
3. と から二項分布で種子数の生存確率を求める
4. 1500サンプル・100個体分の生存確率の平均値を計算する

こうして得られた生存確率から、100個体を調べた時の生存種子数の頻度をプロットしたのが以下となります。過分散を予測できていることがわかります。

10.4 ベイズモデルで使うさまざまな事前分布

事前分布の決め方は3パターン。

• ある分布を信じて、主観的に決める (基本的には使わない)
  • 他のデータセットや統計モデルで得られた分布を使う場合など
• 無情報事前分布
• 階層事前分布

無情報事前分布は、個体ごとに自由に決めてしまうと、フルモデルと同等になるデメリットがある。そのため10.2のように、個体ごとに異なるパラメータなどは階層事前分布とし、階層事前分布の(超)パラメータのみを無情報事前分布とするのが基本戦略となる。
したがって、パラメータの影響が局所的な場合は階層事前分布、大域的な場合は無情報事前分布と使い分ける。

10.5 個体差＋場所差の階層ベイズモデル

個体差に加えて、場所差を組み込んだモデルを見ていきます。

架空植物100個体の種子数を計測し、施肥処理と種子数の影響について調べたい、というケースを想定します。
架空植物は10個の植木鉢A〜Jで10個体ずつ植えられてます。植木鉢A〜Eの50個体は施肥処理無し('C')、F〜Jの50個体は施肥処理有り('T')で育てられています。

個体ごとの種子数をプロットしてみますと、施肥処理無し(左側50個体)の方が若干種子数が多いように見えます。

しかし、植木鉢ごとの種子数を箱ひげ図でプロットしてみますと、種子数の分布が異なっており、植木鉢による差が少なからずあるものと思われます。

個体差と植木鉢差を考慮したモデルを定義していきます。

種子数は離散値・上限なしのため、ポアソン分布を使います。
個体 におけるパラメータ は、線形予測子と対数リンク関数でモデルを定義します。

$$ \log \lambda_i = \beta_1 + \beta_2 f_i + r_i + r_{j(i)} $$

• は、平均0・標準偏差100のの正規分布
  • は個体 の施肥処理の有無(有りの場合に1となる)
• は個体差で、平均0・標準偏差 の正規分布 (階層事前分布)
• は植木鉢差で、平均0・標準偏差 の正規分布 (階層事前分布)
  • パラメータ数は10個(は0〜9)
• は0〜10000の連続一様分布

PyMC3での実装は以下のとおりです。

with pm.Model() as model:
    # β1とβ2の事前分布をN(0, 100)の正規分布で設定(無情報事前分布)
    beta1 = pm.Normal('beta1', mu=0, sd=100)
    beta2 = pm.Normal('beta2', mu=0, sd=100)
    
    # sとspの(超)事前分布をU(0, 10000)の連続一様分布で設定(無情報事前分布)
    s = pm.Uniform('s', lower=0, upper=10000)
    sp = pm.Uniform('sp', lower=0, upper=10000)
    
    # 個体差パラメータrの事前分布をN(0, s)の正規分布で設定(階層事前分布)
    r = pm.Normal('r', mu=0, sd=s, shape=len(data.y))
    
    # 場所差パラメータrpの事前分布をN(0, sp)の正規分布で設定(階層事前分布)
    rp = pm.Normal('rp', mu=0, sd=sp, shape=len(data.POT.unique()))
    
    # ログ関数を設定し、ポアソン分布で推定する
    y = pm.Poisson('y', mu=np.exp(beta1 + beta2 * data.F + r + rp[data.POT]), observed=data.y.values)
    
# HMCでサンプリング
with model:
    # 101個目からサンプルを取得するチェインを3つ作る
    trace = pm.sample(1500, step=pm.HamiltonianMC(), tune=100, njobs=3, random_seed=0)

pm.summary(trace).loc[['beta1', 'beta2', 's', 'sp']]

に着目すると、95%信用区間で-2.243〜0.714となっており、0を跨いでいるので施肥処理による効果は見られないようです。