Yoshi Nishikawa Blog

医学となにかのインタラクティブ

MCMC

ベイズ的＜ポストp値時代＞の統計学 R 統計学 計算機科学

はじめての統計データ分析

「はじめての統計データ分析」 に関して覚書。
1章の記事はこちら。

2章のハイライト

1) MCMC(Markov chain Monte Carlo method)により事後分布のデータを用い、
2) MCMCには、Metropolis-Hasting method, Gibbs sampling methods, Hamiltonian Monte Carlo methodなどが提案されています。
3) MCMCの欠点は、事後確率が定まらないことです。小数点xx桁まで行くと、施行毎に異なる結果を生じます。

この章をマスターすれば、あとは同じ概念の繰り返しという印象。

survival analysisをRのggplot2で行う

R ggplot2 医学

ggplot2でsurvival analysisのプロットをしたいときは、 GGally パッケージに含まれる、ggsurv関数がある。

https://cran.r-project.org/web/packages/GGally/GGally.pdf

自分でggsurv関数を入れ込みたい場合は以下を参照。

www.r-statistics.com

上記から転載。

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で文献検索

医学 研究 疫学 R以前の問題

文献検索には型がある

ネット上で色々と教材が有るので、貼り付けておきます。

20170424_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]

  • 東大 小出先生

文献検索とCritical Reading

  • 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より引用)

調べた方法を検索式として提示できるようにし、網羅的で再現性のある検索を心がけましょう。

すべては再現性です！

グラフ理論

R以前の問題 数学

読んだ。

1章　グラフの基礎
2章　最小全域木
3章　最短経路問題
4章　オイラー回路とハミルトン閉路
5章　グラフの彩色
6章　最大流問題
7章　マッチング

グラフ理論の入門編としておすすめできる。

私が入門編の書物に求める要素

1. 読みやすい

具体例とともに簡潔に示されている。

2. 用語が明記されている

次の段階も学ぶことを考えると、英語でも併記されていることがベスト。用語がわかれば検索できる。

3. 何となく見渡せる

体系的、というと怒られるかもしれないが、簡潔に、かつ、その分野の雰囲気がわかるものがベスト。

本書は三拍子揃っていた。

大腸がんスクリーニング (Colorectal Cancer Screening)

医学 臨床 疫学

Colorectal Cancer Screening に関する書籍をざっと読んだ。大腸がんのスクリーニングが焦点。分子生物学的な観点も含めて記述してあって知識の確認になった。

  1. Colorectal Cancer Pathways
    Petr Protiva
  2. Risk Factors and Screening for Colorectal Cancer
    Joseph C. Anderson
  3. Hereditary Adenomatous Colorectal Cancer Syndromes
    Maqsood Khan and Carol A. Burke
  4. Screening and Surveillance Guidelines
    Robert J. Chehade and Douglas J. Robertson
  5. Barriers to Colorectal Cancer Screening: Patient, Physician, and System Factors
    Catherine R. Messina
  6. Screening for Colorectal Cancer Using Colonoscopy
    Douglas K. Rex
  7. New Colonoscopic Technologies for Colorectal Cancer Screening
    Douglas K. Rex
  8. Screening for CRC Using CT Colonography
    Brooks D. Cash
  9. Noninvasive Screening Tests
    Nabil Fayad and Thomas F. Imperiale
  10. Removal of Difficult Colon Polyps
    Jerome Waye
  11. Screening for Colorectal Cancer in the Elderly
    Charles J. Kahi
  12. Chemoprevention
    Jeffrey Singerman and Petr Protiva

Rで21世紀の相関係数を算出する(MIC)

R R package 統計学 ネ申エクセル

相関

2変数とも連続データで，正規分布に従っているならばPearsonの相関係数を用いることができる。
少なくとも1変数が非連続データの時にはノンパラメトリック検定のSpearmanやKendallの相関係数を用いる。これらは、実データでなく、順位付け（大小関係）で判定する。
これらの相関係数は、線形相関を見ている。 （訂正：Pearsonは線形の相関を調べるが、SpearmanやKendallは順位だけで見るので線形の仮定は要らない）

非線形相関を実データで調べる

Rにminervaというpackageがある。 MIC(Maximal information coefficient)は実データによる非線形相関係数も応用可能。
直線的な相関はピアソンで良いのですが、そうでないものは、順位検定になります。
直線ではない視覚的には明らかな相関を実データを用いて証明したい場合には、MICは有効かも。

www.r-bloggers.com

R二乗値

エクセルでのR二乗値についての質問を受けたので、覚書程度に。
Microsoft officeのサイトによると、R二乗値（決定係数といいます）が近似曲線と並んで出る。近似曲線を引いた場合、その近似曲線がどのくらい実データを近似できているか、ということ。0.8を超えていたらよく、1に近いほどgood。

ここで、決定係数は、直線近似では相関係数の二乗ですが、非線形回帰では相関を示すものではない。 http://note.chiebukuro.yahoo.co.jp/detail/n168494

なので、「指数関数にこのくらい近似できる！！」とはいえても、「これくらい相関している」と言うためには注意を要する。

MICは、Science誌でも「21世紀の相関」として取り上げられている手法。

今度相関を調べる機会に使ってみようと思います。