統計学といくつかのよしなしごと

つれづれなるままに、日くらし統計学にむかひて、心にうつりゆくよしなし事を、そこはかとなく書きつくれば、あやしうこそものぐるほしけれ。

時系列分析の理論と実装を概観する①ARIMA過程まで

1.時系列分析における個人的諸混乱

 実務ではしょっちゅう時系列データに出くわすが、これを使って動向の予測やメカニズムの解釈をしたいと思ったときのデータの扱いが難しくて混乱しがちだ。ただ正しく使えば実務上大いに役に立つし、改めて回帰の枠組みを見直し理解するうえでも良い教材になる。ということでこの記事では時系列分析の標準的な手法であるARIMAモデルの利用*1を前提に、分析の基礎的な理論部分と実装手順を概観する。主な参考文献は以下。

経済・ファイナンスデータの計量時系列分析 (統計ライブラリー)

経済・ファイナンスデータの計量時系列分析 (統計ライブラリー)

時系列分析と状態空間モデルの基礎: RとStanで学ぶ理論と実装

時系列分析と状態空間モデルの基礎: RとStanで学ぶ理論と実装

現場ですぐ使える時系列データ分析 ~データサイエンティストのための基礎知識~

現場ですぐ使える時系列データ分析 ~データサイエンティストのための基礎知識~

2.1変量時系列データの分析概要

 時系列データに眺めた際には、「今日の株価は下落したから明日は買いが増えて株価は上昇するだろう*2」とか「直近5年間で売上は順調に伸びているから来年も伸びるだろう」とか言った具合に、過去の実現値が未来の値に影響する(時間依存する)ことを期待する。これを『自己相関』と言って、通常の回帰と枠組みが異なってくる*3。分析においてはこの自己相関をうまくモデル化してやる必要があるが、基本的なアプローチに過去の自分自身に回帰させる方法、 {\displaystyle y_t } {\displaystyle y_{t-1} }の間に共通の成分を含めてやる方法の2通りがある。前者を『自己回帰(AR)過程(autoregressive process)』、後者を『移動平均(MA)過程(moving average process)』という。それぞれの式の形を以下に示しつつ、基本になるいくつかのモデルと押さえるべき諸条件を確認する。

AR(p)過程


