NEJMに載るようなRCTでも再解析したら結果が異なるのですが??

R Rpackage Rを使いこなす 医学 数学 数学いらずの医科統計学 統計

面白い話題を見つけたので考えてみようと思ったらデータもあったのでやってみた。
ことの発端はここで
産科編 妊娠中期 妊娠中、虫歯の治療はできますか? (周産期医学 49巻13号) | 医書.jp
妊婦に歯科治療を行うと37週未満の早産が減らせるかどうかというRCTについて、
2006年に報告された（最初）NEJMの論文（OPT試験）では、早産が減らなかった（negative result）が、
Treatment of periodontal disease and the risk of preterm birth

対照群には自然流産や死産が多く、無作為化されていたとしてもフォロー終了直前では介入群よりリスクが低い患者たちが対照群に残ってしまうことで、実際の治療介入効果がないように見えてしまうバイアスがかかってしまうため、survivor average causal effect (SACE) もしくは principal strata effect という調整を行うと、ある程度のリスクがある患者では介入によって流産が減るかも、と2018年では言っている。
Periodontal treatment among mothers with mild to moderate periodontal disease and preterm birth: reanalysis of OPT trial data accounting for selective survival

OPT試験の結果を再現すると（出来てないのだが）、37週未満で早産したかのほかに、人工流産（Elective abortion）か死産（stillbirth）かというのがあって、これが介入群で5例、対照群で14例あるからこれがp=0.04で偏っているし、糖尿病もしくは高血圧がある患者が介入群で対照群の2倍くらいいてこれも偏っている、と。
そもそも無作為割付でこのように偏るのが正直言ってよくわからない。そして、「死産したあとの残りの対照群の患者たちは、介入群の患者と比べて早産のリスクが低い」と述べられていてこれも正直理解できていないが、興味の対象が早産ではあるがそもそも生産であるので、死産するような患者はもともと背景が高リスクであり、死産と早産はたいていリスクが共通している。となれば、死産した患者を「早産（生産）ではない」として除外して解析する以上、死産でたくさん除外される対象群は、高リスク群が除かれがち、ということのようである。
ただ、こういうバイアスがあるので、これを考慮して解析する必要があるらしい。
Y：37週未満の早産があった場合1、そうでない場合0
A：介入群なら1、対照群なら0
S：児が生存していれば1、死亡なら0
という基本の表記に対して、反事実（counterfactual）の変数として
Y_1：介入された母で早産が生じたとき、Y_0：介入されていない母で早産が生じたとき
S_1：介入された母で生産が生じたとき、S_1=1 なら介入された母で生産、S_1=0 なら介入された母で死産
S_0：介入されてない母で生産が生じたとき、S_0=1 なら介入されてない母で生産、S_0=0 なら介入されてない母で死産
実際のところ、S_1 は介入された（A=1）場合でしか観測できないし、S_0 は介入されてない（A=0）場合でしか観測できない

真の効果はE [Y|A =1, S=1]- E[ Y|A =0, S=1] で推定されるが、バイアスがあると正しくないのでSACEを以下で定義すると
E[Y_1-Y_0|S_1=S_0=1]
となるが、counterfactual な部分が含まれるので、その部分を\alpha（\propto と書いてあったり本文で突然\alpha になったりするがSACEの元論文を見る限り\alpha である）として適当にいじれるようにすると
E [Y|A =1, S=1]-E [Y|A =0, S=1]-\underbrace{E [Y_1|A =1, S=1]-E [Y_1|A =0, S=1]}_\alpha
となる。\alpha の設定というかSACEが機能する前提というのがあるがここでは省く。
\alpha だが、例えばいろいろな報告から早産に対するリスクの変化が1の前後で適当に増えたり減ったりする、とする。また、早産の一般的な頻度は12%であり、今回のOPT試験もほぼ同等である。いま、対照群に死産が多くて、結果として残った対照群患者は早産リスクにおいて介入群より低い人たちが残っている、という状況の時、介入群（E [Y_1|A =1, S=1]）は対象群（E [Y_1|A =0, S=1]=12%）に比べてリスクが1.33倍、と仮定すると、
1.33*12\%-12\%=4=\alpha となる。ここで論文の例をいじってわざと\alpha=4 にしたのは、結果からの天下りである。

SACEの解析にはrstanもやってみたが普通に確率分布の和を使って
E(X+Y)=E(X)+E(Y)
V(X+Y)=\frac{\hat{p}_1(1-\hat{p}_1)}{n_1}+\frac{\hat{p}_2(1-\hat{p}_2)}{n_2}
とすれば論文でいっている値と一致した。

