傾向スコア法はサンプルサイズが小さくても大丈夫なのか：IPTWの例

【最終更新　2022/5/13】

シミュレーションでのバイアス評価で、推定したい真値を周辺オッズ比とし、モンテカルロ法で算出するようにしました。
条件を変えたシミュレーションをfurrrパッケージで行うようにコード例を修正しました。
その他、細かな表現を改めた箇所があります。

はじめに

観察データを用いて処置の効果を評価したい場合、処置の割り当てがランダムに行われていないことが原因で、単純にアウトカム変数を比較したのでは正しい結果が得られないことがあります。
例えば、セールのお知らせを送る場合と送らない場合とでセール期間中に来店する会員の割合がどれだけ違うかを調べたいとすると、お知らせを送る相手がランダムに決定されていないのであれば、割合の差を単純計算するとミスリーディングな結果を導く可能性があります。

このような場合に有効な解析方法としてよく知られているのが傾向スコア法(Propensity score method)です。*1

ところで、実務上、数十〜数百ほどのサンプルサイズのデータで処置効果の評価をしたい場合がしばしばあります。
ロジスティック回帰などの統計モデルでは、サンプルサイズ（正確にはイベント数）に対する説明変数の数の目安がよく知られており、サンプルサイズが小さい場合には気をつけないといけないことがあります。*2
では、傾向スコア法についてはサンプルサイズがどれくらいあれば大丈夫なのでしょうか？

この点を調べた研究はそれほど多くないようなのですが、例えば、

Pirracchioらによる報告：
Evaluation of the Propensity score methods for estimating marginal odds ratios in case of small sample size | BMC Medical Research Methodology | Full Text
特にサンプルサイズが小さい場合の、傾向スコアマッチング(Propensity score matching、以下マッチングとします)と逆確率重みづけ(Inverse probability of treatment weighting、以下IPTWとします)の性能を調べている

Bottigliengoらによる報告：
Oversampling and replacement strategies in propensity score matching: a critical review focused on small sample size in clinical settings | BMC Medical Research Methodology | Full Text
同じく、サンプルサイズが小さい場合に、1対1マッチングと1対多マッチング、復元マッチングと非復元マッチング*3との性能の違いを調べている

などがあります。

今回はPirracchioらの報告（以下、Pirracchio論文と書きます）について紹介するとともに、Rを使って再現シミュレーションを行なってみます。

サンプルサイズがかなり小さくても問題はないとのことだが...

Pirracchio論文では、処置群と対照群の2グループ間での比較を行う状況を想定し、グループ間での背景因子の偏りをマッチングやIPTWで補正することで2値アウトカムに関するオッズ比（対照群のオッズに対する、処置群のオッズの比）を正しく推定できるかどうかをシミュレーション実験により検討しています。
検討された各グループのサンプルサイズは最大で1000、最小で40でした。*4

全体的な結論として、傾向スコアのモデルが適切に設定されていれば、サンプルサイズが小さい場合でも（検討された範囲では、各グループ40でも）処置効果を正しく推定できるだろう、とされています。
具体的には、さまざまな設定のもとでも、処置効果の推定のバイアスはせいぜい10%以内であり、また第一種の過誤(Type 1 error)の生じる確率も名目値である5%を大きく超えることはない*5という結果でした。

シミュレーション

設定

Pirracchio論文ではデータの分布についてさまざまなシナリオを検討していますが、今回はTable1の一番左の列に記載の「共変量の効果は中程度（条件付きオッズ比1.5）、処置効果は強い（条件付きオッズ比2.5）」というシナリオのみ再現してみます。*6
設定について、Rでのコード例とともに説明します。

＜背景因子（共変量）＞
$X_1, X_2, X_3, X_4$ の4つがあるとします。それぞれ、互いに独立に平均 $0$ 、分散 $0.5^2$ の正規分布に従います。
例えば以下のようにMASSパッケージのmvrnormでサンプリングができます。なおここでは各グループのサンプルサイズを100としています。

library(MASS)

n <- 100
X <- mvrnorm(n = 2*n, mu = rep(0, 4), Sigma = diag(rep(0.5^2, 4)))