{\displaystyle
y_t = c + φ_1 y_{t-1} + φ_2 y_{t-2} + φ_3 y_{t-3} + ... + φ_p y_{t-p} + ε_t, \ \ \ ε_t \ ~ \ W.N.(σ^2)
}
 
 過去p時点までの自分自身の線形和で表現するため、p次のAR過程と言ったりする。ここで {\displaystyle ε_t }は時刻t時点における撹乱項(誤差項)であり、ホワイトノイズ( {\displaystyle } W.N. )と呼ばれる以下を満たす確率分布である。
 {\displaystyle 
E(ε_t) = 0
\\
\\
γ_k = E(ε_t ε_{t-k}) = \left\{ \begin{array}{ll}
    σ^2, & k = 0 \\
    0, & k ≠ 0
\end{array} \right.
}
 一見通常の回帰で誤差項に対して置く仮定と同様のようだが、独立で同一の分布に従うことを要求していない点がポイントになっている。これは実務上この仮定を必ずしも必要としないことや、実際のデータがこれを満たすと考えにくいことが背景にあるよう。
 AR(p)過程の性質についてここでは詳しく述べないが、ポイントとしては自己相関が指数的に減衰する、つまりラグが大きければ大きい程自己相関が弱くなることと、偏自己相関*4がp+1以降(理論上は)0になることである*5

MA(q)過程

 {\displaystyle 
y_t = μ + ε_t + θ_1 ε_{t-1} + θ_2 ε_{t-2} + θ_3 ε_{t-3} + ... + θ_q ε_{t-q}, \ \ \ ε_t \ ~ \ W.N.(σ^2)
}

 こちらのアプローチでは {\displaystyle y_t } を過去の撹乱項の線形和で表現する。時点に限らず撹乱項がホワイトノイズであることに注意すると、q+1以降の自己相関が0になることが分かる。また、後述するがMA(q)過程は特定条件を満たすときAR(∞)過程で書き直すことができ、これにより偏自己相関はラグの大きさに応じて減衰していく。

ARMA(p, q)過程

 AR(p)過程とMA(q)過程を組み合わせたものを『自己回帰移動平均(ARMA)過程(autoregressive moving average process)』という。
 {\displaystyle 
y_t = c + φ_1 y_{t-1} + ... + φ_p y_{t-p} + θ_1 ε_{t-1} + ... + θ_q ε_{t-q}, \ \ \ ε_t \ ~ \ W.N.(σ^2)
}

 ARMA(p, q)はその名の通り両方の過程の性質を併せ持つ過程で、いずれか強いほうの性質に従う。そのため自己相関、偏自己相関はともにラグの大きさに応じて減衰していく。

定常性の仮定

 ところで、これらのモデルは『定常性』の仮定を満たす必要がある。定常性には『弱定常性』と『強定常性』があるが、ここでは弱定常性のみ記載する。強定常性は弱定常性により強い条件を課した概念である(断りなく定常性という場合は普通弱定常性を指す)。

 {\displaystyle 
弱定常性:任意のtとkに対して、
\\
E(y_t) = μ
\\
Cov(y_t, \ y_{t-k}) = E [(y_t - μ)(y_{t-k} - μ) ] = γ_k
 }

 要は弱定常仮定の下では時点を問わず {\displaystyle y_t } の期待値が一定で、かつ自己共分散(よって自己相関も)時点に依存せず、時間差 {\displaystyle k } のみに依存することを課す条件である。この定常性を満たした時系列モデルを基礎として構築したうえで、季節性とかトレンドとか各種の性質をもった実際の時系列データにアプローチしていくことになる。
 定常性を条件に課す背景には、時系列データの特殊性がある。通常の場合ではある確率変数について複数の実現値があるため、それらを用いて平均なり分散なりといった母数を推定するが、時系列データの場合、1つの時点に対して1つの実現値しか存在しない。どういうことかというと、「2018年8月31日の店舗Aの売上は100万円から300万円の可能性があったが(確率変数)、実現値としてはただ1つ120万円が得られた」というようなイメージ。1つの確率変数に対して1つの値しかないため、分散からもう推定できなくなる。そこで、この定常性を課すことで得られた実現値全体からなる一つの過程を考えパラメータの推定を可能にしようというわけである。

定常性と反転可能性の条件

 時系列モデルではまずこの定常性を満たすかどうかの確認から始まるが、①MA過程は常に定常になるため、②AR過程(またはARMA過程のAR過程部分)について定常性を満たすかどうかを検討すればよい。
 ①について。MA(q)過程のモデル式は以下であった。
 {\displaystyle 
y_t = μ + ε_t + θ_1 ε_{t-1} + θ_2 ε_{t-2} + θ_3 ε_{t-3} + ... + θ_q ε_{t-q}, \ \ \ ε_t \ ~ \ W.N.(σ^2)
}

 これについて、
 {\displaystyle 
\begin{eqnarray}
E(y_t) &=& E(μ + ε_t + θ_1 ε_{t-1} + θ_2 ε_{t-2} + θ_3 ε_{t-3} + ... + θ_q ε_{t-q})
\\
&=& E(μ) + E(ε_t) + θ_1 E(ε_{t-1}) + ... + θ_q E(ε_{t-q})
\\
&=& E(μ)
\\
&=& μ
\end{eqnarray}
 }
 {\displaystyle 
\begin{eqnarray}
γ_k = Cov(y_t, \ y_{t-k}) &=& E [ (y_t - μ)(y_{t-k} - μ) ]
\\
&=& E [ (ε_t + θ_1 ε_{t-1} + ... + θ_q ε_{t-q})(ε_{t-k} + θ_1 ε_{t-1-k} + ... + θ_q ε_{t-q-k}) ]
\\
&=& E [ (ε_t + θ_1 ε_{t-1} + ... + θ_k ε_{t-k} + ... + θ_q ε_{t-q})(ε_{t-k} + θ_1 ε_{t-1-k} + ... + θ_{q-k} ε_{t-q} + ... + θ_q ε_{t-q-k}) ]
\\
&=& (θ_k + θ_1 θ_{k+1} + ... + θ_{q-k} θ_q)σ^2
\\
\end{eqnarray}
 }
 {\displaystyle (ただし1 \leq k \leq q, \ \ \ k \geq q+1 ではγ_k = 0) }

 となって、定常性の条件を満たしていることが分かる。

 ②について。AR(p)過程のモデル式は以下であった。

{\displaystyle
y_t = c + φ_1 y_{t-1} + φ_2 y_{t-2} + φ_3 y_{t-3} + ... + φ_p y_{t-p} + ε_t, \ \ \ ε_t \ ~ \ W.N.(σ^2)
}

 これについて、以下で表される『AR特性方程式』のすべての解の絶対値が1より大きくなる( {\displaystyle |z| > 1 } )ときAR過程は定常となる。ここの数理はしっかりとは追っていないが、 {\displaystyle y_t } の分散を有限にするために必要な制約っぽい。

 {\displaystyle 1 - φ_1 z - ... - φ_p z^p = 0 }

 この下でAR(p)過程の期待値と自己共分散を確認しておく。定常であれば {\displaystyle μ = E(y_t) = E(y_{t-1} = ... = E(y_{t-p})) } であることに注意して、期待値は
 {\displaystyle 
\begin{eqnarray}
μ = E(y_t) &=& E(c + φ_1 y_{t-1} + φ_2 y_{t-2} + φ_3 y_{t-3} + ... + φ_p y_{t-p} + ε_t)
\\
&=& E(c) + φ_1 E(y_{t-1}) + φ_2 E(y_{t-2}) + ... + φ_p E(y_{t-p}) + E(ε_t)
\\
&=& c + (φ_1 + φ_2 + ... + φ_p)μ
\\
\end{eqnarray}
 }
 {\displaystyle \therefore μ = E(y_t) = \frac {c}{1 - φ_1 - φ_2 - ... - φ_p} }

一方自己共分散は、
 {\displaystyle 
\begin{eqnarray}
γ_k &=& Cov(y_t, \ y_{t-k})
\\
&=& Cov(c + φ_1 y_{t-1} + φ_2 y_{t-2} + φ_3 y_{t-3} + ... + φ_p y_{t-p} + ε_t, \ y_{t-k})
\\
&=& Cov(c, \ y_{t-k}) + Cov(φ_1 y_{t-1}, \ y_{t-k}) + Cov(φ_2 y_{t-2}, \ y_{t-k}) + ... + Cov(φ_p y_{t-p}, \ y_{t-k}) + Cov(ε_t, \ y_{t-k})
\\
&=& φ_1 γ_1 + φ_2 γ_2 + ... + φ_p γ_p
\end{eqnarray}
 }

 ちなみに、この両辺を分散 {\displaystyle γ_0 = Var(y_t) } で割って自己相関の方程式に変換したもの*6を『ユール・ウォーカー方程式』という。グラフなど描くと自己相関の絶対値が減衰していくことが確認できる。

 ここまで定常性について見てきたが、パラメータの推定という観点でもう一つ気にしなければならないことにMA過程の『反転可能性』というものがある。これはMA過程をAR(∞)に書き直せるかどうかという話で、MA過程においては同じ期待値と自己相関構造を持つパラメータの組が複数ありうるので、その同定のために課す条件である。反転可能、つまり撹乱項 {\displaystyle ε_t }  {\displaystyle y_t } と関数として表現できるならば、過去の {\displaystyle y }  {\displaystyle y_t } を予測したときの予測誤差と解釈できて自然なので、反転可能性を満たすパラメータの組を選択しようという判断基準が与えられる。反転可能条件は以下のMA特性方程式のすべての解の絶対値が1より大きい( {\displaystyle |z| > 1 } )ことである。

 {\displaystyle 1 + θ_1 z + ... + θ_q z^q = 0 }

 このあたりを満たすようなパラメータを導いてやろうというのが、時系列データに対するARMAモデルの推定の考え方となる。

非定常過程と単位根とARIMA過程

 ARMAはデータ系列の定常性の下で推定や予測を行うが、現実のデータは非定常な過程が多く、その場合定常過程とは性質が大きく異なってくる。有名なものに『単位根過程』があり、以下で定義される(「経済・ファイナンスデータの計量時系列分析」P105)。

原系列*7 {\displaystyle y_t } が非定常過程であり、差分系列 {\displaystyle Δy_t = y_t - y_{t-1} } が定常過程であるとき、過程は単位根過程(unit root process)といわれる。

 今回は単位根過程の性質については詳しく述べないが、予測のMSE(平均二乗誤差)が長期的に増大していったり、二変数間のモデル化において『見せかけの回帰*8』を引き起こしたりといろいろ厄介な挙動を示すので追加の処理を検討する必要が出てくる。単位根過程は『1次和分過程』とも言われ、和分過程が何かというと(「経済・ファイナンスデータの計量時系列分析」P106)、

d-1階差分をとった系列は非定常過程であるが、d階差分をとった系列が定常過程に従う過程は、d次和分過程もしくはI(d)過程と呼ばれる。また、I(0)は定常過程で定義される。

 差分系列が定常になるならそれはARMA過程として表現できるので、何階差分をとれば定常になるかを付与したものが『自己回帰和分移動平均(ARIMA)過程』と呼ばれる過程である。定義は以下の通り(「経済・ファイナンスデータの計量時系列分析」P106)。

d階差分をとった系列が定常かつ反転可能なARMA(p, q)過程に従う過程は次数(p, d, q)の自己回帰和分移動平均過程もしくはARIMA(p, d, q)過程と呼ばれる。

 次に、こうした基礎を頭に入れつつ分析の手順とそれに伴う結果を確認してみる。

3.Rによる分析の手順と結果の確認

 Rに用意されているNileデータを使って分析の手順を確認する。Nileは1871年から1970年までの100年間におけるナイル川の流量推移のデータ。原系列のプロットはこんな感じ。

> plot(Nile)

f:id:yanbow221:20180830161535j:plain

 実際の予測時に見立てるために最後の20年分は予測期間として、最初80年分を手元にあるデータとして以降分析していく。

> library(tidyverse)
> library(tseries)
> library(forecast)

> Nile_train <- window(Nile, end = 1950)     #学習データ
> Nile_test <- window(Nile, start = 1951)    #テストデータ
モデルの推定

 はじめに単位根検定を行い原系列が単位根過程かそうでないかを確認する。単位根を持っている場合、例えば差分系列をとってやって定常過程にするといった追加処理を検討する。単位根検定には結構種類があってその詳細は調べていないが、オーソドックスな方法としてここではADF検定を用いる。ADF検定は真のモデルをAR(p)過程として、帰無仮説に「単位根を持つ」、対立仮説に「単位根を持たない」を置いた検定である。p値は0.087で有意水準5%下で帰無仮説を棄却できない、すなわち単位根を持ってないとはいえないなぁという結果になる。

> adf.test(Nile_train) #ADF検定

	Augmented Dickey-Fuller Test

data:  Nile_train
Dickey-Fuller = -3.2443, Lag order = 4, p-value = 0.08675
alternative hypothesis: stationary

 ということで差分系列を作成し、これに対して再び単位根検定を行ってみると、

> Nile_train_diff <- diff(Nile_train, lag = 1)  #差分系列をとる
> adf.test(Nile_train_diff)                     #差分系列に対してADF検定

	Augmented Dickey-Fuller Test

data:  Nile_train_diff
Dickey-Fuller = -5.8527, Lag order = 4, p-value = 0.01
alternative hypothesis: stationary

Warning message:
In adf.test(Nile_train_diff) : p-value smaller than printed p-value

 今度は有意な結果を示したので、この1階差分の系列に対して分析を進める。真の過程が不明な状態なので、各過程を真としたときの自己相関、偏自己相関の性質から仮定するモデルを選択する。「経済・ファイナンスデータの計量時系列分析」P49にある表を抜粋しておく。
f:id:yanbow221:20180830165637p:plain

 これを判断するために差分系列の自己相関、偏自己相関をプロットするとこんな感じ。

> par(mfrow = c(2, 1))
> acf(Nile_train_diff)                     #自己相関プロット
> acf(Nile_train_diff, type = "partial")   #偏自己相関プロット

f:id:yanbow221:20180830175425j:plain

 自己相関(ACF)については1次の自己相関が有意であり、2次以降は減衰というよりは切断されている(0とみなす)気もするが、とりあえずどちらの可能性も持たせておく。偏自己相関(PACF)については1次が有意で2次が微妙だが帰無仮説を棄却できておらず、3次以降は(10次時点が有意になっているが)減衰ないし切断ととらえて、1次までをモデルに盛り込むことを検討する。というわけで過程の候補をAR(1)、MA(1)、ARMA(1, 1)とする。これは差分系列に対してのものなので、それぞれARIMA(1, 1, 0)、ARIMA(0, 1, 1)、ARIMA(1, 1, 1)と同等となる。
 上で選択した候補においてパラメータの推定を行っていく。それぞれ定数項のありなし2種類において推定するので、計6つの過程において推定を行う。どのモデルが優れているかについて今回はAIC(赤池情報量基準)を用いる。

> # パラメータ推定
> model_ar1 <- arima(Nile_train_diff, order = c(1, 0, 0), include.mean = F)      #定数項なしAR(1)
> model_ar1_intercept <- arima(Nile_train_diff, order = c(1, 0, 0))              #定数項ありAR(1)
> 
> model_ma1 <- arima(Nile_train_diff, order = c(0, 0, 1), include.mean = F)      #定数項なしMA(1)
> model_ma1_intercept <- arima(Nile_train_diff, order = c(0, 0, 1))              #定数項ありMA(1)
> 
> model_arma1.1 <- arima(Nile_train_diff, order = c(1, 0, 1), include.mean = F)  #定数項なしARMA(1, 1)
> model_arma1.1_intercept <- arima(Nile_train_diff, order = c(1, 0, 1))          #定数項ありARMA(1, 1)

> # AICの比較
> result_aic <- as.matrix(c(model_ar1$aic, model_ar1_intercept$aic,
+                           model_ma1$aic, model_ma1_intercept$aic,
+                           model_arma1.1$aic, model_arma1.1_intercept$aic), ncol = 1)
> rownames(result_aic) <- c("AR(1)定数項なし", "AR(1)定数項あり",
+                           "MA(1)定数項なし", "MA(1)定数項あり",
+                           "ARMA(1,1)定数項なし", "ARMA(1,1)定数項あり")
> colnames(result_aic) <- "AIC"
> result_aic
                         AIC
AR(1)定数項なし     1027.300
AR(1)定数項あり     1029.236
MA(1)定数項なし     1018.066
MA(1)定数項あり     1019.439
ARMA(1,1)定数項なし 1016.574
ARMA(1,1)定数項あり 1015.936

 AICが最小になるモデルはARMA(1, 1)過程となった。このモデルの下での推定パラメータは以下の通り。

> model_arma1.1_intercept

Call:
arima(x = Nile_train_diff, order = c(1, 0, 1))

Coefficients:
         ar1      ma1  intercept
      0.3396  -0.9997    -4.2092
s.e.  0.1074   0.0664     1.0029

sigma^2 estimated as 19426:  log likelihood = -503.97,  aic = 1015.94

 実はここまでの流れは、Rにおいてはパッケージforecastに入っているauto.arimaが一発でやってくれる。auto.arimaを使うと総当たり的にモデルを推定し、指定した情報量基準の下で最適なモデルを選択してくれる。そのまま原系列に対してやってみるとこんな感じ。

> model_nile <- auto.arima(Nile_train, ic = "aic", trace = T, stepwise = F, approximation = F)

 ARIMA(0,1,0)                    : 1038.336
 ARIMA(0,1,0) with drift         : 1040.313
 ARIMA(0,1,1)                    : 1018.066
 ARIMA(0,1,1) with drift         : 1019.439
 ARIMA(0,1,2)                    : 1017.266
 ARIMA(0,1,2) with drift         : 1017.486
 ARIMA(0,1,3)                    : 1018.316
 ARIMA(0,1,3) with drift         : Inf
 ARIMA(0,1,4)                    : 1020.194
 ARIMA(0,1,4) with drift         : Inf
 ARIMA(0,1,5)                    : 1022.184
 ARIMA(0,1,5) with drift         : Inf
 ARIMA(1,1,0)                    : 1027.3
 ARIMA(1,1,0) with drift         : 1029.236
 ARIMA(1,1,1)                    : 1016.574
 ARIMA(1,1,1) with drift         : Inf
 ARIMA(1,1,2)                    : 1018.371
 ARIMA(1,1,2) with drift         : Inf
 ARIMA(1,1,3)                    : 1020.233
 ARIMA(1,1,3) with drift         : Inf
 ARIMA(1,1,4)                    : Inf
 ARIMA(1,1,4) with drift         : Inf
 ARIMA(2,1,0)                    : 1025.79
 ARIMA(2,1,0) with drift         : 1027.705
 ARIMA(2,1,1)                    : 1018.308
 ARIMA(2,1,1) with drift         : Inf
 ARIMA(2,1,2)                    : 1020.333
 ARIMA(2,1,2) with drift         : Inf
 ARIMA(2,1,3)                    : Inf
 ARIMA(2,1,3) with drift         : Inf
 ARIMA(3,1,0)                    : 1025.755
 ARIMA(3,1,0) with drift         : 1027.631
 ARIMA(3,1,1)                    : 1020.253
 ARIMA(3,1,1) with drift         : Inf
 ARIMA(3,1,2)                    : Inf
 ARIMA(3,1,2) with drift         : Inf
 ARIMA(4,1,0)                    : 1025.598
 ARIMA(4,1,0) with drift         : 1027.434
 ARIMA(4,1,1)                    : 1022.102
 ARIMA(4,1,1) with drift         : Inf
 ARIMA(5,1,0)                    : 1025.951
 ARIMA(5,1,0) with drift         : 1027.72



 Best model: ARIMA(1,1,1)

 こちらではドリフト率(先ほどまで述べていた定数項に当たるもの)なしのARIMA(1, 1, 1)が提案された。先ほどは定数項ありのARMA(1, 1)、すなわちドリフト率ありのARIMA(1, 1, 1)を選択したが、auto.arimaでは結果が「Inf」となっている。これはどうもAR過程の定常性やMA過程の反転可能性が危うい、つまり特性方程式の解が1より大きくなるかどうか怪しいときに出るようで、実際過程における特性方程式の解を確認してみると、反転可能性がちょっと怪しいという結果が確認できる。

> abs(polyroot(c(1, -coef(model_arma1.1_intercept)["ar1"])))  #差分系列における定数項ありARMA(1,1)のAR特性方程式の解
[1] 2.944355
> abs(polyroot(c(1, coef(model_arma1.1_intercept)["ma1"])))   #差分系列における定数項ありARMA(1,1)のMA特性方程式の解
[1] 1.000263
> abs(polyroot(c(1, -coef(model_nile)["ar1"])))               #原系列におけるドリフト率なしARIMA(1,1,1)のAR特性方程式の解
[1] 3.768194
> abs(polyroot(c(1, coef(model_nile)["ma1"])))                #原系列におけるドリフト率なしARIMA(1,1,1)のMA特性方程式の解
[1] 1.150311

 実務ではauto.arimaに頼ることが多いと思うので、こちらで選択されたモデル、原系列におけるARIMA(1, 1, 1)過程について以降のプロセスを進めていく。モデルの定常性、反転可能性は既に確認したので、次に残差について確認する。選択したモデルで原系列の自己相関が十分に説明されているか?という意味で残差 {\displaystyle e_t = y_t - \hat{y_t} } の自己相関が0になっていることを確認する。また、モデルの撹乱項に置いたホワイトノイズ自体は正規分布を仮定しないが、一般に撹乱項は正規性を仮定することが多い(「時系列分析と状態空間モデルの基礎」p30, 68)ので、併せて {\displaystyle e_t } の正規性も検定する。

> checkresiduals(model_nile)             #残差の自己相関の検定

	Ljung-Box test

data:  Residuals from ARIMA(1,1,1)
Q* = 9.6295, df = 8, p-value = 0.292

Model df: 2.   Total lags used: 10

> jarque.bera.test(resid(model_nile))    #残差の正規性の検定

	Jarque Bera Test

data:  resid(model_nile)
X-squared = 0.42784, df = 2, p-value = 0.8074

f:id:yanbow221:20180830182732j:plain

 10期の自己相関に怪しい値が出ているが、 {\displaystyle p > 0.05 } となって有意な自己相関はないことからひとまずOKとする。同様に正規性も問題なさそうなので、今回はこのモデルを採用する。

モデルによる予測

 作成したモデルでテスト期間20期分について予測すると、こんな感じになる。

> # テストデータ期間の予測
> forecast_nile <- forecast(
+   model_nile,
+   h = 20,
+   level = c(95, 80)
+ )
> autoplot(forecast_nile) +
+   theme_bw()

f:id:yanbow221:20180831132507j:plain

 差分系列が定常過程であり、AR過程、MA過程ともに次数が1であることもあって、3期後あたりから横ばいの予測値になっている。定常過程は平均回帰性を持っており、差分系列が定常過程であり定数項がないことによって、差分は徐々に0に近づいていく。そのため原系列も一定値で推移するようになる。一方信頼区間は期が増えるに従い線形に増大しており、これは原系列が単位根を持っていることによる特徴である(定常過程では分散は一定値に収束していく)。この2点から、素朴なARIMAが長期的な予測には向かないと言える。非定常過程のその他の手法として状態空間モデル等あるが、ARIMAを基礎にした工夫としては例えば以下が参考になる。先に紹介した「時系列分析と状態空間モデルの基礎」の著者の方のブログで、時系列予測に関して各記事とても勉強になる。本記事を書くにあたっても書籍と併せて大いに参考にした。今回は分析の大まかな流れを追うため、以降の工夫は省略する。

logics-of-blue.com

 さて、次に予測の良し悪しを評価するが、そのためには比較対象が必要になる。上回るべき最初の比較対象としては「学習データの平均値」や「学習データの最終時点の値」などがあり、こうした原系列からの素朴な予測値を『ナイーブな予測』という(「時系列分析と状態空間モデルの基礎」P113)。評価指標としてはMSE(平均二乗誤差)がよく用いられる。パッケージforecastに入っているaccuracyを使えばこれの平方根をとったRMSE(平均平方二乗誤差)をはじめいくつかの指標を計算してくれる。可視化と併せてこんな感じ。

> plot(Nile_train, xlim = c(1871, 1970))
> lines(Nile_test, lty = "dashed")
> lines(forecast_nile$mean, col = "steelblue", lwd = 2)
> lines(naive_f_mean$mean, col = "purple", lwd = 2)
> lines(naive_f_latest$mean, col = "rosybrown", lwd = 2)
> legend("topright", 
+        legend = c("ARIMA(1,1,1)", "学習データ平均", "学習データ最終点"),
+        col = c("steelblue", "purple", "rosybrown"),
+        lty = "solid", lwd = 2)

f:id:yanbow221:20180831140332j:plain

> accuracy(forecast_nile, x = Nile_test)    #ARIMA(1,1)による予測の評価
                    ME     RMSE      MAE        MPE     MAPE      MASE        ACF1 Theil's U
Training set -18.24538 143.2728 111.4261 -4.3359424 12.97822 0.8310673 -0.03829202        NA
Test set      22.29061 125.0452 106.1632  0.6760674 11.96723 0.7918141  0.19619479  0.869674
> accuracy(naive_f_mean, x = Nile_test)     #ナイーブ予測①(学習データの平均)の評価
                        ME     RMSE      MAE       MPE     MAPE      MASE      ACF1 Theil's U
Training set  4.546749e-14 176.4544 148.3694 -3.840484 16.70686 1.1066069 0.5301332        NA
Test set     -5.287500e+01 133.3132 108.0125 -8.046549 13.15615 0.8056068 0.1922875 0.9354761
> accuracy(naive_f_latest, x = Nile_test)   #ナイーブ予測②(学習データの最終点)の評価
                     ME     RMSE      MAE       MPE     MAPE      MASE       ACF1 Theil's U
Training set  -2.911392 170.7412 134.0759 -2.198271 15.14585 1.0000000 -0.3925129        NA
Test set     -12.950000 123.0624 101.9500 -3.407725 11.96265 0.7603899  0.1922875 0.8623788

 RMSEについて、ナイーブ予測①よりは小さいもののナイーブ予測②よりも大きい、つまり素朴な予測よりもイマイチな予測になっていることが分かった。ここからは前述したような工夫もしつつ、精度向上の手を考えていくことになる。

4.終わりに

 今回は最もよく用いられるだろうARIMAによる時系列分析について概観した。その他季節性や外因性を考慮したモデルや多変量への拡張、状態空間モデル等、整理してそのうちまとめたい。

*1:AR、MA、ARMAを内包する

*2:株のことはよく知らないのでこのような判断があるのかどうかわからないです。完全にたとえとしての妄想

*3:通常の回帰では同一変量内の値の間に独立を仮定する

*4:時点tと時点t-kにおいて、その間の時点の影響を取り除いた自己相関のこと

*5:p時点前までの自分自身で回帰するのでp+1時点以降影響がなくなるのは直感的に分かりやすいと思う

*6:自己相関 {\displaystyle ρ_k = γ_k / γ_0 }

*7:関心を持っている時系列データそのもののこと

*8:明らかに関係のない変数間に有意な関係を示してしまうこと

コレポンは何をやっているのか②検討の流れ

前回

 前回の記事でマーケティングリサーチにおいて多用されるコレスポンデンス分析について、その可視化までの数理を概観した。今回は出力されたマップやその他の統計量を用いた検討のポイントを見ていく。前回同様主な参考文献は以下。

対応分析入門 原理から応用まで 解説◆Rで検算しながら理解する

対応分析入門 原理から応用まで 解説◆Rで検算しながら理解する

「対応分析法・数量化Ⅲ類の考え方」
http://www.wordminer.org/wp-content/uploads/2013/04/63_18.pdf

f:id:yanbow221:20180717203905j:plain
前回出力した同時布置図

固有値と寄与率、累積寄与率

 コレスポンデンス分析によってサンプル、銘柄ともに5次元の成分スコアが得られ、上記の同時布置図ではそのうち第1成分、第2成分を軸に取った。このように取っ掛かりとしては出力された成分の順に沿って2軸を取ることが多いが、これは対応する固有値の大きさに由来している。5成分の固有値を確認してみると、

> result_Xi$values[1:5] # 行方向に主成分分析を行って得た固有値
[1] 0.62603096 0.18766678 0.13451054 0.04519841 0.02603776
> result_Xj$values[1:5] # 列行方向に主成分分析を行って得た固有値
[1] 0.62603096 0.18766678 0.13451054 0.04519841 0.02603776

行方向、列方向それぞれで得た固有値の値が一致している(6成分以降の固有値はほとんど0になっている)。解いた固有値問題が元データに対していくつか操作を加えて得たデータの分散共分散行列に対するものだったことを思い出すと、この固有値の総和はデータが持っていた分散の総和に等しいことが分かる。実際それぞれ計算してみると、

> sum(diag(Vxi)) # 元のデータの分散の総和(行方向)
[1] 1.019444
> sum(diag(Vxj)) # 元のデータの分散の総和(列方向)
[1] 1.019444
> sum(result_Xi$values) # 固有値の総和(行方向)
[1] 1.019444
> sum(result_Xj$values) # 固有値の総和(列方向)
[1] 1.019444

となっていてこれが確認できる。なお、分散の総和と固有値の総和が等しくなるメカニズムについては、例によって以下で丁寧に解説してある。

多変量解析法入門 (ライブラリ新数学大系)

多変量解析法入門 (ライブラリ新数学大系)

 このように、元のデータが持つ分散をできるだけ表現できる軸を順に求めなおしてやって作った空間上に再配置するのが主成分分析である。というわけでこの固有値を使って、元のデータをどのくらい表現できているかを評価するための統計量が定義される。


{\displaystyle 
v_k = \frac {λ_k}{\sum^K_{k=1} λ_k} × 100(\%) \ \ \ k = 1,2,...K \ \ \ (第k成分の寄与率)
\\
\\
v_{k'} = \frac {\sum^{k'}_{k=1} λ_k}{\sum^K_{k=1} λ_k} × 100(\%) \ \ \ 1 \leq k' \leq K \ \ \ (第k'成分までの累積寄与率)
}

 今回は第1,2成分までで分散の約80%を説明できていることになるので、これらで軸をとってやってデータの意味付けを検討していくのが基本になる。このように元のデータの持つデータ項目が多く考えにくい時に、データの持つ情報をなるべく担保しつつ考えやすいより少ない次元に落とし込んでやろう*1というのが主成分分析とかコレスポンデンス分析の意義である。軸の選択については「累積寄与率がX%以上(領域によって慣例的に囁かれている)」とか「スクリープロット*2」の利用等いろいろと成分の採否に関する考え方があるが、コレスポンデンス分析の場合はマッピングに基づく解釈がより重要になるので、寄与率の大きな軸だけで見ずにいろいろと見てみるのが良いと思う*3