論文では\alpha=4 でSACEの95%信頼区間上限が0を下回りますよ、とか言っている。


library(survival)
library(prodlim)
library(stringr)

load(url("https://www.causeweb.org/tshs/datasets/OPT_Study_PersonLevel_Data.RData"))
opt$Birth.outcome <- str_remove(opt$Birth.outcome, " +$")
# opt$censor <- str_detect(opt$Birth.outcome, "Lost")+0
opt$censor <- str_detect(opt$Birth.outcome, "Lost|Elective")+0

sfit <- prodlim(Hist(GA.at.outcome/7, censor) ~ Group, data=opt, reverse=TRUE)
spv <- survdiff(Surv(GA.at.outcome/7, censor) ~ Group, data=opt)

cols <- c(Control="red", Treatment="blue")
xl <- c(12, 38)
yl <- c(0, 0.14)
at.t <- seq(2, 38, by=2)
at.t <- at.t[at.t >= min(xl)]
axis2.at <- seq(0, 1, by=0.02)
axis2.labels <- sprintf("%s", axis2.at*100)

par(mar=c(3.5, 4, 2.5, 1.5), cex.lab=2, cex.axis=1.5, las=1)
plot(sfit, atrisk.labels=sprintf("%s group: ", names(cols)), xlim=xl, ylim=yl,
     xlab="Gestational Age [week]", ylab="Cumulative Percent of Patients", col=cols, atrisk.title="",
     axis1.at=at.t, axis1.labels=at.t,
     axis2.at=axis2.at, axis2.labels=axis2.labels,
     atrisk.at=at.t, atrisk.font=2, atrisk.cex=1.3,
     legend.x="topleft", legend.cex=1, legend=FALSE,
     confint=FALSE, marktime=TRUE, background=TRUE,
     background.horizontal=NA,
     logrank=TRUE, lwd=3, atrisk.line=0:1*1.2+1)
paxy <- par()$usr
legend("topleft", legend=names(cols), col=cols, pch=15, cex=2)
abline(v=37, lty=3)
lpv <- pchisq(spv$chisq, df=length(spv$n)-1, lower.tail=FALSE)
txt <- list(
  as.expression(substitute(atop(x1, phantom(0)), list(x1="log-rank test"))),
  as.expression(substitute(atop(phantom(0), italic(p)~"="~x2), list(x2=sprintf("%.2f", lpv))))
  )
for(j in seq(txt)) legend(paxy[1], paxy[3], legend=txt[[j]], bty="n", cex=2, yjust=-1)

# SACE
library(vioplot)
library(rstan)
rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())

code <- "
opta{
  int N;
  int B[N]; // preterm birth 1 
  int I[N]; // index T or C
}
parameters{
  real<lower=0, upper=1> p[2]; // preterm probability T and C
}
model{
  for(i in 1:N){
    if(I[i] == 1){
      B[i] ~ bernoulli(p[1]);
    } else {
      B[i] ~ bernoulli(p[2]);
    }
  }
}
"

m0 <- stan_model(model_code=code)
standata <- list(N=nrow(opt), B=(opt$GA.at.outcome/7 < 37 & opt$Birth.outcome!="Elective abortion" & opt$Birth.outcome!="Lost to FU")+0, I=ifelse(opt$Group=="T", 1, 0))

fit <- sampling(m0, stanopta, iter=25000, warmup=5000, chain=4)
ex <- extract(fit, pars=head(fit@model_pars, -1))
p <- ex$p

alpha <- 0.05
cia <- c(alpha/2, 0.5, 1-alpha/2)
rr <- quantile(-apply(p, 1, diff), cia)*100

# binomial distribution
m <- tapply(standata$B, standata$I, mean)
v <- sum(m*(1-m)/table(standata$I))
rr <- 100*(diff(m) + qnorm(cia)*sqrt(v))

prop <- seq(-2.5, 6, length=1000)
sace <- t(outer(rr, prop, "-"))
opt_alpha <- prop[which(apply(sace<0, 1, all))[1]]
yl <- c(-15, 10)

par(mar=c(5, 5, 2, 1), las=1, cex.axis=1.5, cex.lab=2)
matplot(prop, sace, type="l", lwd=3, col=1, lty=c(3,1,3), xlab="alpha", ylab="SACE", ylim=yl)
abline(h=0, lty=3)
abline(v=opt_alpha, lty=3)
text(opt_alpha, par()$usr[4], sprintf("%.2f", opt_alpha), pos=3, xpd=TRUE, cex=1.5)