観察研究で治療効果を推定するとき、重症度の高い患者ほど治療を受けやすい集団では、未調整の群間差がその治療の効果をそのまま表しているとは限りません。解析用データセットに並んだ共変量のうち、どの共変量をモデルへ入れるかで推定値は目に見えて動きます。
ただし、治療とも結果とも関連する変数を片端から入れると、交絡を取り除くつもりで別の偏りを持ち込むことになります。この回では、交絡が起きている状態の定義から始めて、手元の共変量のどれを交絡因子と呼ぶのかまでを扱います。
交絡が起きている状態の定義
交絡という言葉は、日常的には他の要因が混ざっているという程度の意味で使われますが、因果推論では、比較する治療群と対照群の交換可能性が崩れている状態として捉えます。Greenland・Robins（1986）は、治療群と対照群が反事実的な結果について交換可能か、という条件から交絡を整理しています。
$A$は治療の有無、$Y$は観測される結果、$Y^{a}$は$A=a$を与えたときの潜在結果、$\mathbf{L}$は調整に使う共変量の集合を表すとした場合、以下の式が成り立ちます。
$$Y^{a} \perp A \mid \mathbf{L} \quad (a = 0,\ 1) \tag{1}$$
式(1)は、共変量$\mathbf{L}$の値を揃えた集団の中では、実際にどちらの治療を受けたかと、それぞれの治療を受けた場合の潜在結果とが独立になるという条件を表しています。必要な共変量で条件づけてもこの交換可能性を満たせない状態では、交絡が残ります。
この定義で目を引くのは、交絡が個々の変数へ貼られた札としてではなく、比較しようとしている2つの集団の関係として書かれている点です。ある変数を交絡因子と呼べるかどうかは、その変数の分布や相関係数だけからは決まりません。
古典的な基準が拾いすぎる変数
教科書でよく見る基準は、治療とも結果とも関連し、かつ因果経路の途中にない変数を交絡因子と呼ぶ、というものです。この基準がどこまで働くかを見るために、治療の割り付けが重症度と施設の処方傾向から決まり、結果が治療、中間変数、重症度、予後因子から決まる6,000例のランダムに生成したデータセットを用意しました。
このデータでは治療群が2,691例、対照群が3,309例で、構造式に置いた真の値は、中間変数を経由しない直接効果が1.0、治療から中間変数への効果が0.8、中間変数から結果への効果が1.5です。
$$\tau = \beta_A + \alpha_M \beta_M = 1.0 + 0.8 \times 1.5 = 2.2 \tag{2}$$
式(2)の$\tau$は治療が結果へ与える総効果で、$\beta_A$は中間変数を経由しない直接効果、$\alpha_M$は治療から中間変数への効果、$\beta_M$は中間変数から結果への効果を表します。中間変数を通る分が1.2あるため、推定の目標になる総効果は2.2です。
この2.2に対して、共変量を何も入れない比較では3.693（95%信頼区間3.550〜3.836）となり、真の値より1.49大きい側へ出ました。重症度が高いほど治療されやすく、かつ結果の値も高いため、未調整の差には治療の効果と重症度の差が重なっています。
ところが、治療とも結果とも関連するという基準だけで候補変数を選別すると、重症度（治療との相関0.354、結果との相関0.718）のほかに、施設の処方傾向、中間変数、追跡中の追加検査も候補に入ります。結果だけの原因になる予後因子は、治療との相関が0.008、治療を説明変数にした回帰のp値が0.545で、この基準では落ちています。
つまり、関連の有無だけで選ぶと、本当に調整が必要な変数と、調整してはいけない変数、調整しても交絡の除去には寄与しない変数を適切に区別できません。Hernánら（2002）が述べているように、因果構造の仮定を先に置かない限り、同じ関連の統計量は交絡にも選択による偏りにも対応しえます。
治療だけの原因を調整したときの代償
施設の処方傾向を重症度に足して調整すると、推定値は2.257から2.263へほとんど動かず、標準誤差だけが0.058から0.064へ上がりました。同じ生成過程で300回作り直しても、平均の偏りは0.000のままで、推定値の標準偏差が0.058から0.065へ広がっています。
この変数は治療の受けやすさを左右しますが、結果へ向かう矢印を持ちません。そのため調整しても偏りは減らず、治療の説明されない変動だけが削られて、推定の精度が落ちます。
VanderWeele（2019）は、治療の原因または結果の原因となる共変量を調整候補とし、治療だけの原因になる変数は原則として外す方針を勧めています。私も、処方の決定に強く効くという理由だけで施設要因を調整集合へ入れることは適切ではありません。
その理由は、偏りが減らないまま区間が広がるうえに、未測定の交絡が残っている場合には、残った偏りをかえって強めることがあるためです。
この判断は、傾向スコアのモデルに入れる変数を選ぶ場面でも使えます。治療の受けやすさを最もよく予測する変数を集めるという発想だけでモデルを組むと、傾向スコアマッチングでも同様に精度を損なうことがあります。
中間変数の調整による過剰調整
治療開始後に測った検査値が結果とよく相関していると、説明力が上がるという理由で回帰へ入れたくなります。この場合の推定値は0.990（95%信頼区間0.907〜1.073）で、総効果2.2からは1.21小さい側へ動きました。
この0.99という値は、精度を欠いた推定値というよりも別の量に当たっており、300回の反復平均1.000は構造式に置いた直接効果1.0とほぼ一致します。中間変数を固定すると、治療がその変数を通して結果へ与える1.2の分が推定から外れるため、残った値が直接効果の側へ寄ります。
Schistermanら（2009）は、この種の調整を過剰調整として整理し、総効果を推定したい解析では避ける操作として扱っています。治療の後に測られた変数を共変量へ入れるかどうかは、測定時期だけではなく、その変数が治療の影響を受ける下流の変数かどうかで決まります。
もう一つ、実務で難しくなるのは、投与開始後の検査値が、次の治療の決定に効く交絡因子の更新値でもある場面です。治療で動く変数が次の治療を決めているなら、単純な回帰調整では目的とする総効果を一般には適切に推定できず、周辺構造モデルのような重み付けの方法へ移ることになります。
合流点の調整が作る見かけの関連
治療と結果の両方から矢印を受ける変数を調整すると、推定値は-0.207（95%信頼区間-0.270〜-0.144）まで落ち、符号まで逆転しました。ここで入れた追加検査の実施は、治療を受けた人ほど行われやすく、結果の値が高い人ほど行われやすいという、両者に共通する下流の変数です。
この符号の逆転は、追加検査の値で条件づけると、治療を受けていることと結果が高いことの間に見かけ上の負の関連が生じるために起こります。300回の反復でも平均は-0.281で、標本を取り直しても中心は0付近へ戻りません。
関連が強いことは、調整してよい理由になりません。
なお、この推定値の95%信頼区間の幅は、今回比べたどの調整集合よりも狭く、300回の標準偏差も0.033で最小でした。区間が狭いことは偏りが小さいことを保証しないという、分かりやすい例になっています。
調整集合として定義する交絡因子
ここまで見た5つの候補には、治療と結果の両方に関連する変数だけでなく、どちらか一方だけの原因になる変数や、治療後の下流にある変数も含まれていました。VanderWeele・Shpitser（2013）は、交絡因子を1つ1つの変数の性質として定義する試みを整理したうえで、識別に十分な調整集合との関係で定義する立場を示しています。
この立場に立つと、どの変数を交絡因子と呼ぶかという言い方は、変数を1つ取り出しただけでは答えが決まらない形になります。今回の生成過程では、重症度を調整することで式(1)の条件を満たす調整集合を作れますが、別の因果構造のもとでは同じ変数がその役目を持つとは限りません。
もう一点、調整の前後で推定値が動いたことを交絡の証拠として読む習慣にも注意が要ります。Greenlandら（1999）が整理したように、オッズ比やハザード比では、交絡がなくても共変量を入れると値が動く非縮約可能性が働くため、値の変化だけからは交絡の有無を判定できません。
今回の線形モデルでは差の尺度を使っているため、オッズ比のような非縮約可能性による変化ではありません。推定値の移動は、この生成過程で置いた交絡、過剰調整、合流点への調整の影響として読めます。
もっとも、推定の枠組みを変えても調整集合を決める作業は残り、重み付けに移ったときに何がどう変わるかは、回帰調整とIPTWの違いを扱った回で数値と一緒に説明しています。
報告と研究での位置づけ
報告の面では、調整した共変量の一覧だけでなく、調整集合を選んだ根拠も示しておく必要があります。Hernánら（2002）は先天異常の疫学研究を例に、同じ共変量が、仮定する因果構造によって交絡因子にも選択による偏りの源にもなることを扱っています。
そのため、調整した共変量の一覧を表へ並べるだけでなく、どの変数が治療の原因で、どの変数が結果の原因かを、投稿前に文章か図で書き残しておくと、査読での説明が短く済みます。私が解析を担当する研究では、変数選択をp値の大小からではなく、この因果の向きの整理から始めます。
もう一つ書き残す価値があるのは、未測定の共変量が残る可能性です。観察研究では式(1)の条件をデータだけから検証できないため、どの程度の未測定交絡があれば結論が覆るかを示す感度分析を報告に加える方法があります。
Rによる調整集合ごとの比較
ここまでの数値を出したコードを、実行できる形で載せます。使うのはbase Rとstatsだけで、外部パッケージは読み込みません。
set.seed(42)
n_obs <- 6000
generate_cohort <- function(n) {
severity <- rnorm(n) # 治療と結果の共通原因
pattern <- rnorm(n) # 治療だけの原因
prognostic <- rnorm(n) # 結果だけの原因
treatment <- rbinom(n, size = 1, prob = plogis(-0.2 + severity + 1.2 * pattern))
mediator <- 0.8 * treatment + rnorm(n) # 治療の下流、結果の原因
outcome <- treatment + 1.5 * mediator + 2 * severity + prognostic + rnorm(n)
detection <- 0.9 * treatment + 0.9 * outcome + rnorm(n) # 治療と結果の共通の下流
data.frame(treatment, outcome, severity, pattern, prognostic, mediator, detection)
}
cohort <- generate_cohort(n_obs)ここで用意したのは、治療と結果の共通原因に当たる重症度、治療だけの原因になる施設の処方傾向、結果だけの原因になる予後因子、治療の下流にあって結果の原因にもなる中間変数、治療と結果の両方から影響を受ける追加検査を持つ、6,000例のデータセットです。乱数の種は42に固定してあります。
adjustment_sets <- c(crude = "treatment", severity = "treatment + severity",
pattern = "treatment + severity + pattern",
prognostic = "treatment + severity + prognostic",
mediator = "treatment + severity + mediator",
detection = "treatment + severity + detection")
fit_effect <- function(rhs, data) {
unname(coef(lm(as.formula(paste("outcome ~", rhs)), data = data))["treatment"])
}
estimates <- vapply(adjustment_sets, fit_effect, numeric(1), data = cohort)
round(estimates, 3)#> crude severity pattern prognostic mediator detection
#> 3.693 2.257 2.263 2.233 0.990 -0.207この出力で真の値2.2に近いのは、重症度だけを入れた2.257と、そこへ施設要因や予後因子を足した2.263、2.233です。読み違えやすいのは後ろの2つで、中間変数を入れた0.990は総効果ではなく直接効果に近い量を捉えており、追加検査を入れた-0.207は合流点への条件づけによって偏りが生じた結果です。
n_rep <- 300
rep_matrix <- matrix(NA_real_, nrow = n_rep, ncol = length(adjustment_sets),
dimnames = list(NULL, names(adjustment_sets)))
for (i in seq_len(n_rep)) {
rep_data <- generate_cohort(n_obs)
rep_matrix[i, ] <- vapply(adjustment_sets, fit_effect, numeric(1), data = rep_data)
}
round(rbind(bias = colMeans(rep_matrix) - 2.2, sd = apply(rep_matrix, 2, sd)), 3)#> crude severity pattern prognostic mediator detection
#> bias 1.389 -0.001 0.000 -0.003 -1.200 -2.481
#> sd 0.071 0.058 0.065 0.050 0.042 0.033反復の出力では、未調整の偏りが1.389、重症度だけを調整したときが-0.001で、共通原因を1つ入れるかどうかの差がはっきり出ています。同時に、施設要因を足すと標準偏差が0.058から0.065へ広がり、予後因子を足すと0.050へ縮んでおり、偏りに効く調整と精度に効く調整が別々に見えます。
ここで注意したいのは、中間変数と追加検査の列では、総効果2.2からの差が-1.200と-2.481のまま、標準偏差が最も小さくなっている点です。中間変数では推定対象が総効果から直接効果の側へ変わり、追加検査では合流点への調整による偏りが残っているため、ばらつきが小さくても総効果2.2へ近づきません。
plot(estimates, 6:1, xlim = c(-0.8, 4.1), pch = 19, yaxt = "n", ylab = "",
xlab = "Estimated treatment effect")
abline(v = 2.2, lty = 2)
axis(2, at = 6:1, labels = names(estimates), las = 1)
boxplot(rep_matrix, horizontal = TRUE, las = 1, outline = FALSE,
xlab = "Treatment coefficient over 300 replications")
abline(v = 2.2, lty = 2)
解析計画で先に決めること
解析を始める前に決めておく価値があるのは、共変量の一覧そのものではなく、各変数が治療より前にあるのか後にあるのか、そして治療や結果とどのような因果関係にあるのかです。この向きを整理できていれば、調整集合は式(1)の条件を満たすために必要な変数から考えられます。
今回のデータで言えば、重症度だけを入れた推定が真の値に近く、そこへ結果だけの原因になる予後因子を足す調整は、偏りをほとんど変えずに区間を狭める方向へ働きました。一方、中間変数を入れると推定対象が総効果から直接効果の側へ変わり、追加検査を入れると合流点への調整による偏りが生じました。
実際のデータでは、測定の時期、欠測の入り方、治療方針の変更の記録まで見たうえで調整集合を決めることになり、Dr.データサイエンスでは、こうした解析方針の相談を承っています。
参考文献
- Greenland S, Robins JM. Identifiability, exchangeability, and epidemiological confounding. International Journal of Epidemiology. 1986;15(3):413-419.
- Greenland S, Pearl J, Robins JM. Confounding and collapsibility in causal inference. Statistical Science. 1999;14(1):29-46.
- Hernán MA, Hernández-Díaz S, Werler MM, Mitchell AA. Causal knowledge as a prerequisite for confounding evaluation: an application to birth defects epidemiology. American Journal of Epidemiology. 2002;155(2):176-184.
- Schisterman EF, Cole SR, Platt RW. Overadjustment bias and unnecessary adjustment in epidemiologic studies. Epidemiology. 2009;20(4):488-495.
- VanderWeele TJ, Shpitser I. On the definition of a confounder. The Annals of Statistics. 2013;41(1):196-220.
- VanderWeele TJ. Principles of confounder selection. European Journal of Epidemiology. 2019;34(3):211-219.
