先日、ある変数の2次元分布をヒートマップを使って可視化する方法をまとめました。
ヒートマップで変数の区切りを細かくしていくと、滑らかなグラデーションにすることができます。 見た目が綺麗になるのはいいんですが、色の差がわかりにくくなってしまい、どことどこが同じくらいなのか読み取れません。 こうなると等高線が欲しくなってきます。
モデルの当てはめ
使用するパッケージとデータを準備する
今回もgtsummaryパッケージのtrialデータを使います。一度ヒートマップの回でまとめた内容なので、説明は短めにします。
library(tidyverse) data(trial, package="gtsummary")
> trial %>% glimpse Rows: 200 Columns: 8 $ trt <chr> "Drug A", "Drug B", "Drug A", "Drug A", "Drug A", "Drug B", "D… $ age <dbl> 23, 9, 31, NA, 51, 39, 37, 32, 31, 34, 42, 63, 54, 21, 48, 71,… $ marker <dbl> 0.160, 1.107, 0.277, 2.067, 2.767, 0.613, 0.354, 1.739, 0.144,… $ stage <fct> T1, T2, T1, T3, T4, T4, T1, T1, T1, T3, T1, T3, T4, T4, T1, T4… $ grade <fct> II, I, II, III, III, I, II, I, II, I, III, I, III, I, I, III, … $ response <int> 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0,… $ death <int> 0, 0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, 0, 0, 0,… $ ttdeath <dbl> 24.00, 24.00, 24.00, 17.64, 16.43, 15.64, 24.00, 18.43, 24.00,…
ロジスティック回帰モデルを当てはめる
deathをアウトカム変数、ageとmarkerを説明変数としてロジスティック回帰モデルを当てはめます。
fit <- glm(death ~ age*marker, data = trial, family = binomial(link = "logit"))
作図用データセットの作成
仮想的に説明変数を作成する
ageは0歳から80歳まで1歳刻みに、markerは0.5から3.5まで0.1刻みのカテゴリーにします。
age_img <- seq(from=10, to=80, by=1) marker_img <- seq(from=0.5, to=3.5, by=0.1)
expand.grid関数を使って全ての組み合わせを含むデータフレームを作る
71×31=2201通りの組み合わせを手書きで書くのは無理なので、expand.grid( )を使います。tidyrパッケージにはexpand_grid( )という同じような関数があるようです。
変数名を変えて、あとでモデルから予測確率を計算できるようにしておきます。
data_img <- as.data.frame(expand.grid(age_img, marker_img)) %>% rename(age = Var1, marker = Var2)
predict関数を使って予測確率を計算する
data_img$pred <- predict(fit, data_img, type="response")
作図
geom_tile関数を使ってヒートマップを描く
前回同様、まずはgeom_tile( )でヒートマップを描いてみます。
q <- ggplot(data=data_img, aes(x=age, y=marker, fill=pred)) + geom_tile() print(q)
geom_contour関数を使って等高線を描く
先程描いたヒートマップに等高線を足してみます。
z軸方向の変数としてpredを指定し、線を描く高さをbreakで指定します。下のコードでは予測確率が0.5になる線*1を白で描きます。
q + geom_contour(aes(x=age, y=marker, z=pred), breaks = 0.5, colour = "white")
breakは複数の値を一度に指定することができますが、等高線のラベル（=標高）を振ってくれないので、どの線がどの高さなのか分からなくななってしまいます（解決方法勉強中）。
zの値が粗いと等高線も粗くなるので、滑らかに見える程度のカテゴリー数が必要です。
おわりに
- geom_contour( )のもう少し詳しい使い方や、geom_density_2d( )、geom_countour_filled( )など3次元データを2次元で表現するための関数についてはまた別の機会にまとめてみようと思います。
*1:予測確率が>0.5となる点と<0.5となる点を分ける線という言い方の方が正しい？
