一日一膳(当社比)

RとJavaと時々数学

Rで始める傾向スコア分析による効果検証

最近因果推論の勉強を始めた管理人です. 因果効果の推定の基礎となるポテンシャルアウトカムフレームワーク, 選択バイアス, 傾向スコア(特に逆確率重み付き推定)についてRでの実装つきでまとめたので記事にしました. 大まかな内容については、以下の本を参考にしています.

参考 - 効果検証入門(安井)

1. 因果効果

対象への働きかけによって, 結果にどの程度影響が及ぶかについて推定したい状況を考えます.
例えば, ある新薬を患者に投与することで症状の改善が見られるかを測りたい状況をイメージするとよいでしょう.
因果推論の分野において, 対象への働きかけを介入処置と呼び, 介入によって生じた効果を因果効果と呼びます.

1.1 ポテンシャルアウトカムフレームワーク

ポテンシャルアウトカムフレームワーク(potential outcome framework)は, 因果効果を統計的に取り扱う考え方です. 潜在的な観測量(potential outcome)とは, 介入が行われていない場合に観測される結果を指します.

潜在的な観測量を表す確率変数をY^{(0)}, 介入時の観測量を表す確率変数をY^{(1)}とします.
この時, 介入の効果は次の量で表されることになります(下図参照).

\tau = Y^{(1)} - Y^{(0)}

f:id:kimigayoseishou:20200715225206p:plain
因果効果

この時, 平均的な因果効果は次のように表されます,

\hat{\tau} = E[Y^{(1)} - Y^{(0)}]

これは, 平均処置効果(ATE/Average Treatment Effect)と呼ばれる量です.

平均処置効果を推定するには, 十分な数のサンプルを用意し各サンプルに対してY^{(1)} - Y^{(0)}を観測しそれらの平均値を取るという方法が考えられるでしょう.

しかし, 実際のケースとして各サンプルに対して介入時の結果と介入なしの場合のどちらかしか 観測できないことが多く, 上記のような平均処置効果の推定は現実的ではありません(因果推論の根本問題).

そこで, 同一の対象に対して介入時の結果と介入なしの場合の測定をするのではなく,
介入が行われたグループと介入が行われなかったグループに分た各グループ間の比較による因果効果の推定(対照実験)が応用上重要となります.

1.2 選択バイアスと無作為化実験

対照実験によって得られたデータを用いて因果効果を分析する場合, 意図せず介入の仕方が偏ってしまう場合があります.
例えば, 血圧を下げる効果が期待される新薬の治験において, 新薬を投与した患者のグループの投薬前の血圧が他方のグループより高い数値であった場合などがこれにあたります.

このような状況において, 結果へ影響を与えるような性質の偏りが比較したいグループ間に存在することによって因果効果が正しく推定できないことを選択バイアスといいます.
以下, どのようにして選択バイアスが発生するのか説明します.

記号

  • N : サンプルサイズ
  • Z_i (i = 1, \ldots, N) : 介入の有無を表す変数(介入ありの場合1, そうでない場合0).
  • Y_i (i = 1, \ldots, N) : 観測結果
  • Y^{(1)}_i (i = 1, \ldots, N) : 介入があった場合の観測結果
  • Y^{(0)}_i (i = 1, \ldots, N) : 介入がない場合の観測結果

介入の効果を推定するのに, 介入を行ったグループの平均値とそうでないグループの平均値の差もしばしば用いられます,

\hat{\tau}_{naive} = \frac{1}{\sum_{i=1}^N Z_i}\sum_{i = 1}^N Z_iY_i - \frac{1}{\sum_{i=1}^N (1- Z_i)}\sum_{i = 1}^N (1-Z_i)Y_i

サンプルサイズNを増やす時, グループ間の平均値の差は以下の量に等しくなります.

\tau_{naive} = E[Y^{(1)} | Z = 1 ] - E[Y^{(0)} | Z = 0]

\tau_{naive}\tauを比較するために式変形を行います,


\tau_{naive}  \\
 = E[Y^{(1)} | Z = 1] - E[Y^{(0)} | Z = 0] \\
 = E[Y^{(1)} | Z = 1] - E[Y^{(0)} | Z = 1] + E[Y^{(0)} | Z = 1] - E[Y^{(0)} | Z = 0] \\
= E[Y^{(1)} - Y^{(0)} | Z = 1] + E[Y^{(0)} | Z = 1] - E[Y^{(0)} | Z = 0]

ここで, 因果効果Y^{(1)} - Y^{(0)}と介入の有無が独立だと仮定することにより E[Y^{(1)} - Y^{(0)} | Z = 1] = E[Y^{(1)} - Y^{(0)}] = \hat{\tau}となるため, 平均処置効果とグループ間の平均差はE[Y^{(0)} | Z = 1] - E[Y^{(0)} | Z = 0]だけずれることがわかります. この差は偏った介入の仕方に起因するもので選択バイアス(selection bias)と呼ばれています.

例えば, Z_iY^{(0)}_iに正の相関があるような場合, Y^{(0)}_iの大きい部分ではZ=1の傾向が強くなりますが, そのようなサンプルは観測されません. その結果, 介入なしのグループのサンプルの平均値は真の平均値E[Y^{(0)}]よりも 小さく推定されてしまいます(下図).

f:id:kimigayoseishou:20200715225333p:plain
選択バイアス

選択バイアスによる影響をなるべく減らすために対象への介入を無作為に割り当てて行う対照実験は 無作為化比較実験(RCT : Randomized Controlled Trial)と呼ばれています.

1.3 Rによる選択バイアスの数値実験

ここでは介入に偏りのあるデータを作成し, 選択バイアスを擬似的に発生させる数値実験をRを用いて行ってみます.

(a)介入が無作為である場合

まずは介入が無作為である場合です. 平均0, 標準偏差0.5の正規分布から200個のサンプルを生成し, これを介入が行われていない場合のデータとします(y0). また介入によって結果は1増えると仮定し, 介入が行われた場合のデータを作成します(y1). これらのデータy0, y1に対して, 無作為に介入を行い, 介入の有無によって2つのグループに分けてその平均値を比較します.

library(dplyr)
library(ggplot2)
set.seed(1)
# 介入なし
y0 <- rnorm(200, 0, 0.5)
# 介入が行われた場合, 観測結果は1増加
y1 <- y0 + 1
# 無作為な介入
z <- sample(c(0, 1), 200, TRUE)

# 介入は無作為になっている
data <- tidyr::tibble(Y = y1 * z + y0 * (1-z), Z = z)

data %>%
    group_by(Z) %>%
    summarise(mean_Y = mean(Y))

# 赤点は各グループの平均, 赤縦線は平均 ± 標準偏差の範囲を示す
ggplot(data, aes(x = as.factor(Z), y = Y)) + 
    geom_point(alpha=.4) +
    stat_summary(
        fun.y = mean,
        fun.ymin = function(x) mean(x) - sd(x), 
        fun.ymax = function(x) mean(x) + sd(x),
        color = "red"
        ) +
    xlab("介入の有無(介入あり Z = 1, 介入なし Z = 0)") + 
    ggtitle('介入ありグループと介入なしグループの平均値の比較(RCT)')

f:id:kimigayoseishou:20200715225417p:plain
RCT グループごとの平均値

出力結果を確認すると, 介入なしグループの平均値は約-0.02, 介入ありグループの平均値は約1.06となっており, その平均値の差=1.06 - (-0.02) =1.08は真の因果効果=1と近い値になります.

f:id:kimigayoseishou:20200715225446p:plain
RCTプロット

(b)介入に偏りがある場合

次に介入に偏りがある場合を考えます. 以下のコードでは, 「介入なしの場合の値が0より大きい場合は必ず介入」 というルールを追加して意図的に介入が偏ったデータを作成します.

# 意図的に選択バイアスを発生させる
data_biased <- tidyr::tibble(Z = ifelse(y0 > 0, 1, z)) %>%
    mutate(Y = y1 * Z + y0 * (1-Z))

data_biased %>%
    group_by(Z) %>%
    summarise(mean_Y = mean(Y))

# 赤点は各グループの平均, 赤縦線は平均 ± 標準偏差の範囲を示す
ggplot(data_biased, aes(x = as.factor(Z), y = Y)) + 
    geom_point(alpha=.4) +
    stat_summary(
        fun.y = mean,
        fun.ymin = function(x) mean(x) - sd(x), 
        fun.ymax = function(x) mean(x) + sd(x),
        color = "red"
        ) +
    xlab("介入の有無(介入あり Z = 1, 介入なし Z = 0)") + 
    ggtitle('介入ありグループと介入なしグループの平均値の比較(非RCT)')

f:id:kimigayoseishou:20200715225522p:plain
非RCT 平均値

