シミュレーションで理解するランダム化比較試験:(2)処置の効果を正しく推定するために

はじめに

前回は、異なる処置をランダムに割り当てれば、処置のグループ間で対象者の属性が偏らずにバランスするというお話をしました。
mstour.hatenablog.com
具体的には、内容の異なるクーポンAとクーポンBのどちらがより売り上げの増加につながるかを検討するという事例を考えました。ランダム化をしない場合・した場合それぞれでクーポンAを配布されたグループ(以下、グループAとします)とクーポンBを配布されたグループ(同じくグループBとします)の属性を集計し、ランダム化をしない場合には属性がグループ間でアンバランスになりうること、対してランダム化によってうまく属性のバランスが取れそうだということを確認しました。

今回は、ユーザーのある属性がクーポン配布後の購入金額に強く影響を及ぼすようなケースを考え、ランダム化を行わないと処置(クーポン)の効果を誤って推定してしまう恐れがあることをシミュレーションで見ていきたいと思います。


設定

前回と同じく、異なる内容のクーポンAとクーポンBのどちらがより購入金額を増加させるかを調べたいという架空の事例で説明していきます*1。状況設定を再掲すると以下の通りです。

* * *
とあるアパレルショップでは、新規にアプリ会員登録をしたお客さんを対象に、次回の来店時に使用できるクーポンを配布することにしました。ただし、クーポンの内容をどうするかについては社内で議論があり、以下の2つのどちらが次回の購入金額アップにつながるかを確かめたいということになりました。

(クーポンA):5,000円以上の買い物をした場合に1,000円を割引
(クーポンB):10,000円以上の買い物をした場合に購入金額の10%を割引

さらに、一度の購入単価はお客さんの属性(年齢、性別、収入、など・・・)によってかなり違いそうだという意見が出たため、どちらのクーポンを配布するかをランダムに決めることでクーポンA/Bそれぞれの配布されるお客さんの属性が偏らないよう配慮し、平均的にどちらのクーポンを配布したグループの売り上げが高くなるかを評価しよう、と会議で決定しました。
* * *

今回は、グループAとグループBの次回の購入金額を比較するという状況を考えます。

さて、もう少し詳しく設定を作っていきましょう。
まず、購入金額は女性の方が男性よりも平均で5,000円ほど大きい傾向にあるものとします。したがって、今回のシミュレーションでも購入金額のデータ生成にはそのことを反映させます。また、話を単純にするために、それ以外の属性は全く購入金額には影響していないものとします。
もう一つ重要な点として、クーポンAとクーポンBには効果の違いはないものとします。したがって、正しく推定を行うことができれば、グループAとグループBの購入金額の平均は同じくらいになるはずです。

そして、前回と同じく「ランダム化を行わないケース」「ランダム化を行うケース」それぞれで結果がどうなるかを見ていくことにします。
割り当て方法1として前半の500人にクーポンA、後半の500人にクーポンBを配布するという安直な方法を、割り当て方法2としてランダムにどちらのクーポンを配布するかを決める方法を設定します。

今回のポイントは、「処置(A/B異なるクーポンの配布)はアウトカム(購入金額)には影響を与えないが」、「処置以外の因子(性別)はアウトカム(購入金額)に影響を与える」という点です。
したがって、もしクーポンの割り当てと性別とに関連性があった場合には「交絡」が生じ、クーポンの違いと購入金額との間に見せかけの関連性が現れてしまいます。

f:id:mstour:20210411211416j:plain

この点を、ランダム化をしない場合(割り当て方法1)とランダム化した場合(割り当て方法2)との対比で確認していきましょう。


シミュレーション

前回と同じ要領で乱数を利用して模擬データを作っていきますが、今回は「同じ模擬データ生成プロセスを1000回繰り返す」ということを行います。
つまり、もしもクーポン配布の実験を同じ条件で1000回繰り返したとするとどうなるのか、を見ていきます。1回だけの結果は偶然に左右される恐れがありますが、大量の反復結果を平均することで真の状況に近い結果を得ることが期待できます。


