RMSTの算出

生存時間分析でRMSTを出したい。 こちらではsurvRM2パッケージを使うと良いとある。

nshi.jp

survRM2パッケージのvignetteはこちら。 https://cran.r-project.org/web/packages/survRM2/vignettes/survRM2-vignette3-2.html

こんな感じでRMSTを算出してくれる。

> smp <- pbc %>% filter(status %in% c(0,2), !is.na(trt)) %>% 
+   mutate(status = if_else(status == 2, 1, 0),
+          trt = trt-1 # 後述のrmst2()のamr引数が0,1しか受け付けないため変換
+          ) # サンプルデータ

>library(survRM2)
> with(smp, rmst2(time, status, trt, tau = 3650))

The truncation time: tau = 3650  was specified. 

Restricted Mean Survival Time (RMST) by arm 
                 Est.      se lower .95 upper .95
RMST (arm=1) 2614.267 111.805  2395.133  2833.401
RMST (arm=0) 2574.660 105.944  2367.014  2782.305


Restricted Mean Time Lost (RMTL) by arm 
                 Est.      se lower .95 upper .95
RMTL (arm=1) 1035.733 111.805   816.599  1254.867
RMTL (arm=0) 1075.340 105.944   867.695  1282.986


Between-group contrast 
                       Est. lower .95 upper .95     p
RMST (arm=1)-(arm=0) 39.608  -262.281   341.496 0.797
RMST (arm=1)/(arm=0)  1.015     0.904     1.141 0.797
RMTL (arm=1)/(arm=0)  0.963     0.723     1.283 0.797

しかし、adjustされていないシンプルなRMSTであれば、survfitの結果をprint()すれば表示してくれる。

> fit <- survfit( Surv(time, status) ~ trt, data =  smp)
> print(fit, print.rmean=TRUE, digits = 5, rmean = 3650) 
Call: survfit(formula = Surv(time, status) ~ trt, data = smp)

        n events *rmean *se(rmean) median 0.95LCL 0.95UCL
trt=0 148     65 2574.7     105.94   3222    2540      NA
trt=1 145     60 2614.3     111.81   3395    3090      NA
    * restricted mean with upper limit =  3650 

IPW+Cox回帰

ATEとしてハザード比を求める必要があったので、IPWを使った。 ダミーデータで備忘録を残しておく。

library(survminer)
library(survival)
library(cobalt)
library(tidyverse)

# データの準備
diabetic_mod <- diabetic %>% 
  filter(eye == "left") %>% 
  mutate(flag_highrisk = if_else(risk >= 10, 1, 0))

# KM曲線の描画
surv <- survfit(Surv(time, status) ~ trt, data = diabetic_mod)
ggsurvplot(surv)

# Cox回帰
fit <- coxph(Surv(time, status) ~ trt, data = diabetic_mod)
summary(fit)

# IPW + Cox回帰
# weightsにIPWの重みを与える
library(WeightIt)
model_ipw <- weightit(trt ~ flag_highrisk,
                  data = diabetic_mod,
                  estimand = "ATE", 
                  method = "ps"
)

fit_ipw <- coxph(Surv(time, status) ~ trt, 
                 data = diabetic_mod, 
                 weights = get.w(model_ipw))
summary(fit_ipw)

まあ、専門外の人に説明するときは傾向スコアマッチングの方が説明しやすいんだけど、それだとATTとなってしまうので...

重み付けCox回帰だとcoxphwパッケージを用いる方が適切かもしれないが、今回はナイーブにcoxph()で実施している。 結果はcoxphw()のtemplate="PH"とした場合とほぼ一致する。 stats.stackexchange.com

メモリに乗らないデータを読み込むパッケージの備忘録

メモリに乗らないデータをRでどう処理するかという問題は定期的に話題になる。

そのたびに「DBにつっこんでSQLで処理するじゃろjk」「せっかくだから俺はawkを使うぜ」「私のメモリは64GBです」などと喧々諤々しており既視感にとらわれるのだが、我々が古い道具でさばいている間にも世界は進歩しており、この問題に対応するパッケージなどが日々開発されている。

 

ということでざっと調べた結果を記録しようと思ったら、私自身も5年前にこういう記事を書いていた。何も覚えていない。ここで言及していたreadrパッケージもchunkedパッケージもまだメンテナンスされている。

dichika.hateblo.jp

 

私は知らなかったが、data.table::freadを内部的に用いるパッケージとしてbigreadrパッケージというものもあるようだ。

github.com

 

とりあえず備忘録としてこの記事を残しておく。

timezoneだけを変換する場合はforce_tzで時刻も変換する場合はwith_tz

Rのlubridateパッケージにおけるtimezone変換をいつも忘れてヘルプを見ているのでメモ。 timezoneに合わせて時刻も変換したい場合はwith_tz、 timezoneだけを置き換えたい場合はforce_tz 。

library(lubridate)

> (tmp <- Sys.time())
[1] "2021-04-23 10:54:58.933 JST"

> with_tz(tmp, "UTC") # timezoneも時刻も変換
[1] "2021-04-23 01:54:58.933 UTC"

> force_tz(tmp, "UTC") # timezoneのみ変換
[1] "2021-04-23 10:54:58.933 UTC"

ggplot2の描画データを取り出したいときはggplot_build()

geom_line()で折れ線プロットを描いた後、geom_smooth()して平滑化することなどが多々ある。

で、その平滑化した結果から例えばピーク検出したい場合もあったりするが、その際、描画結果のデータが欲しくなる。

もちろん、geom_smooth()で用いた平滑化の手法(gamとかloessとか)を使えばいいんだが面倒くさい。

そんなときはggplot_build()を用いる。

使い方については↓を参照。

stackoverflow.com