PCI手術は有効か？

The Lancet に掲載された論文 Percutaneous coronary intervention in stable angina (ORBITA): a double-blind, randomised controlled trial（安定狭心症における経皮的冠動脈インターベンション（ORBITA）: 二重盲検ランダム化比較試験）が注目を浴びている。同誌のコメント Last nail in the coffin for PCI in stable angina? や New York Times の記事 ‘Unbelievable’: Heart Stents Fail to Ease Chest Pain（「信じられない」：心臓のステントは胸の痛みを和らげない），日本語の解説 愛し野塾　第144回　経皮冠動脈インターベンションのターニングポイント も参照されたい。

狭心症の治療のためにカテーテルを挿入してステントを留置する経皮的冠動脈インターベンション（PCI）という手術が全世界で毎年50万回行われている。その有効性を調べるため，200人をランダムにほぼ半々に分け，一方（PCI群）には本物のPCI手術を，他方（プラセボ群）には「偽の手術」を行い，効果を比較した。これだけの規模で偽の手術を行うのは異例である。

このような試験では，p ハッキングを避けるため，どういう項目を調べるかを，あらかじめ登録しておく。この ORBITA 試験は ClinicalTrials.gov に Objective Randomised Blinded Investigation With Optimal Medical Therapy of Angioplasty in Stable Angina (ORBITA) として登録されている。それによれば，主要評価項目（the primary endpoint）はPCI群・プラセボ群のトレッドミルでの運動時間の増加の差である。これ以外にも，副次的評価項目（secondary endpoints）としていろいろ調べるが，たくさん調べて一つでも統計的に有意（p < 0.05）な項目があれば，事後的に最初からそれを調べたかったのだということにして「有意な差が見られた」と報告するのが p ハッキングである。こういうことを避けるため，事前に何を主要評価項目とするかを公的な機関に登録しておくのが正しい流儀である。

結果は次の通りであった（時間の単位は秒）：

さて，この論文の統計的手法に，有名な統計学者 Gelman 先生たちが，いちゃもんをつけている：ORBITA and coronary stents: A case study in the analysis and reporting of clinical trials。ポイントは，統計的に有意な差がないことと効果がないこととは違うだろうということだが，その中で Gelman 先生たちはこの p 値を 0.200 から 0.09 に「改良」する方法を示している。おもしろいので，試してみた。

生データは与えられていないので，上の結果に合致するデータを作ってみる。四捨五入の関係でPCI群の増加が 28.4 にならないので，処置後 556.3 は 556.4 と改変した：

set.seed(180210)
d = 0.5423874  # 後述
x1 = rnorm(104)
y1 = x1 + rnorm(104, 0, d)
x2 = rnorm(90)
y2 = x2 + rnorm(90, 0, d)
x1 = (x1 - mean(x1)) / sd(x1) * 178.7 + 528.0
x2 = (x2 - mean(x2)) / sd(x2) * 195.0 + 490.0
y1 = (y1 - mean(y1)) / sd(y1) * 178.7 + 556.4
y2 = (y2 - mean(y2)) / sd(y2) * 190.9 + 501.8

d の値は p 値がちょうど 0.200 になるように決めた：

f = function(d) {
    set.seed(180210)
    x1 = rnorm(104)
    y1 = x1 + rnorm(104, 0, d)
    x2 = rnorm(90)
    y2 = x2 + rnorm(90, 0, d)
    x1 = (x1 - mean(x1)) / sd(x1) * 178.7 + 528.0
    x2 = (x2 - mean(x2)) / sd(x2) * 195.0 + 490.0
    y1 = (y1 - mean(y1)) / sd(y1) * 178.7 + 556.4
    y2 = (y2 - mean(y2)) / sd(y2) * 190.9 + 501.8
    t.test(y1-x1, y2-x2)$p.value - 0.2
}
uniroot(f, c(0,1))

これで t.test(y1-x1, y2-x2) を計算してみれば，確かに p 値は 0.2 になる：

t.test(y1-x1, y2-x2)

	Welch Two Sample t-test

data:  y1 - x1 and y2 - x2
t = 1.286, df = 189.46, p-value = 0.2
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -8.861419 42.061419
sample estimates:
mean of x mean of y 
     28.4      11.8 

さて，Gelman 先生たちは，この p 値を 0.09 に改良する。おそらく次のようにやったのではないかと思われる：

x = c(x1, x2)
y = c(y1, y2)
z = c(rep(1,104), rep(0,90))
r = lm(y ~ x + z)

結果は次のようになる：

summary(r)

Call:
lm(formula = y ~ x + z)

Residuals:
     Min       1Q   Median       3Q      Max 
-238.326  -58.680   -3.507   57.322  231.587 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 73.87740   18.86142   3.917 0.000125 ***
x            0.87331    0.03365  25.956  < 2e-16 ***
z           21.41417   12.57851   1.702 0.090300 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 86.92 on 191 degrees of freedom
Multiple R-squared:  0.7839,	Adjusted R-squared:  0.7816 
F-statistic: 346.4 on 2 and 191 DF,  p-value: < 2.2e-16

確かに z の係数の p 値は 0.09 になる。そして，こちらのほうがもっともらしい。

では後者のほうが検出力が大きい（p 値が小さくなる傾向がある）といえるか？ 同様な問題を何度も解いて試してみる：

f = function() {
    x1 = rnorm(100)
    y1 = (x1 + rnorm(100, 0, 0.5)) / sqrt(1.25) + 0.05
    x2 = rnorm(100)
    y2 = (x2 + rnorm(100, 0, 0.5)) / sqrt(1.25) + 0.15
    p1 = t.test(y1-x1, y2-x2)$p.value  # 前者のp値
    x = c(x1, x2)
    y = c(y1, y2)
    z = c(rep(1,100), rep(0,100))
    r = lm(y ~ x + z)
    p2 = summary(r)$coefficients[3,4]  # 後者のp値
    c(p1, p2)
}
a = replicate(1000, f())
median(a[1,])
median(a[2,])
mean(a[1,] > a[2,])

どちらかといえば後者の p 値が小さい傾向が見られるが，顕著ではない。

今回はたまたま処置群の処置前の値が高い方にゆらいでいたので，それが「平均への回帰」で低い方に動くことで，処置群の変化が大きくなり，p 値が減ったのであろう。