ランダム化比較試験

View slides in full screen

1 多重比較

回帰分析から得られた推定値のp値を比較回数でかける

pairwise.t.test(y, treatment, p.adjust.method = "bonferroni")

もしtreatmentが5グループなら比較回数は10通りであり、回帰分析のp値はpairwise.t.test(y, treatment, p.adjust.method = "bonferroni")の10分の1

mc_df <- tibble(treatment = rep(LETTERS[1:5], each = 100),
                y         = rep(1:100, 5))

mc_df |> 
  reframe(mean   = mean(y),
          median = median(y),
          sd     = sd(y),
          min    = min(y),
          max    = max(y),
          .by    = treatment)

mc_sim <- function(n = 100) {
  temp_df <- mc_df |> 
    slice_sample(n = n)
  
  temp_t <- pairwise.t.test(temp_df$y, temp_df$treatment, p.adjust.method = "none")
  temp_p <- temp_t$p.value
  temp_p <- as.vector(temp_p)
  
  sum(temp_p < 0.05, na.rm = TRUE)
}

mc_result <- replicate(10000, mc_sim())

table(mc_result)

mc_result |> sum()

mc_sim_bf <- function(n = 100) {
  temp_df <- mc_df |> 
    slice_sample(n = n)
  
  temp_t <- pairwise.t.test(temp_df$y, temp_df$treatment, p.adjust.method = "bonferroni")
  temp_p <- temp_t$p.value
  temp_p <- as.vector(temp_p)
  
  sum(temp_p < 0.05, na.rm = TRUE)
}

mc_result_bf <- replicate(10000, mc_sim_bf())

table(mc_result_bf)

mc_result_bf |> sum()

1. 母集団pからサンプルsを抽出
2. sを3つのグループに分割（a、b、c）
3. グループ間の差は理論上0
4. グループ間の差を回帰分析で推定（ベースラインはa）
5. a vs. b、a vs. cの計2回の比較が行われる
6. 以上の作業を複数回繰り返す
7. 一つでもp値が0.05を下回るケースは全施行の中、約1割
8. 閾値を比較回数である2で割ったら0.25となり、0.25を基準とすると約0.5割が統計的有意に

set.seed(19861008)
p <- rnorm(10000, 0, 1)
result <- data.frame(id = NA, b = NA, c = NA)

for (i in 1:10000) {
  s <- sample(p, 1000)
  names(s) <- sample(letters[1:3], 1000, replace = TRUE)
  pvalue <- enframe(s) |> 
    lm(value ~ name, data = _) |> 
    broom::tidy() |> 
    filter(term != "(Intercept)") |> 
    pull(p.value)
  
  result[i, ] <- c(1, pvalue)
  
}

table(result$b < 0.05)
table(result$c < 0.05)

result |> 
  mutate(b2 = if_else(b < 0.05, 1, 0),
         c2 = if_else(c < 0.05, 1, 0)) |> 
  mutate(p = if_else(b2 + c2 >= 1, 1, 0)) |> 
  pull(p) |> table()
  
result |> 
  mutate(b2 = if_else(b < 0.05 / 2, 1, 0),
         c2 = if_else(c < 0.05 / 2, 1, 0)) |> 
  mutate(p = if_else(b2 + c2 >= 1, 1, 0)) |> 
  pull(p) |> table()