同時布置図の検討

 そんな形で作成した同時布置図を見て、あれこれ「軸」と「点」について解釈していくのが次の段階になる。検討に際し有用な統計量は大きく2つで、「ある軸の形成に対して各々の行項目*4や列項目*5がどの程度寄与しているか」を示す絶対寄与度と、「各々の行項目/列項目が各軸(成分)によってどの程度表現されているか」を示す平方相関(相対寄与度)がある。

(1)絶対寄与度
 各行項目への絶対寄与度 {\displaystyle C_k(i) }、各列項目への絶対寄与度 {\displaystyle C_k(j) }はそれぞれ以下で計算できる。

{\displaystyle 
C_k(i) = \frac {p_{i+}(z_{ik})^2}{λ_k}, \ \ \ \sum^{m}_{i=1} C_k(i) = 1
}


{\displaystyle 
C_k(j) = \frac {p_{+j}(z_{jk}^*)^2}{λ_k}, \ \ \ \sum^{n}_{j=1} C_k(j) = 1
}


{\displaystyle 
\left( i \in I,\ \  j \in J,\ \  k = 1,2,...,K,\ \ K = min\{m, n\} - 1 \right)
}

 それぞれ求めてみるとこんな感じ。

> round((Pi %*% (rscore2 ^ 2) %*% diag(1 / result_Xj$values[1:5])) * 100, 2)
       [,1]  [,2]  [,3]  [,4]  [,5]
 [1,]  0.47  4.85 34.84  3.34 10.87
 [2,]  9.26 12.15 30.48 18.59  7.57
 [3,] 14.67  5.67  8.17  1.56  0.85
 [4,]  3.16  3.87  0.11  8.16 23.86
 [5,]  8.58  1.36  1.06 17.63 33.29
 [6,]  1.30 13.29 15.36  9.15  0.05
 [7,] 18.25  4.65  0.21 14.78 11.34
 [8,] 14.67  5.67  8.17  1.56  0.85
 [9,]  7.22 17.15  0.99 24.64  4.36
[10,] 22.43 31.35  0.60  0.61  6.95
> round((Pj %*% (cscore2 ^ 2) %*% diag(1 / result_Xj$values[1:5])) * 100, 2)
      [,1]  [,2]  [,3]  [,4]  [,5]
[1,] 19.53  0.01  4.77 46.71 15.20
[2,]  0.44 32.65  4.12 22.15 19.95
[3,] 17.61  0.01 61.52  1.47  2.14
[4,] 18.94 49.50  0.18 23.97  0.51
[5,]  8.68  4.14  4.82  0.66 61.02
[6,] 34.80 13.69 24.60  5.04  1.18

 例えば第1成分の軸については、縦方向に見れば順にサンプル10, 7, 3等が、横方向に見れば順に銘柄F, D, C等が軸の形成に寄与していることが分かる。同様にして第2成分の軸については、縦方向に見れば順にサンプル10, 9, 6, 2等が、横方向に見れば順に銘柄D, B, F等が寄与している。一方でサンプル1や2は第1成分、2成分には寄与しておらず、むしろ第3成分の形成に寄与していることが分かる。銘柄Cも第3成分の形成において大きく寄与している。こういったことを確認しながら、必要に応じて軸の解釈も行いつつデータ項目の関係性を確認していく。

(2)平方相関(相対寄与度)
 各行項目に対する平方相関 {\displaystyle C^*_k(i) }、各列項目に対する平方相関 {\displaystyle C^*_k(j) }はそれぞれ以下で計算できる。

{\displaystyle 
C^*_k(i) = \frac {z_{ik}^2}{ \sum^{n}_{j=1} p_{+j} \left( \frac {p_{ij} - p_{i+} p_{+j}}{p_{i+}p_{+j}} \right) }
}


{\displaystyle 
C^*_k(j) = \frac {(z_{ik}^*)^2}{ \sum^{m}_{i=1} p_{i+} \left( \frac {p_{ij} - p_{i+} p_{+j}}{p_{i+}p_{+j}} \right) }
}


{\displaystyle 
\left( i \in I,\ \  j \in J,\ \  k = 1,2,...,K,\ \ K = min\{m, n\} - 1 \right)
}

 こんな感じで求まる。

> round(diag(1 / apply((rscore2 ^ 2), 1, sum)) %*% (rscore2 ^ 2) * 100, 2)
       成分スコア
        第1成分 第2成分 第3成分 第4成分 第5成分
   [1,]    4.61   14.40   74.13    2.39    4.48
   [2,]   43.88   17.26   31.02    6.36    1.49
   [3,]   80.28    9.30    9.61    0.61    0.19
   [4,]   53.31   19.59    0.39    9.95   16.76
   [5,]   72.27    3.43    1.92   10.72   11.66
   [6,]   14.08   43.08   35.68    7.14    0.02
   [7,]   85.97    6.57    0.21    5.02    2.22
   [8,]   80.28    9.30    9.61    0.61    0.19
   [9,]   49.68   35.36    1.47   12.24    1.25
  [10,]   69.46   29.11    0.40    0.14    0.89
> round(diag(1 / apply((cscore2 ^ 2), 1, sum)) %*% (cscore2 ^ 2) * 100, 2)
      成分スコア
       第1成分 第2成分 第3成分 第4成分 第5成分
  [1,]   79.51    0.02    4.17   13.73    2.57
  [2,]    3.25   72.28    6.54   11.81    6.13
  [3,]   56.76    0.01   42.60    0.34    0.29
  [4,]   53.25   41.72    0.11    4.87    0.06
  [5,]   64.11    9.16    7.64    0.35   18.74
  [6,]   78.02    9.20   11.85    0.82    0.11

 サンプル3, 7, 8、銘柄1や6等は第1成分によってよく表現されていることがわかる。ここで注意しないといけないのは、例えばサンプル1を同時布置上で見てみると真ん中あたりにプロットされているが、これを第1成分、第2成分で定義された空間上での「平均」と解釈するのは誤りである。サンプル1の第1成分、第2成分への絶対寄与度は0.47%, 4.85%となっておりそもそも軸の形成に寄与していない。また平方相関を見ても第1成分、第2成分それぞれ4.61%, 14.40%でありこれらの軸はサンプル1にそれ程強く寄与していないことが分かる。サンプル1においてはむしろ第3成分との関連が強くあり、これはつまり第1成分、第2成分で描かれた空間ではサンプル1は「表現されていない」と解釈した方が良さそうということになる。
 このように、同時布置図でプロットされた全ての点の位置関係をもって素朴に解釈すると思わぬ誤りを犯してしまうことがあるので、求めた成分軸それぞれで同時布置図と統計量を確認することが重要となる。以下は「対応分析法・数量化Ⅲ類の考え方」から引用。

ここで留意すべきことは、散布図の情報は"そこでいま眺めている成分軸の組み合わせの中での射影図"であり多次元データとしての全情報ではない、ということである。よって、絶対寄与度、相対寄与度も図と併せて観察し、どの成分でどう寄与するかを多次元的に知ることが必要である。

その他いくつかの補足

(1)サプリメンタリーポイント
 対応分析は行方向の分析と列方向の分析を遷移方程式を使って関係づけるが、この式によって追加処理を行いデータの解釈を助けることができる。例えば先の10サンプルについて実は「性別」という属性があったとする。つまりデータセットがこんな感じであったとする。

> sex_sample <- matrix(c(1, 0, 1, 1, 1, 0, 0, 1, 0, 0,
+                        0, 1, 0, 0, 0, 1, 1, 0, 1, 1), nrow = 10, ncol = 2, byrow = F)
> colnames(sex_sample) <- c("男性", "女性")
> (F_col_added <- cbind(F, sex_sample))
           銘柄A 銘柄B 銘柄C 銘柄D 銘柄E 銘柄F 男性 女性
サンプル1      0     1     0     0     1     1    1    0
サンプル2      0     0     0     0     0     1    0    1
サンプル3      0     0     1     0     0     1    1    0
サンプル4      0     1     1     0     1     1    1    0
サンプル5      0     1     1     0     0     1    1    0
サンプル6      1     1     1     0     1     0    0    1
サンプル7      1     1     0     1     1     0    0    1
サンプル8      0     0     1     0     0     1    1    0
サンプル9      1     1     0     0     1     0    0    1
サンプル10     1     0     0     1     1     0    0    1

 遷移方程式を使うことで、銘柄項目に対して算出した成分スコアを担保した状態で新たに性別項目の成分スコアを求めることができる。図と併せてこんな感じ。

> # 追加処理(サプリメンタリーポイント)
> P_col_added <- F_col_added / sum(F_col_added)
> pj_col_added <- matrix(apply(P_col_added, 2, sum))
> C_col_added <- diag(1 / pj_col_added[, 1]) %*% t(P_col_added)
> Dm1 <- diag(1 / sqrt(result_Xi$values[1:5]))
> cscore_col_added <- C_col_added %*% rscore2 %*% Dm1 # 遷移方程式
> 
> # 銘柄の成分スコアが担保されていることの確認
> cscore2 # 元々算出した成分スコア
       成分スコア
銘柄       第1成分     第2成分     第3成分     第4成分     第5成分
  銘柄A -0.9413812 -0.01306953 -0.21557326  0.39122449  0.16936819
  銘柄B -0.1153204  0.54418081  0.16369085 -0.21996465  0.15846620
  銘柄C  0.7996653  0.01204477 -0.69280014 -0.06201891 -0.05690933
  銘柄D -1.3113409 -1.16064830 -0.05848639 -0.39637626  0.04373046
  銘柄E -0.5125064  0.19370276  0.17693374  0.03799750 -0.27710881
  銘柄F  1.0261401 -0.35232510  0.39991984  0.10495867  0.03857811
> cscore_col_added # 性別を付与したデータセットについて算出した成分スコア
           [,1]        [,2]        [,3]        [,4]        [,5]