＜処置の割り当て＞
$Z$ は処置群の場合1を、対照群の場合0をとる変数です。以下のロジスティックモデルを通して、背景因子が割り当て確率に関連しているものとします。つまり、背景因子のうち $X_3, X_4$ が処置の割り当てに影響する状況を考えます。処置の割り当て比率が平均的に1:1になるよう、切片は0に設定されています。

$\log \left ( \frac{P(Z=1|X)}{1-P(Z=1|X)} \right ) = 0.41 X_3 + 0.41 X_4$

以下では、まずロジスティックモデルのパラメータalphaを定義し、上の式の右辺を計算しています。次にその結果を確率 $P(Z=1|X)$ に変換し、二項分布から $Z$ をサンプリングします。

alpha_0 <- 0
alpha_1to4 <- c(0, 0, 0.41, 0.41)

logit_p_trt <- alpha_0 + X %*% alpha_1to4
p_trt <- 1 / (1 + exp(-logit_p_trt))

Z <- rbinom(n = 2*n, size = 1, prob = p_trt)

＜アウトカム＞
$Y$ は関心ある2値アウトカムを表し、処置および共変量との関連はロジスティックモデルで表されるとします。 $X_2$ はアウトカムのみに影響し、 $X_4$ は割り当てとアウトカムの双方に影響する、いわゆる交絡因子となる状況です。
また $Z$ の回帰係数0.92は、処置群の対照群に対する（条件付き）オッズ比が約2.5であることを意味しています。

$\log \left ( \frac{P(Y=1|Z, X)}{1-P(Y=1|Z, X)} \right ) = \beta_0 + 0.92 Z + 0.41 X_2 + 0.41 X_4$

なお $\beta_0$ は、1,000,000サンプルを用いて、アウトカムの発現割合が平均で50%になるように最適化計算で求めることとされています。計算方法は色々あるかと思いますが、ここでは最小二乗法と同じ考え方で以下の誤差、つまり「ロジスティックモデルによる予測確率と50%との誤差（=差の二乗和）」が最小になるような値を $\beta_0$ としました。

$\sum_{i=1}^{1000000} \left( \frac{1}{1 + \exp( - (\beta_0 + 0.92 Z_i + 0.41 X_{2i} + 0.41 X_{4i}) )} - 0.5 \right)^2$

beta_1to4 <- c(0, 0.41, 0, 0.41)
beta_T <- 0.92

# beta_0の決定
## 1,000,000サンプルの生成
n_p <- 1000000
X_p <- mvrnorm(n = n_p, mu = rep(0, 4), Sigma = diag(rep(0.5^2, 4)))
logit_p_trt_p <- alpha_0 + X_p %*% alpha_1to4
p_trt_p <- 1 / (1 + exp(-logit_p_trt_p))
Z_p <- rbinom(n = n_p, size = 1, prob = p_trt_p)

## 誤差関数
targetfunc <- function(b) {
  sum((1/(1+exp(-(b + Z_p*beta_T + X_p%*%beta_1to4))) - rep(0.5, n_p))^2)
}
## 誤差の最小化
opt <- optimize(targetfunc, interval = c(-5, 5))
beta_0 <- opt$minimum

# Yのサンプリング
logit_p_out <- rep(beta_0, 2*n) + beta_T * Z + X %*% beta_1to4
p_out <- 1 / (1 + exp(-logit_p_out))
Y <- rbinom(n = 2*n, size = 1, prob = p_out)

# データフレームに整理
simdf <- data.frame(ID = 1:(2*n), Z = Z, Y = Y,
                    X_1 = X[,1], X_2 = X[,2], X_3 = X[,3], X_4 = X[,4])

このようにして生成されたデータは、共変量が処置群・対照群の間でバランスしておらず、処置効果について交絡が生じています。そのため、傾向スコア法によって共変量のアンバランスを補正した解析を行います。
Pirracchio論文ではマッチングとIPTWの両方が検討されていますが、今回の再現ではIPTWのみを扱います。

今回は各グループのサンプルサイズが40、60、100、500の場合のIPTWの推定性能を確認しています。