まずはシミュレーションのコード全体を掲載します:

# 使用するパッケージ
library(tidyverse)
library(TruncatedNormal)
library(gt)

# 設定
n = 1000
set.seed(1234)
nsim = 1000

# 各回の結果を格納する入れ物
result1 = NULL
result2 = NULL
test = NULL

# シミュレーション
for (i in 1:nsim){
  # 割り当て方法1
  alloc1 = c(rep("A", n/2), rep("B", n/2))
  # 割り当て方法2
  alloc2 = sample(alloc1)
  # 性別
  gender = c(rbinom(n/2, 1, 0.7), rbinom(n/2, 1, 0.3))
  # 購入金額
  buy = round(rtnorm(n = n, mu = 15000, sd = 5000, lb = 3000, ub = 50000))
  buy[gender == 1] = buy[gender == 1] + 5000 # 女性の場合平均5,000円高い
  # データフレーム
  simdata = data.frame(alloc1 = factor(alloc1),
                       alloc2 = factor(alloc2),
                       gender = gender,
                       buy = buy) %>%
    mutate(genderc = factor(case_when(
      gender == 1 ~ "F",
      gender == 0 ~ "M"
    )))
  # 割り当て方法1の集計
  summary1 = group_by(simdata, alloc1) %>%
    summarize(buy_mean = mean(buy), female_prop = mean(gender)) %>%
    mutate(simno = i)
  # 割り当て方法2の集計
  summary2 = group_by(simdata, alloc2) %>%
    summarize(buy_mean = mean(buy), female_prop = mean(gender)) %>%
    mutate(simno = i)
  # t検定と95%信頼区間
  test1 = t.test(buy ~ alloc1, data = simdata, var.equal = TRUE)
  test2 = t.test(buy ~ alloc2, data = simdata, var.equal = TRUE)
  testres = data.frame(allocation = c(1, 2),
                       pvalue = c(test1$p.value, test2$p.value),
                       CIlower = c(test1$conf.int[1], test2$conf.int[1]),
                       CIupper = c(test1$conf.int[2], test2$conf.int[2]))
  # 結果を格納
  result1 = rbind(result1, summary1)
  result2 = rbind(result2, summary2)
  test = rbind(test, testres)
}

ざっくりと、中身を部分に分けて説明していきます。

library(tidyverse)
library(TruncatedNormal)
library(gt)

n = 1000
set.seed(1234)
nsim = 1000

今回使用するパッケージの読み込みと、シミュレーションの設定です。クーポンを配布するユーザーは合計1000人、シミュレーションの反復回数は1000回にしています。また乱数のシード値を指定することで同じ結果が再現できるようにしています。

result1 = NULL
result2 = NULL
test = NULL

これらにシミュレーションの各回の結果を格納していきます。

for (i in 1:nsim){

}

シミュレーションのためのfor文(反復処理)です。Rでfor文は処理速度がとても遅いらしく推奨されていないようなのですが、直感的に理解しやすいので今回はこれで実施しています*2
以下は、シミュレーションの各ステップの内容(for文の中身)です。

alloc1 = c(rep("A", n/2), rep("B", n/2))
alloc2 = sample(alloc1)
gender = c(rbinom(n/2, 1, 0.7), rbinom(n/2, 1, 0.3))
buy = round(rtnorm(n = n, mu = 15000, sd = 5000, lb = 3000, ub = 50000))
buy[gender == 1] = buy[gender == 1] + 5000

割り当て方法は、「設定」のところで記述した通りです。性別は前回と同じ方法でサンプリングします。つまり、前半の500人には女性が、後半の500人には男性が多くなります。
クーポン配布後の購入金額については最小3000円・最大50000円となるように切断正規分布からサンプリングしますが、女性のほうが平均的に5000円購入金額が高いという設定を行いました。
さらに、購入金額は配布されたクーポンがAかBかには関係なくサンプリングされています。クーポンAとBとの間には差がないという状況にするため、このように設定しています。