[1,] -0.9413812 -0.01306953 -0.21557326  0.39122449  0.16936819
[2,] -0.1153204  0.54418081  0.16369085 -0.21996465  0.15846620
[3,]  0.7996653  0.01204477 -0.69280014 -0.06201891 -0.05690933
[4,] -1.3113409 -1.16064830 -0.05848639 -0.39637626  0.04373046
[5,] -0.5125064  0.19370276  0.17693374  0.03799750 -0.27710881
[6,]  1.0261401 -0.35232510  0.39991984  0.10495867  0.03857811
[7,]  0.9035441 -0.04731066 -0.11472595 -0.33848483 -0.25003063
[8,] -0.4252808 -0.38593508  0.42217115  0.77741483  0.43181891
> 
> # 追加処理後の同時布置図
> library(gridExtra)
> 
> dimnames(cscore_col_added) <- list(銘柄 = c("銘柄A", "銘柄B", "銘柄C", "銘柄D", "銘柄E", "銘柄F", "男性", "女性"),
+                             成分スコア = c("第1成分", "第2成分", "第3成分", "第4成分", "第5成分"))
> dat_map_col_added <- rbind(rscore2, cscore_col_added)
> dat_map_col_added <- transform(dat_map_col_added, label = c(rep("サンプル", 10), rep("銘柄", 6), rep("性別", 2)))
> 
> map2 <- ggplot(dat_map_col_added, aes(x =1成分, y =2成分, color = label)) + 
+   xlim(-1.5, 1.5) + ylim(-1.5, 1.5) +
+   geom_hline(yintercept = 0, color = "black") +
+   geom_vline(xintercept = 0, color = "black") +
+   geom_point(size = 5) +
+   scale_color_manual(
+     values = c("サンプル" = "tomato", "銘柄" = "skyblue", "性別" = "mediumseagreen")
+   ) +
+   geom_text(aes(label = rownames(dat_map_col_added)), vjust = -1, size = 5) +
+   theme_bw()
> 
> grid.arrange(map1, map2, ncol = 2) # 左図:元の同時布置図/右図:性別を追加処理した同時布置図

f:id:yanbow221:20180805122244j:plain
左図:元の同時布置図/右図:性別を追加処理した同時布置図

 これはとても便利で、今のように属性を追加して図の解釈を深堀りたいときや、外れ値があった際にそれを一時的に除外したり再配置したりして外れ値の影響を見てみるなどの使い方が考えられる。

(2)馬蹄形効果
 ところで同時布置図をよく見ると上に凸な放物線上に点がプロットされていることに気づく。これは対応分析でよく見られる現象で、馬蹄形効果と言われる。「対応分析入門」のP30によると、

この効果は第1軸が非常に支配的なときによく起こり、第2軸は、第1軸の線形変換になっている。

 とある。馬蹄形効果自体はそれ程問題としないことが多いようだけれど*6、例えば年代のような順序関係のあるデータが馬蹄形に同時布置上を推移しているとき、その変化は直線的でなく曲線的であるということを気に留めておく、というような注意はした方がよさそう。

終わりに

 生物系で発展した手法として、対応分析と回帰分析を結び付けた正準対応分析(CCA;Canonical Correspondence Analysis)というものがある。いろいろ使い道がありそうな気がするのでそのうち検討してみたい。

*1:縮約という

*2:各軸の持つ寄与率を可視化したもの。寄与率の大きさに加え変化量(軸間の崖の大きさ)の観点も含めて検討するための図

*3:実際単純に寄与率の大きい軸だけで考えると思わぬ間違った解釈をしてしまうことがある

*4:今回で言えばサンプル

*5:今回で言えば銘柄

*6:馬蹄形を修正するような手法も提案されているみたい

コレポンは何をやっているのか①数理の概要

告解

 ここ数か月とても忙しく仕事をしていて*1、前の投稿から半年以上経ってしまった。とはいえ何よりも自分の知識や経験の定着を一番の目的にブログを書いているので、これだけ空いたのは怠慢でしかないんだよなぁ。。もう少し頑張ろうと思います。

マーケティングリサーチでよくやるいわゆるコレポン

 気を取り直して。マーケティングリサーチではコレスポンデンス分析*2、いわゆるコレポンがよく用いられる。マーケティングはもちろん、戦略や戦術を考えるといった際に僕たちは二軸に切ったマップを頻繁に利用する。会社であれ製品・サービスであれ顧客ニーズであれ、切れ味の良さそうな軸を定義してその上にそれらを置いてやれば、何となく洞察が得られた気になる。実際関心のある状況を簡略化、可視化することは思考の大きな助けになるし、その点マッピングはとても有用である。しかし経験的直感だけでマッピングするのは心許ないから、データに基づいて確からしマッピングしたい。これに応える手法の一つがコレスポンデンス分析である。
 例として、あるモデル事務所が目新しいモデルを打ち出したいと思っており、目新しさとは何か検討するために従来のモデルの外見的特徴をマッピングする、という架空のストーリーを考える。外見の9割は目と髪の色で決まると考え、各モデルにこれらの値を付与したデータセットを作成した*3。これに対しコレスポンデンス分析は以下のようなマップを出力する。図の黒文字が目の色、赤文字が髪の色である。

> caith
       fair red medium dark black
blue    326  38    241  110     3
light   688 116    584  188     4
medium  343  84    909  412    26
dark     98  48    403  681    85
> # 表側:目の色
> # 表頭:髪の色
> # 各セルの値は該当する人数

f:id:yanbow221:20180716230923j:plain

 これを見ると目の色と髪の色の組み合わせとして[dark , black][blue or light , fair]などのモデルが多いことが分かる。プロットされた点から何となくモデルの外見的特徴は「暗め」「明るめ」「標準(medium)」の3タイプに分かれそうな気がするし、少なくとも横軸については「明るい-暗い」といった意味付けをしたくなる。このようにコレスポンデンス分析は関心のあるデータ項目を数理的にプロットしてくれるため、ビジネスを検討する上での重要かつ使い勝手の良いツールとなるのである。
 しかし、そもそもこの横軸縦軸は何を表しているのか、なぜその2軸が選ばれたのか、なぜ異なる意味合いを持つ項目を同じ図に布置できるのか、軸に意味付けをすることは正当か、正当ならば意味付けに際して注意するべきことはあるか、などなどの疑問が湧く。分かりやすい可視化で容易に利用できる手法であるからこそ、分析者は手法に対する理解を深め、その後のビジネス上の検討を誤った方向にリードしてしまわないように注意する必要があるよなぁ。と思って以下の本で一通り勉強しなおしたので備忘としてブログに残しておくことにする。

対応分析入門 原理から応用まで 解説◆Rで検算しながら理解する

対応分析入門 原理から応用まで 解説◆Rで検算しながら理解する

 Rでの検算の章もあって理解の助けになるが、文章に癖があるためかややとっつきにくいので以下リンクの「対応分析法・数量化Ⅲ類の考え方」等で補完するのが良いかも。
http://www.wordminer.org/wp-content/uploads/2013/04/63_18.pdf

コレスポンデンス分析の数理の概要

 コレスポンデンス分析は対応分析とも言われるように、先の例で見たように表側の項目と表頭の項目を対応させ一つの図に同時に布置する分析手法である。イメージはこんな感じ。
f:id:yanbow221:20180717102910j:plain

  1. 行、列それぞれの方向に対してデータを変換する
  2. 変換したデータについて分析してそれぞれプロットする
  3. 重ね合わせる

 3ステップで直感的にはとても分かりやすい。では具体的にどのようにこれを行っているのか実際に計算しながら見ていくことにする。

データセット

 上でリンクしたPDF「対応分析法・数量化Ⅲ類の考え方」に記載されている以下の架空データを使う。10サンプルについてある商品の「好きな銘柄」を選んでもらったという設定。なので「好き=1, 好きでない=0」の二値データになっている。

## データセット
F <- matrix(c(0, 0, 0, 0, 0, 1, 1, 0, 1, 1,
              1, 0, 0, 1, 1, 1, 1, 0, 1, 0,
              0, 0, 1, 1, 1, 1, 0, 1, 0, 0,
              0, 0, 0, 0, 0, 0, 1, 0, 0, 1,
              1, 0, 0, 1, 0, 1, 1, 0, 1, 1,
              1, 1, 1, 1, 1, 0, 0, 1, 0, 0),
            nrow = 10, ncol = 6, byrow = F)
dimnames(F) <- list(サンプル = c("サンプル1", "サンプル2", "サンプル3", "サンプル4", "サンプル5",
                             "サンプル6", "サンプル7", "サンプル8", "サンプル9", "サンプル10"),
                        銘柄 = c("銘柄A", "銘柄B", "銘柄C", "銘柄D", "銘柄E", "銘柄F"))

> F
            銘柄
サンプル     銘柄A 銘柄B 銘柄C 銘柄D 銘柄E 銘柄F
  サンプル1      0     1     0     0     1     1
  サンプル2      0     0     0     0     0     1
  サンプル3      0     0     1     0     0     1
  サンプル4      0     1     1     0     1     1
  サンプル5      0     1     1     0     0     1
  サンプル6      1     1     1     0     1     0
  サンプル7      1     1     0     1     1     0
  サンプル8      0     0     1     0     0     1
  サンプル9      1     1     0     0     1     0
  サンプル10     1     0     0     1     1     0
ステップ1:行、列それぞれの方向に対してデータを変換する

 まずはマッピングするためのデータを行、列それぞれに作っていくステップ。ここでもいくつか踏む手順があるので、順番に見ていく。
1.元のデータセットをの確率行列を作成する
 第1段階として同時確率分布{\displaystyle \bf P_{ij}}、行方向の周辺確率分布を対角要素にした対角行列{\displaystyle \bf P_I}、列方向の周辺確率分布を対角要素にした対角行列{\displaystyle \bf P_J}を用意する。

{\displaystyle 
\begin{eqnarray}
{\bf P_{ij}} &=& \left(p_{ij}\right) = \left(\frac{f_{ij}}{N}\right),\ \ \ N = \sum^m_{i=1}\sum^n_{j=1}f_{ij}
\\
{\bf P_I} &=& diag \left( p_{i+} \right) = \left(\frac{f_{i+}}{N}\right) = \left(\frac{\sum^n_{j=1} f_{ij}}{N}\right)
\\
{\bf P_J} &=& diag \left( p_{+j} \right) = \left(\frac{f_{+j}}{N}\right) = \left(\frac{\sum^m_{i=1} f_{ij}}{N}\right)
\end{eqnarray}
}

2.確率行列を使って行、列それぞれのプロファイルを作成する
 ここで行、列それぞれに対応する行列が作られる。要は行方向に見たときの条件付き確率分布、列方向に見たときの条件付き確率分布のことだが、対応分析においてはそれぞれ行プロファイル列プロファイルと呼ばれる。定式化するとこんな感じ。

{\displaystyle 
\begin{eqnarray}
{\bf N_I} &=& \left( q_{ij} \right) &=& \left( \frac{p_{ij}}{p_{i+}} \right) \ \ \ (行プロファイル)
\\
{\bf N_J} &=& \left( q_{ij}^* \right) &=& \left( \frac{p_{ij}}{p_{+j}} \right) \ \ \ (列プロファイル)
\end{eqnarray}
}

3.ストレッチプロファイルを作成し平均値で中心化する
 上記で作成した行列に対して距離を求めた場合はユークリッド距離となるが、独立性の検定に用いられるカイ二乗統計量に結び付けるためカイ二乗距離になるようデータを変換したものをストレッチプロファイルと呼ぶ(右辺第1項)。こうして作った行列から平均値(右辺第2項)を引いて中心化した式が以下で、これが対応分析に利用する行列になる。

{\displaystyle 
\begin{eqnarray}
{\bf X} &=& \left( x_{ij} \right) &=& \left( \frac{q_{ij}}{\sqrt{p_{+j}}} - \sqrt{p_{+j}} \right)
\\
{\bf X^*} &=& \left( x_{ij}^* \right) &=& \left( \frac{q_{ij}^*}{\sqrt{p_{i+}}} - \sqrt{p_{i+}} \right)
\end{eqnarray}
}

 なお、この式に関する詳細な説明は以下の「対応分析法の基本的な考え方(数理の要点)」が詳しい。
http://wordminer.org/wp-content/uploads/2013/04/2d3ed444f1290ae0e0992a05b64a5817.pdf

 ここまでの作業を書いたコードは以下。

# 確率行列
P <- F / sum(F)

# 行の確率要素(P_iの対角要素)と対角行列化
pi <- matrix(apply(P, 1, sum))
Pi <- diag(pi[, 1])

# 列の確率要素(P_jの対角要素)と対角行列化
pj <- matrix(apply(P, 2, sum))
Pj <- diag(pj[, 1])

# 行プロファイル
Ni <- diag(1 / pi[, 1]) %*% P

# カイ二乗距離を求めるためのストレッチプロファイル
Ni2 <- Ni %*% diag(1 / sqrt(pj[, 1]))
dist(Ni2) # ストレッチしているためdist関数で求める距離はカイ二乗距離になる

# 平均値(列質量の平方根)による中心化
Xi <- sweep(Ni2, 2, sqrt(pj), FUN = "-") # 対応分析における基本のデータ行列①(行)


# 列プロファイル
Nj <- diag(1 / pj[, 1]) %*% t(P)

# カイ二乗距離を求めるためのストレッチプロファイル
Nj2 <- Nj %*% diag(1 / sqrt(pi[, 1]))
dist(Nj2) # ストレッチしているためdist関数で求める距離はカイ二乗距離になる

# 平均値(行質量の平方根)による中心化
Xj <- sweep(Nj2, 2, sqrt(pi), FUN = "-") # 対応分析における基本のデータ行列②(列)
ステップ2:変換したデータについて分析してそれぞれプロットする

 こうして作られた行列{\displaystyle \bf X}{\displaystyle \bf X^*}に対して以下の操作を加えてやるとそれぞれの分散共分散行列になる。

{\displaystyle 
\begin{eqnarray}
{\bf V} &=& {}^t\! {\bf X} {\bf P_I} {\bf X}
\\
{\bf V^*} &=& {}^t\! {\bf X^*} {\bf P_J} {\bf X^*}
\end{eqnarray}
}

 これらの行列についてそれぞれ固有値問題を解けば固有値固有ベクトルが求まるが、これはまさに主成分分析である。つまり対応分析は行方向と列方向のそれぞれに対して主成分分析を行っていることになる。
 ちなみに主成分分析が固有値問題に帰着する理由については毎度のことながら永田先生の本がとても分かりやすい。

多変量解析法入門 (ライブラリ新数学大系)

多変量解析法入門 (ライブラリ新数学大系)

 ブログだと以下が丁寧に書かれてある。
blog.aidemy.net

 こうして求めた固有ベクトルで元のデータ(ストレッチプロファイル)を線形変換してやることで目当ての座標系が得られる。コードはこんな感じ。

### 固有値問題に持ち込む
## 行方向
# 分散共分散行列
Vxi <- t(Xi) %*% Pi %*% Xi

# 固有値問題を解く
result_Xi <- eigen(Vxi)

# 結果の確認
result_Xi$values # 固有値
result_Xi$vectors # 固有ベクトル

# 固有ベクトルによる座標変換
rscore <- Ni2 %*% result_Xi$vectors

## 列方向
# 分散共分散行列
Vxj <- t(Xj) %*% Pj %*% Xj

# 固有値問題を解く
result_Xj <- eigen(Vxj)

# 結果の確認
result_Xj$values # 固有値
result_Xj$vectors # 固有ベクトル

