行列を1行ずつリストに変換する
行列を1行ずつリストに変換したいという非常にニッチなニーズがあり、少し苦労したのでコードをメモ。
lapply(seq_len(nrow(mat)), function(i) mat[i,])
RMSTの算出
生存時間分析でRMSTを出したい。 こちらではsurvRM2パッケージを使うと良いとある。
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パッケージもまだメンテナンスされている。
私は知らなかったが、data.table::freadを内部的に用いるパッケージとしてbigreadrパッケージというものもあるようだ。
とりあえず備忘録としてこの記事を残しておく。
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"
rowSumsはrowwise+acrossで代替する
特定の条件のもとカウント回数を算出したいときは以下のようにやる。
iris %>% rowwise() %>% mutate(hoge = sum(across(Sepal.Length:Sepal.Width) >= 1))
ggplot2の描画データを取り出したいときはggplot_build()
geom_line()で折れ線プロットを描いた後、geom_smooth()して平滑化することなどが多々ある。
で、その平滑化した結果から例えばピーク検出したい場合もあったりするが、その際、描画結果のデータが欲しくなる。
もちろん、geom_smooth()で用いた平滑化の手法(gamとかloessとか)を使えばいいんだが面倒くさい。
そんなときはggplot_build()を用いる。
使い方については↓を参照。