続いて上記で作成した変数をデータフレームの形にまとめます。

simdata = data.frame(alloc1 = factor(alloc1),
                       alloc2 = factor(alloc2),
                       gender = gender,
                       buy = buy) %>%
mutate(genderc = factor(case_when(
      gender == 1 ~ "F",
      gender == 0 ~ "M"
)))


以下のコードでは、割り当て方法1の場合・2の場合それぞれで、グループA/Bそれぞれの女性の割合と購入金額の平均値を集計します。
また、各割り当て方法でグループA/B間の購入金額の差に関するt検定のp値と、95%信頼区間を計算します*3
最後にこれらの結果を格納して、次の回へ進みます。

  # 割り当て方法1の集計
  summary1 = group_by(simdata, alloc1) %>%
    summarize(buy_mean = mean(buy), female_prop = mean(gender)) %>%
    mutate(simno = i)
  # 割り当て方法2の集計
  summary2 = group_by(simdata, alloc2) %>%
    summarize(buy_mean = mean(buy), female_prop = mean(gender)) %>%
    mutate(simno = i)
  # t検定と95%信頼区間
  test1 = t.test(buy ~ alloc1, data = simdata, var.equal = TRUE)
  test2 = t.test(buy ~ alloc2, data = simdata, var.equal = TRUE)
  testres = data.frame(allocation = c(1, 2),
                       pvalue = c(test1$p.value, test2$p.value),
                       CIlower = c(test1$conf.int[1], test2$conf.int[1]),
                       CIupper = c(test1$conf.int[2], test2$conf.int[2]))
  # 結果を格納
  result1 = rbind(result1, summary1)
  result2 = rbind(result2, summary2)
  test = rbind(test, testres)

結果

結果の出力にはgtパッケージを使ってみることにしました。
例えば以下のページが非常に参考になりました。
onoshima.github.io


まず、割り当て方法1を行った場合の集計結果は次のようになりました。

f:id:mstour:20210411211441j:plain

一番右の列「Female(Percent)」は、女性の割合の1000回のシミュレーション結果を平均したものです*4
グループAには女性が、グループBには男性が多く含まれていることがわかります(そのようにデータを作っているのですが)。
真ん中の列「Buy(Yen)」が、クーポン配布後の購入金額の全シミュレーションを通した平均を示しています。グループAの方がグループBよりも2000円ほど高くなっています。
本当はクーポンAとクーポンBとで購入金額に差はないのですが、テキトウにクーポンの配布を行うと間違ってクーポンAの方が売り上げアップにつながると結論してしまう可能性があることがわかります。


次に、割り当て方法2(ランダム化)の場合の集計結果は以下のようになりました。

f:id:mstour:20210411211453j:plain

同じく1000回のシミュレーションを通した平均を示していますが、グループAとBの間には違いはほとんどありません。
もちろん、偶然グループAとBの間に差が生じる可能性はありますが、ランダムにクーポンを配布すれば、クーポンの効果を間違って結論づけてしまう危険はかなり避けられると言えるでしょう。
この点は、t検定のp値と95%信頼区間の結果を見るとより理解できます。

f:id:mstour:20210411211505j:plain

左の列「Method」は、割り当て方法1(ランダム化なし)と2(ランダム化)を示しています。

真ん中の列「Type1 error rate」には、1000回のシミュレーションのうち、t検定のp値が0.05未満になった割合を計算したものです。
割り当て方法1のType1 error rateは100%であり、これは1000回中1000回でp<0.05となったことを意味します。本当はグループ間に差はないはずなので、有意水準5%でt検定を行った場合に確実に間違って有意差ありと判定してしまいます*5
一方、割り当て方法2のType1 error rateは5.5%であり、有意水準5%に近い値となっています。この結果は、有意水準を5%とした場合の「本当は差がないのに、間違って差があると結論づける確率が5%」という検定の前提条件に整合しており、割り当て方法1のような重大な問題が起こっていないことがわかります。

