サブグループ解析とは、全体からある特性を持った一部の集団（例：男性）を取り出して解析をすること。ランダム化が崩れたり、解析対象集団のサイズが小さくなるから検出力が落ちたり、多重検定の問題が生じたり、と色々注意はある。
でもこの記事では、そこらへんの問題は置いといて、「対象集団抽出 → 解析 → 結果表示」という単純だが面倒なプロセスを繰り返さないといけなくなるので、tidyverse系のpurrrパッケージで楽に行えるようになりたい。
パッケージとデータセットの準備
gtsummaryパッケージに含まれるtrialというシミュレーションデータを使う。tidyverseパッケージに加え、簡単・キレイに結果を表示するためのbroomパッケージを入れておく（これまでも下の記事などで、少し使ったことあり）。
> library(tidyverse) > library(broom) > data(trial, package = "gtsummary") > glimpse(trial) Rows: 200 Columns: 8 $ trt <chr> "Drug A", "Drug B", "D… $ age <dbl> 23, 9, 31, NA, 51, 39,… $ marker <dbl> 0.160, 1.107, 0.277, 2… $ stage <fct> T1, T2, T1, T3, T4, T4… $ grade <fct> II, I, II, III, III, I… $ response <int> 0, 1, 0, 1, 1, 0, 0, 0… $ death <int> 0, 0, 0, 1, 1, 1, 0, 1… $ ttdeath <dbl> 24.00, 24.00, 24.00, 1…
trt：化学療法age：年齢marker：腫瘍マーカーstage：Tステージgrade：グレードresponse：腫瘍の反応性death：死亡death_cr：腫瘍関連死亡、腫瘍以外による死亡、打ち切りttdeath：死亡または打ち切りまでの時間
ネストしたデータ
「ネストした（nested）」は「入れ子構造になった」ということ。Rで慣れ親しんだデータ集合の形式は「data frame型」と言われるもので、各セルには、列毎に同じデータ形式（文字列なら文字列）の要素が1つだけ入るようになっている。
これに対してtidyverse系で推奨されるデータ集合形式は「tibble型」と言われるもので、各セルの要素としてデータの集合体（データフレームや解析結果など）を取ることができる。「要素としてデータの集合体を含んだデータの集合体」なので、nested dataと呼ばれる。
tibble型はdata frame型の上位互換で、要素に普通の数字・文字列が入るときは違いを意識しなくても困らない。
group_by( )とnest( )を使ってネストデータを作成する
例えばtrialデータにおいて、stageに基づいてデータを分けてサブグループ解析を行いたいときは、次のようにstageでグループ化してnest( )に渡す。これはよく使われるので、group_nest( )というまとめて行える関数も用意されている。
nested_data <- trial %>% group_by(stage) %>% nest()
> nested_data # A tibble: 4 × 2 # Groups: stage [4] stage data <fct> <list> 1 T1 <tibble [53 × 7]> 2 T2 <tibble [54 × 7]> 3 T3 <tibble [43 × 7]> 4 T4 <tibble [50 × 7]>
dataという変数の中に、43-53症例分の7つの変数を持ったtibble型データが格納されている。dataを開いてみると下のような感じ。
> nested_data$data [[1]] # A tibble: 53 × 7 trt age marker grade response death ttdeath <chr> <dbl> <dbl> <fct> <int> <int> <dbl> 1 Drug A 23 0.16 II 0 0 24 2 Drug A 31 0.277 II 0 0 24 3 Drug A 37 0.354 II 0 0 24 4 Drug A 32 1.74 I 0 1 18.4 5 Drug A 31 0.144 II 0 0 24 ~~~~~（以下省略）~~~~~
同じ処理を一括で行う
先程作成したstage別の4つのデータそれぞれに対して、死亡deathをアウトカム、治療trtと年齢ageを説明変数としたロジスティック回帰を行ってみる。最終的にはtrtの効果（=調整OR）を取り出したい。
function( )で処理関数を書く
ここでいう「処理関数」とは、データを引数として受け取って、前述のロジスティック回帰を行い、結果を出力する関数。次のように定義する。
fun <- function(d){ fit <- glm(death ~ trt + age, family = binomial(link="logit"), data = d) return(fit) }
map( )で各データに処理を適用する
map(引数, 処理関数)の要領で、さっき定義した関数を各データ（dataの要素）に適用する。結果を、fitという変数名で元のtibbleに追加する。処理関数が2つ、あるはそれ以上の引数を必要とする場合は、map( )の代わりにmap2( ), pmap( )を使う。
nested_data %>% mutate(fit = map(data, fun))
この段階で以下が出力されたはず。
# A tibble: 4 × 3 stage data fit <fct> <list<tibble[,7]>> <list> 1 T1 [53 × 7] <glm> 2 T2 [54 × 7] <glm> 3 T3 [43 × 7] <glm> 4 T4 [50 × 7] <glm>
右端にfitという名前でglm( )の結果が格納されている。
ついでに、broomパッケージのtidy( )を使って回帰の結果をキレイにした形で保存するコードも追加する（estという名前を付けた）。
nested_data %>% mutate(fit = map(data,fun), est = map(fit, tidy, conf.int = TRUE, exponentiate = TRUE))
tidy( )で指定する条件（conf.int = TRUEやexponatiate = TRUE）は、関数名の後に書く。tidy( )のカッコの中に書くとエラーになる。
unnest( )で入れ子構造を解く
上記のままだと結果の中身がわからないので、unnest( )で入れ子構造を解除する。解除したい列をcolsで指定する。
nested_data %>% mutate(fit = map(data,fun), est = map(fit, tidy, conf.int = TRUE, exponentiate = TRUE)) %>% unnest(est)
# A tibble: 12 × 10 # Groups: stage [4] stage data fit term estimate std.error statistic p.value conf.low conf.high <fct> <list> <list> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> 1 T1 <tibble> <glm> (Inte… 0.122 1.05 -2.00 0.0454 0.0131 0.869 2 T1 <tibble> <glm> trtDr… 1.24 0.592 0.369 0.712 0.385 4.01 3 T1 <tibble> <glm> age 1.04 0.0213 1.71 0.0865 0.996 1.08 4 T2 <tibble> <glm> (Inte… 0.235 1.23 -1.17 0.240 0.0185 2.50 5 T2 <tibble> <glm> trtDr… 1.61 0.569 0.839 0.402 0.533 5.05
変数毎（Intercept, trtDrug B, age）の結果が縦長に繋がっているので、trtの結果のみ取り出しておく。
nested_data %>% mutate(fit = map(data,fun), est = map(fit, tidy, conf.int = TRUE, exponentiate = TRUE)) %>% unnest(est) %>% filter(term == "trtDrug B")
# A tibble: 4 × 10 # Groups: stage [4] stage data fit term estimate std.error statistic p.value conf.low conf.high <fct> <list> <list> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> 1 T1 <tibble> <glm> trtDru… 1.24 0.592 0.369 0.712 0.385 4.01 2 T2 <tibble> <glm> trtDru… 1.61 0.569 0.839 0.402 0.533 5.05 3 T3 <tibble> <glm> trtDru… 0.623 0.639 -0.740 0.459 0.174 2.17 4 T4 <tibble> <glm> trtDru… 3.46 0.783 1.59 0.112 0.798 18.7
見栄えはまだイマイチだが、概ね欲しい結果は取り出せた。
おわりに
- 自分の確認のため、一つ一つ分割した説明になってしまいました。
- 以下のことはもう少し勉強してからまとめてみたいです（これだけあって、今回は何をまとめたんだろう...）
- チルダ式
- 引数が複数の場合
- unnest( )の挙動
- purrrパッケージのリスト操作関数
- 私は「猫が膝に乗る → 気づくと寝ている → 何も進んでない」の繰り返しです。
参考文献
- 「RStudio 導入」でググると上位に出てくる詳しい解説スライドの作成者、矢内勇生先生が宋財泫 (Jaehyun Song)先生と執筆されているeBookです。私のようにデータサイエンスに縁遠い者としては、「なぜそうなるのか？」みたいな細かい説明があって、非常に重宝しています。