# 固有ベクトルによる座標変換
cscore <- Nj2 %*% result_Xj$vectors

 なお固有値の個数は{\displaystyle min\{m, n\} - 1}で、今回の場合は5個となる。1を引いているのは元々プロファイルが確率行列で和が行和ないし列和が1となる制約条件による。解いた固有値固有ベクトルを確認してみると確かに5成分で表現できていそう。
 ところでこれまでは行と列個別に見てきたが、これら2つには双対性という性質があり、列成分スコアは行成分スコアを用いて、行成分スコアは列成分スコアを用いて求めることができる。算出式を遷移方程式と言って、行成分スコア(rscore)を{\displaystyle \bf Z_{ik}}、列成分スコア(cscore)を{\displaystyle \bf Z_{jk}^*}とすれば(kは固有値/成分に該当)、

{\displaystyle 
\begin{eqnarray}
{\bf Z_{ik}} &=& \frac{1}{\sqrt{λ_k}} \sum^n_{j=1} \left( \frac{p_{ij}}{p_{i+}} \right) z_{jk}^*
\\
{\bf Z_{ik}^*} &=& \frac{1}{\sqrt{λ_k}} \sum^m_{i=1} \left( \frac{p_{ij}}{p_{+j}} \right) z_{ik}^*
\end{eqnarray}
}

で定式化される。コードは以下。とりあえずマッピングに際しては列成分スコアを遷移方程式で求めなおして利用することにする。

# rscoreによる遷移方程式を用いてcscoreを求めなおす
rscore2 <- rscore[, 1:5] # 対象の5成分を抽出
C <- diag(1 / pj[, 1]) %*% t(P) # 列プロファイルの再算出
Dm1 <- diag(1 / sqrt(result_Xi$values[1:5]))
cscore2 <- C %*% rscore2 %*% Dm1 # 遷移方程式

# 成分スコアの確認
dimnames(rscore2) <- list(サンプル = c("サンプル1", "サンプル2", "サンプル3", "サンプル4", "サンプル5",
                                   "サンプル6", "サンプル7", "サンプル8", "サンプル9", "サンプル10"),
                              成分スコア = c("第1成分", "第2成分", "第3成分", "第4成分", "第5成分"))

dimnames(cscore2) <- list(銘柄 = c("銘柄A", "銘柄B", "銘柄C", "銘柄D", "銘柄E", "銘柄F"),
                            成分スコア = c("第1成分", "第2成分", "第3成分", "第4成分", "第5成分"))
rscore2 # 行成分スコア
cscore2 # 列成分スコア

 

ステップ3:重ね合わせる

 ひとまず横軸に第1成分、縦軸に第2成分をとって可視化してみると、こんな感じ。

# 第1成分、第2成分による同時布置図
library(ggplot2)

dat_map <- rbind(rscore2, cscore2)
dat_map <- transform(dat_map, label = c(rep("サンプル", 10), rep("銘柄", 6)))

map1 <- ggplot(dat_map, aes(x =1成分, y =2成分, color = label)) + 
  xlim(-1.5, 1.5) + ylim(-1.5, 1.5) +
  geom_hline(yintercept = 0, color = "black") +
  geom_vline(xintercept = 0, color = "black") +
  geom_point(size = 5) +
  scale_color_manual(
    values = c("サンプル" = "tomato", "銘柄" = "skyblue")
  ) +
  geom_text(aes(label = rownames(dat_map)), vjust = -1, size = 5) +
  theme_bw()

map1

f:id:yanbow221:20180717203905j:plain

これで対応分析の主要アウトプットである同時布置図を描くことができた。

釈然としない

 対応分析の流れを追ってみたが、行成分スコアと列成分スコアが同じ次元にマッピングされる正当性について明瞭に示されている感じはしない。実は特殊なケースを除いてこれら2つのスコアが同次元にはならず、かい離があることが指摘されていて、同じ空間にマッピングされるのは解釈に便利だから、ということらしい。さらに、実際計算してみると分かるが、上で紹介した「対応分析法・数量化Ⅲ類の考え方」に記載されている行成分スコア・列成分スコアとRで計算したスコアは一部符号が反転している成分がある。これは固有ベクトルの持つ性質*4とそれに対する利用ソフトウェアのアルゴリズムが異なるためで、つまり同じデータセットに対して対応分析を行ってもソフトウェアによってマップが視覚的に異なる可能性があることを示している。次元が異なるならうかつに軸に意味付けすることができないし、やはり対応分析の結果の解釈については注意が必要であることを裏付けている。

 長くなってしまったで対応分析の結果をどのように評価していけば良いかについては次のブログでまとめることにします。

*1:お客さんとの雑談の中でPJメンバーの一人がリアルに「欲しがりません勝つまでは」と言うくらい大変だった

*2:対応分析とも言われる

*3:Rのパッケージ「MASS」で提供されているデータセット「caith」を使うためこのようなあるかどうかも分からない設定にした。実際にはcaithはイギリス在住者5,387名に対する調査結果らしい

*4:固有ベクトルvに対してスカラーkを掛けたkvもまた固有ベクトルであるという性質。このkをどう定めるかがソフトウェアによって異なるだろうため、符号の反転が起こる

顧客のセグメント移動を概観するはなし

優良顧客はいつまで優良か

 企業に対する顧客の貢献度を何等か数値化したうえでセグメントを設定し、優良顧客にはこれこれ、離脱顧客にはそれそれといったアプローチをとるのは王道中の王道になっている。差し当たっての問題は「いかにして顧客セグメントを設定するか」だが、ひとまず大枠をとらえて方針を決める程度ならRFMやパレートの法則などのコンセプトに則って設定するというので当面は十分に機能する場合も少なくないと思う。一方でこれだけでは一地点で顧客のステータスを決め打っているに過ぎないため、定期的な更新が必要になるのだが、そこで「そもそも自社の顧客はどの程度セグメントを移ろいやすいのか」という疑問が出てくる。優良顧客はずっと優良顧客でいるのか、あるいはしばらくすると優良でなくなってしまうのかを企業単位、あるいはブランド単位で知っておくことは施策のサイクルや方向性を検討するうえで有意義だろうから、そのために定点観測する指標を持っておきたい。今回はあまり複雑なモデリングをせずに、社会移動データの分析で使われる手法を用いてそれを概観できるのではないかというはなし。

移動指標について

 ここでは移動指標の考え方について『人文・社会科学の統計学』の第9章に則って簡単に説明する。以下のクロス表を用いる。

f:id:yanbow221:20180102154038p:plain

 この表は父から息子への階層移動表で、例えば676名いたホワイトカラーにおいてその息子もホワイトカラーであったのは517名であるという具合の表になっている。移動研究においては息子の階層は父の階層からどの程度影響しているかに焦点を当てていて、それを「開放的」だとか「閉鎖的」だと言う風に言うらしい。ぱっと見たところホワイトカラーの息子はホワイトカラーで居続けていて、農業の息子は各階層に散らばっているように印象を持つ。

①構造移動と純粋移動

 移動性を考える際にはその移動を「構造移動」と「純粋移動」に区別する。構造移動とは対角セルについて、その行と列の周辺分布の違いによって必然的に起きる移動のことを言う。例えば(1, 1)セル(父子ともにホワイトカラー)に注目すると、{\displaystyle n_{1・} = 676}である一方{\displaystyle n_{・1} = 1011}なので、{\displaystyle 1011-676=335}名は必然的にホワイトカラー以外の父階層から入ってこなければならず、これを構造移動と定義する。実際の移動量からこの構造移動量を除いたものは周辺分布の違いによらない移動と考えられるため、これを純粋移動と定義する。先のケースで言えば実際の他階層からホワイトカラーへの移動は{\displaystyle 244 + 250 = 494}なので、この差{\displaystyle 494 - 335 = 159}名がホワイトカラー階層における純粋移動となる*1。定式化すると

{\displaystyle
構造移動\ \ \ \ \ |n_{i・} - n_{・i}| \\
純粋移動\ \ \ \ \ min\{n_{i・},\ n_{・i}\} - n_{ii}
}

となる。全ケース数を{\displaystyle N}とすれば、これらを用いて

{\displaystyle
粗移動率\ \ \ \ \ \ \ \frac{N-\Sigma_{i} n_{ii}}{N} \\
構造移動率\ \ \ \ \ \frac{\Sigma_{i} |n_{i・}-n_{・i}|}{2N} \\
純粋移動率\ \ \ \ \ \frac{\Sigma_{i} (min\{n_{i・}, \ n_{・i}\} - n_{ii})}{N}
}

が定義される。ここで粗移動率=構造移動率+純粋移動率である。

②安田係数

 上記の考え方は階層移動を考えるうえで有用だが、移動率は移動表全体に関するものであり、一方構造移動/純粋移動は各階層における絶対数を扱っているため階層間で単純には比較できない。そのため各階層の移動性を比較する指標が必要となるが、その一つが安田係数である*2。以下で定義される。

{\displaystyle
\begin{eqnarray}
y_{i} &\equiv& \frac{階層iの純粋移動の実現値}{独立移動における階層iの純粋移動の期待値}
\\ 
&=& \frac{min\{n_{i・}, \ n_{・i}\} - n_{ii}}{min\{n_{i・}, \ n_{・i}\} - n_{i・}n_{・i}/N}
\end{eqnarray}
}

 {\displaystyle y }の値は実現値が独立移動と等しければ({\displaystyle n_{ii}=n_{i・}n_{・i}/N })1であり、純粋移動が少なく非移動が多くなると0に近づく。これをもって各階層の移動性を比較、評価することができるようになった。

架空の顧客セグメントデータに適用する

 これらの指標を使ってセグメント間の移動を評価するシミュレーションを行う。想定としては以下の通り。

  1. 特定企業の扱うブランドA, Bについて、同時点でそれぞれの購入顧客を何等かの基準によりhigh(h), medium(m), low(l)の3セグメントに設定した
  2. 各ブランドで対象となった顧客はどちらも1,500名で、最初の時点における各セグメントのサイズはそれぞれ500名ずつとした
  3. 一定期間経った後、各ブランドの対象となった1,500名それぞれを先の基準に照らしセグメントを振り分けた*3

 生成したデータはこんな感じ。どちらも同数の顧客を対象にしているので、一見してブランドAのほうが同セグメントに留まりやすく、ブランドBはセグメント間の移動がより激しくなっていることが分かると思う。
f:id:yanbow221:20180102171522p:plain
f:id:yanbow221:20180102171535p:plain

 これについて紹介した各指標を算出してみる。

#データの生成
set.seed(221)
dat.start <- c(500, 500, 500)
prob.a <- sample(c("h", "m", "l"), rep = T, 500, prob = c(0.7, 0.2, 0.1))
prob.b <- sample(c("h", "m", "l"), rep = T, 500, prob = c(0.2, 0.4, 0.3))
prob.c <- sample(c("h", "m", "l"), rep = T, 500, prob = c(0.05, 0.15, 0.8))
prob.d <- sample(c("h", "m", "l"), rep = T, 500, prob = c(0.35, 0.35, 0.3))
prob.e <- sample(c("h", "m", "l"), rep = T, 500, prob = c(0.3, 0.4, 0.3))
prob.f <- sample(c("h", "m", "l"), rep = T, 500, prob = c(0.33, 0.32, 0.35))

a <- table(prob.a)
b <- table(prob.b)
c <- table(prob.c)
d <- table(prob.d)
e <- table(prob.e)
f <- table(prob.f)

dat.A <- rbind(a, b, c)
row.names(dat.A) <- c("h", "m", "l")

dat.B <- rbind(d, e, f)
row.names(dat.B) <- c("h", "m", "l")

dat.A <- dat.A[, c(1,3,2)]
dat.B <- dat.B[, c(1,3,2)]

#移動率と安田係数を求める関数
ido <- function(data){
  if(nrow(data) != ncol(data)){stop("移動前後のセグメント定義が一致しません")}
  
  ido.kozo <- matrix(NA, nrow = 1, ncol = ncol(data))
  ido.jun <- matrix(NA, nrow = 1, ncol = ncol(data))
  yasuda <- matrix(NA, nrow = 1, ncol = ncol(data))
  for(i in 1:nrow(data)){
    ido.kozo[i] <- max(rowSums(data)[i], colSums(data)[i]) - min(rowSums(data)[i], colSums(data)[i]) #構造移動
    ido.jun[i] <- min(rowSums(data)[i], colSums(data)[i]) - data[i, i] #純粋移動
    yasuda[i] <- ido.jun[i] / (min(rowSums(data)[i], colSums(data)[i]) - (rowSums(data)[i] * colSums(data)[i] / sum(data))) #安田係数
  }
  colnames(ido.kozo) <- colnames(ido.jun) <- colnames(yasuda) <- colnames(data)

  ratio.ara <- (sum(data) - sum(diag(data))) / sum(data) #粗移動率
  ratio.kozo <- sum(ido.kozo) / (2 * sum(data)) #構造移動率
  ratio.jun <- sum(ido.jun) / sum(data) #純粋移動率
  ratio.ido <- rbind(ratio.ara, ratio.kozo, ratio.jun)
  
  kiyo.ara <- ratio.ara / ratio.ara
  kiyo.kozo <- ratio.kozo / ratio.ara #構造移動率 / 粗移動率
  kiyo.jun <- ratio.jun / ratio.ara #純粋移動率 / 粗移動率
  kiyo.ido <- rbind(kiyo.ara, kiyo.kozo, kiyo.jun)

  matome <- cbind(ratio.ido, kiyo.ido)
  dimnames(matome) <- list(c("粗移動率", "構造移動率", "純粋移動率"), c("移動率", "寄与度"))
  
  result <- list(ido.kozo, ido.jun, matome, yasuda)
  names(result) <- c("構造移動", "純粋移動", "クロス表全体の移動率", "安田係数")
  return(result)
}

結果はこんな感じ。

> ido(dat.A) #ブランドAにおける各指標の値
$構造移動
      h   m   l
[1,] 10 112 122

$純粋移動
       h   m  l
[1,] 132 172 94

$クロス表全体の移動率
               移動率    寄与度
粗移動率   0.34666667 1.0000000
構造移動率 0.08133333 0.2346154
純粋移動率 0.26533333 0.7653846

$安田係数
             h         m         l
[1,] 0.4040816 0.6649485 0.3211845

> ido(dat.B) #ブランドBにおける各指標の値
$構造移動
      h  m  l
[1,] 25 44 19

$純粋移動
       h   m   l
[1,] 298 298 289

$クロス表全体の移動率
               移動率     寄与度
粗移動率   0.61933333 1.00000000
構造移動率 0.02933333 0.04736276
純粋移動率 0.59000000 0.95263724

$安田係数
             h         m         l