ここで、IPTWで推定する対象は処置に関する周辺オッズ比となるのですが、周辺オッズ比は条件付きオッズ比のように回帰係数を指数変換するだけでは求められず、共変量についての積分計算が必要になります。
今回は切片 $\beta_0$ を計算する際に使用した1,000,000サンプルを利用して、モンテカルロ法に基づく周辺オッズ比を算出しました（参考文献[2]や[4]を参照）。この値を各シミュレーションにおける真値として参照します。

## 1,000,000サンプル個々のアウトカム発現確率
logit_p_out_e <- rep(beta_0, n_p) + beta_T * Z_p + X_p %*% beta_1to4
p_out_e <- 1 / (1 + exp(-logit_p_out_e))

## 群ごとにアウトカム発現確率を平均
p_out_mean <- data.frame(p_out = p_out_e, Z = Z_p) %>%
group_by(Z) %>% summarise(p = mean(p_out))

## 周辺オッズ比
true_marginal_OR <- (p_out_mean[2, ]$p/(1-p_out_mean[2, ]$p))/(p_out_mean[1, ]$p/(1-p_out_mean[1, ]$p))

＜IPTW＞
ロジスティック回帰モデルなどによって推定した、処置に関する個人 $i$ の傾向スコアを $p_i$ とするとき、各個人に以下の重みをつけて解析する方法がIPTWです。

$IPTW_i = \frac{Z_i}{p_i} + \frac{1-Z_i}{1-p_i}$

つまり処置群に属する個人には $\frac{1}{p_i}$ の重みが、対照群に属する個人には $\frac{1}{1-p_i}$ の重みが与えられます*7。

傾向スコアを推定するためのモデルが正しく設定されているならば（交絡因子がすべて含まれているなど）、IPTWによって処置効果を正しく（ばらつきはあるが、バイアスなく）推定することができます。

Rで傾向スコアを用いたIPTWを実施する方法はいくつかありますが、今回は重みの計算にipwパッケージを使用し、重みづけ解析の部分はPirracchio論文と同様にsurveyパッケージで行いました。*8
ipwパッケージとsurveyパッケージを用いたIPTWの実行については以下の論文でも解説されているので、詳細はこちらをご参照ください。
ipw: An R Package for Inverse Probability Weighting | Journal of Statistical Software

以下のコード例では、最初にipwパッケージのipwpoint関数で重みの計算を行います。重みにもさまざまな種類がありますが、ここではPirracchio論文と同様の基本的なものを使用しました。
重みは元のデータフレームへ格納します。
次にsurveyパッケージのsvyglmを使用し、IPTWによるロジスティック回帰を行います。
最後に、回帰の出力結果をbroomパッケージのtidy関数を利用してデータフレームにまとめます。

library(ipw)
iptw <- ipwpoint(exposure = Z, family = "binomial", link = "logit",
                 numerator = ~ 1, denominator = ~ X_1 + X_2 + X_3 + X_4,
                 data = simdf, trunc = NULL)
simdf_w <- simdf %>% mutate(IPTW = iptw$ipw.weights)

library(survey)
fit <- svyglm(Y ~ Z, design = svydesign(id = ~ 1, weights = ~ IPTW, data = simdf_w),
              family = quasibinomial())

library(broom)
res <- fit %>% tidy() %>% filter(term == "Z")

次項でシミュレーション全体のコードと結果を報告します。

シミュレーション結果

ここまでの流れをシミュレーション用にまとめると、以下のようになります。

library(tidyverse)
library(broom)
library(MASS)
library(ipw)
library(survey)