出力結果を確認すると, 介入なしグループの平均値は約-0.34, 介入ありグループの平均値は約1.17となっており, その平均値の差=1.17 - (-0.34) = 1.51は真の因果効果=1と比較して大きく推定されてしまいます. このデータではy0 > 0のデータについて必ず介入を行ったことでy0 > 0のサンプルが欠測となり, その結果介入なしのグループ平均値は真の平均値E[Y^{(0)}=0]よりも小さく推定されていることがわかります.

f:id:kimigayoseishou:20200715225608p:plain
非RCT プロット

上記の例をまとめると以下のようになります.

  • RCTの場合は, グループ毎の平均値の差によって平均処置効果を推定できる.
  • 非RCTの場合は, 選択バイアスによりグループ毎の平均値の差は平均処置効果と解離することがある.

特にグループ毎の平均値を用いて因果効果の推定を行う場合は, 介入の仕方に偏りがないか注意を払う必要があります.


2 傾向スコア

傾向スコア(propensity score)は, 無作為でない介入が行われたサンプルからでも因果効果を推定できる手法として広く用いられています. 以下, 傾向スコアについて最小限の事項をまとめています.

詳細については,

などなどを参照されたし....

2.1 バランシングスコアと傾向スコア

以下傾向スコアを導入するにあたり, いくつか用語の準備をします.

記号

  • X : 説明変数
  • Y : 応答変数
  • Y^{(1)} : 介入が行われた場合の結果
  • Y^{(0)} : 介入が行われな買った場合の結果
  • Z : 介入の有無を表す変数(介入ありの場合1, そうでない場合0)

定義-バランシングスコア(balancing score)

ベクトル値関数b(X)バランシングスコア(balancing score)であるとは, 次の条件が成り立つことをいう

X  \perp Z | b(X).

定義-傾向スコア(propensity score)

xの関数P:x \mapsto p(Z = 1|X =x)傾向スコア(propensity score)と定義する.

定義-強く無視できる割り当て(strongly ignorable treatment assignment)

介入の割り当てが確率変数Wのもとで強く無視できる割り当て(strongly ignorable treatment assignment)であるとは, 次の条件が満たされることをいう

(Y^{(1)}, Y^{(0)}) \perp Z | W.

さらにXに関して強く無視できる割り当てである場合, 単に強く無視できる割り当てと呼ぶ.

以上の定義のもとで次が成り立ちます.

定理
1. 傾向スコアはバランシングスコアである.
2. 介入の割り当てが強く無視できる割り当て => 任意のバランシングスコアのもとで強く無視できる割り当てになっている.

証明

こちらのスライドを一部参考にさせていただきました.

