データベース研究で薬の効果を推定しようとすると、その薬を使っている患者と使っていない患者では治療が始まった理由そのものが異なることがあり、観察されたイベント率の差には治療の効果と病状の差が混ざって入ってきます。この混ざり方は解析手法を選ぶ前の段階、つまり誰を対照群に置き、追跡の起点をいつに定めるかで、大きさも向きも変わります。
そこで今回は、乱数で生成した12,000例のデータセットを使い、薬剤Aの効果を真のハザード比0.80と置いたうえで、非使用者を対照にした場合、同じ適応に使われる薬剤Bを対照にした場合、そして薬剤Aの既存使用者まで取り込んだ場合で推定値がどこへ動くかを比べます。
この比較を通して判断できるようになるのは、対照群をどう取るかという設計上の選択が、交絡の調整よりも前に効果の推定値を左右するという関係です。新規使用者デザインと実薬対照デザインがそれぞれ何を揃えるための設計なのか、そしてどこから先は調整や感度分析に任せるしかないのかまで扱います。
薬が選ばれる理由が結果に混ざる
重症な患者ほど薬物治療を受けやすい状況では、治療群と非治療群のイベント率をそのまま比べると、薬の効果と重症度の差が同じ数字の中に入り込み、薬が有害に見えることさえあります。
この現象をSalasら（1999）は適応による交絡（confounding by indication）という用語を取り上げ、疫学研究でその使われ方に幅があることを整理しました。同じ論文が指摘しているのは、用語が指す範囲が論文ごとに異なるため、この語を使うときには、どの要因を介した交絡を問題にしているのかを具体的に示す必要があるという点です。
ここからの推定には、Cox比例ハザードモデルを固定して使います。時点$t$のハザードを$\lambda(t)$、薬剤Aの使用を$A$、測定できた共変量ベクトルを$\mathbf{x}$と書くと、次の形になります。
$$\lambda(t \mid A, \mathbf{x}) = \lambda_0(t)\exp(\beta A + \boldsymbol{\gamma}^\top \mathbf{x}) \tag{1}$$
式(1)は、薬剤Aを使ったかどうかと共変量の値によって、基準となるハザード$\lambda_0(t)$が何倍になるかを表しています。$\beta$は薬剤Aのログハザード比、$\boldsymbol{\gamma}$は共変量の係数ベクトルで、$\exp(\beta)$が共変量で条件づけた薬剤Aのハザード比です。
なお、同じデータに当てはめても回帰調整と重み付けでは推定する量が異なり得るという論点は、生存時間解析における因果推論で扱っています。
適応による交絡はどちら向きに効くか
今回のデータセットでは、重症度スコアが高いほど薬物治療を受ける確率が上がり、同時に重症度そのものがハザードを押し上げるように設定しました。生成された12,000例のうち非使用者が6,709例、薬剤Aが2,698例、薬剤Bが2,593例で、3年間の追跡で2,429件のイベントが発生しています。
この設定のもとで群ごとの重症度の平均を見ると、非使用者が-0.448、薬剤Bが0.502、薬剤Aが0.579でした。群間の均衡は平均差そのものより標準化差で読む方が扱いやすく、2群の平均を$\bar{x}_1$、$\bar{x}_0$、分散を$s_1^2$、$s_0^2$として次のように計算します。
$$d = \frac{\bar{x}_1 – \bar{x}_0}{\sqrt{(s_1^2 + s_0^2)/2}} \tag{2}$$
式(2)は、2群の平均の差を、両群の分散から求めた標準偏差で割った値です。この定義で計算すると、薬剤A群と非使用者群の標準化差は1.167、薬剤A群と薬剤B群では0.088にとどまりました。
この2つの標準化差を並べると、今回のデータでは薬を使う人と使わない人の間に重症度で1標準偏差を超える開きがあるのに対し、2つの薬剤の間では差が小さいことが分かります。
新規使用者デザインが揃えるもの
新規使用者デザイン（new-user design）は、対象者を薬剤の使用開始時点で捉え、その時点を追跡起点に置く設計です。Ray（2003）は臨床試験の外で薬効を評価する枠組みとしてこれを整理し、既に使っている患者を集めた解析では、治療開始前の共変量を適切に捉えにくいことと、開始直後のリスク期間が観察から落ちることが同時に起こり得ると述べています。
この2点のうち、実務で見落とされやすいのは後者です。使用開始から最初の数か月に有害事象が集中する薬剤では、その期間を通過した患者だけを集めることで、開始直後にイベントを起こした患者が解析対象から除かれます。
もっとも、追跡起点の置き方が推定値を動かす仕組みは、ターゲット試験エミュレーションで扱った時間ゼロの議論にもつながります。Suissa（2008）が整理したイモータルタイムバイアスは、曝露分類のために将来の生存を必要とする期間を誤って曝露時間へ含めることで生じるため、既存使用者を組み入れることで生じる選択とは区別して考える必要があります。
既存使用者を入れると誰が消えるか
解析用のデータセットでは、薬剤Aについてだけ研究開始より前から使っている患者を取り込み、薬剤Bは新規開始者に限る、という非対称な集め方を再現しました。使用歴の長さは0年から8年の一様分布とし、研究開始時点までにイベントを起こした薬剤A使用者はデータに入らない扱いにしています。
この操作で薬剤A群2,698例のうち観察対象として残ったのは1,911例、割合にして70.8%でした。落ちた787例は研究開始時点までにイベントを起こした患者なので、残った集団の重症度の平均は0.579から0.382へ下がっています。
落ちた患者がいるという点だけを見れば、傾向スコアマッチングで対照が絞られる状況と似ていますが、こちらは解析ではなく組み入れの段階で、しかも研究開始時点までイベントを起こさず生存したことを条件にして選別が起きています。マッチングなら誰が外れたかを解析対象者の中で確認できますが、既存使用者だけを研究開始時点で抽出したデータでは、それ以前にイベントを起こした患者の情報を解析データから直接確認できない場合があります。
実薬対照が減らせる偏りと残る偏り
実薬対照デザイン（active-comparator design）は、非使用者ではなく同じ適応に使われる別の薬剤を対照に置く設計です。Lundら（2015）はこの設計と新規使用者デザインを組み合わせた枠組みを薬剤疫学で広く用いられる設計として整理し、Yoshidaら（2015）も、非使用者との比較では治療を受けるかどうかに関わる要因が交絡として入りやすいことから、両者を組み合わせる利点を説明しています。
同じ適応に使われる実薬対照を選べる場合には、2剤のどちらを選ぶかに関わる要因の分布が、薬物治療を受ける人と受けない人を比べる場合より近くなることがあります。今回の設定では重症度が薬剤の選択にほとんど効かないため、図1のとおり2群の分布がほぼ重なりました。
もっとも、実薬対照に替えれば交絡が消えるわけではありません。腎機能や併用薬の制約によって2剤のどちらかが選ばれる場面では、その要因が薬剤選択とアウトカムの両方に効いていれば偏りは残り、群間差が小さいために見つけにくいこともあります。
調整で足りる場面と足りない場面
今回の設定では重症度を変数として持ち、さらにデータ生成に対応したCoxモデルを当てはめているため、非使用者を対照にした場合でも重症度を調整すると推定値は真値の0.80とほぼ一致しました。
とはいえ、この一致から対照の取り方は調整で置き換えられると読むのは早すぎます。今回は交絡因子を重症度1つだけに限り、しかもその値を完全に観測している設定にしてあり、実際の診療記録では治療方針を決めた要因のうち記録に残らない部分が存在することがあります。
そこで私なら、研究目的に合う実薬対照が取れる状況では非使用者対照を主解析にしません。未測定の適応要因が残る可能性を考えると、調整前から患者背景が近い比較を作る方が、未測定交絡の余地を小さくする方向に働くためです。
Rで対照と時間ゼロを入れ替える
設定の全体は、重症度が治療を受けるかどうかとハザードの両方に強く効き、薬剤の選択には弱くしか効かない構造にしてあります。survivalパッケージのcoxph関数を使い、対照群と組み入れ方を変えた5通りの推定値を比べます。
library(survival)
set.seed(42)
n_patients <- 12000
true_log_hr <- log(0.80) # 薬剤Aの真のログハザード比
admin_censor_years <- 3.0 # 追跡上限
severity <- rnorm(n_patients, mean = 0, sd = 1)
treated <- rbinom(n_patients, size = 1, prob = plogis(-0.30 + 1.30 * severity))
drug_a <- ifelse(treated == 1,
rbinom(n_patients, size = 1, prob = plogis(0.00 + 0.10 * severity)), 0)
group <- ifelse(treated == 0, "nonuser", ifelse(drug_a == 1, "drug_a", "drug_b"))このデータ生成では、基礎発生率0.06に重症度と薬剤Aの効果を掛けたハザードを置き、イベント時間を指数分布から発生させています。3年を超えた分は打ち切りとして扱いました。
hazard_rate <- 0.06 * exp(0.90 * severity + true_log_hr * drug_a)
event_time <- rexp(n_patients, rate = hazard_rate)
followup <- pmin(event_time, admin_censor_years)
event <- as.integer(event_time <= admin_censor_years)
cohort <- data.frame(severity = severity, group = group, drug_a = drug_a,
event_time = event_time, followup = followup, event = event,
stringsAsFactors = FALSE)続いて、非使用者を対照にした集合と薬剤Bを対照にした集合をそれぞれ作り、調整なしと重症度調整ありの両方を当てはめます。
nonuser_set <- cohort[cohort$group %in% c("nonuser", "drug_a"), ]
active_set <- cohort[cohort$group %in% c("drug_b", "drug_a"), ]
fit_nonuser <- coxph(Surv(followup, event) ~ drug_a, data = nonuser_set)
fit_nonuser_adj <- coxph(Surv(followup, event) ~ drug_a + severity, data = nonuser_set)
fit_active <- coxph(Surv(followup, event) ~ drug_a, data = active_set)
fit_active_adj <- coxph(Surv(followup, event) ~ drug_a + severity, data = active_set)さらに、既存使用者を取り込んだ集合は、薬剤A使用者に0年から8年の使用歴を与え、研究開始時点までイベントを起こさずに到達した患者だけを残して作ります。
prior_use_years <- runif(n_patients, min = 0, max = 8.0)
survived_to_entry <- cohort$group == "drug_a" & cohort$event_time > prior_use_years
prevalent_a <- cohort[survived_to_entry, ]
prevalent_a$followup <- pmin(prevalent_a$event_time - prior_use_years[survived_to_entry],
admin_censor_years)
prevalent_a$event <- as.integer(prevalent_a$event_time - prior_use_years[survived_to_entry] <=
admin_censor_years)
prevalent_set <- rbind(prevalent_a, cohort[cohort$group == "drug_b", ])
fit_prevalent <- coxph(Surv(followup, event) ~ drug_a, data = prevalent_set)#> HR non-user comparator : 1.924
#> HR non-user adjusted : 0.798
#> HR active comparator : 0.868
#> HR active adjusted : 0.786
#> HR prevalent user : 0.766
#> True HR : 0.800この出力でまず目を引くのは、非使用者を対照にした未調整の推定値1.924で、95%信頼区間は1.746から2.121でした。真値を0.80と置いた設定から1を大きく超える値が出ているので、この生成過程では薬剤Aの効果とは反対方向の関連が観察されています。
読み違えやすいのは、この1.924という値が信頼区間の狭い、いかにも安定した推定値に見える点です。同じ生成過程のもとで標本数を増やしても、系統的な群間差そのものは解消されないため、精度の高さだけから偏りの有無を判断することはできません。
そこで対照を薬剤Bに替えると、調整を何もしないまま推定値は0.868（0.783から0.962）まで真値へ近づきました。重症度で調整すると0.786（0.709から0.871）、非使用者対照でも調整すれば0.798（0.713から0.892）となり、今回の設定ではどちらも真値の近くに落ち着きます。
一方、既存使用者を取り込んだ推定値は0.766（0.682から0.861）で、新規使用者に限った実薬対照の0.868より小さい側へ動きました。この生成過程では対照群を同じ薬剤B群に固定しているため、この差には薬剤A群の組み入れ方を変えた影響が表れています。
これらの図は、上のコードで作ったcohortと当てはめ結果を使って描いています。図1の箱ひげ図は次の指定で書き出しました。
png("severity-by-group.png", width = 1000, height = 620, res = 130)
boxplot(severity ~ factor(group, levels = c("nonuser", "drug_b", "drug_a")), data = cohort,
names = c("Non-user", "Drug B", "Drug A"), xlab = "Treatment group",
ylab = "Disease severity score",
main = "Severity differs by treatment, not between the two drugs")
dev.off()公表研究と報告での扱われ方
設計の違いが実データで推定値を変えた例として、Hernánら（2008）が閉経後ホルモン療法と冠動脈疾患の関係を再解析した研究があります。観察データを無作為化試験に近い形で組み立て直し、使用開始時点を追跡起点に置いて解析したところ、それ以前の観察研究で報告されていた強い予防的な関連とは異なる推定値が得られたと報告しています。
この再解析から読み取れるのは、組み入れ条件と追跡起点を臨床試験に近い形へ揃えることで、それ以前の観察解析とは推定値が変わり得るという点です。論文の方法欄では、対象者をいつ組み入れ、どの時点を追跡起点としたかを一文で書けるかどうかが、そのまま読者にとっての判断材料になります。
報告の実務でもう一点、私が投稿前の原稿で確認するのは、対照群の定義に非使用者という語が出てくるかどうかです。この語が出てきた場合には、なぜ実薬対照を取らなかったのかという理由が本文のどこかに書かれている方が、比較集団を選んだ根拠が伝わります。
研究計画で先に決めておくこと
解析前に決まっているべきものは、共変量の一覧より先に、誰を組み入れ、いつから数え始めるかの2点です。今回のデータセットでは、この2点と調整条件を変えることで推定値が1.924から0.766まで動きました。
ただし、実薬対照が取れるかどうかは疾患領域や研究目的によって決まるため、常にこの設計を選べるわけではありません。それでも、研究目的に合う同じ適応の薬剤が実薬対照として利用できるのに非使用者を対照に選ぶ計画は、未測定の適応要因による群間差を大きく残す可能性があるため、適切ではありません。
実際のデータベース研究では、過去の処方記録がどれだけ遡れるか、洗い出し期間を何か月に置くか、実薬対照の候補が同じ適応と言えるか、といった条件を、データの構造と研究目的の両方から判断することになります。
そのため、Dr.データサイエンスでは、こうした設計と解析方針のご相談を承っています。
参考文献
- Ray WA. Evaluating medication effects outside of clinical trials: new-user designs. American Journal of Epidemiology. 2003;158(9):915-920.
- Lund JL, Richardson DB, Stürmer T. The active comparator, new user study design in pharmacoepidemiology: historical foundations and contemporary application. Current Epidemiology Reports. 2015;2(4):221-228.
- Yoshida K, Solomon DH, Kim SC. Active-comparator design and new-user design in observational studies. Nature Reviews Rheumatology. 2015;11(7):437-441.
- Salas M, Hofman A, Stricker BH. Confounding by indication: an example of variation in the use of epidemiologic terminology. American Journal of Epidemiology. 1999;149(11):981-983.
- Hernán MA, Alonso A, Logan R, Grodstein F, Michels KB, Willett WC, Manson JE, Robins JM. Observational studies analyzed like randomized experiments: an application to postmenopausal hormone therapy and coronary heart disease. Epidemiology. 2008;19(6):766-779.
- Suissa S. Immortal time bias in pharmaco-epidemiology. American Journal of Epidemiology. 2008;167(4):492-499.
