降圧薬の観察研究で、年齢と治療前の血圧を調整して12週後の収縮期血圧を比較した解析が、方法欄に「多変量解析を行った」とだけ書かれていることがあります。説明変数は複数ありますが、アウトカムは収縮期血圧の1つです。英語論文では、このような回帰はmultivariable regressionと呼ばれます。これに対して、収縮期血圧と拡張期血圧のように複数のアウトカムを同時に扱う回帰はmultivariate regressionです。
日本語では両者が「多変量解析」と書かれることがあり、方法欄だけではモデルの左辺に何本のアウトカムが置かれているのか分からなくなります。回帰の文脈では、この区別は用語だけの問題ではありません。何を推定できるのか、アウトカム間の共分散がどこで必要になるのかまで変わってきます。
この違いを確かめるため、年齢によって治療を受ける確率が異なる600例のデータセットをランダムに生成しました。収縮期血圧と拡張期血圧を別々に解析した場合と、2つを同時に解析した場合を比べます。先に押さえておきたいのは、同じ対象者、同じ説明変数を使う通常の多変量線形モデルでは、各アウトカムの回帰係数とその標準誤差は、アウトカムごとに当てはめた線形回帰と一致することです。違いが現れるのは、2つのアウトカムをまたいだ検定や対比を求める場面です。
左辺が1本か複数かで意味が変わる
回帰分析の目的と単回帰分析と重回帰分析の違いでは、説明変数とアウトカムの関係を扱いました。ここでは、その回帰式の左辺に注目します。対象者$i$の説明変数を並べたベクトルを$\mathbf{x}_i$、係数を$\boldsymbol{\beta}$とすると、1つのアウトカムを扱う線形回帰は次のように書けます。
$$ y_i = \beta_0 + \boldsymbol{\beta}^\top \mathbf{x}_i + \varepsilon_i \tag{1} $$
式(1)では、対象者$i$について観測されるアウトカム$y_i$は1つです。右辺に治療、年齢、ベースライン血圧などを何個置いても、左辺は1本のままです。このように、1つのアウトカムを複数の説明変数で説明する回帰がmultivariable regressionです。この記事では「多変数回帰」と呼びます。
これに対して、同じ対象者から収縮期血圧と拡張期血圧を測定し、2つを同時にモデルへ入れると左辺がベクトルになります。
$$ \mathbf{y}_i = \mathbf{B}^\top \mathbf{x}_i + \boldsymbol{\varepsilon}_i, \qquad \boldsymbol{\varepsilon}_i \sim N(\mathbf{0}, \boldsymbol{\Sigma}) \tag{2} $$
式(2)の$\mathbf{y}_i$には複数のアウトカムが並び、$\mathbf{B}$には説明変数ごとの係数がアウトカム別に入ります。$\boldsymbol{\Sigma}$は残差の共分散行列で、対角成分がそれぞれのアウトカムの残差分散、非対角成分がアウトカム間の残差共分散です。
回帰に限って両者を見分けるなら、私なら最初に左辺を確認します。説明変数が何個あるかではなく、同時にモデル化しているアウトカムが1つなのか複数なのかを見る方が混乱しません。
600例では調整後の推定値が0.76mmHg動いた
今回用意した600例では、治療群が289例、対照群が311例でした。年齢が高いほど治療を受ける確率が低くなるように生成したため、治療群の平均年齢は56.8歳、対照群は62.5歳となり、両群には5.8歳の差があります。12週後の血圧に置いた真の治療効果は、収縮期血圧で-8mmHg、拡張期血圧で-4mmHgです。
まず収縮期血圧だけを見ます。年齢を調整せず、治療の有無だけを説明変数にした回帰では、治療効果は-9.46mmHgで、95%信頼区間は-11.07から-7.84mmHgでした。年齢とベースライン収縮期血圧を加えると、推定値は-8.70mmHg、95%信頼区間は-10.05から-7.34mmHgへ移動しています。
この0.76mmHgの変化は、多変数回帰そのものが推定値を0へ近づける性質を持つからではありません。今回のデータでは、治療群の方が若く、年齢が高いほど12週後の収縮期血圧も高くなるように生成しています。そのため、年齢を無視した解析では、治療群がもともと低い血圧になりやすいことまで治療効果へ混ざりました。年齢を調整すると、その分が取り除かれて真値の-8mmHgに近い推定値になっています。
交絡の向きが変われば、調整後の推定値が0から遠ざかることもあります。「調整すると効果が小さくなる」と読むのではなく、どの変数が治療とアウトカムの双方にどう関係しているかを見る必要があります。
2つを同時に扱っても各係数は動かなかった
次に、説明変数はそのままにして、12週後の収縮期血圧と拡張期血圧を同時に式(2)へ入れます。生成時には2つの誤差に0.6の相関を置き、当てはめたモデルの残差相関は0.580でした。同じ対象者から測定した2種類の血圧が連動している状況です。
ここで、2つのアウトカムを同時に扱えば、それぞれの治療効果も精密になると考えたくなるかもしれません。今回のように全対象者で両アウトカムが観測され、各アウトカムに同じ説明変数を置く通常の多変量線形モデルでは、そのような変化は起こりません。
治療係数は収縮期血圧で-8.70mmHg、拡張期血圧で-4.18mmHgでした。これは、それぞれのアウトカムへ別々に多変数回帰を当てはめた係数と一致します。収縮期血圧の治療係数の標準誤差も0.690mmHgのままでした。同じ説明変数行列を使う最小二乗推定では、各アウトカムの係数推定をアウトカムごとに分けて計算できるためです。
同時に扱うことで新しく得られるのは、アウトカムをまたいだ情報です。治療効果が2つの血圧について同時に0であるという帰無仮説を検定すると、Pillaiのトレースは0.2623、近似F値は105.78、自由度は2と595でした。F分布から求めたp値は約4.96×10-40です。
ただし、この非常に小さいp値から治療効果の臨床的な大きさは判断できません。ここで検定しているのは2つの治療係数をまとめた帰無仮説です。Pocockら（1987）は、複数のアウトカムを扱う臨床試験では、多重性とアウトカム間の相関を考える必要があることを論じ、複数の評価項目に対する統合的な検定も選択肢として扱っています。
効果差を比べると共分散が必要になる
多変量モデルの意味がはっきりするのは、収縮期血圧と拡張期血圧の治療効果を互いに比較するときです。今回の推定値なら、収縮期血圧の治療効果-8.70mmHgと拡張期血圧の-4.18mmHgとの差は、単純な引き算で-4.51mmHgになります。
点推定値だけなら引き算で足ります。難しいのは、その差の標準誤差です。
$$ \mathrm{Var}(\hat{\beta}_{1} – \hat{\beta}_{2}) = \left( \Sigma_{11} + \Sigma_{22} – 2\Sigma_{12} \right) \left[ (\mathbf{X}^\top \mathbf{X})^{-1} \right]_{tt} \tag{3} $$
式(3)の$\hat{\beta}_{1}$と$\hat{\beta}_{2}$は2つのアウトカムに対する同じ説明変数の係数です。$\Sigma_{11}$と$\Sigma_{22}$は各アウトカムの残差分散、$\Sigma_{12}$は残差共分散、$[(\mathbf{X}^\top\mathbf{X})^{-1}]_{tt}$は比較している説明変数に対応する成分を表します。
ここで効いてくるのが$-2\Sigma_{12}$です。2本の単変量回帰について通常の係数表だけを保存しても、このアウトカム間の共分散は含まれていません。同じ対象者の対応を保った残差から別途推定するか、多変量モデルとして共分散まで扱う必要があります。
今回の収縮期血圧の標準誤差は0.690mmHg、拡張期血圧は0.426mmHgでした。この2つの推定値を独立と仮定して二乗和の平方根を取ると、差の標準誤差は0.811mmHgになります。実際には残差に正の相関があるため、式(3)で共分散を含めると0.563mmHgまで小さくなりました。
その結果、治療効果の差は-4.51mmHg、95%信頼区間は-5.62から-3.41mmHg、t値は-8.02でした。独立と仮定した0.811mmHgは、共分散を含めた0.563mmHgの約1.44倍です。2つの推定値の差を評価するときに相関を無視すると、今回のデータでは区間を必要以上に広く見積もる方向へ働きました。標準誤差の作り方で信頼区間が変わるという点は、ロバスト標準誤差で扱った話ともつながります。
Teixeira-PintoとNormand（2009）は、相関する2つのアウトカムを扱う方法を論じており、アウトカム間の関連を無視できない研究目的では同時モデル化が必要になることを説明しています。逆に、今回と同じ条件で各アウトカムの係数と信頼区間を別々に報告するだけなら、同時モデルへ変えても各係数の推定精度は上がりません。私なら、統合検定や式(3)のようなアウトカム間の対比が研究目的に含まれているかを、多変量モデルを選ぶ基準にします。
報告ではmultivariableと調整変数を明記する
HidalgoとGoodman（2013）は、医学・公衆衛生分野でmultivariateとmultivariableが混同されていることを指摘し、両者を区別しています。Katz（2003）も、複数の説明変数を用いて1つのアウトカムを解析する方法をmultivariable analysisとして解説しています。
報告ガイドラインにも同じ語法が見られます。2015年のTRIPOD声明は、表題から「multivariable prediction model」という語を使っています。複数の予測因子から診断や予後に関するアウトカムを予測するモデルを扱うためです。
STROBE声明の解説文書でも、交絡因子を調整する回帰はmultivariable regressionと書かれています。また、STROBEの項目16では、粗推定値と交絡因子で調整した推定値、その精度を報告し、何を調整したのかを明らかにするよう求めています。今回なら、未調整の-9.46mmHgと調整後の-8.70mmHgを示し、後者では年齢とベースライン収縮期血圧を調整したと書く形です。
そのため、日本語論文で「多変量解析を行った」とだけ書く方法は適切ではありません。この一文からは、アウトカムが複数なのか、1つのアウトカムに複数の共変量を入れただけなのか判断できないからです。「12週後の収縮期血圧を目的変数とし、治療群、年齢、ベースライン収縮期血圧を説明変数とする多変数線形回帰を行った」のように、モデルの左右に何を置いたのかまで書く方が正確です。
Rで600例の結果を確かめる
ここまでの結果は、base Rとstatsで再現できます。まず、年齢によって治療を受ける確率が異なる600例を生成し、収縮期血圧と拡張期血圧には相関0.6の誤差を加えています。
set.seed(42)
n_total <- 600
age <- rnorm(n_total, mean = 60, sd = 10)
sbp_baseline <- rnorm(n_total, mean = 150, sd = 15)
treat_prob <- plogis(-0.06 * (age - 60))
treatment <- rbinom(n_total, size = 1, prob = treat_prob)
sd_sbp <- 8
sd_dbp <- 5
rho_error <- 0.6
sigma_true <- matrix(c(sd_sbp^2, rho_error * sd_sbp * sd_dbp,
rho_error * sd_sbp * sd_dbp, sd_dbp^2), nrow = 2)
errors <- matrix(rnorm(2 * n_total), ncol = 2) %*% chol(sigma_true)ここで用意した誤差には収縮期血圧と拡張期血圧の相関が入っています。続いて、真の治療効果を収縮期血圧-8mmHg、拡張期血圧-4mmHgとして12週後の値を生成します。
sbp_week12 <- 120 + 0.35 * (sbp_baseline - 150) +
0.25 * (age - 60) - 8 * treatment + errors[, 1]
dbp_week12 <- 75 + 0.20 * (sbp_baseline - 150) +
0.10 * (age - 60) - 4 * treatment + errors[, 2]
trial <- data.frame(treatment, age, sbp_baseline,
sbp_week12, dbp_week12)収縮期血圧だけをアウトカムにした未調整モデルと調整済みモデルでは、治療係数が次のように変わります。
fit_crude <- lm(sbp_week12 ~ treatment, data = trial)
fit_multivariable <- lm(sbp_week12 ~ treatment + age + sbp_baseline,
data = trial)
coef(fit_crude)["treatment"]
coef(fit_multivariable)["treatment"]#> treatment
#> -9.456380
#> treatment
#> -8.698122この出力で読むべきなのは、調整後の値が小さくなったという事実だけではありません。今回のデータ生成では、年齢による交絡を取り除いた結果、-9.46mmHgから真値-8mmHgに近い-8.70mmHgへ移動しています。交絡の構造が異なるデータで同じ方向へ動くとは限りません。
次に、2つの血圧を同時にアウトカムへ置きます。
fit_multivariate <- lm(cbind(sbp_week12, dbp_week12) ~
treatment + age + sbp_baseline, data = trial)
coef(fit_multivariate)["treatment", ]
resid_cov <- crossprod(residuals(fit_multivariate)) /
fit_multivariate$df.residual
resid_cov[1, 2] / sqrt(resid_cov[1, 1] * resid_cov[2, 2])
anova(fit_multivariate)#> sbp_week12 dbp_week12
#> -8.698122 -4.184924
#> [1] 0.5798800
#> Df Pillai approx F num Df den Df Pr(>F)
#> treatment 1 0.26229 105.776 2 595 < 2.2e-16収縮期血圧の治療係数-8.698122は、単独で当てはめた多変数回帰と一致しています。一方、Pillaiの検定は収縮期血圧だけの係数を調べているのではなく、2つの治療係数をまとめた仮説を扱っています。Rの表示はp<2.2e-16で打ち切られていますが、F値105.776、自由度2と595から尾確率を求めると約4.96×10-40です。
最後に、収縮期血圧と拡張期血圧の差そのものをアウトカムにすると、式(3)と同じ治療効果差を確認できます。
fit_difference <- lm(I(sbp_week12 - dbp_week12) ~
treatment + age + sbp_baseline, data = trial)
summary(fit_difference)$coefficients["treatment", ]
confint(fit_difference)["treatment", ]#> Estimate Std. Error t value Pr(>|t|)
#> -4.513199e+00 5.626183e-01 -8.021777e+00 5.548278e-15
#> 2.5 % 97.5 %
#> -5.618154 -3.408243標準誤差は0.5626mmHgで、2つの治療係数を独立と仮定した0.811mmHgより小さくなっています。今回は同じ説明変数を置いた線形モデルなので、2つのアウトカムの差を直接回帰する方法と、残差共分散を使って式(3)から求める方法が一致します。別のモデル構造でも常にこの置き換えができる、と一般化しない方がよいでしょう。
自分の研究では推定したい量から決める
12週後の収縮期血圧という1つのアウトカムについて、治療、年齢、ベースライン血圧を説明変数に入れたのであれば、英語ではmultivariable regressionです。日本語で「多変量解析」とだけ書くより、アウトカムと調整変数を文章で明記した方が解析内容は正確に伝わります。
複数のアウトカムを同時に扱う必要があるかどうかは、アウトカムの数だけでは決まりません。今回の条件では、収縮期血圧と拡張期血圧を同時にモデル化しても、それぞれの治療係数や標準誤差は変わりませんでした。統合検定を行いたい、2つの治療効果の差を推定したい、アウトカム間の共分散そのものを推論に使いたいといった研究目的があって、初めて同時モデル化の意味が出てきます。
実際の研究では、主要アウトカムをどれに定めているか、複数アウトカムに対する多重性をどう扱うか、アウトカム間の関連をどの推論に使うのかまで含めて解析方針を決めます。Dr.データサイエンスでは、こうした研究目的と解析モデルの対応についてご相談を承っています。
参考文献
- Hidalgo B, Goodman M. Multivariate or Multivariable Regression? American Journal of Public Health. 2013;103(1):39-40.
- Katz MH. Multivariable analysis: a primer for readers of medical research. Annals of Internal Medicine. 2003;138(8):644-650.
- Pocock SJ, Geller NL, Tsiatis AA. The analysis of multiple endpoints in clinical trials. Biometrics. 1987;43(3):487-498.
- Teixeira-Pinto A, Normand SL. Correlated bivariate continuous and binary outcomes: issues and applications. Statistics in Medicine. 2009;28(13):1753-1773.
- von Elm E, Altman DG, Egger M, Pocock SJ, Gøtzsche PC, Vandenbroucke JP. The Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) statement: guidelines for reporting observational studies. Journal of Clinical Epidemiology. 2008;61(4):344-349.
- Collins GS, Reitsma JB, Altman DG, Moons KGM. Transparent Reporting of a multivariable prediction model for Individual Prognosis Or Diagnosis (TRIPOD): the TRIPOD statement. Annals of Internal Medicine. 2015;162(1):55-63.