右の列「Coverage probability」は、シミュレーションで計算した1000個の95%信頼区間のうち、本当のグループ間の差である「0」を含んでいるものの割合を算出したものです(日本語では、「被覆確率」と呼ばれます)。
割り当て方法1では、0を含む信頼区間は一つもありませんでした。正しく推定が行われていたならば、計算した信頼区間のうち95%は0を含むことが期待されるので、やはり割り当て方法1には欠陥があることがわかります。
割り当て方法2では、Coverage probabilityは94.5%であり、期待される95%にほぼ等しい結果になりました。

集計に使用したRのコードは以下の通りです。長くなりますので、記法の詳細については前述の参考資料などをご覧いただけると幸いです。

simall_1 %>%
  dplyr::select(-female_mean_all) %>%
  gt() %>%
  fmt_number(columns = vars("buy_mean_all", "female_pct"), decimals = 1) %>%
  cols_label(alloc1 = "Group", buy_mean_all = "Buy(Yen)", female_pct = "Female(Percent)")
simall_2 %>%
  dplyr::select(-female_mean_all) %>%
  gt() %>%
  fmt_number(columns = vars("buy_mean_all", "female_pct"), decimals = 1) %>%
  cols_label(alloc2 = "Group", buy_mean_all = "Buy(Yen)", female_pct = "Female(Percent)")
testall_sum %>%
  gt() %>%
  cols_label(allocation = "Method", alpha_error = "Type1 error rate", coverage = "Coverage probability")

まとめ

関心のある結果(購入金額)に影響を与える因子(性別)によって処置(異なるクーポンの配布)に見せかけの効果が現れるような状況を想定し、ランダム化をしない場合・した場合それぞれでどのようなことが起こるかをシミュレーションで見ていきました。
ランダム化を行うことにより、結果に影響する因子と処置との間の関係性を取り除くことができるため、見せかけの処置効果は現れず設定通りグループ間に差はないという結果となりました。一方、ランダム化をしなかった場合には、見事に見せかけの処置効果が現れました。

このようにランダム化は処置の効果を正しく推定するための非常に強力な方法なのですが、それでも、偶然によって属性のアンバランスが生じることはありえます。例えば今回の例でも、シミュレーションのうち何回かは性別に無視できないほどの偏りが生じた可能性があります。
特に結果に大きく影響するような属性に関しては、グループ間でバランスが取れるように祈るだけでなく、もっと積極的にバランスを取りにいった方が望ましい場合もあります。そのような方法としては、例えば層別ランダム化が広く利用されています。
次回以降は、そのような方法をまたシミュレーションでやっていこうかと考えています。

今回の記事作成では以下の資料を参考にしました。

[1] 岩崎学(2015), "統計的因果推論", 朝倉書店.
[2] J.L.Fleiss(2004), "臨床試験のデザインと解析", アーム.
[3] gtパッケージについて - TAKAHIRO ONOSHIMA'S WEBSITE
https://onoshima.github.io/stat/gt/

*1:どちらのクーポンも何もしない時と大して結果は変わらない、という可能性は今回は無視することにします。が、このようなケースで一つネタが作れそうなので、今後また出てくるかもしれません。

*2:というよりは、それ以外にシミュレーションをやるための方法が思いつきませんでした。

*3:以降の結果は、p値、信頼区間ともに両側で算出しています。

*4:シミュレーションを通した平均を計算しているので、グループAとBの合計は1にはなっていません。

*5:このような結果になるのは、サンプルサイズが1000とかなり大きいため、些細な差も統計的有意となってしまうからかもしれません。しかしながら、先に述べた単純集計の結果からもこの割り当て方法がまずいことは明らかです。