dichika.hateblo.jp Open in urlscan Pro
35.75.255.9  Public Scan

Form analysis
1 forms found in the DOM

GET https://dichika.hateblo.jp/search

<form class="search-form" role="search" action="https://dichika.hateblo.jp/search" method="get">
  <input type="text" name="q" class="search-module-input" value="" placeholder="記事を検索" required="">
  <input type="submit" value="検索" class="search-module-button">
</form>

Text Content

読者になる



盆栽日記

この広告は、90日以上更新していないブログに表示しています。

2021-10-25


行列を１行ずつリストに変換する

R

行列を１行ずつリストに変換したいという非常にニッチなニーズがあり、少し苦労したのでコードをメモ。

lapply(seq_len(nrow(mat)), function(i) mat[i,]) 


dichika (id:dichika) 2年前




広告を非表示にする

 * もっと読む

コメントを書く
2021-07-09


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 


dichika (id:dichika) 2年前




広告を非表示にする

 * もっと読む

コメントを書く
2021-06-30


IPW+COX回帰

R

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

dichika (id:dichika) 2年前




広告を非表示にする

 * もっと読む

コメントを書く
2021-06-04


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

R

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

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

 

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

dichika.hateblo.jp

 

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

github.com

 

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

dichika (id:dichika) 2年前




広告を非表示にする

 * もっと読む

コメントを書く
2021-04-23


TIMEZONEだけを変換する場合はFORCE_TZで時刻も変換する場合はWITH_TZ

R

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"


dichika (id:dichika) 2年前




広告を非表示にする

 * もっと読む

コメントを書く
2021-03-12


ROWSUMSはROWWISE＋ACROSSで代替する

特定の条件のもとカウント回数を算出したいときは以下のようにやる。

iris %>% 
rowwise() %>% 
mutate(hoge = sum(across(Sepal.Length:Sepal.Width) >= 1)) 


dichika (id:dichika) 2年前




広告を非表示にする

 * もっと読む

コメントを書く
2020-10-20


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

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

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

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

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

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

stackoverflow.com

 

dichika (id:dichika) 3年前




広告を非表示にする

 * もっと読む

コメントを書く
次のページ

プロフィール
dichika (id:dichika)
読者です 読者をやめる 読者になる 読者になる
29
このブログについて
検索

リンク
 * はてなブログ
 * ブログをはじめる
 * 週刊はてなブログ
 * はてなブログPro

最新記事
 * 行列を１行ずつリストに変換する
 * RMSTの算出
 * IPW+Cox回帰
 * メモリに乗らないデータを読み込むパッケージの備忘録
 * timezoneだけを変換する場合はforce_tzで時刻も変換する場合はwith_tz

月別アーカイブ
 * ▼ ▶
   2021 (6)
   * 2021 / 10 (1)
   * 2021 / 7 (1)
   * 2021 / 6 (2)
   * 2021 / 4 (1)
   * 2021 / 3 (1)
 * ▼ ▶
   2020 (2)
   * 2020 / 10 (1)
   * 2020 / 9 (1)
 * ▼ ▶
   2019 (18)
   * 2019 / 11 (4)
   * 2019 / 10 (2)
   * 2019 / 7 (1)
   * 2019 / 6 (1)
   * 2019 / 4 (1)
   * 2019 / 2 (8)
   * 2019 / 1 (1)
 * ▼ ▶
   2018 (6)
   * 2018 / 12 (1)
   * 2018 / 7 (1)
   * 2018 / 6 (2)
   * 2018 / 2 (1)
   * 2018 / 1 (1)
 * ▼ ▶
   2017 (15)
   * 2017 / 9 (1)
   * 2017 / 8 (1)
   * 2017 / 7 (2)
   * 2017 / 5 (3)
   * 2017 / 4 (6)
   * 2017 / 3 (2)
 * ▼ ▶
   2016 (31)
   * 2016 / 11 (3)
   * 2016 / 10 (4)
   * 2016 / 9 (4)
   * 2016 / 8 (3)
   * 2016 / 6 (2)
   * 2016 / 4 (2)
   * 2016 / 3 (6)
   * 2016 / 2 (5)
   * 2016 / 1 (2)
 * ▼ ▶
   2015 (67)
   * 2015 / 12 (24)
   * 2015 / 11 (6)
   * 2015 / 10 (1)
   * 2015 / 9 (2)
   * 2015 / 8 (2)
   * 2015 / 7 (1)
   * 2015 / 6 (2)
   * 2015 / 5 (3)
   * 2015 / 4 (2)
   * 2015 / 3 (8)
   * 2015 / 2 (6)
   * 2015 / 1 (10)
 * ▼ ▶
   2014 (107)
   * 2014 / 12 (11)
   * 2014 / 11 (6)
   * 2014 / 10 (7)
   * 2014 / 9 (4)
   * 2014 / 8 (5)
   * 2014 / 7 (9)
   * 2014 / 6 (5)
   * 2014 / 5 (5)
   * 2014 / 4 (5)
   * 2014 / 3 (9)
   * 2014 / 2 (21)
   * 2014 / 1 (20)
 * ▼ ▶
   2013 (25)
   * 2013 / 12 (1)
   * 2013 / 11 (2)
   * 2013 / 10 (1)
   * 2013 / 9 (2)
   * 2013 / 8 (2)
   * 2013 / 7 (3)
   * 2013 / 5 (4)
   * 2013 / 4 (4)
   * 2013 / 3 (1)
   * 2013 / 2 (2)
   * 2013 / 1 (3)
 * ▼ ▶
   2012 (42)
   * 2012 / 12 (2)
   * 2012 / 11 (2)
   * 2012 / 10 (1)
   * 2012 / 9 (2)
   * 2012 / 8 (13)
   * 2012 / 7 (4)
   * 2012 / 6 (2)
   * 2012 / 5 (1)
   * 2012 / 4 (3)
   * 2012 / 3 (2)
   * 2012 / 2 (8)
   * 2012 / 1 (2)
 * ▼ ▶
   2011 (67)
   * 2011 / 12 (4)
   * 2011 / 11 (5)
   * 2011 / 10 (5)
   * 2011 / 9 (4)
   * 2011 / 8 (4)
   * 2011 / 7 (4)
   * 2011 / 6 (4)
   * 2011 / 5 (4)
   * 2011 / 4 (6)
   * 2011 / 3 (5)
   * 2011 / 2 (15)
   * 2011 / 1 (7)
 * ▼ ▶
   2010 (29)
   * 2010 / 12 (3)
   * 2010 / 11 (1)
   * 2010 / 8 (1)
   * 2010 / 7 (1)
   * 2010 / 6 (6)
   * 2010 / 5 (5)
   * 2010 / 3 (1)
   * 2010 / 2 (11)

盆栽日記

Powered by Hatena Blog | ブログを報告する




引用をストックしました

ストック一覧を見る 閉じる

引用するにはまずログインしてください

ログイン 閉じる

引用をストックできませんでした。再度お試しください

閉じる

限定公開記事のため引用できません。

読者です 読者をやめる 読者になる 読者になる
29