1 p(Z = 1, X=x'|P = e) = p(Z = 1|P=e)p(X=x'|P=e)を示せば良い. P(x') \neq eの場合は両辺とも0となる. 以下, P(x') = eと仮定する. このとき, 示すべき等式の左辺は次のようになる. 
p(Z = 1, X = x'|P = e) \\
= p(Z = 1|X = x', P = e)p(X=x'|P = e) \\
= p(Z = 1|X = x')p(X=x'|P = e) \\
= P(x')p(X=x'|P = e) \\
= ep(X=x'|P = e)

また, p(Z=1 | P=e)Xに関して周辺化し計算を行う. 
p(Z = 1|P = e) \\
 = \int p(Z = 1, X=x|P = e) dx \\
= \int p(Z = 1| X=x, P = e) p(X = x|P=e)dx \\
= \int_{P(x) = e} p(Z = 1, | X=x, P = e) p(X = x|P=e)dx \\
= \int_{P(x) = e} p(Z = 1, | X=x) p(X = x|P=e)dx \\
= \int_{P(x) = e} P(x) p(X = x|P=P(x_2))dx \\
= e \int_{P(x) = e} p(X = x|P=e)dx \\
= e

よって, 示すべき等式の左辺=ep(X=x'|P = e) = p(Z = 1|P = e)p(X=x'|P = e)=示すべき等式の右辺となり結果が従う.

2 p(Z=z, (Y^{(1)}, Y^{(0)}) = (y_1, y_0)|b = b') = p(Z=z|b = b')p((Y^{(1)}, Y^{(0)}) = (y_1, y_0)|b = b')を示せば良い. 以下のようにXについて周辺化を行うことで, 示すことができる. 
p(Z=z, (Y^{(1)}, Y^{(0)}) = (y_1, y_0)|b = b') \\
= \int p(Z=z, (Y^{(1)}, Y^{(0)}) = (y_1, y_0), X = x|b = b') dx \\
= \int p(Z=z, (Y^{(1)}, Y^{(0)}) = (y_1, y_0)|X = x, b = b') p(X = x|b=b')dx \\
= \int_{b(x) = b'} p(Z=z, (Y^{(1)}, Y^{(0)}) = (y_1, y_0)|X = x) p(X = x|b=b')dx \\
=  \int_{b(x) = b'} p(Z=z|X=x)p((Y^{(1)}, Y^{(0)}) = (y_1, y_0)|X = x) p(X = x|b=b')dx \quad (\because (Y^{(1)}, Y^{(0)}) \perp Z | X)\\
= \int_{b(x) = b'} p((Y^{(1)}, Y^{(0)}) = (y_1, y_0)|X = x) p(Z = z, X = x|b=b')dx \\
= \int_{b(x) = b'} p((Y^{(1)}, Y^{(0)}) = (y_1, y_0)|X = x) p(Z = z|b=b')p(X=x|b=b')dx \quad (\because X \perp Z |b)\\
= p(Z = z|b=b') \int_{b(x) = b'} p((Y^{(1)}, Y^{(0)}) = (y_1, y_0)|X = x) p(X=x|b=b')dx \\
= p(Z = z|b=b') \int p((Y^{(1)}, Y^{(0)}) = (y_1, y_0), X = x|b = b')dx \\
= p(Z = z|b=b')p((Y^{(1)}, Y^{(0)}) = (y_1, y_0)|b = b').

以上の結果より, 介入の割り当てが強く無視できる割り当てであるとき, E[Y^{(i)}|Z=i, P]\ (i=0, 1)は真の平均値E[Y^{(i)}]の不偏推定量になっていることが次のようにして導かれます,

E[Y^{(i)}] = E[E[Y^{(i)}|P]] = E[E[Y^{(i)}|Z=i, P]] \quad (\because Z \perp Y^{(i)} | P).

注目すべきはY^{(i)}|Z=i, Pは欠測とはならないため, (傾向スコアの推定を行った上で)サンプルから推定が行える点です.

以上をまとめると, 傾向スコアを用いて平均処置効果を推定する大まかな流れは以下のようになります.

  1. 介入の割り当てが強く無視できる割り当てとみなせるような説明変数Xの選定を行う.
  2. (GLM等のモデルを用いた)傾向スコアP(x) = p(Z = 1 | X =x)の推定.
  3. 各傾向スコアPに対して, E[Y^{(i)}|Z=i, P]をサンプルから推定し, 真の平均値E[Y^{(i)}]を推定する.
  4. 3で推定した値の差をとった量を平均処置効果の推定量とする.

2.2 逆確率重み付き推定

逆確率重み付き推定(IPW : Inverse Probability Weighting)は傾向スコアを利用した真の平均値E[Y^{(i)}]の推定法の一つです.

Nをサンプルサイズ, X_iをサンプルiにおける説明変数, Y_i(i = 1, \ldots, N)をサンプルiにおける観測値, Z_i(i = 1. \ldots, N)をサンプルiにおける介入の有無とする. サンプルiの傾向スコアをP_i = p(Z = 1|X = X_i)とするとき, \hat{Y}^{(i)}_N(i=0, 1)

\hat{Y}^{(1)}_N = \sum_{j=1} ^N \frac{Y_iZ_i}{P_i} \bigr/ \sum_{j=1}^N \frac{Z_i}{P_i} \\
\hat{Y}^{(0)}_N = \sum_{j=1} ^N \frac{Y_i(1-Z_i)}{1-P_i}  \bigr/  \sum_{j=1}^N \frac{1-Z_i}{1-P_i}

と定義します.

命題-(逆確率重み付き推定)
\hat{Y}^{(i)}_N(i=0, 1)E[Y^{(i)}](i=0, 1)の一致推定量である.

証明

\hat{Y}^{(i)}_N \rightarrow E[Y^{(i)}](N \rightarrow \infty)を示せば良い. i=1の場合のみ示す.


E[YZ|P] \\
= E[(ZY^{(1)} + (1-Z)Y^{(0)})Z | P] \\
= E[ZY^{(1)}| P] \\
= E[Z| P]E[Y^{(1)}|P] \quad ( \because Y^{(1)} \perp Z | P) \\
= P E[Y^{(1)}|P]

であるので, 大数の法則より


\frac{1}{N}\sum_{j=1}^N\frac{Y_iZ_i}{P_i} \\
\rightarrow E \biggl[ \frac{YZ}{P} \biggr] \quad (N \rightarrow \infty)\\
= E \biggl[ \frac{E[YZ | P]}{P} \biggr] \\
= E \biggl[ \frac{PE[Y^{(1)}|P] }{P} \biggr] \\
= E [E[Y^{(1)}|P] ] \\
= E[Y^{(1)}] \\

となる. 再び大数の法則より,


\frac{1}{N}\sum_{j=1}^N \frac{Z_i}{P_i} \\
\rightarrow E \biggl[ \frac{Z}{P} \biggr] \quad (N \rightarrow \infty)\\
= E \biggl[ \frac{E[Z|P]}{P} \biggr] \\
=  E \biggl[ \frac{P}{P} \biggr] = 1

以上より, 
\hat{Y}^{(i)}_N = \frac{ \frac{1}{N}\sum_{j=1}^N\frac{Y_iZ_i}{P_i}}{\frac{1}{N}\sum_{j=1}^N \frac{Z_i}{P_i}} \rightarrow \frac{E[Y^{(1)}]}{1} = E[Y^{(1)}] \quad (N \rightarrow \infty).
となるから, 確かに\hat{Y}^{(1)}E[Y^{(1)}]の一致推定量である.

逆確率重み付き推定の数値実験

ここでは, 擬似データを用いて逆確率重み付き推定による因果効果の推定の数値実験を行います.

あるアパレル小売店Aでは, 販売促進のためにクーポン券の配布を行いました. そこでクーポンの利用によってクーポン利用可能期間中のユーザーの購買額が増加したかを分析したいとします.

  • X_i : ユーザーiの衣料品に対する平均支出額
  • Y^{(0)}_i : ユーザーiのクーポンがない場合の購入金額
  • Y^{(1)}_i : ユーザーiのクーポンを利用した場合の購入金額
  • Y_i : ユーザーiの実際の購入金額
  • Z_i : クーポン利用の有無(クーポンを利用した場合1, そうでない場合0)

クーポンがない場合の購入金額は衣料品に対する平均支出額を説明変数とした線形モデルで表せるとし、 さらにクーポン利用時は\lambdaだけ購入金額が増加すると仮定します.

Y^{(0)}_i \sim {\sf Normal}(\alpha_0 + \alpha_1 X_i, \tau^{2}) \quad (\tau  {\text >}0) \\
Y^{(1)}_i = Y^{(0)}_i + \lambda.

次にクーポン利用についてです. ここでは, クーポンの利用者は非利用者よりも衣料品支出額が多いという仮定をおきます.

クーポン非利用者の衣料品に対する平均支出額は平均\mu_0, 分散\sigma^{2}正規分布に, クーポン利用者の衣料品に対する平均支出額は平均\mu_1, 分散\sigma^{2}正規分布に, それぞれ従うとします

X|Z=j \sim {\sf Normal}(\mu_j, \sigma^{2})\quad (\mu_j  {\text >} 0, \sigma^{2}  {\text >} 0, j= 0, 1),

以上のモデルにしたがって200人分(うちクーポン利用者100人)の擬似データを作成します. 各パラメータは次のようにしています,

\alpha_0=0.1, \alpha_1 = 0.05, \tau = 0.05, \lambda = 0.2, \mu_1 = 10, \mu_0 = 6, \sigma = 2

# 擬似データの作成
library(dplyr)
set.seed(1)
alpha0 <- .1
alpha1 <- .05
tau <- .1
lambda <- .2
N <- 200L
mu1 <- 10
mu0 <- 6
sigma <- 2
# クーポン利用者の衣料品に対する平均支出額
x1 <- rnorm(N/2, mu1, sigma)
# クーポン非利用者の衣料品に対する平均支出額
x0 <- rnorm(N/2, mu0, sigma)
# 全ユーザーの衣料品に対する平均支出額
x <- c(x0, x1)

# クーポン非利用時の衣料品に対する平均支出額
y0 <- alpha0 + alpha1 * x + rnorm(N, 0, tau)
y1 <- y0 + lambda
z <- c(rep(0L, N/2), rep(1L, N/2))

# 200件の購買データ
purchase <- tidyr::tibble(Z = z, Y = z * y1 + (1-z) * y0, X = x)

purchase %>%
    ggplot(aes(x = Y, fill = Z)) +
    geom_histogram(alpha = .6, position = "identity") +
    ggtitle('実際の購入金額の分布')

f:id:kimigayoseishou:20200715225733p:plain
購買金額のヒストグラム

まず, 介入ありグループと介入なしグループで平均値の差を確認します.

lm(data = purchase, formula = Y ~ Z) %>%
    broom::tidy()

f:id:kimigayoseishou:20200715225656p:plain
グループ毎平均値

出力結果のterm=Zの行のestimateの値(=0.42)は, グループ毎の平均値による平均処置効果の推定値になります. この値は真の平均処置効果\lambda = 0.2の約2.1倍の数値で, 選択バイアスを含んだ値となっています.

そこで選択バイアスの影響を除去するため, 逆確率重み付き推定による平均処置効果の推定値を行います.

まずクーポンの利用率p(Z = 1) = 0.5として, 傾向スコアを計算すると,


P(x) \\
= p(Z = 1|X=x) \\
= \frac{p(X = x|Z =1)p(Z=1)}{p(X = x|Z =0)p(Z=0) + p(X = x|Z =1)p(Z=1)} \\
= \frac{p(X = x|Z =1)/p(X = x|Z =0)}{1 + p(X = x|Z =1)/p(X = x|Z =0)} \\
= \frac{e^{\beta_0 + \beta_1x}}{1 + e^{\beta_0 + \beta_1x}}

となり傾向スコアはロジスティック関数で表されます. ただし, ここで

\beta_0 = - \frac{\mu_1 ^{2}-\mu_0 ^{2}}{2\sigma ^{2}}, \beta_1 = \frac{\mu_1-\mu_0}{\sigma^{2}
}

としています.

注意 今の状況下では, \beta_0, \beta_1\mu_0, \mu_1, \sigmaからカンニングできてしまいますが, 実際のデータ分析時にはデータから推定する必要があります.

この式から傾向スコアの算出を行います.

# 傾向スコアの計算
beta0 <- -.5 * (mu1^2-mu0^2)/sigma^2
beta1 <- (mu1 - mu0) / sigma^2
prop_score <- 1 / (1 + exp(-(beta0 + beta1 * purchase$X)))

# 逆確率重み付き推定
# クーポン利用時の購買金額推定値
y1_estim  <- sum(purchase$Y * purchase$Z / prop_score) / sum(purchase$Z / prop_score)

# クーポン非利用時の購買金額推定値
y0_estim <- sum(purchase$Y * (1 - purchase$Z) / (1 - prop_score)) / sum((1 - purchase$Z) / (1 - prop_score))

# 逆確率重み付き推定によるクーポン利用時/非利用時の購買金額推定値の差が, 
# 平均的な因果効果の推定値となる
y1_estim - y0_estim

f:id:kimigayoseishou:20200715225816p:plain

逆確率重み付き推定による因果効果の推定値は0.27となり, 真の平均処置効果\lambda = 0.2の約1.4倍の数値です. 選択バイアスの影響が取り除かれた分, グループ毎の平均値による推定値よりも小さく推定されていることが確認できます.


3. まとめ

かなり急ぎ足で, ポテンシャルアウトカムフレームワークにおける因果効果の推定の方法について 見てきました. 今回Rでの分析で用いたデータは擬似データのみでしたので面白みは感じられないと思いますが, 傾向スコアを用いた実証研究などについてさらに調べて見たいと思います.

【R】Search Consoleレポートからクリック率を高めるキーワードを推定する

1. 検索キーワードとクリック率

Search Consoleで取得できる検索パフォーマンスレポートでは, Search Consoleを導入済みウェブサイトがどのような検索キーワードからアクセスされているかを確認できます.

このデータを元に「どんなキーワードで自分のウェブサイトがクリックされやすいか?」を意識して サイト内容を見直すことで, 検索エンジンからのアクセス数を改善することができます.

本記事では検索パフォーマンスレポートのデータを基に, 検索に使われる高クリック率なキーワードの分析をRを使って試してみます.

f:id:kimigayoseishou:20200610132907p:plain

2. モデル化

2.1 ロジスティック回帰

検索クエリは, 複数キーワードの組み合わせ(例えば, 「おいしいカレー / 作り方」)からなります. キーワードごとのクリック率の傾向を見るためには検索クエリをキーワード単位で分解する必要があります.
そこでまず, 検索クエリ内に現れるキーワードを予めリストアップしておき, w1, \ldots, w_Lとします. 次に検索クエリQに対して,

q_j = \begin{cases} 1 ,  \quad ({\rm キーワード}w_j{\rm が}Q{\rm に含まれる}) \\ 0 , \quad ({\rm  上記以外})\end{cases}

とすることで, 検索クエリをキーワードの有無を各成分にもつベクトルq = (q_j) \in {0,1}^Lに変換できます.

検索クエリQで検索を行ったユーザーが自分のサイトにアクセスした場合1, そうでない場合0とするbinaryな確率変数をYとします.

検索クエリと自分のサイトのクリック率の関係を探るためモデルを導入します. 今回は最も単純なモデルとして, Yを目的変数, クエリQのベクトル化qを説明変数としたロジスティック回帰(logistic regression) :

\begin{align} Y | Q \sim {\sf Binomial}(\pi), \quad \pi = \frac{\exp(\beta_0 + \beta ^\top q)}{1 + \exp(\beta_0 + \beta ^\top q)} \end{align}

を採用します, ただし\beta_0 \in {\mathbb R}, \beta = (\beta_1, \ldots, \beta_L) \in {\mathbb R}^Lはパラメータです.

\pi = p(Y = 1 |Q)はクリック率に対応するため, \beta_jが大きいキーワードw_jほど, 高クリック率になりやすいキーワードだと考えられます.

2.2 正則化

N個の検索クエリ\underline{Q} = (Q_1, \ldots, Q_N)と実際にクリックされたかどうかの結果\underline{y} = (y_1, \ldots, y_N) \in {0, 1} ^ Nを観測値としてロジスティック回帰モデルのパラメータの推定を行いたい.

最尤推定では対数尤度l(\beta_0, \ \beta) = \log p(\underline{y}|\underline{Q})の最大化によりパラメータ推定を行いますが, 今回の観測データのように検索キーワードを説明変数とするモデルでは, 説明変数が多くなり過学習に陥ってしまう恐れがあります.

このような場合, 適当な罰則項を考慮した最適化問題

\begin{align} \underset{(\beta_0, \ \beta) \in {\mathbb R}^{L+1}}{{\sf argmin}}\ - \frac{1}{N} l(\beta_0, \ \beta) + \lambda ({\text 罰則項}) \end{align}

によるパラメータ推定を行うことで過学習を回避できることがあります(正則化). 代表的な罰則項としては次のようなものがあります


3. glmnetを利用したLasso推定

以下, Ridge/Lasso(及びその線型和であるElastic Net)正則化回帰モデルを扱うRパッケージ glmnetパッケージを利用して, 検索パフォーマンスレポートのデータの分析を行います.

使用するデータは当ブログの過去3ヶ月のデータですが, Search Consoleを利用している方は 自サイトのデータで分析するのもいいかもしれません.

3.1. データの確認

# require(tidyverse)
# データ読み込み
search_console_dat <- readxl::read_excel('自分のブログ_Performance-on-Search-2020-06-06.xlsx')

search_console_dat %>%
    mutate(検索キーワード = reorder(検索キーワード, CTR)) %>%
    top_n(10, CTR) %>%
    ggplot(aes(x = as.factor(検索キーワード), y = CTR)) +
    geom_bar(stat = "identity") +
    xlab("検索クエリ") +
    ylab("クリック率") +
    coord_flip() + 
    ggtitle("クリック率上位10件の検索クエリ")

3.2. 前処理

次に前処理として, 検索クエリをキーワード単位でベクトル化します.

# 検索クエリ内に現れるキーワードをリスト化
words_in_query <- search_console_dat$検索キーワード %>%
    strsplit(" ") %>%
    unlist() %>%
    unique()

wordlist_in_query <- search_console_dat$検索キーワード %>%
    lapply(function(str) as.numeric(str_detect(str, words_in_query)))

# を列とし, 各検索キーワードを行とする行列
# 各行の検索キーワードがある列に対応する単語を含む場合1, そうでなければ0
word_in_query <- rlang::exec(rbind, !!!wordlist_in_query) 
colnames(word_in_query) <- words_in_query

count_word_in_query <- as_tibble(word_in_query) %>%
    add_column(Imp = search_console_dat$表示回数, .before = 1) %>%
    add_column(click = search_console_dat$クリック数, .before = 1)

3.3 パラメータ推定

それでは実際にパラメータ推定を行います. 今回はLasso推定を実際に試してみます. Lasso回帰を行うには, glmnet::glmnet関数, もしくはglmnet::cv.glmnet関数を利用します. 特に, glmnet::cv.glmnet関数では罰則項の大きさを指定するパラメータまで交差検証により 自動で推定が行われます.

# see https://cran.r-project.org/web/packages/glmnet/vignettes/glmnet.pdf
require(glmnet)
imp <- search_console_dat$表示回数
click <- search_console_dat$クリック数

# 交差検証により自動でlambdaを調整してくれる
fit_lasso.cv <- cv.glmnet(x = as.matrix(count_word_in_query[,-c(1:2)]), 
                          y = cbind(imp-click, click), 
                          nfolds = 4, family = "binomial", alpha=1)

glmnet::cv.glmnetの引数の意味は以下の通りです

  • x:説明変数(matrixで指定する)

  • y:目的変数

  • nfolds:交差検証の分割数

  • family:回帰モデルの指定("gaussian", "binomial", "poisson", "multinomial", "cox", "mgaussian"のいずれかが指定できる)

  • alpha:罰則項の形を決定する. 具体的な形は\alpha ||\beta||_1 + (1-\alpha)||\beta||^2_2 / 2となる.

3.4 推定結果の確認

最後に推定結果の確認を行います.

plot(fit_lasso.cv)

coef(fit_lasso.cv, s="lambda.min") %>% 
    broom::tidy() %>%
    mutate(row = reorder(row, value)) %>%
    ggplot(aes(x = as.factor(row), y = value)) +
    geom_bar(stat = "identity") +
    xlab("検索キーワード") +
    ylab("回帰係数の推定値") +
    coord_flip() 

注目すべきキーワードとして「mac」や「catalina」といったPCのOSに関するキーワードが 回帰係数が大きく, 高クリック率に寄与するという点が挙げられます.

このブログのようなプログラミングに関する内容のサイトでは, OSキーワードを記事に含めることで競合サイトが減り, 結果クリック率が上がるのかもしれません.

f:id:kimigayoseishou:20200610133119p:plain
Lasso推定による回帰係数の推定結果

Power Queryの操作をpowershellで自動化する

Power Queryはデータの取得・加工を行うことができるExcelアドインツールで、 Excel2016以降のExcelであれば標準搭載されています.

普段から利用している人も多いのではないでしょうか.

さてPower QueryはGUIベースでも十分使いやすい代物ですが, 自動化を行うことでより強力なツールとなります.

以下では, powershellを使ったpower queryの自動化について紹介します.

power queryをpowershellから実行する

以下はファイルシステム内の複数のExcelデータを一つの表にまとめるpowershellスクリプトになります.

####################################################################################
# Power Queryを利用して複数Excel表データを結合するスクリプト
####################################################################################

# 設定 #############################################################################

# Excelのパスとシート名を記述
$DATA_SOURCE_CONFIG = @(
        @{"path"="<PATH_TO_XLSX_FILE1>"; "sheetName"="<SHEETNAME>"},
        @{"path"="<PATH_TO_XLSX_FILE2>"; "sheetName"="<SHEETNAME>"}
        # etc...
    )

# Excelの出力先
$OUTPUT = "$($PSScriptRoot)\report.xlsx"

# 設定ここまで ######################################################################
function mFormulaBuilder {
    param (
        $dataSourceList
    )

    $dataSourceObjList = $dataSourceList | % { "Excel.Workbook(File.Contents(""$($_.path)""), true){[Name=""$($_.sheetName)""]}[Data]" }

    return "Table.Combine({$($dataSourceObjList -join ",")})"
}

function main {
    
    $excel = New-Object -ComObject Excel.Application
    $book = $excel.Workbooks.Add()

    # クエリの登録
    $QUERY_NAME = "query1"
    $mFormulaString = mFormulaBuilder($DATA_SOURCE_CONFIG)
    $book.Queries.Add($QUERY_NAME, $mFormulaString)

    $sheet = $book.Worksheets.Item(1)
    # Sheet1のA1にQueryTableを作成
    $qt = $sheet.QueryTables.Add(
        # 接続文字列
        "OLEDB;Provider=Microsoft.Mashup.OleDb.1;Data Source=`$Workbook`$;Location=${QUERY_NAME};Extended Properties=""""",
        # 表範囲
        $sheet.Range('$A$1')
    )
    # sql文
    $qt.CommandText = "SELECT * FROM [${QUERY_NAME}]"
    # sql文実行
    $qt.Refresh($false)

    # ワークブックを保存する
    $book.SaveAs($OUTPUT)

    $excel.Quit()
    $book = $null
    $excel = $null
    
    [GC]::Collect()
}

main

powershell(>= 6.0)であれば, 以下のコマンドを実行することで複数のExcelファイルのデータを1つにまとめたExcelファイルが出力されます.

pwsh ./bind_excel_data.ps1

One drive上のExcelファイルへのアクセス

前節では, ファイルシステム管理下のExcelファイルデータのアクセスをおこないました.

以下の記事のよう法人向けOne driveはHTTPSプロトコルによるファイルアクセスに対応しています.

Power Query から https 経由で SharePoint の Excel ブックを開く

前節のスクリプトファイルに対し以下のようにFile.ContentsWeb.Contentsと修正を加えることで, (法人向け)One drive上の複数Excelファイルを統合できるようになります.

修正前

function mFormulaBuilder {
    param (
        $dataSourceList
    )

    $dataSourceObjList = $dataSourceList | % { "Excel.Workbook(File.Contents(""$($_.path)""), true){[Name=""$($_.sheetName)""]}[Data]" }

    return "Table.Combine({$($dataSourceObjList -join ",")})"
}

修正後

function mFormulaBuilder {
    param (
        $dataSourceList
    )

    $dataSourceObjList = $dataSourceList | % { "Excel.Workbook(Web.Contents(""$($_.path)""), true){[Name=""$($_.sheetName)""]}[Data]" }

    return "Table.Combine({$($dataSourceObjList -join ",")})"
}

参考リンク

【R】エーレンフェストモデル

1.エーレンフェストモデル

エーレンフェストモデル(Ehrenfest model)は, 外界から閉ざされた容器内での気体の拡散を説明する場面に使われる確率モデルの一つで、 マルコフ連鎖の例としてしばしば取り上げられます.

以下、Rでのシミュレーションと合わせてこのモデルの基本性質について取り上げます.

1.1 定式化

2つの箱(それぞれA, Bとします)内に合計r個のボールが入っており, ボールには1~rまでの番号が書かれているとします. 次のルールに従ってボールの遷移を考えます.

  • ルール: r個のボールから1つをランダムに選び, もう他方の箱へ移す.

この遷移がt回行われたときの箱A内のボールの個数X(t)は、 S = \{ 0, 1, \ldots, r \} 上のマルコフ連鎖になります.

この遷移確率は次で与えられることがわかります,

p_{ij}:=p(X(t+1)=j \ | \ X(t)= i) = \begin{cases} \frac{i}{r}  \quad (j=i-1) \\ \frac{r-i}{r}  \quad (j=i+1) \\ 0 \quad  (\text{otherwise}). \end{cases}

1.2 定常分布

十分に長い時間が経過したとき, X(t)の分布は定常分布\piに近づきます. \pi_i = \pi(X=i)としたとき, 定常分布(equilibrium distribution)は 次の詳細釣り合いの式から求めることが可能です,

\pi_i p_{ij} = \pi_j p_{ji} \quad (0 \leq i,j \leq r)

特にj=i-1とすることで次が得られます,

\pi_i = \frac{p_{i-1i}}{p_{ii-1}}\pi_{i-1} = \frac{r-i+1}{i}\pi_{i-1}=\cdots = \binom{r}{i}\pi_0.

分布の条件1 = \pi_0+\cdots+\pi_rから \pi_0 = 1 / 2^{r} ともとまる.

エーレンフェストモデルの定常分布 エーレンフェストモデルの定常分布は二項分布  \begin{equation} {\sf Binom}\bigl( r, \frac{1}{2} \bigr) \end{equation} である.

1.3 Rを使った定常分布のシミュレーション

require(tidyverse)
ehrenfest.path <- function(r, len, init) {
    path <- c(init)
    for (t in 1:(len-1)) {
        x_next <- ifelse(rbinom(1, 1, path[t]/r) == 1, path[t]-1, path[t]+1)
        path <- c(path, x_next)
    }
    return(path)
}

# 長さ4000のパス
sample.path <- ehrenfest.path(r=100, len=4000, init=0)

data.frame(value=sample.path) %>%
    ggplot(aes(x=value)) + 
    geom_density(fill="blue", alpha=.4) +
    geom_col(aes(x=x, y=y), data.frame(x=c(0:100), y=dbinom(c(0:100), size=100, prob=.5)), alpha=.7) +
    xlim(c(0,100)) +
    ggtitle('histgram of path v.s. equilibrium distribution')

f:id:kimigayoseishou:20200427200246p:plain
サンプルの密度分布(薄青)と, 定常分布(試行回数=100, 成功確率=1/2の二項分布)の比較


2.1 ギブスエントロピー

有限個の状態をとりうる系について, ギブスエントロピー(Gibbs entropy)をつぎの式で定義されます.

H = - \sum_i p_i \ln p_i

熱力学第二法則によれば, 不可逆な変化のもとでギブスエントロピーは増大します. エーレンフェストモデルでは箱A内のボール個数の分布は定常分布へ近づくことになり, 不可逆な変化が発生します. 以下では, エーレンフェストモデルのギブスエントロピーが時間と共に増加することをシミュレーションにより観察します.

2.1 エーレンフェストモデルのシミュレーション

エーレンフェストモデルはマルコフ連鎖に従うため,  X(t)の確率分布から X(t+1)の確率分布を次のように明示的に表すことができます,

p(X(t+1) = i) \\ = \sum_{j=0}^r p(X(t+1)=i,\ X(t) = j) \\ = \sum_{j=0}^r p(X(t+1)=i | \ X(t) = j) p(X(t) = j) \\ =  p(X(t)=i-1)p_{i-1i} + p(X(t)=i+1)p_{i+1i} \\ = \frac{r-(i-1)}{r}p(X(t)=i-1) + \frac{i+1}{r}p(X(t)=i+1)

この式による, エーレンフェストモデルの時間変化のシミュレーションのRコード以下のようになります.

#' @param T_max 
#' @param N 時刻を刻む個数
#' @param r ボールの個数
#' @param init r-vetor 初期値
Ehrenfest.sim <- function(T_max, N, r, init) {
    # 横ベクトルの要素を左にずらす
    shift.p <- function(x) {
        return(c(x, 0)[c(-1)])
    }
    # 横ベクトルの要素を右にずらす
    shift.m <- function(x) {
        return(c(0, x)[c(-length(x)-1)])
    }

    sim <- matrix(0, nrow = T_max, ncol = r+1)
    sim[1,] <- init
    for (i in 2:T_max) {
        sim[i,] <- (1 + 1 / r - c(0:r)/r) * shift.m(sim[i-1,]) + (1 / r + c(0:r) / r) * shift.p(sim[i-1,])
    }
    return(sim)
}

ehrenfest.sim <- Ehrenfest.sim(T_max=200, N=100, r=100, init=c(c(1, 1)/2, rep(0, 99)))

# シミュレーション結果をプロット
plt.sim <- function(sim, time) {
    x.range <- c(0:(ncol(sim)-1))
    plot(x=x.range, y=sim[1,], xlim = c(0, max(x.range)), ylim = c(0, 0.6), type="l", 
         ylab = 'probability')
    
    for (t in time) {
        lines(x=x.range, y=sim[t,], lty="solid")
    }
}

# 10回の遷移づつ確率分布をプロット
# (横軸 - ボールの個数, 縦軸 - 確率)
plt.sim(ehrenfest.sim, c(1:10) * 20)

f:id:kimigayoseishou:20200427200154p:plain
確率質量分布の時間変化

プロットの結果より, 左端に局在していた確率分布が時間の経過と共に定常分布に近づく様子が確認できます.

次にギブスエントロピーのシミュレーションを行います.

gibbs_entropy <- function(x) - sum(x * log(x), na.rm = TRUE)
sim.entropy <- apply(ehrenfest.sim, MARGIN = 1, gibbs_entropy)
plot(sim.entropy, type="l", xlab = 'time', ylab='gibbs entropy')

f:id:kimigayoseishou:20200427200052p:plain

プロット結果から, エーレンフェストモデルのギブスエントロピーが時間と共に増大することが 確認できます.

【参考】

マルコフ連鎖から格子確率モデルへ(Amazonリンク)

【R】StanでMAKE AMERICA GREAT AGAINの変化点検知

2020年はアメリカで大統領選挙が行われる年ですが, 前回の選挙戦でのトランプ氏は「MAKE AMERICA GREAT AGAIN」というスローガンで注目を集めていましたね。 今回の記事では、トランプ氏の「MAKE AMERICA GREAT AGAIN」ツイート投稿数の推移を調べてみました。

トランプ大統領の「MAKE AMERICA GREAT AGAIN」ツイート数の可視化

  • データ
    トランプ氏のツイッターの2015年3月 ~ 2020年2月までの5年間の投稿ツイートのうち, 「MAKE AMERICA GREAT AGAIN」を含む436ツイートを取得した.

「MAKE AMERICA GREAT AGAIN」の可視化

require(tidyverse)
trump_tweets <- read_csv('./trump_tweets.csv')

# 月ごとにツイートを集計
trump_tweets.count <- trump_tweets %>%
  count(date_to_month = as.POSIXct(投稿日, format="%Y-%m-%d") %>% round(units = 'month') %>% as.POSIXct())

# 多い順
trump_tweets.count %>% arrange(-n)
# 可視化
ggplot(trump_tweets.count, aes(x=date_to_month, y=n)) +
  geom_col() +
  xlab('year_to_month') + ylab('tweets per month') +
  ggtitle('# "MAKE AMERICA GREAT AGAIN" per month in tweets of @realDonaldTrump')

f:id:kimigayoseishou:20200315232505p:plain
「MAKE AMERICA GREAT AGAIN」ツイート月間投稿数の過去5年推移

グラフをみてみると, 2015年~2016年にかけツイート数のピークがみられます.
トランプの大統領選出馬立候補日が2015/6/16なので, 選挙活動開始タイミング〜大統領当選までの期間とおおよそ重なります.


変化点検出

ツイート数について2016年付近と比較して2020年付近ではツイート数が減少しているように見受けられます. 以下ツイート数変化の境がいつになるのか、変化点検知を試します.

モデル

\tauを変化点, その前後でt期における投稿数tweets_tがパラメータがそれぞれ\lambda_1,\lambda_2 のPoisson分布で与えられるモデルを立てる。

 tweets_t \sim {\sf Poisson}(t \geq  \tau  \ ? \ \lambda_2 : \lambda_1)

\lambda_1, \lambda_2, \tauには次のような事前分布を課します


\lambda_1 \sim {\sf Exponential}(a),
\lambda_2 \sim {\sf Exponential}(b),
\tau \sim {\sf Uniform}(T).

このモデルの同時分布は次のようになります,

 p(tweets_{1:T}, \lambda_1, \lambda_2, \tau)

 = p(\lambda_1)p(\lambda_2)p(\tau)p(tweets_{1:T} | \lambda_1, \lambda_2, \tau)

 =p(\lambda_1)p(\lambda_2)p(\tau) \prod_{t=1}^T p(tweets_t | \lambda_1, \lambda_2)

 = {\sf Exponential}(\lambda_1 | a){\sf Exponential} (\lambda_2 | b){\sf Uniform}(\tau | T) \prod_{t=1}^T {\sf Poisson}(tweets_t | t  \geq  \tau \ ? \  \lambda_2 : \lambda_1).

Stanでのサンプリング

Stan Users Guideを参考に、 Stanを用いたパラメータ推定を行います。

Poisson分布のパラメータのサンプリング

ハイパーパラメータに離散確率変数を含むモデルでは, Stanを用いて直接サンプリングは行えないため, 離散確率変数に関して周辺化を行う必要がある.

 \tauに関する周辺化を行った上で\lambda_1, \lambda_2からサンプリングを行う.

\lambda_1, \lambda_2からサンプリングを行うため, p(tweets_{1:T}|\lambda_1, \lambda_2)を計算する.

 p(tweets_{1:T} | \lambda_1, \lambda_2)

 = \sum_{t = 1}^T p(tweets_{1:T}, t | \lambda_1, \lambda_2)

 =\sum_{t = 1}^T p(\tau = t)p(tweets_{1:T} | \lambda_1, \lambda_2, t)

 =\sum_{t = 1}^T \exp \biggl( \log \ p(\tau = t) + \sum_{t'=1}^T \log \ {\sf Poisson}(tweets_{t'} | t'  \geq  t \ ? \  \lambda_2 : \lambda_1) \biggr)

stanコードで実装したものが次になります.

data {
  int<lower=0> T_max;
  int<lower=0> tweets[T_max];
  real<lower=0> a;
  real<lower=0> b;
}

transformed data {
  real log_unif;
  log_unif = -log(T_max);
}

parameters {
  real<lower=0> lambda_1;
  real<lower=0> lambda_2;
}

transformed parameters {
  vector[T_max] lp[T_max];// lp[t, s] = log p(tweets[t] | lambda_1, lambda_2, s)
  vector[T_max] lp_marginal;// lp[t] = log p(tweets[t] | lambda_1, lambda_2)
  
  for (s in 1:T_max)
    for (t in 1:T_max)
      lp[t, s] = poisson_lpmf(tweets[t] | t < s ? lambda_1 : lambda_2);
  
  lp_marginal = rep_vector(log_unif, T_max);
  for (s in 1:T_max)
    for (t in 1:T_max)
      lp_marginal[s] = lp_marginal[s] + lp[t, s];
}

model {
  lambda_1 ~ exponential(a);
  lambda_2 ~ exponential(b);
  target+= log_sum_exp(lp_marginal);
}

R側で上のコードをキックしてサンプリングを行います.

.model = rstan::stan_model(file = './stan/change_point.stan')

N <- 4000
tweets <- trump_tweets.count[['n']]
T_max <- length(tweets)

# 各chain 3,000の長さのサンプル(前半) を4つ作成
.fit = rstan::sampling(.model, data = list(T_max=T_max, tweets=tweets, a=0.1, b=0.1), 
                       iter=N, chains=4, warmup=1000)

rstan::stan_trace(.fit, c('lambda_1', 'lambda_2'))

rstan::stan_hist(.fit, pars = c('lambda_1', 'lambda_2'))

.fit %>% 
    rstan::extract(pars = c('lambda_1', 'lambda_2')) %>%
    lapply(summary)

f:id:kimigayoseishou:20200315235713p:plain
事後分布からのサンプリング

中央値では, \lambda_1=16.2, \lambda_2 = 4.5と推定されました. 特に, 変化点の前半では平均して約3.6倍多く「MAKE AMERICA GREAT AGAIN」ツイートが投稿されていたと考えられます.

変化点の事後分布とその可視化

\tauの事後分布はMonte Carlo法を用いて次のように計算できる.

p(\tau | tweets_{1:T})

 \propto p(tweets_{1:T} | \tau)p(\tau)

 \underset{\text{Monte Carlo}}{\sim} p(\tau) \frac{1}{N_{\text{sample}}} \sum_{i=1}^{N_{\text{sample}}} p(tweets_{1:T} | \lambda_1^{(i)}, \lambda_2^{(i)}, \tau)

 = p(\tau) \frac{1}{N_{\text{sample}}}\sum_{i=1}^{N_{\text{sample}}} \exp \biggl( \sum_{t = 1}^T \log \  p(tweets_t | \lambda_1^{(i)}, \lambda_2^{(i)}, \tau) \biggr)

先ほどのサンプリング結果から\tauの事後確率をR側で計算するには次のようにします,

# p(tau | tweets)
tau_post <- rep(0, T_max)

for (tau in c(1:T_max)) {
  tau_post[tau] <- .fit %>% 
    rstan::extract(pars = str_glue('lp[{c(1:T_max)}, {tau}]')) %>%
    as_tibble() %>%
    apply(MARGIN = 1, function(y) exp(sum(y))) %>%
    mean
}

tau_post <- tau_post / sum(tau_post)

tau_post.df <- data.frame(date_to_month = trump_tweets.count$date_to_month, prob = tau_post)

tau_post.df %>%
  ggplot(aes(x=date_to_month, y=prob, label=if_else(prob > 0.1, date_to_month, NULL))) +
  geom_col() +
  geom_text(vjust = 0, nudge_y = 0.05) +
  ylab('posterior probability of tau') +
  ggtitle('posterior probability of tau')

f:id:kimigayoseishou:20200316000450p:plain
変化点の事後確率分布

変化点の事後確率のピークは2016年6月ということがみて取れます. 2016年7月に開かれた共和党全国大会の時期にかけて、 集中的に「MAKE AMERICA GREAT AGAIN」と呟いていたということがわかります.

【Node.js】 Googleスプレッドシートを簡易データベースとして使う

目標

Googleスプレッドシートをデータベースと見立て操作を行うプログラムを作成します.

ソースコードは一応gitに上げているので参考までに.

【Node.js】 Googleスプレッドシートを簡易データベースとして使うGitレポジトリ

Google Sheets APIクライアントアカウントを作成する.

1. プロジェクトの作成

Google開発者向けサービス の上部タブより[プロジェクトの選択] > [新しいプロジェクト]へ進みプロジェクトを作成する.

2. Google Sheets APIの有効化

1で作成したプロジェクト [ダッシュボード]上で[APIとサービスを有効化]をクリックしたあと, Google Sheets APIを選択し, これを有効化.

f:id:kimigayoseishou:20200301221747p:plain
APIとサービスを有効化

f:id:kimigayoseishou:20200301221825p:plain
Google Sheets APIを選択

f:id:kimigayoseishou:20200301221915p:plain
Google Sheets API有効化

3. サービスアカウントの作成

1で作成したプロジェクトの[認証情報]上で[認証情報を作成]>[サービス アカウントの作成]へ進み, サービスアカウントを作成.

f:id:kimigayoseishou:20200301221946p:plain
サービスアカウント作成

4. 認証キーの取得

プロジェクトの[認証情報]画面上で, 3で作成したアカウントの詳細に進み, [キーを作成]ボタンをクリック. 認証キーが作成されるのでjson形式で保存する.

f:id:kimigayoseishou:20200301222009p:plain
認証キー作成


set up

npmでjsプロジェクを作成しましょう.

また今回の開発ではスプレッドシート操作のため, Google Sheets API v4のjs wrappergoogle-spreadsheetを同時にinstallします.

npm init -y
npm i google-spreadsheet

スプレッドシートの情報にアクセスする

サービスアカウントにスプレッドシートの閲覧権限を与えておくことにより, サービスアカウントによる認証を行った上でスプレッドシートの情報にアクセスすることができます. (これに対して, 一般ユーザーに公開されているスプレッドシートへのアクセスにはAPIキー認証で十分です.)

以下は, 既存のスプレッドシートへの編集権限がサービスアカウントに与えられていることを前提とします. まず手始めに, スプレッドシートのキーからタイトルを取得してみましょう.

  • getSpreadsheetTitle.js
const { GoogleSpreadsheet } = require('google-spreadsheet');
// 認証情報jsonファイルを読み込む
const CREDIT = require('<認証情報jsonファイルへのパス>')
// スプレッドシートキー
const SPREADSHEET_KEY = '<spreadsheetのキー>'

const getSpreadsheetTitleByKey = async (spreasheetKey) => {
    // 一般ユーザーに公開していないスプレッドシートへアクセスしたい場合, 作成したサービスアカウントに対し
    // 閲覧権限を与えておく.
    const doc = new GoogleSpreadsheet(spreasheetKey);
    
    // サービスアカウントによる認証
    await doc.useServiceAccountAuth({
        client_email: CREDIT.client_email,
        private_key: CREDIT.private_key,
    });

    // スプレッドシートの情報を読み込みを行い, タイトルを取得
    await doc.loadInfo(); 
    console.log(doc.title);
}

getSpreadsheetTitleByKey(SPREADSHEET_KEY)
// <スプレッドシートのタイトル>

Googleスプレッドシートをデータベースとして使う

以下, Googleスプレッドシートをデータベースとして使うことを想定して次のようなデモを行います.

  1. 空白のスプレッドシートに列名を設定する

  2. スプレッドシートに対する更新・読み取り・更新・削除処理(CRUD処理)

1. 空白のスプレッドシートに列名を設定する

setHeaderRowメソッドはスプレッドシートに1行目に列名を書き込みます.

const { GoogleSpreadsheet } = require('google-spreadsheet');
// 認証情報jsonファイルを読み込む
const CREDIT = require('<認証情報jsonファイルへのパス>')
// スプレッドシートキー
const SPREADSHEET_KEY = '<spreadsheetのキー>'

const setHeaderToSpreadsheet = async (spreasheetKey, sheetIndex, headerValues) => {
    
    const doc = new GoogleSpreadsheet(spreasheetKey);
    
    await doc.useServiceAccountAuth({
        client_email: CREDIT.client_email,
        private_key: CREDIT.private_key,
    });

    // スプレッドシートの情報を読み込みを行い, タイトルを取得
    await doc.loadInfo(); 
    const sheet = doc.sheetsByIndex[sheetIndex]

    // ヘッダー行を作成する
    await sheet.setHeaderRow(headerValues)
}

setHeaderToSpreadsheet(SPREADSHEET_KEY, 0, ['id', 'name', 'age'])

2. 作成・読み取り・更新・削除処理の実装

Googleスプレッドシートをデータベースとして利用するため, リソースに対する更新・読み取り・更新・削除(CRUD処理)を実装したものが以下になります.

  • spreadSheetService.js
class SpreadSheetService {
    /**
     * コンストラクター
     * @param {*} spreadsheetKey スプレッドシートキー
     */
    constructor(spreadsheetKey) {
        this.doc = new GoogleSpreadsheet(spreadsheetKey);
    }
    /**
     * サービスアカウントを用いて認証を行う
     * @param {*} credit 
     */
    async authorize(credit) {
        await this.doc.useServiceAccountAuth({
            client_email: credit.client_email,
            private_key: credit.private_key,
        });
    }
    /**
     * 行データを返す
     * @param {*} index 
     */
    async getRows(index) {
        await this.doc.loadInfo(); 
        const sheet = this.doc.sheetsByIndex[index]
        return sheet.getRows();
    }
    /**
     * 行を追加する
     * @param {*} value 
     */
    async insert(value) {
        await this.doc.loadInfo(); 
        const sheet = this.doc.sheetsByIndex[0]
        return await sheet.addRow(value);
    }
    /**
     * データを取得する
     * @param {*} callBack 
     */
    async select(callBack) {
        const rows = await this.getRows(0)
        const data = []
        for (const row of rows) {
            if (callBack(row)) {
                data.push({id: row.id, name: row.name, age:row.age})
            }
        }
        return data
    }
    /** 
     * idに紐づくユーザーの情報を更新する
    */
    async updateById(id, value) {
        const rows = await this.getRows(0);
        for (const row of rows) {
            if (row.id == id) {
                for (const attr in value) {
                    row[attr] = value[attr]
                    await row.save()
                }
            }
        }
    }
    /**
     * idに紐づくユーザーを削除する
     * @param {*} id 
     */
    async deleteById(id) {
        const rows = await this.getRows(0);
        for (const row of rows) {
            if (row.id == id) {
                await row.delete()
            }
        }
    }
}

module.exports = SpreadSheetService

では実際にスプレッドシートの操作を行ってみます.

リソースの追加

  • demo.js
const SpreadSheetService = require('./spreadSheetService')
// 認証情報jsonファイルを読み込む
const CREDIT = require('<認証情報jsonファイルへのパス>')
// スプレッドシートキー
const SPREADSHEET_KEY = '<spreadsheetのキー>'

// データを4件追加
const insertMany = async () => {
    await spreadSheetService.insert({id:1, name:'John Doe', age:40})
    await spreadSheetService.insert({id:2, name:'Jane Doe', age:30})
    await spreadSheetService.insert({id:3, name:'山田太郎', age:20})
    await spreadSheetService.insert({id:4, name:'山田花子', age:30})
}

insertMany()

スプレッドシートの内容が更新されます.

id name age
1 John Doe 40
2 Jane Doe 30
3 山田太郎 20
4 山田花子 30

リソースの読み取り

  • demo.js
// ageが30であるユーザーの情報を取得
spreadSheetService.select(row => row.age == 30)
.then(data => console.log(data))

出力結果

[{ id: '2', name: 'Jane Doe', age: '30' },
{ id: '4', name: '山田花子', age: '30' }]

リソースの更新

  • demo.js
// id=1のユーザーのnameを「Tom Doe」に更新
spreadSheetService.updateById(1, {name: 'Tom Doe'})

スプレッドシートのidが1のユーザーのnameが「Tom Doe」に更新されます.

id name age
1 Tom Doe 40
2 Jane Doe 30
3 山田太郎 20
4 山田花子 30

リソースの削除

  • demo.js
// id=4のユーザーを削除
spreadSheetService.deleteById(4)

結果をスプレッドシートで確認すると,

id name age
1 Tom Doe 40
2 Jane Doe 30
3 山田太郎 20

と, idが4のユーザーが削除されていることが確認できます.

さいごに

今回はGoogleスプレッドシートを簡易データベースと見立ててスプレッドシート操作を行うプログラムを作成しました. Googleスプレッドシートをデータベースとして用いてアプリを作成することにより, SQLなどの知識がないユーザーでも簡単にデータの取得・レポート作成ができるというメリットがあります. そのため分析を目的としたデータ蓄積用のアプリケーション開発に用途があるのではと考えています.

インタラクティブな表をネイティブjavascriptで実装

javascriptの練習をしようと思い, ブラウザで検索・ソートが行える表を実装しました.

サンプルはこちらのURLで確認できます.

またソースコードgitにおいてます.

準備

まず今回用にjsプロジェクトを作成します.

mkdir interactive-table
cd interactive-table
npm init 

またコンテンツ格納用のフォルダを作成します.

mkdir public
mkdir public/style
mkdir public/js

HTML/CSS

最低限のhtmlとcssを記述

public/index.html

<!DOCTYPE html>
<html>
<head>
<title>intaractibe table</title>
<meta charset="utf-8">
<!-- bootstrap css読み込み -->
<link rel="stylesheet" href="https://stackpath.bootstrapcdn.com/bootstrap/4.3.1/css/bootstrap.min.css" integrity="sha384-ggOyR0iXCbMQv3Xipma34MD+dH/1fQ784/j6cY/iJTQUOhcWr7x9JvoRxT2MZw1T" crossorigin="anonymous">
<link rel="stylesheet" href="style/default.css">
<link rel="shortcut icon" href="">
</head>
<body>
<div class="contents-wrapper">
    <div class="header">
      <p>在庫管理</p></br>
    </div>
    <div class="form-group">
      <label for="item_name">名称</label><button type="button" class="btn btn-primary">検索</button>
      <input type="text" id="item_name" class="form-control">
    </div>
    <div class="table-wrapper">
        <table class="table table-striped">
            <thead class="table-dark">
              <tr class="row-header">
                <th scope="col" data-colName="id">id</th>
                <th scope="col" data-colName="name">名称</th>
                <th scope="col" data-colName="unit_price">定価(単位 : 円)</th>
                <th scope="col" data-colName="stock">在庫</th>
              </tr>
            </thead>
            <tbody class="items">
            </tbody>
          </table>
    </div>
</div>
<!-- bootstrap読み込み -->
<script src="https://code.jquery.com/jquery-3.3.1.slim.min.js" integrity="sha384-q8i/X+965DzO0rT7abK41JStQIAqVgRVzpbzo5smXKp4YfRvH+8abtTE1Pi6jizo" crossorigin="anonymous"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/popper.js/1.14.7/umd/popper.min.js" integrity="sha384-UO2eT0CpHqdSJQ6hJty5KVphtPhzWj9WO1clHTMGa3JDZwrnQq4sF86dIHNDz0W1" crossorigin="anonymous"></script>
<script src="https://stackpath.bootstrapcdn.com/bootstrap/4.3.1/js/bootstrap.min.js" integrity="sha384-JjSmVgyd0p3pXB1rRibZUAYoIIy6OrQ6VrjIEaFf/nJGzIxFDsf4x0xIM+B07jRM" crossorigin="anonymous"></script>
<script src="js/app.js"></script>
</body>
</html>

public/style/default.css

body {
    font-family: Arial, Helvetica, sans-serif;
    font-size: 16px;
}

.header {
    margin-top: 40px;
    font-size: 24px;
}

.contents-wrapper {
    padding-left: 120px;
    padding-right: 120px;
}

button {
    margin: 10px;
}

.form-control {
    width: 300px;
}

静的な表を出力

テーブルの全レコードは配列としてTableModeldataとして管理します. また表示するテーブルのレコードはtableItemsで管理するようにし, TableView側でtableItemsを元にHTML文字列を組み立てを行い描画すれば表の出力ができます.

public/js/app.js

const items = [
    {id:1, name:'apple', unit_price:100, stock:10},
    {id:2, name:'orange', unit_price:40, stock:30},
    {id:3, name:'grape', unit_price:1000, stock:60},
    {id:4, name:'banana', unit_price:300, stock:100},
]

class TableModel {
    constructor(data) {
        this.data = data;
        this.tableItems = data;
        this._listeneres = {'search':[], 'sort':[]};
    }
    /**
     * イベントの登録を行う
     * @param {*} type イベント
     * @param {*} callback イベントハンドラ
     */
    on(type, callback) {
        this._listeneres[type].push(callback);
    }
    /**
     * 登録イベントを実行
     * @param {*} type イベント
     */
    trigger(type) {
        this._listeneres[type].forEach(callback => {
            callback();
        });
    }
}

class TableView {
    constructor(data) {
        this.tableItems = data;
        this.model = new TableModel(data);
        this.tableElement = document.querySelector('.items');
        this.init()
    }
    createTableItemsHTML(tableItems) {
        return(tableItems.
            map((e) => `<tr><th scope="row">${e.id}</th><td>${e.name}</td><td>${e.unit_price}</td><td>${e.stock}</td></tr>`).
            join('\n'))
    }
    init() {
        this.tableElement.innerHTML = this.createTableItemsHTML(this.tableItems);
    }
    
}

const app = new TableView(items);

表が正しく出力されると, つぎのようになります.

f:id:kimigayoseishou:20200216224055p:plain


検索ボタン押下時の動作を実装する

検索ボタンを押したときに指定した名称のレコードのみ表示する動作を実装します. この動作のフローは以下になります.

  1. 検索ボタン押下
  2. TableModelのtableItemsを更新
  3. 更新後のtableItemsに基づいてHTMLの表を再描画

TableModelの修正

名前をキーにしてtableItemsを抽出するTableModelのメソッドselectByNameを実装します.

public/js/app.js

    /**
     * tableItemsに対しnameをキーにして抽出を行う
     * @param {*} name 名称
     */
    selectByName(name) {
        this.tableItems = this.data.filter((e) => e.name===name);
        this.trigger('search');
    }

TableViewの修正

TableViewのinitメソッド内部に以下を実装します

  1. search時にHTMLの表を再描画するメソッドをイベントリスナ-に追加
  2. 検索ボタン押下時にTableModel.selectByNameを呼び出すイベントの登録

public/js/app.js

class TableView {
    constructor(data) {
        this.tableItems = data;
        this.model = new TableModel(data);
        this.tableElement = document.querySelector('.items'); 
        this.searchButton = document.getElementsByTagName('button')[0];
        this.init()
    }
    createTableItemsHTML(tableItems) {
        return(tableItems.
            map((e) => `<tr><th scope="row">${e.id}</th><td>${e.name}</td><td>${e.unit_price}</td><td>${e.stock}</td></tr>`).
            join('\n'))
    }
    init() {
        this.tableElement.innerHTML = this.createTableItemsHTML(this.tableItems);
        
        // search時に対して表を再描画するメソッド
        this.model.on('search', () => {
            this.tableElement.innerHTML = this.createTableItemsHTML(this.model.tableItems);
        })

        // 検索ボタン押下時にmodelのtableItemsを更新する
        this.searchButton.addEventListener('click', (e) => {
            this.model.selectByName(document.getElementById('item_name').value);
        })
    }
    
}

テキストボックスに「apple」と入力して、検索ボタンを押すと1件だけ出力されます.

f:id:kimigayoseishou:20200216223947p:plain


ソート動作を実装する.

ソート動作の実装は検索機能実装とほとんど同じ流れなので、 修正箇所だけ列挙します.

TableModelの修正

列名に関してtableItemsをソートするTableModelのメソッドsortByを実装します.

    /**
     * 引数で指定した列名に関してtableItemsをsortを行う
     * @param {*} columnName 列名
     */
    sortBy(columnName) {
        this.tableItems = this.data.sort((e1, e2) => e1[columnName] - e2[columnName]);
        this.trigger('sort');
    }

TableViewを修正する

TableViewのinitメソッド内部に以下を実装します

  1. search時にHTMLの表を再描画するメソッドをイベントリスナ-に追加
  2. 列名押下時にTableModel.sortByを呼び出すイベントの登録
class TableView {
    constructor(data) {
        this.tableItems = data;
        this.model = new TableModel(data);
        this.tableElement = document.querySelector('.items');
        this.searchButton = document.getElementsByTagName('button')[0];
        this.tableColumnElements = document.querySelectorAll('.row-header > th');
        this.init()
    }
    createTableItemsHTML(tableItems) {
        return(tableItems.
            map((e) => `<tr><th scope="row">${e.id}</th><td>${e.name}</td><td>${e.unit_price}</td><td>${e.stock}</td></tr>`).
            join('\n'))
    }
    init() {
        this.tableElement.innerHTML = this.createTableItemsHTML(this.tableItems);
        
        // search時に対して再描画するメソッド
        this.model.on('search', () => {
            this.tableElement.innerHTML = this.createTableItemsHTML(this.model.tableItems);
        })

        // 検索ボタン押下時にmodelのtableItemsを更新する
        this.searchButton.addEventListener('click', (e) => {
            this.model.selectByName(document.getElementById('item_name').value);
        })

        // sort時にmodelのtableItemsを更新する
        this.model.on('sort', () => {
            this.tableElement.innerHTML = this.createTableItemsHTML(this.model.tableItems);
        })

        // 列名押下時にmodelのtableItemsを更新する
        this.tableColumnElements.forEach((colElment) => {
            colElment.addEventListener('click', (e) => {
                this.model.sortBy(e.currentTarget.getAttribute('data-colName'));
            })
        })
    }
}

f:id:kimigayoseishou:20200216224031p:plain