RCTで多変量解析#2

http://r-statistics-fan.hatenablog.com/entry/2015/02/19/000248
前回の続き

RCTで多変量解析をするとPower検出力が上がる。
しかし、弱い予測因子を投入するとPowerは下がるはずだというのを
前回実証しようとして失敗した。

今回は、投入する因子を増やして、シミュレーション回数も増やしてみる。

結果、かえってPowerが下がる現象を証明し得た。

#回数増やしたので時間がかかる

library(tcltk)

##各変数の相関は今回は考えない
set.seed(1)
nsim <- 30000 #sim回数
nn <- 200   #各studyのn
b0 <- 0     #x0の真の係数 差がない場合
b0d <- 1    ##x0の真の係数 差がある場合
b1 <- 5     #大きな影響因子
b2b <- b2a <- b2 <- 0.01  #小さな影響因子
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)
      x2a <- rnorm(nn)
      x2b <- rnorm(nn)
      x3 <- rnorm(nn)
      x4 <- rnorm(nn)
      temp <- b1 * x1  + b2 * x2 +b2a * x2a + b2b * x2b + 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 + x2a + x2b))$coefficients[2,4]
      glm_p1_b2[i] <- summary(glm(y1 ~ x0 + x2 + x2a + x2b))$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

pow2 <- (glm_p1_b2 < 0.05)
pow1 <- (p1 < 0.05)
table(pow1, pow2)
mcnemar.test(pow1, pow2, correct = FALSE)
binom.test(367, 819)

> binom.test(367, 819)

Exact binomial test

data: 367 and 819
number of successes = 367, number of trials = 819, p-value = 0.003309
alternative hypothesis: true probability of success is not equal to 0.5

ということで、RCTで弱い予測因子３つを入れて多変量解析をすると
かえって有意に検出力が下がった。

> table(pow1, pow2)
pow2
pow1 FALSE TRUE
FALSE 23118 367
TRUE 452 6063

> binom.test(367, 819)

Exact binomial test

data: 367 and 819
number of successes = 367, number of trials = 819, p-value = 0.003309
alternative hypothesis: true probability of success is not equal to 0.5

> mcnemar.test(pow1, pow2, correct = FALSE)

McNemar's Chi-squared test

data: pow1 and pow2
McNemar's chi-squared = 8.8217, df = 1, p-value = 0.002977

と、いっても、かなり微妙。３万回でたった85回の差でしかない。