[1,] 0.9410526 0.9351464 0.9012474

 まず移動率における構造移動率と純粋移動率の比を見てみると、ブランドBはブランドAよりも純粋移動率の比が高く、95%程度になっている。このことからブランドBは顧客の入れ替わりが激しく、少なくともこの2地点間隔ではセグメントに基づく施策の打ち分けが機能しない懸念があるか、あるいはそもそもの商品特性として(ブランド単体においては)企業への貢献度に基づくセグメント設定に意味がないなどの仮説が立つだろう。
 安田係数でも同様で、ブランドBについてはセグメントh, m, lすべて0.9を超えておりセグメントの移動性が独立に近い。一方でブランドAは全体的にブランドBほど係数が大きくなく、最初の時点で設定したセグメントがある程度継続されていることがわかる。その中でセグメントmはやや係数が高めで他セグメントへの移動性が強いことを示しており、これが実データであれば「程々な顧客は上位にも下位にも移ろいやすい*4」ということが言えるだろう。

おわりに

 今回は顧客のセグメントの移動を社会移動データとみなして移動パターンを概観する方法を書いた。先の本には、

社会移動の研究の主要関心の一つは、異なる時点や社会の間で、階層移動のパターンにちがいがあるかどうかという問題である
(P275)

という一節がある。購買データについても同様に

  • 一定間隔で指標を追い、どの程度のスパンで顧客が移り変わっていくか
  • 商品やブランド間で指標を比較し、セグメントの移動パターンがどのように異なっているか

を見ることは施策の更新頻度や顧客セグメントに基づく施策方針(その実施有無を含め)を検討するうえで一つの参考になるんじゃないかと思う。さらに考察を深めるための方法論もあるので、機会があったらまた検討してみたい。

*1:行周辺度数と列周辺度数のどちらが大きいかで行に考えるか列に考えるか変わるのでそこだけ注意

*2:前述の『人文・社会科学の統計学』では結合指数/分離指数という指標の欠点を克服した指標であることや安田係数にも問題があることなど詳しく書いてある

*3:各セグメント同数としたのは最初の時点のみで、その後のセグメント設定では基準に設けた値(例えば累積の購買金額)のセグメント閾値を用いて振り分けを行う

*4:浮動票的な

ビジネス現場における統計的仮説検定のはなし

弁明

 「何をいまさら」感のあるタイトルです。ビジネス現場における仮説検定のはなしは各種ブログで散々書かれているし、機械学習全盛の時代にあまり興味を持たれるはなしでもなさそうな気もする。それでもやはり現場では「有意差」を気にされる場面は相変わらずあるし、またそれなりに正しい知識の下で使わないと害になったり、期待した役割を仮説検定が果たさないような場面もやっぱり相変わらずあるので、改めて実務でありがちなケースを書き留めておきたい。

ビジネス現場で極めて起こりがちな仮説検定にまつわるやりとり

 今回はビジネス現場、なかでもマーケティング活動における仮説検定のはなしを書く。これまで僕が仕事をしてきた中で経験した『危うい仮説検定の扱い』はざっくり3つのケースに分類される。

1.「これって有意差はありますか?」

 顧客セグメント間の施策への反応率でも、異なる施策の効果の差でもなんでも良い。何らか複数対象の数値差を比較するときによく聞かれる質問である。この問い合わせ自体は自然なことだが、有意差のあるなしだけで物事を判断するのは実は必ずしも正しくない。後ほど解説するように、実質的な差がないにも関わらず有意差が出ることは多々ある。有意差だけに頼ってしまうと、実質価値のない仮説を掘り下げてアサッテの方向に施策を展開していくということになりかねない。

2.「そこまでしなくていいよ」

 統計の枠組みを持ち出そうとするとき、「そこまでする必要はない」という判断が下されることがある。もちろん仮説検定をする必要のない状況はある。特にデジタルマーケティングのような環境変化が激しい領域では、実質的に仮説検定に価値がない、あるいはその枠組みに基づく知見の蓄積が困難な状況はあるように思える。ただ上記のセリフが出るとき、往々にしてその発言者はデータから見い出した何らかの差について深く考察を掘り下げていくため、なぜ仮説検定をする必要がないのか不明であることが多い。

3.「『CHITEST』で一発です」

 Excelの関数である。ほかにもいくつか仮説検定を行うための関数はあり、データさえあれば手軽に検定結果を出力することができる。「ツールに突っ込めばOK」という認識はツール開発側からすれば目指すべき最高のユースケースだとは思うが、やはりある程度正しく使わなければミスリードを招く危険性が統計学にはある。

 いずれのケースでも、根本には共通して「謎の権威化」があると思う。仮説検定は万能の方法論ではないし、だからこそ出力結果を鵜呑みにするのは危険である。なにより、仮説検定を「科学的方法論」としてでなく単に「権威」として用いるのは、あまり誠実とは言えない*1統計学の枠組みを用いたデータ分析をするなら、少なくともその分析担当者はそれに関する正しい知識を積み重ねていく姿勢について責任を持つべきだと思う*2
 今回は、そうした有意差だけを判断材料にした分析が期待した成果に貢献しない可能性について、仮説検定の基本的な考え方と実際のシミュレーションを通して検討したいと思う。

 

仮説検定の基礎的な解説

 仮説検定は何らかの重要な仮説の正否を検証することを期待して用いられる。が、実はそれ自体は実際には仮説の正しさを証明することには寄与しない*3。仮説検定のプロセスで定量化された値を人間が見て、実用上「仮説は正しい」や「仮説は間違っている」などと決めているだけである。どういうことか。実務でよく出てきそうなケースを用いて検定の枠組みを説明する。以下のようなケースは、現場では仮説検定を行うモチベーションが高まることが多いだろう。

新しいキャンペーン施策がa、bの2案検討されている。当然より有効な方を打っていきたいので、テストマーケティングを行い効果の良かった方をマーケット全体に対し打つことにした。テストの枠組みはざっくり言うと、

  • 顧客の中から適当な数を抽出し、それをランダムに2グループに分ける
  • 一方には施策aを、もう一方には施策bを打ち、その成果(購買率など、何等かのコンバージョン率としよう)を比較する

なお、テスト顧客の決定や2グループへの分割においてはランダム化が担保されており*4、2グループのテストサイズは同数とする。
ついでに、知見を蓄積したいので、施策aとbの結果の差に何らか違いがあると認められるなら、なぜその違いが生まれたかについても考えたい*5

このとき、仮説検定は以下の手順で行うことになる。

1. 期待している結果を考える。ここでは担当者は施策aと施策bは明らかに違う切り口のアプローチであるため効果は異なる(が実際どちらが良い悪いかはよく分からない)と考えているため、期待している結果は「施策aと施策bのコンバージョン率(CVR)は異なる」となる。

2. 仮説検定の枠組みに則って、2つの仮説を設定する。ここでは手順1で設定した期待している結果を『対立仮説』として設定し、その対立仮説と否定の関係にある仮説を『帰無仮説』として設定する。つまりこう。


{\displaystyle
\begin{eqnarray} 
帰無仮説\ H_0\ :\ \ 施策aのCVR\,&=&\,施策bのCVR \\
対立仮説\ H_1\ :\ \ 施策aのCVR\,&\neq&\,施策bのCVR
\end{eqnarray}
}

 仮説検定は基本的に、「{\displaystyle H_0}が正しいと仮定した下で」検討が行われる。これが極めて重要な部分で、これを押さえていないことが誤用乱用を招く原因になっている、と個人的に思う。
 検定においては、以下に説明する『有意水準』と『p値』の値をもって帰無仮説を棄却するか否かを判定する。いずれもある確率を表しており(これについては後で説明する)、「予め決めておいた有意水準をp値が下回っていた場合、帰無仮説{\displaystyle H_0}を棄却する」、つまり「有意である」と結論することになる。

3. 有意水準を決める
 この確率を下回っていたら帰無仮説を棄却するという基準が『有意水準』である。領域によって異なるが、一般的には5%や1%と設定されることが多い*6

4. p値を算出する
 実際のデータからp値を算出する。検定においてはp値を算出するために『検定統計量』という量が用いられる。少しややこしい話になるが、検定統計量とは、仮説検定の枠組みに当てはめるために、得られたデータ(今回で言えばテスト結果のデータ)を変換した量、と考えれば良い。p値とは「{\displaystyle H_0}が正しいと仮定した下で期待される検定統計量の確率分布において、今回のデータから得られた検定統計量以上に離れた量が得られる確率」のことになる。

5. 帰無仮説の棄却可否を判定する
 ここまでの手順で帰無仮説の棄却可否を判定する材料は揃っているが、文章だと我ながら分かりにくいので、図示するとこんな感じになる。

・帰無仮説を棄却する
f:id:yanbow221:20171022123425p:plain
・帰無仮説を棄却しない
f:id:yanbow221:20171009175649p:plain

 つまり、確率密度関数において、起きることが稀である両端から一定の累積確率(有意水準)を指定すると、それに対応する検定統計量の領域が得られる(『棄却域』という)。今回得たデータから計算した検定統計量が棄却領域内に入っていれば、「{\displaystyle H_0}が正しいと仮定した下では得られる可能性が低い結果が得られた、ということは元の仮説が間違ってと考える方が自然だ」と考えて帰無仮説を棄却し、検定統計量が棄却領域内に入っていなければ、「{\displaystyle H_0}が正しいと仮定した下で得られることが珍しくない結果が得られた、ということは元の仮説を否定する材料は今のところ得られていない」と考えて帰無仮説を棄却しないという判定を行う。これはまさに背理法の考え方になっている。
 一点注意しなければいけないことは、仮説検定においては「帰無仮説を棄却する」か「帰無仮説を棄却しない」かの二択であり、「帰無仮説を採択する」という結論はないという点である。これは以下の仮説について考えるとわかりやすいだろう。


{\displaystyle
帰無仮説\ H_0\ :\ \ 白鳥は白い \\
対立仮説\ H_1\ :\ \ 白鳥は白くない
}

 もし例えば黒い白鳥を見つけた場合、帰無仮説は明確に否定(棄却)される一方で、白くない白鳥を見つけられない(白い白鳥しか見つからなかった)としても、それは「白鳥は白くない」ことを示す根拠がないだけで、積極的に「白鳥は白い」と結論づけることは難しい*7。こうした非対称性が仮説検定には存在する*8。ちなみに『統計学を拓いた異才たち』では、フィッシャーの言った言葉としてこんな一説が紹介されている(文庫版P168)。

ある仮説が扱うことのできる事実に矛盾していないという理由だけで、その仮説の正しさが証明されたと信じるような論理的誤謬は、統計的な論拠においても、その他の科学的な論拠においても、受け入れられない。……それゆえ、有意性検定で、データが仮説に矛盾する場合には仮説を棄却したり無効にしたりできるが、矛盾しない場合には仮説が正しいとはいえない、と一般に広く理解されたならば、有意性検定の明快さは大いに増すだろう。

統計学を拓いた異才たち - デイヴィッド・サルツブルグ,竹内惠行,熊谷悦生|日本経済新聞出版社

 こうした手続きを経て、今回で言えば2つの施策aとbは効果に差があるとみなせるのか、効果はどうも差がなさそうだ(それを否定する材料は今回は得られなかった)と結論を下すことになる。

シミュレーション

 ここまで仮説検定の基本的な考え方を紹介した。これを前提に、得られたデータに盲目的に仮説検定を適用するだけは危険、というのがどういうことか、シミュレーションで検討したい。今回は全く同じ結果が得られてもその判断がサンプルサイズによって異なるということをシミュレーションで確認し、p値だけを見てその他の物事(今回で言えばサンプルサイズ)について注意を払わないと、一貫した結論が得られないしその後どうすれば良いのかの指針も立たない、ということを示す。
 想定するケースは前述のテストマーケティングで、2つの施策の効果に差があるかどうかを仮説検定によって検討するが、その際得られるp値がテストに用いたサンプルサイズによってどのように影響を受けるかを確認する。コードはこんな感じ。

N <- 100000 #母集団サイズ。使わないけど念のため
n <- seq(100, 1000, 10) #標本サイズ
cvr_a <- seq(0.01, 0.3, 0.005) #施策aのコンバージョンレート
cvr_b <- seq(0.01, 0.3, 0.005) #施策bのコンバージョンレート
cvr <- as.data.frame(expand.grid(cvr_a, cvr_b)) #コンバージョンレートの組み合わせ
colnames(cvr) <- c("cvr_a", "cvr_b")
size <- as.data.frame(matrix(0, nrow = nrow(cvr), ncol = length(n)))
for(i in 1:ncol(size)){
  colnames(size)[i] <- paste("n=", n[i])
}

#2群の母比率の差の検定(p < .05)
#2群の標本サイズは常に等しく、完全無作為抽出が担保されていると仮定
for(i in 1:length(n)){
  cv <- round(cvr * n[i])
  for(j in 1:nrow(cvr)){
    size[j, i] <- prop.test(c(cv[j, 1], cv[j, 2]), c(n[i], n[i]))$p.value
  }
} #sizeオブジェクトに施策a, bのコンバージョンレート、標本サイズの組み合わせごとのp-valueを格納する

dat <- cbind(cvr, size)

size2 <- size
for(i in 1:ncol(size2)){
  size2[, i] <- ifelse(size[, i] >= 0.05, 0, 1)
}
dat2 <- cbind(cvr, size2)

#ヒートマップアニメーション
library(ggplot2)
library(animation)

ani.options(convert = "C:/ImageMagick-7.0.7-5-portable-Q16-x64/convert.exe")
animation::saveGIF({
  for(i in 1:ncol(size2)){
    print(ggplot(dat2, aes(as.factor(cvr_a), as.factor(cvr_b))) +
            geom_tile(aes(fill = dat2[, i + 2]), color = "grey") +
            scale_fill_gradient2(low="white",high="steelblue") +
            labs(x = "CVR of A", y = "CVR of B", title = colnames(size2[i])) +
            theme(legend.position = "none") +
            geom_abline(intercept = 0, slope = 1, col = "pink2"))
  }
}, interval = 0.1, movie.name = "p-value.gif", ani.width = 1500, ani.height = 1000)

f:id:yanbow221:20171021233348g:plain
 添付のアニメーションは横軸に施策aのCVR、縦軸に施策bのCVRをとっていて、青色のセルは「有意水準5%において帰無仮説が棄却される(有意と判定される)」セルであり、その範囲がサンプルサイズの変化({\displaystyle n=100}から{\displaystyle n=1000}までを10刻みで変化)によって大きくなることを表している*9。つまり、施策a, bのCVRが全く同じ組み合わせであっても、それに用いたサンプルサイズによって帰無仮説を棄却できるか棄却できないかの判断が変わってしまうということになる。
 なぜこのようなことが起きるかというと、仮説を検討するために必要な検定統計量は基本的にその数式の分母に{\displaystyle n }の逆数が入っているからである。今回で言うと、