# シミュレーション実行用関数
## nsim：シミュレーション反復数
## n：各グループのサンプルサイズ
## beta_T：処置に関する回帰係数の真値
sim <- function(nsim, n, beta_T){

  # ロジスティックモデルのパラメータ
  alpha_0 <- 0
  alpha_1to4 <- c(0, 0, 0.41, 0.41)
  beta_1to4 <- c(0, 0.41, 0, 0.41)
  beta_T <- beta_T
  
  # 最適化計算でbeta_0を決めます
  n_p <- 1000000
  X_p <- mvrnorm(n = n_p, mu = rep(0, 4), Sigma = diag(rep(0.5^2, 4)))
  logit_p_trt_p <- alpha_0 + X_p %*% alpha_1to4
  p_trt_p <- 1 / (1 + exp(-logit_p_trt_p))
  Z_p <- rbinom(n = n_p, size = 1, prob = p_trt_p)

  targetfunc <- function(b) {
    sum((1/(1+exp(-(b + Z_p*beta_T + X_p%*%beta_1to4))) - rep(0.5, n_p))^2)
  }
  opt <- optimize(targetfunc, interval = c(-5, 5))
  beta_0 <- opt$minimum

  # 周辺オッズ比の算出
  logit_p_out_e <- rep(beta_0, n_p) + beta_T * Z_p + X_p %*% beta_1to4
  p_out_e <- 1 / (1 + exp(-logit_p_out_e))
  p_out_mean <- data.frame(p_out = p_out_e, Z = Z_p) %>%
    group_by(Z) %>% summarise(p = mean(p_out))
  true_marginal_OR <- (p_out_mean[2, ]$p/(1-p_out_mean[2, ]$p))/(p_out_mean[1, ]$p/(1-p_out_mean[1, ]$p))

  # 結果を格納するデータフレームを作っておきます
  sim_result <- data.frame(n = rep(n, nsim),
                           beta_T = rep(beta_T, nsim),
                           true_or = rep(true_marginal_OR, nsim),
                           estimate = rep(0, nsim),
                           p.value = rep(0, nsim))
  
  # シミュレーション反復開始です
  for (i in 1:nsim){
    ## データを生成します
    X <- mvrnorm(n = 2*n, mu = rep(0, 4), Sigma = diag(rep(0.5^2, 4)))

    logit_p_trt <- alpha_0 + X %*% alpha_1to4
    p_trt <- 1 / (1 + exp(-logit_p_trt))
    Z <- rbinom(n = 2*n, size = 1, prob = p_trt)

    logit_p_out <- rep(beta_0, 2*n) + beta_T * Z + X %*% beta_1to4
    p_out <- 1 / (1 + exp(-logit_p_out))
    Y <- rbinom(n = 2*n, size = 1, prob = p_out)

    simdf <- data.frame(ID = 1:(2*n), Z = Z, Y = Y, 
                        X_1 = X[,1], X_2 = X[,2], X_3 = X[,3], X_4 = X[,4])

    ## IPTWで解析します
    iptw <- ipwpoint(exposure = Z, family = "binomial", link = "logit",
                     numerator = ~ 1, denominator = ~ X_1 + X_2 + X_3 + X_4,
                     data = simdf, trunc = NULL)
    simdf_w <- simdf %>% mutate(IPTW = iptw$ipw.weights)

    fit <- svyglm(Y ~ Z, design = svydesign(id = ~ 1, weights = ~ IPTW, data = simdf_w),
                    family = quasibinomial())
    res <- fit %>% tidy() %>% filter(term == "Z")

    ## 結果を格納します
    sim_result[i,]$estimate <- res$estimate
    sim_result[i,]$p.value <- res$p.value
  }
  return(sim_result)
}

IPTWの性能評価の指標として、Pirracchio論文で検討されているもののうち、以下の項目を確認することにします。

Type1 errorの確率：処置効果がないのに、統計的有意に効果ありと判断してしまう確率。今回は両側5%に設定しているものとします。
Bias：処置効果の推定値と、真の処置効果とに平均的にどれだけズレがあるかを表す。報告では数式は記載されていないが、 $Bias = \hat{OR}_{est} - OR_{true}$ により計算していると思われます。 $\hat{OR}_{est}$ は処置効果（今回はオッズ比で示しています）の推定値のシミュレーションを通した平均、 $OR_{true}$ は真の処置効果で、モンテカルロ法で算出した結果おおむね2.56程度になりました。
Percent bias：バイアスをパーセント表示したものです。こちらも明記はありませんが、 $PercentBias = \frac{\hat{OR}_{est} - OR_{true}}{OR_{true}} \times 100$ で計算していると考えられます。

