survival analysisをRのggplot2で行う
ggplot2でsurvival analysisのプロットをしたいときは、 GGally パッケージに含まれる、ggsurv関数がある。
https://cran.r-project.org/web/packages/GGally/GGally.pdf
自分でggsurv関数を入れ込みたい場合は以下を参照。
上記から転載。
ggsurv <- function(s, CI = 'def', plot.cens = T, surv.col = 'gg.def', cens.col = 'red', lty.est = 1, lty.ci = 2, cens.shape = 3, back.white = F, xlab = 'Time', ylab = 'Survival', main = ''){ library(ggplot2) strata <- ifelse(is.null(s$strata) ==T, 1, length(s$strata)) stopifnot(length(surv.col) == 1 | length(surv.col) == strata) stopifnot(length(lty.est) == 1 | length(lty.est) == strata) ggsurv.s <- function(s, CI = 'def', plot.cens = T, surv.col = 'gg.def', cens.col = 'red', lty.est = 1, lty.ci = 2, cens.shape = 3, back.white = F, xlab = 'Time', ylab = 'Survival', main = ''){ dat <- data.frame(time = c(0, s$time), surv = c(1, s$surv), up = c(1, s$upper), low = c(1, s$lower), cens = c(0, s$n.censor)) dat.cens <- subset(dat, cens != 0) col <- ifelse(surv.col == 'gg.def', 'black', surv.col) pl <- ggplot(dat, aes(x = time, y = surv)) + xlab(xlab) + ylab(ylab) + ggtitle(main) + geom_step(col = col, lty = lty.est) pl <- if(CI == T | CI == 'def') { pl + geom_step(aes(y = up), color = col, lty = lty.ci) + geom_step(aes(y = low), color = col, lty = lty.ci) } else (pl) pl <- if(plot.cens == T & length(dat.cens) > 0){ pl + geom_point(data = dat.cens, aes(y = surv), shape = cens.shape, col = cens.col) } else if (plot.cens == T & length(dat.cens) == 0){ stop ('There are no censored observations') } else(pl) pl <- if(back.white == T) {pl + theme_bw() } else (pl) pl } ggsurv.m <- function(s, CI = 'def', plot.cens = T, surv.col = 'gg.def', cens.col = 'red', lty.est = 1, lty.ci = 2, cens.shape = 3, back.white = F, xlab = 'Time', ylab = 'Survival', main = '') { n <- s$strata groups <- factor(unlist(strsplit(names (s$strata), '='))[seq(2, 2*strata, by = 2)]) gr.name <- unlist(strsplit(names(s$strata), '='))[1] gr.df <- vector('list', strata) ind <- vector('list', strata) n.ind <- c(0,n); n.ind <- cumsum(n.ind) for(i in 1:strata) ind[[i]] <- (n.ind[i]+1):n.ind[i+1] for(i in 1:strata){ gr.df[[i]] <- data.frame( time = c(0, s$time[ ind[[i]] ]), surv = c(1, s$surv[ ind[[i]] ]), up = c(1, s$upper[ ind[[i]] ]), low = c(1, s$lower[ ind[[i]] ]), cens = c(0, s$n.censor[ ind[[i]] ]), group = rep(groups[i], n[i] + 1)) } dat <- do.call(rbind, gr.df) dat.cens <- subset(dat, cens != 0) pl <- ggplot(dat, aes(x = time, y = surv, group = group)) + xlab(xlab) + ylab(ylab) + ggtitle(main) + geom_step(aes(col = group, lty = group)) col <- if(length(surv.col == 1)){ scale_colour_manual(name = gr.name, values = rep(surv.col, strata)) } else{ scale_colour_manual(name = gr.name, values = surv.col) } pl <- if(surv.col[1] != 'gg.def'){ pl + col } else {pl + scale_colour_discrete(name = gr.name)} line <- if(length(lty.est) == 1){ scale_linetype_manual(name = gr.name, values = rep(lty.est, strata)) } else {scale_linetype_manual(name = gr.name, values = lty.est)} pl <- pl + line pl <- if(CI == T) { if(length(surv.col) > 1 && length(lty.est) > 1){ stop('Either surv.col or lty.est should be of length 1 in order to plot 95% CI with multiple strata') }else if((length(surv.col) > 1 | surv.col == 'gg.def')[1]){ pl + geom_step(aes(y = up, color = group), lty = lty.ci) + geom_step(aes(y = low, color = group), lty = lty.ci) } else{pl + geom_step(aes(y = up, lty = group), col = surv.col) + geom_step(aes(y = low,lty = group), col = surv.col)} } else {pl} pl <- if(plot.cens == T & length(dat.cens) > 0){ pl + geom_point(data = dat.cens, aes(y = surv), shape = cens.shape, col = cens.col) } else if (plot.cens == T & length(dat.cens) == 0){ stop ('There are no censored observations') } else(pl) pl <- if(back.white == T) {pl + theme_bw() } else (pl) pl } pl <- if(strata == 1) {ggsurv.s(s, CI , plot.cens, surv.col , cens.col, lty.est, lty.ci, cens.shape, back.white, xlab, ylab, main) } else {ggsurv.m(s, CI, plot.cens, surv.col , cens.col, lty.est, lty.ci, cens.shape, back.white, xlab, ylab, main)} pl }
Pubmedで文献検索
文献検索には型がある
ネット上で色々と教材が有るので、貼り付けておきます。
- 音羽病院廣江先生
- 東京大学図書館
http://www.lib.m.u-tokyo.ac.jp/manual/pubmedmanual.pdf
http://www.med.lib.keio.ac.jp/pdf/ug/ug_pubmed.pdf
- 京都大学図書館(基礎編)
http://www.kulib.kyoto-u.ac.jp/modules/refguide/docs/RefGuide_PubMed_basic.pdf
- The SPELL
EBM資料集−最もシンプルで分かりやすいPubMed 検索法 [The SPELL]
- 東大 小出先生
- editage
若手研究者のためのシステマティックレビューの書き方指南 | Editage Insights
Pubmedはいろいろと気を利かしてくれるのですが、その気配りの仕組を知り、MeSH termを用い、型通りに検索するのが大切です。
(MeSH termとは、Medical Subject Headingsの略)
胃のことをstomachという論文もあれば,gastricという論文もあるでしょう.癌についても同様で,cancerという論文もあるでしょうし,adenocarcinomaという論文もあります.これでは,思いついた単語だけでなく似たような単語を全部入れてみないと肝心な論文が漏れてしまう可能性があります.また逆に,この研究で肺癌の患者は除いた,という記述があるだけで検索されてしまいます.こうした不都合をなくすため,胃癌なら胃癌という語句は何か統一したキーワードを登録して検索を容易にしようという工夫がなされています.このキーワードがMeSH termであり,MeSH termを用いることで効率よく検索ができるのです.実際には胃癌は,Stomach neoplasmsというMeSH termで登録されています.(The SPELLより引用)
調べた方法を検索式として提示できるようにし、網羅的で再現性のある検索を心がけましょう。
すべては再現性です!
ggraphでグラフ作成
ggraph/README.md at master · thomasp85/ggraph · GitHub
ggplotのextensionはこんなところにも。
大腸がんスクリーニング (Colorectal Cancer Screening)
Colorectal Cancer Screening に関する書籍をざっと読んだ。大腸がんのスクリーニングが焦点。分子生物学的な観点も含めて記述してあって知識の確認になった。
- Colorectal Cancer Pathways
Petr Protiva - Risk Factors and Screening for Colorectal Cancer
Joseph C. Anderson - Hereditary Adenomatous Colorectal Cancer Syndromes
Maqsood Khan and Carol A. Burke - Screening and Surveillance Guidelines
Robert J. Chehade and Douglas J. Robertson - Barriers to Colorectal Cancer Screening: Patient, Physician, and System Factors
Catherine R. Messina - Screening for Colorectal Cancer Using Colonoscopy
Douglas K. Rex - New Colonoscopic Technologies for Colorectal Cancer Screening
Douglas K. Rex - Screening for CRC Using CT Colonography
Brooks D. Cash - Noninvasive Screening Tests
Nabil Fayad and Thomas F. Imperiale - Removal of Difficult Colon Polyps
Jerome Waye - Screening for Colorectal Cancer in the Elderly
Charles J. Kahi - Chemoprevention
Jeffrey Singerman and Petr Protiva
Rで21世紀の相関係数を算出する(MIC)
相関
2変数とも連続データで,正規分布に従っているならばPearsonの相関係数を用いることができる。
少なくとも1変数が非連続データの時にはノンパラメトリック検定のSpearmanやKendallの相関係数を用いる。これらは、実データでなく、順位付け(大小関係)で判定する。
これらの相関係数は、線形相関を見ている。
(訂正:Pearsonは線形の相関を調べるが、SpearmanやKendallは順位だけで見るので線形の仮定は要らない)
非線形相関を実データで調べる
Rにminervaというpackageがある。
MIC(Maximal information coefficient)は実データによる非線形相関係数も応用可能。
直線的な相関はピアソンで良いのですが、そうでないものは、順位検定になります。
直線ではない視覚的には明らかな相関を実データを用いて証明したい場合には、MICは有効かも。
R二乗値
エクセルでのR二乗値についての質問を受けたので、覚書程度に。
Microsoft officeのサイトによると、R二乗値(決定係数といいます)が近似曲線と並んで出る。近似曲線を引いた場合、その近似曲線がどのくらい実データを近似できているか、ということ。0.8を超えていたらよく、1に近いほどgood。
ここで、決定係数は、直線近似では相関係数の二乗ですが、非線形回帰では相関を示すものではない。 http://note.chiebukuro.yahoo.co.jp/detail/n168494
なので、「指数関数にこのくらい近似できる!!」とはいえても、「これくらい相関している」と言うためには注意を要する。
MICは、Science誌でも「21世紀の相関」として取り上げられている手法。
今度相関を調べる機会に使ってみようと思います。
データ解析チートシート cheat sheet
チートシート
R Studioのウェブサイトにまとまっており、
以下のように様々なチートシートがある。
Data Wrangling with dplyr and tidyr Cheat Sheet
Data Visualization with ggplot2 Cheat Sheet
これは便利だ。