{\displaystyle
\begin{eqnarray}
T &=& \frac{CVR_a - CVR_b}{\sqrt{(\frac{1}{n_a}-\frac{1}{n_b})CVR(1-CVR)}}
\\
ただし、
\\
CVR &=& \frac{n_a CVR_a + n_b CVR_b}{n_a + n_b}
\end{eqnarray}
}

が検定統計量の算出式で*10、サンプルサイズが大きくなるほど分母が小さくなり検定統計量の値が大きくなる、つまり帰無仮説を棄却しやすくなることがわかる。
 これから、とりあえずテストしたり手元にあるデータに対して仮説検定を適用してp値だけで判断を下すということが、危なっかしいということがわかると思う。今回の記事では扱わないが、仮説検定においてはp値の他に検出力、効果量、そして今回扱ったサンプルサイズなど、いくつかの量を検討したうえでのテスト計画が必要になってくる*11

おわりに

 今回はビジネス現場における仮説検定について、p値乱用の危険性について見てきた。とはいえ仮説検定にまつわる議論と乱用はこの枠組みが生まれた当初からずっとあるみたいで、掘り下げれば掘り下げるほど何が何だかよくわからなくなってくる*12。うまくバランスを取りながら、「ここを押さえてないとマズイ」ポイントはしっかり押さえたうえで、現場に適用していきたいなぁ。

統計学を拓いた異才たち(日経ビジネス人文庫)

統計学を拓いた異才たち(日経ビジネス人文庫)

統計学入門 (基礎統計学?)

統計学入門 (基礎統計学?)

*1:とはいえどこまでも突き詰めると際限はないし、結局は人それぞれの線引きラインの違いでしかないのかもしれない。僕もこの後疑わずにRのパッケージを使うし…

*2:本当は「正しい知識を用いる」とか書きたかったけど、自分自身が正しい知識を持っていると確信できる日は少なくとも僕には永遠に来ないと思うので、「姿勢」に対して責任を持つというニュアンスにしてます

*3:さらに言うと、先に書いた「何らか複数対象の数値差」の差の大きさや重要性についても仮説検定自体は示唆しない。

*4:要は変な偏りなくテストできてるよという理解で概ね問題ない

*5:単に効果が高そうな方を選ぶだけなら、テスト結果が良かった方を選べば基本的には良いだろうからで、問題は結果の差に意味があると認めてそれを考察したいからである

*6:この水準は基本的には「慣例」であること以上の根拠はない

*7:白鳥を例にしたのはナシム・ニコラス・タレブ著『ブラック・スワン』より。あの類の本ってホント面白いですよね

*8:とはいえ、実際上は帰無仮説が正しい、というかそれを受け入れてその後に展開するというのはやある程度むを得ないとも思う

*9:2つの施策のCVRをどちらもCVR=0.01~0.3としているので、この図はy=xに対して対称であることに留意。後目盛かぶりの都合で横軸と縦軸の幅を同じにしてない点はご容赦ください…

*10:シミュレーションで使ったprop.testはカイ二乗近似で出してるっぽい

*11:後は信頼区間も重要な枠組みで、これもそのうち記事にして検討したい

*12:僕はそうでした、てかそうなってます…

ガンマ分布のはなし

いまいちイメージがつかみにくい*1

 「ガンマ分布!」と言われてもイマイチなんのこっちゃわかりにくい気がする。二項分布やポアソン分布のように直感的なイメージを持ちにくく、教科書でも「ガンマ分布は指数分布の一般化です」などとどことなく味気ないからだと思う。なので今回は、ガンマ分布の成り立ちとこれに従う現実の現象をみることで、なんとなくイメージを持ちやすくしてもらうことを目指して記事を書いてみる。

ガンマ分布

 先に書いた通りガンマ分布は指数分布を一般化したもので、以下の確率密度関数で定義される。


{\displaystyle
f(x) = \frac{λ^α}{Γ(α)}  x^{α-1} e^{-λx} ~~ (x\geq0), ~~~~~ 0~~(x<0)
}

このガンマ分布に対して与えられる意味合いは、「一定期間に1回起きると期待されるランダムな事象が複数回起きるまでの時間の分布」というものだが、なぜこの式がその意味を持つのかを理解するには、幾何分布や指数分布との関連を見ていくのが良い。図示するとこんな感じ。赤枠の確率分布について、以下順に追ってみる。
f:id:yanbow221:20170902095612p:plain

幾何分布について

 幾何分布は成功(確率p)か失敗(確率q)のいずれかをとる試行(ベルヌーイ試行という)において、初めて成功するまでの試行回数の分布と定義される。その確率は、


{\displaystyle
f(x) = pq^{x-1}, ~~x=1,2,3,...
}

初めて成功するまでの試行回数がx回なのだから(x-1)回の失敗と最後1回の成功ということでこれはわかりやすいと思う。この幾何分布は『離散的な待ち時間分布』とも呼ばれる。1回の成功を得る確率がpなのだから、別の見方をすれば1回の成功を得るためには(平均的に)1/pだけの時間を要するということからこの名前がついている。

指数分布について

 先ほどの幾何分布は試行回数と説明したように離散的な分布になっている。これを連続値について扱ったものが指数分布である。確率密度関数は、


{\displaystyle
f(x) = λe^{-λx} ~~(x\geq0), ~~~~~ 0~~(x<0)
}

で定義される。幾何分布が離散的な待ち時間分布と呼ばれたのに対し、この指数分布は『連続的な待ち時間分布』と呼ばれる。幾何分布と指数分布は一見したときの式の形が異なっているため、どちらも待ち時間分布という同じ性質を持っていることを理解するには多少の証明が必要になる。以下がその解説。

<指数分布が連続的な待ち時間分布と呼ばれることのワケ>

実数xに対して、y=[x]をxを超えない最大の整数と定義する。例えばx=3.2ならy=3、x=5.84ならy=5、といった感じ。指数分布に従う確率変数Xに対して確率変数Yはどのような分布に従うかを考える。nを正の整数として、


{\displaystyle
\begin{eqnarray} 
P(Y = n) &=&  P([X ] = n)
\\
&=&
P(n \leq X < n+1)
\end{eqnarray}
}

これより、


{\displaystyle
\begin{eqnarray} 
P(Y \geq n) &=& P(Y = n) + P(Y = n+1) + ~~...
\\
&=&
P(n \leq X < n+1) + P(n+1 \leq X < n+2) + ~~...
\\
&=&
P(X \geq n)
\\
&=&
\int_n^∞ λe^{-λx} dx
\\
&=&
e^{-λn}
\end{eqnarray}
}

よって、

{\displaystyle
\begin{eqnarray}
P(Y=n) &=& P(Y \geq n) - P(Y \geq n+1)
\\
&=&
e^{-λn} - e^{-λ(n+1)}
\\
&=&
e^{-λn} (1-e^{-λ})
\\
&=&
(e^{-λ})^n (1-e^{-λ})
\end{eqnarray}
}


{\displaystyle
q = e^{-λ}
}
とすれば、幾何分布の式と同じ形をしていることがわかる。この指数分布は、「期間1/λあたりに1回起きると期待される事象が1回起きるまでの時間の分布」という意味合いを持つことになる*2

ガンマ分布について

 ガンマ分布は指数分布の一般化である。改めてガンマ分布の確率密度関数は以下。


{\displaystyle
f(x) = \frac{λ^α}{Γ(α)}  x^{α-1} e^{-λx} ~~ (x\geq0), ~~~~~ 0~~(x<0)
}

指数分布と見比べてみるとよく似た形をしていることがわかる。実際α=1としてやると指数分布に一致する。
※補足すると、分母のΓ(α)は階乗n!を一般化したもので、


{\displaystyle
Γ(1)=1
}

となる。このガンマ関数が式内に入っている理由は、確率密度関数


{\displaystyle
f(x) = \int_0^∞ \left \{ \frac{λ^α}{Γ(α)}  x^{α-1} e^{-λx}\right \} dx = 1
}

と1にするため(確率なので定義域全体で積分すると1になる必要がある)に必要な係数だと理解しておいていい。なお、詳しい導出方法は以下が参考になる。

ガンマ関数の導出

 さて、指数分布を一般化したものであるガンマ分布が持つ意味合いとしては、「期間1/λあたりに1回起こると期待されるランダムな事象がα回起こるまでの時間の分布」となる。ガンマ分布に従う現象としては、製品部品の寿命やウイルスの潜伏期間のほか、人の体重なんかも従うらしい*3

 ここまで、ガンマ分布とその他の確率分布の関係性やその意味合いについて何となく見てきた。

ガンマ分布の例

 ここからは実際にデータを使ってガンマ分布がどんな形かを確認してみる。ただ少し探したところ適当なデータセットが見つからなかったので、今回はガンマ分布の意味合いから考えて闇のゲーム*4を題材にしてこれを考えることにする。

<闇のゲーム「モンスター・ワールド」とは>

 今回取り扱うのは、僕世代の男性なら知らない人はいないであろう『遊戯王』で、千年パズルを持つ主人公遊戯と千年リングを持つライバル獏良が戦ったボードゲームである。詳しくは単行本を読んでいただくとして、ここで知っておくべきルールは、

・獏良は1体のボス『ゾーク』を操り、遊戯は5体のキャラクターを操る
・遊戯は盤上でキャラクターを進めていってゾークを倒せば勝ち。逆にキャラクター全員が倒れると負け
・キャラクターとゾークは互いに10面ダイス2個を振って出た目の数で勝負をする。赤いダイスは十の位、白いダイスは一の位の数字を表す(つまり目は00~99の値をとる)。
・ダイスの目の数が「小さい」ほうが良い目で、良い目を出したほうが相手へ攻撃することができる(同じ目の時はやり直し)
・相手より良い目を出しかつその目が「00~09」だった場合、クリティカルヒットになる
・1度クリティカルヒットを出したキャラクターは勝ち抜けとなり、次のキャラクターとゾークのバトルに進む
・遊戯はクリティカルヒットを全キャラクターで出せば(つまり累計5回のクリティカルヒット)、ゾークを倒すことができ勝ちとなる
・獏良は『洗脳ダイス』により常に「10」の目を出してくる
・獏良のほうが良い目であっても、遊戯側のキャラクター『杏子(アンズ)』が即座に回復魔法をかけるため実質遊戯のキャラクターを倒すことはない
・単位時間(10秒)あたりのダイスロールは平均的に0.7回

である。図化するとこんな感じ。なお、ガンマ分布のシミュレーションのため原作のルールを多少(というか大幅に)変更している。
f:id:yanbow221:20170905235554p:plain

 このルールのもと、遊戯は獏良に勝つまでに大体どの程度の時間がかかるかを見てみる。まずキャラクターがゾークにクリティカルヒットを出すには「00, 01, 02, 03, ..., 09」のいずれかの目を出さなければならないので、その確率は1/10である。単位時間10秒あたりダイスロールは0.7回のペースで行われるので、単位時間当たりのクリティカルヒットの可能性は0.7×(1/10)=0.07、すなわちλ=0.07である。これは言い換えればキャラクター1人がクリティカルヒットを出すまでにかかる平均的な時間は(1/0.07)×10≒143秒ということである。これが各キャラクターにおいて出れば遊戯の勝ちだから、その時間の分布はλ=0.07、α=5のガンマ分布に従うと予想される。その確率密度関数は以下のようにして描かれる。

#闇のゲームのシミュレーション
set.seed(221)

#λ=0.035、α=5のガンマ分布
x <- seq(0, 600, 1)
λ <- 0.07
α <- 5

curve(dgamma(x, shape = α, scale = 1/λ),
      xlim = c(0, 600), ylim = c(0, 0.02), xlab = "時間(×10秒)", ylab = "確率密度")
par(new = T)
abline(v = 71.4, col = "red") #平均値

f:id:yanbow221:20170905232317j:plain

これに対して10,000回シミュレーションを回してみてその分布を比較してみるとこんな感じ。

#シミュレーション
hako <- matrix(nrow = 10000, ncol = 1) #結果の格納箱
dice <- seq(0, 99, 1)

for(i in 1:10000){
  count <- 0 #クリティカルヒット数カウント用
  total <- 0 #全ダイスロール回数カウント用
  while(count <= 5){
    x <- sample(dice, 1)
    if(x < 10){
      count <- count + 1
      total <- total + 1
    }else{
      total <- total + 1
    }
  }
  hako[i, 1] <- total
}

hako2 <- (hako - 5) / 0.7 #実質離散分布で、5単位未満で遊戯が勝つ可能性はないため調整

dens <- density(hako2) #確率密度の推定
par(new = T)
lines(dens, col = "blue", xlim = c(0, 600)) #シミュレーションに基づく分布

f:id:yanbow221:20170905233156j:plain

今回はシミュレーションの方法的に実質離散型の分布になっているため(だと思う)多少ずれているが、概ねフィットしているといえそう。遊戯は獏良から勝利を得るのに大体12分ほどかかっている。ちなみに得られたシミュレーションデータについてガンマ分布を仮定したパラメータ推定を行うとα'=5.6、λ'=0.07となる。

> est <- fitdistr(hako2, "gamma")
> est
      shape         rate    
  5.615397633   0.071035304 
 (0.077115818) (0.001020378)

終わりに

 今回は式が複雑でややイメージを掴みにくいガンマ分布についてみてきた。このガンマ分布はベイズ統計でも活躍する分布なので、こちらの勉強がある程度形になったらまた取り上げたい。
 ところで、冒頭で分布の関連を図示した中で指数分布とポアソン分布の間に双方向の矢印を引いた。実はこの2つの分布は同じ事象を別の観点から見たものになる。これについては以下の資料がとても分かりやすい。

www.slideshare.net

 ガンマ分布について見たので次はベータ分布についてなんか書いてみよう。

*1:個人の感想です。

*2:確率変数Xの期待値E(X)を求めることで期間1/λというのが導かれる

*3:体重については待ち時間分布としての性質というより、(1)一般に身長が正規分布に従うこと、(2)ガンマ分布の三乗根は正規分布になること、(3)人の縦横高さの積で体重を近似できそうなこと、から理解できるようなのだけど、この(2)について詳しく確認できていない。詳しい方いらっしゃれば、教えていただけると嬉しいです。m(_ _)m

*4:小学生の頃はこれとMTGというカードゲームになけなしのお小遣いをつぎ込んだ

多重共線性の仕組みとシミュレーション