シミュレーション結果を出力するために、今回は以下のコードを用いました。結果を表形式で出力するためにgtパッケージを使用しています。
サンプルサイズ（40, 60, 100, 500）と回帰係数の値（Type1 errorの評価の場合0、バイアス評価の場合0.92）について複数の条件で実験を行うので、furrrパッケージにより並列処理を行うことにしました（参考文献[5]や[6]を参照）。最後のあたりは汚いコードになってしまいすみません。
furrrパッケージ、purrrパッケージで何をやっているのかについてはまた別の機会に書ければと思います。

library(furrr)
library(tictoc) #このパッケージで処理にかかった時間を測れます

# シミュレーション設定をデータフレームにする
## tidyrパッケージのexpand_gridを使用
sim_params <- expand_grid(
  nsim = 7300,
  n = c(40, 60, 100, 500),
  beta_T = c(0, 0.92)
)

# 各設定でのシミュレーション実行
plan(multisession) # 並列処理
tic() # tic() toc()で時間を測る
result <- future_pmap_dfr(.l = sim_params, .f = sim,
                          .options = furrr_options(seed = 123L),
                          .progress = TRUE)
toc()

# シミュレーション結果を表にする
## Pirracchio論文の結果
bmc_result <- data.frame(N = rep(c("N = 40", "N = 60", "N = 100", "N = 500"), each = 3),
                         Measure = rep(c("Type I error", "Bias", "Percent bias"), 4),
                         IPTW_BMC = c(0.048, 0.057, 6.5, 0.046, 0.036, 4,
                                   0.046, 0.018, 2, 0.046, 0.0007, 0.8))
## 今回のシミュレーションの結果
result_sum <- result %>%
  mutate(est_or = exp(estimate),
         reject = case_when(p.value < 0.05 ~ 1, TRUE ~ 0)) %>%
  group_by(n, beta_T) %>%
  summarise(true_or = mean(true_or),
            est_or = mean(est_or),
            alphaerror_power = mean(reject)) %>%
  mutate(bias = est_or - true_or,
         relbias = (bias/true_or)*100)
sim_null <- result_sum %>% filter(beta_T == 0.00) %>%
  dplyr::select(n, alphaerror_power)
sim_alt <- result_sum %>% filter(beta_T == 0.92) %>%
  dplyr::select(n, bias, relbias)
sim_result <- inner_join(sim_null, sim_alt, by = "n") %>%
  pivot_longer(cols = c(alphaerror_power, bias, relbias),
               names_to = "Measure", values_to = "IPTW_SIM")

all_result <- bmc_result %>% mutate(IPTW_SIM = sim_result$IPTW_SIM)

## 表の作成と保存
library(gt)
res_tbl <- gt(all_result, groupname_col = "N", rowname_col = "Measure") %>%
  fmt_number(columns = c("IPTW_BMC", "IPTW_SIM"), decimals = 3)
gtsave(res_tbl, "sim_restbl_20220411.png")

シミュレーション結果は以下の通りです。
左側（IPTW_BMC）にPirracchio論文での数値を、右側（IPTW_SIM）に今回のシミュレーションで得られた数値を表示しました。

Type1 errorの確率については、今回のシミュレーションでも全てのサンプルサイズで指定した水準の5%を下回っていました。
サンプルサイズが小さい時にType1 errorの確率が増加してしまうのは他の統計手法ではありがちなのですが、IPTWに関しては、サンプルサイズが各グループ40程度でも（今回の設定の範囲、例えば共変量の数が4つだったり、共変量どうしの相関がほぼないような場合においてではありますが）その危険性は小さいように思います。

