UPDATE: 2022-10-23 12:03:16
ここでは、複数の列が存在する場合にt検定を一括で行う方法をまとめておく。使用するパッケージは{tidyverse}と{broom}。
今回想定しているデータは下記のようなテーブルのデータ。segには0と1が入っている。この変数の0と1毎に各v*の値の平均に差があるかどうか、t検定を行いたい。v*は5列なので、今回であれば5行t.test()を書けば済む話だが、傾向スコアマッチングでベースラインの変数が揃ったかどうかを確認するときには、数十回書く必要があったりするので、なんとかしたい。
library(tidyverse)
library(broom)
# 有効桁数と乱数シードの制御
options(pillar.sigfig = 5)
set.seed(1989)
n <- 1000
df <- tibble(seg = sample(c(0,1), n, replace = TRUE),
v1 = rnorm(n, 100, 10) + if_else(seg == 1, 5, 0),
v2 = rnorm(n, 50, 1) + if_else(seg == 1, 5, 0),
v3 = rnorm(n, 200, 5) + if_else(seg == 1, 5, 0),
v4 = rnorm(n, 300, 20) + if_else(seg == 1, 5, 0),
v5 = rnorm(n, 50, 0.5) + if_else(seg == 1, 5, 0))
df## # A tibble: 1,000 × 6
## seg v1 v2 v3 v4 v5
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 88.938 55.294 207.66 327.78 55.020
## 2 0 100.07 48.158 204.41 279.01 49.620
## 3 1 105.98 55.046 207.84 276.06 56.454
## 4 0 114.41 49.544 196.71 287.03 49.906
## 5 0 84.048 49.560 200.85 322.68 50.022
## 6 0 114.60 49.727 206.86 290.33 50.540
## 7 0 102.72 50.094 203.11 270.24 49.980
## 8 1 105.64 55.884 207.95 289.18 55.308
## 9 0 102.91 51.457 205.94 285.13 49.134
## 10 0 103.88 50.974 199.89 263.93 49.470
## # … with 990 more rows
## # ℹ Use `print(n = ...)` to see more rows{tidyverse}と{broom}を組み合わせることで解決する。これは確かに楽ではあるが、おそらくRユーザー以外には著しく可読性が下がる。for-loopのほうが他の言語を使っている人からしてもイメージはしやすいの確かなのかな。
さておき、種明かしを始めるが、やっていることは単純で、データ構造をwide型からlong型に変換し、サブセットを作って、各サブセットにt検定を行ったのち、必要な情報をデータフレームに戻すということをやっている。
df %>%
gather(key = group, val = vals, -seg) %>%
group_by(group) %>%
nest() %>%
mutate(fit = map(data, ~ t.test(.$vals ~ .$seg, var.equal = FALSE)),
glanced = map(fit, glance)) %>%
unnest(glanced, .drop = TRUE) %>%
select(group, estimate1, estimate2, estimate,conf.low, conf.high, p.value) %>%
set_names(c("group", "seg0_mean", "seg1_mean","diff","conf_low", "conf_high", "pval")) %>%
mutate(sig = if_else(pval < 0.05, "sig", "not sig"))## # A tibble: 5 × 8
## # Groups: group [5]
## group seg0_mean seg1_mean diff conf_low conf_high pval sig
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 v1 99.930 105.58 -5.6499 -6.8700 -4.4299 5.3188e-19 sig
## 2 v2 49.931 54.959 -5.0285 -5.1523 -4.9047 0 sig
## 3 v3 200.54 205.11 -4.5743 -5.2022 -3.9464 2.6233e-42 sig
## 4 v4 299.84 304.93 -5.0856 -7.5052 -2.6660 4.0244e- 5 sig
## 5 v5 49.966 54.950 -4.9833 -5.0449 -4.9217 0 sig処理がミスってないかどうか確認しておく。愚直に5回書いてみたが、問題なさそうである。
# t.test(df$v1 ~ df$seg, var.equal = FALSE)
#
# Welch Two Sample t-test
#
# data: df$v1 by df$seg
# t = -9.0878, df = 997.22, p-value < 2.2e-16
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
# -6.869951 -4.429941
# sample estimates:
# mean in group 0 mean in group 1
# 99.92969 105.57964
#
# t.test(df$v2 ~ df$seg, var.equal = FALSE)
#
# Welch Two Sample t-test
#
# data: df$v2 by df$seg
# t = -79.717, df = 997.7, p-value < 2.2e-16
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
# -5.152287 -4.904719
# sample estimates:
# mean in group 0 mean in group 1
# 49.93081 54.95932
#
# t.test(df$v3 ~ df$seg, var.equal = FALSE)
#
# Welch Two Sample t-test
#
# data: df$v3 by df$seg
# t = -14.295, df = 996.98, p-value < 2.2e-16
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
# -5.202228 -3.946372
# sample estimates:
# mean in group 0 mean in group 1
# 200.5395 205.1138
#
# t.test(df$v4 ~ df$seg, var.equal = FALSE)
#
# Welch Two Sample t-test
#
# data: df$v4 by df$seg
# t = -4.1246, df = 994.62, p-value = 4.024e-05
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
# -7.505204 -2.666034
# sample estimates:
# mean in group 0 mean in group 1
# 299.8409 304.9265
#
# t.test(df$v5 ~ df$seg, var.equal = FALSE)
#
# Welch Two Sample t-test
#
# data: df$v5 by df$seg
# t = -158.77, df = 997.84, p-value < 2.2e-16
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
# -5.044925 -4.921737
# sample estimates:
# mean in group 0 mean in group 1
# 49.96640 54.94973