多重共線性の確認

 マーケティングリサーチで重回帰分析を用いる際、各説明変数の係数に関心が行くことは多いと思う。「これこれの説明変数は他に比べて影響度が強いですね、コイツを上げる(下げる)施策を考えましょう*1」などと言って施策の方向性を定めていくという具合である。
 そこで必ず問題になるのがこの多重共線性で、一般に説明変数間に強い相関がある際に偏回帰係数の推定が不安定になるという状況を指す。この数理については『統計的方法のしくみ』の解説が抜群にわかりやすい。そこで多重共線性の仕組みについて、シミュレーションで確かめつつ記事に残しておくことにする。

統計的方法のしくみ―正しく理解するための30の急所

統計的方法のしくみ―正しく理解するための30の急所

推定の流れ

 この記事では2つの説明変数で重回帰分析を行うことを考える。回帰モデルは


{\displaystyle y_i = β_0 + β_1x_{1i} + β_2x_{2i} + ε_i, ε_i ~ N(0, σ^2)
}

で表現される。最小二乗法では残差平方和を最小にするようパラメータを求めるので、残差を


{\displaystyle e_i = y_i - \hat{y_i}
}

としたとき、残差平方和


{\displaystyle S_e = \sum_{i=1}^n{e_i}^{2} = \sum_{i=1}^n{\{y_i - (\hat{β_0} + \hat{β_1}x_{1i} + \hat{β_2}x_{2i})}\}^2
}

を最小にする{\displaystyle \hat{β_0},\, \hat{β_1},\, \hat{β_2}}を求める。各パラメータで偏微分して=0と置いて求めれば良いのでそのまま云々とやっていくと、行列を用いて以下のように表現できる。


{\displaystyle \begin{pmatrix}
S_{11} & S_{12} \\
S_{12} & S_{22} \\
\end{pmatrix}
\begin{pmatrix}
\hat{β_1} \\
\hat{β_2} \\
\end{pmatrix}
=
\begin{pmatrix}
S_{1y} \\
S_{2y} \\
\end{pmatrix}
}

ここで {\displaystyle S_{11}, S_{22}}はそれぞれ {\displaystyle x_1, x_2}の偏差平方和、 {\displaystyle S_{12}} {\displaystyle x_1} {\displaystyle x_2}の偏差積和、 {\displaystyle S_{1y}, S_{2y}}はそれぞれ {\displaystyle x_1} {\displaystyle y} {\displaystyle x_2} {\displaystyle y}の偏差積和。あとは逆行列を左から掛けてやれば良いので、


{\displaystyle
\begin{eqnarray}
\begin{pmatrix}
\hat{β_1} \\
\hat{β_2} \\
\end{pmatrix}
&=&
\begin{pmatrix}
S_{11} & S_{12} \\
S_{12} & S_{22} \\
\end{pmatrix}
^{-1}
\begin{pmatrix}
S_{1y} \\
S_{2y} \\
\end{pmatrix}
\\
&=& \dfrac{1}{S_{11}S_{22} - S_{12}^{2}}
\begin{pmatrix}
S_{22} & {-S}_{12} \\
{-S}_{12} & S_{11} \\
\end{pmatrix}
\begin{pmatrix}
S_{1y} \\
S_{2y} \\
\end{pmatrix}
\\
&=& \dfrac{1}{S_{11}S_{22} - S_{12}^{2}}
\begin{pmatrix}
S_{22}S_{1y}-S_{12}S_{2y} \\
{-S_{12}S_{1y}+S_{11}S_{2y}} \\
\end{pmatrix}
\end{eqnarray}
}

ふぅ*2。これで偏回帰係数が求められた。

多重共線性のメカニズムと相関との関係

 先の式を見ると行列に \displaystyle{1/(S_{11}S_{22}-S_{12}^{2})}が掛かっており、これが0のとき(すなわち偏差平方和・積和行列に逆行列が存在しないとき)は解が無数に存在する(不定)か解が存在しない(不能)ということになる。この状況を多重共線性が存在するという。
 では、この \displaystyle{1/(S_{11}S_{22}-S_{12}^{2})}=0という状況と説明変数間の相関の関係がどうなっているかというと、

\displaystyle{
\begin{eqnarray}
S_{11}S_{22}-S_{12}^{2}=0 &\Longleftrightarrow& \dfrac{S_{12}^{2}}{S_{11}S_{22}}=1
\\
&\Longleftrightarrow& r_{x_{1}x_{2}}^{2} = \left\{\dfrac{S_{12}}{\sqrt{S_{11}S_{22}}}\right\}^{2} = 1
\\
&\Longleftrightarrow& r_{x_{1}x_{2}} = ±1
\end{eqnarray}
}

これにより、説明変数間の相関が±1という状況で多重共線性が発生するということがわかった*3。また \displaystyle{ r_{x_{1}x_{2}}}が±1に極めて近ければ\displaystyle{ S_{11}S_{22}-S_{12}}が0に極めて近いということになるので、解の変化幅が大きくなり安定しないということになる*4

Rでシミュレーション

 ここまでで多重共線性の仕組みが分かったので、実際回帰式の推定にどのように影響するのかをシミュレーションで確かめてみる。シミュレーション内容はこんな感じ。

  • N=1,000で目的変数 \displaystyle{y}と1つ目の説明変数 \displaystyle{x_{1}}の正規乱数を作る(今回は \displaystyle{r_{yx_{1}}=0.4}で固定しておく*5
  • 2つ目の説明変数 \displaystyle{x_{2}}については \displaystyle{r_{x_{1}x_{2}}}を0, 0.05, 0.1, ..., 1と0.05ずつ刻み、各々の値に対し100個の \displaystyle{x_{2}}セットを作る
  •  \displaystyle{r_{x_{1}x_{2}}}毎にパラメータがどのようにぶれるのか確認
#シミュレーション用データの作成
set.seed(221)

#目的変数と説明変数1の作成
y <- rnorm(1000, mean = 0, sd = 1) #目的変数
e <- rnorm(1000, mean = 0, sd = 1) #説明変数作成のために必要なデータ
p=0.4 #yとx1の相関係数を固定する
x1 <- p * y + sqrt(1 - p ^ 2) * e #yとの相関が0.4の説明変数x1を作成

#説明変数2を作成し、シミュレーションを実施
r = seq(0, 1, by = 0.05) #r_12 = 0, 0.05, ..., 1で設定(21段階)

b0 <- matrix(nrow = 100, ncol = 21) #β0を格納する箱。各相関係数に対し100個ずつ
b1 <- matrix(nrow = 100, ncol = 21) #β1を格納する箱。各相関係数に対し100個ずつ
b2 <- matrix(nrow = 100, ncol = 21) #β2を格納する箱。各相関係数に対し100個ずつ

for(i in 1:21){
  stock <- matrix(nrow = 1000, ncol = 100)
  for(j in 1:100){
    e2 = rnorm(1000, mean = 0, sd = 1)
    x2 <- r[i] * x1 + sqrt(1 - r[i] ^ 2) * e2
    stock[, j] <- x2
    
    dat <- as.data.frame(cbind(y, x1, x2 = stock[, j]))
    res <- lm(y ~ ., data = dat)
    b0[j, i] <- res$coefficients[1]
    b1[j, i] <- res$coefficients[2]
    b2[j, i] <- res$coefficients[3]
  }
}

これで \displaystyle{\hatβ_{0}, \hatβ_{1}, \hatβ_{2}}についてそれぞれ、100(ある \displaystyle{r_{x_{1}x_{2}}}下でのデータセット数)×21( \displaystyle{r_{x_{1}x_{2}}}の段階数)の値が得られた。例えば \displaystyle{\hatβ_{1}}の箱にはこんな感じで値が入っている。

> head(b1)
        r=0    r=0.05     r=0.1    r=0.15     r=0.2    r=0.25     r=0.3    r=0.35     r=0.4    r=0.45     r=0.5
1 0.4306679 0.4325078 0.4283519 0.4317270 0.4361706 0.4260468 0.4184288 0.4508224 0.4204306 0.4261973 0.4178097
2 0.4322390 0.4346633 0.4311487 0.4308968 0.4269444 0.4343129 0.4436880 0.4380596 0.4259308 0.4228499 0.4535279
3 0.4313652 0.4297881 0.4267221 0.4289261 0.4415001 0.4260835 0.4101387 0.4224972 0.4297717 0.4291471 0.4383937
4 0.4304963 0.4330273 0.4283483 0.4204201 0.4253410 0.4492216 0.4335296 0.4235167 0.4709619 0.4294702 0.4160259
5 0.4322628 0.4335394 0.4328325 0.4331517 0.4354849 0.4317349 0.4314967 0.4338696 0.4191118 0.4356418 0.3976430
6 0.4334116 0.4313089 0.4324587 0.4292935 0.4345609 0.4290567 0.4176773 0.4438420 0.4313572 0.4149576 0.4371784
     r=0.55     r=0.6    r=0.65     r=0.7    r=0.75     r=0.8    r=0.85     r=0.9    r=0.95       r=1
1 0.4111970 0.4093910 0.4339163 0.4184000 0.3734337 0.4756565 0.3922978 0.4544657 0.4121455 0.4322047
2 0.4391379 0.4516133 0.3873439 0.4027791 0.4100673 0.4379951 0.3401908 0.3585950 0.3741827 0.4322047
3 0.4076027 0.4311014 0.4450212 0.4278456 0.4481529 0.4148706 0.4109423 0.5119509 0.4214167 0.4322047
4 0.4324386 0.4308051 0.4209028 0.4485870 0.4614509 0.3932755 0.4713389 0.5392622 0.4512087 0.4322047
5 0.4360566 0.4087040 0.4727114 0.4316713 0.4444006 0.4488128 0.4205936 0.4277047 0.4467345 0.4322047
6 0.4354719 0.4178111 0.4086015 0.4483683 0.4363689 0.3716729 0.3623556 0.4071001 0.6429823 0.4322047

> nrow(b1)
[1] 100


どのような結果になっているか箱ひげ図で確認してみると、 \displaystyle{r_{x_{1}x_{2}}}が1に近づくほどブレ幅が大きくなることがわかる。

#b1の偏回帰係数の図示
b1 <- as.data.frame(b1)
colnames(b1) <- colnames_r

b1_long <- b1 %>% 
  tidyr::gather(key = r, value = b1)

ggplot(b1_long, aes(x = r, y = b1)) +
  geom_boxplot() +
  ggtitle("説明変数間の相関の強弱に伴うb1のバラツキ")

#b2の偏回帰係数の図示
b2 <- as.data.frame(b2)
colnames(b2) <- colnames_r

b2_long <- b2 %>% 
  tidyr::gather(key = r, value = b2)

ggplot(b2_long, aes(x = r, y = b2)) +
  geom_boxplot() +
  ggtitle("説明変数間の相関の強弱に伴うb2のバラツキ")

 \displaystyle{\hatβ_{1}}の結果図示
f:id:yanbow221:20170805033034j:plain
 \displaystyle{\hatβ_{2}}の結果図示
f:id:yanbow221:20170805033052j:plain

ちなみに \displaystyle{\hatβ_{0}}について同じように図示するとこうなる。
 \displaystyle{\hatβ_{0}}の結果図示
f:id:yanbow221:20170805033808j:plain

これは、 \displaystyle{\hatβ_{0}}については

\displaystyle{
\hat{β_0}=\bar{y}-\hat{β_1}\bar{x_1}-\hat{β_2}\bar{x_2}
}
を解くことで求まるが、\displaystyle{\bar{x_1},\bar{x_2}}についてはどのデータセットでも同じ値であること、そのうえで\displaystyle{(\hat{β_1},\hat{β_2})}の組み合わせ(足し合わせ)が

> head(b1 + b2)
        r=0    r=0.05     r=0.1    r=0.15     r=0.2    r=0.25     r=0.3
1 0.4062166 0.4154052 0.3839790 0.4027693 0.4471804 0.4732325 0.4681899
2 0.4664639 0.4493094 0.4712963 0.4394332 0.4249584 0.4497007 0.4273785
3 0.4267857 0.4824163 0.4557214 0.4181341 0.4777092 0.4289483 0.4333954
4 0.4379343 0.4109607 0.4335442 0.4378974 0.4731251 0.4555382 0.4145050
5 0.4744823 0.4237386 0.4455729 0.4360844 0.4437070 0.4832751 0.4168798
6 0.4561966 0.4234361 0.4466039 0.3999932 0.4673920 0.4392135 0.4219600
     r=0.35     r=0.4    r=0.45     r=0.5    r=0.55     r=0.6    r=0.65
1 0.4214690 0.3901941 0.4546184 0.4290851 0.4504682 0.4186020 0.4276473
2 0.4738807 0.4074830 0.4196819 0.4178540 0.4350354 0.4318700 0.4474259
3 0.4628155 0.4458085 0.4428629 0.4416364 0.4580332 0.4516505 0.4254354
4 0.4202259 0.4383358 0.4085167 0.4249698 0.4341524 0.4211597 0.4240621
5 0.4645901 0.4214547 0.4467117 0.4362892 0.4431544 0.4008554 0.4301498
6 0.4132844 0.4117565 0.4340837 0.4367685 0.4224682 0.4045233 0.4415290
      r=0.7    r=0.75     r=0.8    r=0.85     r=0.9    r=0.95 r=1
1 0.4308617 0.4302392 0.4172367 0.4228712 0.4474329 0.4418713  NA
2 0.4246184 0.4306662 0.4290404 0.4310337 0.4285018 0.4379655  NA
3 0.4430312 0.4253368 0.4116594 0.4210978 0.4446660 0.4302211  NA
4 0.4359602 0.4359433 0.4402329 0.4316962 0.4229854 0.4388992  NA
5 0.4174786 0.4173395 0.4187569 0.4235048 0.4437448 0.4364285  NA
6 0.4327213 0.4251063 0.4441639 0.4309610 0.4294942 0.4313048  NA

とrによらずほぼ一定であることから理解できる。

終わりに

 ここまで多重共線性の仕組みとその影響を確認してきた。今回は値自体に意味はないため結果の実際的な解釈は行っていないが、多重共線性は係数の大きさや符号を不安定なものにする危険があることがよくわかった。オーソドックスな回帰でも多重共線性をはじめいくつか落とし穴や確認すべき点があるので、結果に基づいた意思決定を行う際は十分注意して分析していこうと思う次第です。

*1:説明変数と目的変数は因果性を示すものではないが、実務的にはそこに因果性を認めて設計することが多いと思う

*2:あまり慣れてないのでLaTeX書くの疲れる

*3:相関係数が±1ということは\displaystyle{x_{1}とx_{2}}が一直線上に並んでいるということ。共通の直線上にあるから「共線」らしい

*4:『人文・社会科学の統計学』の例6.11がわかりやすい

*5:設定した相関係数を担保するデータ作成のための処理がコード内に含まれているが、なぜその式でいけるかは相関のある2つの擬似乱数の生成(Pythonサンプルつき) - Qiitaがとても分かりやすい