バイアスについては「Percent bias」を見てみましょう。今回の実験では、サンプルサイズが各グループ40の時にはPirracchio論文と比べ推定のバイアスが2倍程度になっていました。
特にPirracchio論文では「サンプルサイズが40でもバイアスは10%を下回る」ことが述べられていたのですが、シミュレーションでは13%ほどと少し大きめの結果になりました。これはオッズ比のスケールで表しているので、おおざっぱに言うと、真のオッズ比が2.56であるところを2.89（2.56の1.13倍）と推定していると理解できます。*9
一方、各グループのサンプルサイズが100の場合にはPirracchio論文に近い値に、500の場合では過小推定の方向に若干のバイアスが生じる結果になりました。しかしながらそこまで深刻なバイアスは生じないだろうことが伺えます。

おわりに

サンプルサイズが小さい場合に傾向スコア解析を行うとどうなるのかについて、いくつかの研究報告があるようです。
今回は、比較的古い報告であるPirracchioらの結果を紹介し、IPTWの評価結果の一部についてRで再現シミュレーションを行いました。

おおむね近い結果は得られたように思いますが、サンプルサイズが各グループ40〜60の場合には報告されている結果の2倍近いバイアスが確認されました。Pirracchio論文とはシミュレーションコードや設定条件が異なるためかもしれませんし、今回の反復回数では十分な精度が得られなかった可能性もあります。
また、共変量の数や相関の程度、さらにIPTWの方法（用いる重みなど）によっても結果は異なってくると考えられます。

一応、今回の結果を踏まえると、サンプルサイズが非常に小さい場合には多少のバイアスが生じる可能性があると思っておいた方が無難なのかもしれません。

参考文献

[1] Pirracchio, R., Resche-Rigon, M., & Chevret, S. (2012). Evaluation of the propensity score methods for estimating marginal odds ratios in case of small sample size. BMC medical research methodology, 12(1), 1-10.
[2] Bottigliengo, D., Baldi, I., Lanera, C., Lorenzoni, G., Bejko, J., Bottio, T., ... & Gregori, D. (2021). Oversampling and replacement strategies in propensity score matching: a critical review focused on small sample size in clinical settings. BMC medical research methodology, 21(1), 1-16.
[3] van der Wal, W. M., & Geskus, R. B. (2011). ipw: an R package for inverse probability weighting. Journal of Statistical Software, 43, 1-23.
[4] Austin, P. C., & Stafford, J. (2008). The performance of two data-generation processes for data with specified marginal treatment odds ratios. Communications in Statistics—Simulation and Computation®, 37(6), 1039-1051.
[5] furrr パッケージで R で簡単並列処理 | Atusy's blog
[6] Monte Carlo Simulations: Going parallel with the furrr package

*1:傾向スコアについては優れた解説が色々なところで見つかるので、今回は説明を省略します。また機会があればここでも書きたいと思います。

*2:具体的な数は諸説あります。

*3:マッチングは処置群のある個体に対して対照群の似たような個体をマッチング相手として見つける方法ですが、復元マッチングとは、一度マッチング相手となった個体を、他の処置群の個体のマッチング候補として再利用することを言います。非復元マッチングでは、一度マッチングされた個体は以降はマッチング候補として使用されません。

*4:サンプルサイズの数値がグループごとなのかグループ合わせた合計なのかは明記されていませんでしたが、文中の記述からはグループごとの数値と思われますので、再現シミュレーションもそれを前提に行いました。

*5:有限回のシミュレーションによる推定なので、誤差によって5%を少しだけ超えることはあります。

*6:「他の共変量をある特定の値に固定したときの」オッズ比のことです。これに対して、「他の共変量について、全個人の平均をとったときの」オッズ比である周辺オッズ比という概念もあり、両者は必ずしも一致しません。前者はある共変量の値をもつ特定の個人についての効果量、後者は集団全体での効果量を意味しています。

*7:おおざっぱに言えば「処置群に入りそうにないような特徴を持つ個人のデータには、その珍しさに応じた重みをつける」ことによって、そのような特徴を持つ人の人数が対照群と同じくらいになるような補正を行っています（対照群でも同様）。

*8:重みづけ解析を行う方法にも色々と選択肢があり、注意が必要と言われているのですが、今回の本題とは外れるので別の機会に書ければと思います。

*9:実際には2.89は「同じようなデータを得て、同じ解析を行う」ことを何度も繰り返した時の平均なので、あくまで傾向ではあります。