RCTで多変量解析

某社のwebカンファを見てきた。

RCTで多変量解析するのは、by-chanceで出た差を調整して
確認するためなどという感じのことを言っていた。

しかし、自分の認識は異なる。
RCTで多変量解析するのは、原則的にあらかじめpowerの増加を
期待して計画された時だけと認識している。
差がでたら困る変数は、層別ランダム化しておけって話。
それでもなお生じる差を更に調整する感じだ。

ここで、多変量解析に投入される因子はアウトカムに多大に影響する
既知の因子である。

しかし、あまり関係ない因子を投入すると、かえってpowerは
下がるのではないか？
ここのところを確認したかったので、シミュレーションを行う

#RCTで多変量解析

library(mvtnorm)
library(tcltk)

##各変数の相関は今回は考えない
set.seed(1)
nsim <- 3000 #sim回数
nn <- 100   #各studyのn
b0 <- 0     #x0の真の係数 差がない場合
b0d <- 1    ##x0の真の係数 差がある場合
b1 <- 5     #大きな影響因子
b2 <- 0.02  #小さな影響因子
b3 <- 3     #未測定因子1
b4 <- 2     #未測定因子2

glm_p0_b2 <- glm_p1_b2 <- glm_p1_b1 <- glm_p0_b1 <- p0 <- p1 <- numeric(nsim)

pb <- txtProgressBar(min = 1, max = nsim, style = 3)

for (i in seq_len(nsim)){
      setTxtProgressBar(pb, i)
      x0 <- sample.int(2, nn, replace = TRUE) - 1
      x1 <- rnorm(nn)
      x2 <- rnorm(nn)
      x3 <- rnorm(nn)
      x4 <- rnorm(nn)
      temp <- b1 * x1  + b2 * x2 + b3 * x3 + b4 + x4 + 0.1 * rnorm(1)
      y1 <- b0d * x0 + temp
      y0 <- b0 * x0 + temp
      p0[i] <- t.test(y0 ~ x0, alternative = "two.sided", var.equal = FALSE)$p.value
      p1[i] <- t.test(y1 ~ x0, alternative = "two.sided", var.equal = FALSE)$p.value
      glm_p0_b1[i] <- summary(glm(y0 ~ x0 + x1))$coefficients[2,4]
      glm_p1_b1[i] <- summary(glm(y1 ~ x0 + x1))$coefficients[2,4]
      glm_p0_b2[i] <- summary(glm(y0 ~ x0 + x2))$coefficients[2,4]
      glm_p1_b2[i] <- summary(glm(y1 ~ x0 + x2))$coefficients[2,4]
}

hist(p0)
hist(p1)
hist(glm_p0_b1)
hist(glm_p1_b1)
hist(glm_p0_b2)
hist(glm_p1_b2)

plot1 <- function(xx){
plot(sort(xx), ylim = c(0,1), pch=20, type = "p")
abline(a=0, b=1/nsim, col = "red")
}

plot1(p0)
plot1(p1)
plot1(glm_p0_b1)
plot1(glm_p1_b1)
plot1(glm_p0_b2)
plot1(glm_p1_b2)

sum(p0 < 0.05)/nsim
sum(p1 < 0.05)/nsim
sum(glm_p0_b1 < 0.05)/nsim
sum(glm_p1_b1 < 0.05)/nsim
sum(glm_p0_b2 < 0.05)/nsim
sum(glm_p1_b2 < 0.05)/nsim

<||

>||
> sum(p0 < 0.05)/nsim
[1] 0.05233333   ##RCTのt検定　アルファエラーは0.05に保たれている
> sum(p1 < 0.05)/nsim
[1] 0.1336667　　##RCTのt検定　Powerは0.13しかない
> sum(glm_p0_b1 < 0.05)/nsim
[1] 0.05033333　##RCTの多変量解析　アルファエラーは0.05に保たれている
> sum(glm_p1_b1 < 0.05)/nsim
[1] 0.342　　　　##RCTの多変量解析　強い予測因子で調整するとPowerは0.34にUP!
> sum(glm_p0_b2 < 0.05)/nsim
[1] 0.05233333　　##RCTの多変量解析　アルファエラーは0.05に保たれている
> sum(glm_p1_b2 < 0.05)/nsim
[1] 0.1333333　　##RCTの多変量解析　弱い予測因子で調整するとPowerは0.13で変わらない。
予測ではかえってPower下がると思ったけど下がらなかったわ。
||

ということで、今回のRCTのシミュレーションでは

モデルが正しくないとしても、既知の強力な予測因子で多変量解析で
調整すると、単純なRCTよりも、アルファエラーを損なわずにpowerを
高めることが出来る。

弱い予測因子を入れるとかえってPowerが悪化すると思ったけど、今回は
下がらなかった。１個だけじゃ、あまり影響がないのか。多分、沢山弱い因子を
入れたらPowerが下がるんじゃないかな。今日はもう眠いのでやらない。