Wolfeyes Bioinformatics beta

新しいオモチャへの期待

2年くらい前に10万円ほどする自転車を買った．その自転車はロードバイクの中でもエントリーモデルで特に性能が良いわけでもないのだが，今までの僕の自転車に乗るという体験を全く別のものに変えてしまうほどの乗り心地だった．今でも当時の感覚は鮮明に残っているし，自転車に乗るのはとても楽しい．もっとハイクラスな自転車を欲しくなるし，自転車をもっと上手く乗り回したいという気持ちもある．

自転車を買った当初は，これで僕の生活はガラリと変わるだろうと思っていた．前々からやりたいと思っていたツーリングにも行けるし，あっちこっち遊びに行って写真を撮ったり新しい体験ができる，そして適度な運動もできて健康になるかもしれない．そんな希望とも目標ともつかないようなやんわりとしたイメージを自転車に抱きながら，新しく手に入ったオモチャで遊んでいた．

そのような時期から2年が経ち，改めて現在の自分の生活を振り返ってみると，達成できたものもあり達成できなかったものもある．

ツーリングは続かなかった

自転車を買って，まず荒川の河川敷に出かけるところから始めた．そこには葛西臨海公園近くの河口付近から上流に向かって川沿いに整備された道があり，サイクリングロードではないものの，平坦で長い距離をずっと走ることができる．当時は新しい自転車に慣れていなかったり，そもそも自転車に乗るための筋力が全くなかったこともあって，トレーニングとしてその道をよく自転車で走った．信号機などの障害がほとんど無い道なので，走ろうと思えば30分でも1時間でも走っていられる．また，天気のいい日の荒川河川敷から見る日の出は本当に美しく，東京であることを忘れるくらいの清々しい気分にさせてくれる．このように，自転車乗りとしてはとても良い場所であることには間違いないのだが，徐々に出かける回数が減り，最近ではめっきり足が遠のいてしまった．ある程度飽きてしまったからかもしれないし，時間を確保してその日の天気と体調を確認してといったように，意識して準備しなければいけないことが多いからかもしれない．

都内の移動はもっぱら自転車になった

一方で，都内の主な移動手段としての自転車が，僕の生活の中で非常に大きな位置を占めるようになった．電車での移動がほとんど無くなり，代わりに自転車であちこちに出かけるようになった．それまで電車で移動していた距離くらいなら，自転車を使ったほうが早く着くし交通費が掛からない．もちろん遠出は無理だし天気のいい日に限るが，それでも普段の移動の半分以上は自転車に置き換わったと言っていい．電車と自転車を使い分けるときの境界線は，移動時間がだいたい30分を超えるか超えないかくらいだと思う．それくらいの時間でも，僕の家からだと山手線の中はほとんど網羅できる．そうなると，もはや普段からよく行く場所は自然と自転車での移動となった．おまけに，自転車で移動するコストは非常に低いので，普段電車で行くほどでもないような場所にも気軽に足を運ぶことができる．神保町なんかは自転車に乗るようになってから頻繁に通うようになった場所の一つだ．特に欲しい物が無いときでも，ふらっと出かけては店の表に積まれている安売りの本を眺めたりしている．また，深夜に新宿に出掛けては，終電に乗れないくらいの時間帯の空いている映画館で一人映画を見ることも多くなった．といったように，移動手段の変化によって僕の生活の一部は自転車無しでは考えられないくらいにまでに色々な習慣が定着した．

まとめ

結局のところ自転車を買ったことで生活は色々な面で変わったものの，当初予想していたところからはだいぶ違うところに来てしまったしまった．それでも今の自転車のある生活にはかなり満足しているし，僕の性格がよく表れているような気がする．

私信

ということで，たとえ僕が高級デジイチを買ったとしても途中で飽きそうなので，カメラは買いません．

今年も残すところあと1ヶ月ということで，毎年恒例になりつつあるAdvent Calendarが各所で盛り上がっているようです．私もQiitaで開催されているMachine Learning Advent Calendar 2012の2日目で，ギブスサンプリングに関する記事を書きました．こんな動画を撮って遊んでいます．もしよければ本記事の方も御覧ください．

事の発端

これこれ、これですよ / “PRML読書会第一回 - Seeking for my unique color.” htn.to/w8bGen
— Yokoyama T (@wakuteka) November 26, 2012

先日「PRML読書会第一回 - Seeking for my unique color.」のソースコードが一部界隈で局地的に話題にあがり，やれ

takousikiってなんだよロシア人かよ

だとか

evalが多すぎて引いた

だのといった罵詈雑言が（主に僕の口から）飛び交っていたわけだけれども，さすがに言うだけではsyou6162先生に申し訳ないと感じ，改めて彼のコードをRで書きなおしてみることにした．

多項式曲線フィッティングのコードを書き直す

f <- function(x){sin(2*pi*x)}

takousiki <- function(n){
  x <- numeric(0)
  z <- numeric(0)
  for(i in 1:n){
    x <- cbind(x, seq(0,1,by=0.1)^i)
    z <- cbind(z, seq(0,1,by=0.01)^i)
  }
  data <- data.frame(f(seq(0,1,by=0.1))+rnorm(11, mean=0, sd=0.5), x)

  colnames(data) <- c("y",paste("x",as.character(1:n),sep=""))
  colnames(z) <- paste("x",as.character(1:n),sep="")

  # predict
  fit <- lm(y ~ ., data=data)

  plot(data[,2], data[,1], ylab="", ylim=c(-1.5,1.5), xlab="x", main=paste("M=",n))
  curve(f(x), add=T, lwd=2)
  lines(seq(0,1,by=0.01), predict(fit,data.frame(z)), col="red", type="l")
}

使い方としては，

takousiki(1)
takousiki(3)
takousiki(8)

とすると，多項式で曲線フィッティングした結果がプロットされる．ただし切片だけの0次は書けない．

サンプルデータの生成と推定した曲線の描写用のデータの生成を同時に行なっているため，コードがちょっと見にくいかもしれない．基本的に $x^2$ や $x^3$ などの多項式の各項の値をあらかじめデータフレームとして作っておいて，lm()関数で回帰するときのモデルを簡略化している．あと，関数名はオリジナルを踏襲しているので，あしからず．

他の方法としては

lm()関数などに渡すモデルは，as.formula()関数を使うことによって文字列から生成することができる．この方法もかなり簡単に使うことができるので，データフレームの列名を規則的に付けている場合などは使い勝手は良いと思われる．

> n <- 5
> as.formula(paste("y~",paste("x",1:n,sep="",collapse="+")))
y ~ x1 + x2 + x3 + x4 + x5

読書ノートアーカイブ：http://yagays.github.com/blog/2012/10/20/archive-introducing-monte-carlo-methods-with-r/

練習問題 6.13

この問題は例6.3と同じく，線形回帰におけるパラメータの事後分布からメトロポリス・ヘイスティング・アルゴリズムを用いてサンプリングを行うことを目標とする．線形回帰の対象はmcsmパッケージに収められているchallengerデータセットで，過去に打ち上げられたスペースシャトルにおけるOリングの破損の有無と，その時の気温に関するものである．

a. glm()関数を使ってロジスティック回帰を行う

打ち上げられたスペースシャトルのOリングは破損しているかしていないかの2値で表されるため，今回はロジスティック回帰を用いた当てはめを行う．Oリングが破損している場合は1を，破損していない場合には0が割り当てられており，ある気温 $x_i$ におけるOリングの破損確率は

$P(x_i) = \frac{\exp(\alpha+\beta_i)}{1+\exp(\alpha+\beta_i)}$

で表される．今回はこの $\alpha$ と $\beta$ をglm()関数を使って推定する．

なお，以下のコードではロジスティック関数を式変形して $\frac{1}{1+\exp(-(\alpha+\beta_i))}$ として計算している．

library(mcsm)
data(challenger)

fit <- glm(oring ~ temp, data=challenger, family=binomial)
plot(challenger$temp, challenger$oring, pch=19)
curve(1/(1+exp(-(fit$coef[1]+fit$coef[2]*x))), from=53, to=81, add=T, col="orange", lwd=4)

この図は，横軸に気温（華氏），縦軸にOリングの破損の有無をプロットしたものである．また，glm()関数で推定したロジスティック曲線をオレンジ色で示しており，任意の気温における破損する確率と対応している．この図を見ると，華氏65前後で破損確率が50%となり，気温が高い場合は破損確率が低く，逆に気温が低い時には破損確率が高いことがわかる．

glm()関数のパラメータ推定の詳細は以下の通り．

> summary(fit)

Call:
glm(formula = oring ~ temp, family = binomial, data = challenger)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-1.0611  -0.7613  -0.3783   0.4524   2.2175

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  15.0429     7.3786   2.039   0.0415 *
temp         -0.2322     0.1082  -2.145   0.0320 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 28.267  on 22  degrees of freedom
Residual deviance: 20.315  on 21  degrees of freedom
AIC: 24.315

Number of Fisher Scoring iterations: 5

b. ロジスティック回帰で推定したパラメータの値を用いて，メトロポリス・ヘイスティング・アルゴリズムでサンプリングを行う

これは前回と殆ど変わらず．

y <- challenger$oring
x <- challenger$temp

loglik <- function(a,b){sum(y*(a+b*x) - log(1+exp(a+b*x)))}
logq <- function(a,b){
  dnorm(a,mean=a.est[1],sd=a.std,log=T)
  + dnorm(b,mean=b.est[1],sd=b.std,log=T)
}

Niter <- 5000
n <- length(challenger$temp)
fit.param <- summary(fit)$coef
a.est <- fit.param[1,1]
a.std <- fit.param[1,2]
b.est <- fit.param[2,1]
b.std <- fit.param[2,2]


for(i in 2:Niter){
  candidate <- c(rnorm(1, mean=a.est[1], sd=a.std), rnorm(1, mean=b.est[1], sd=b.std))
  fy <- loglik(candidate[1], candidate[2])
  fx <- loglik(a.est[i-1], b.est[i-1])
  qxy <- logq(a.est[i-1], b.est[i-1])
  qyx <- logq(candidate[1], candidate[2])
  rho <- min(1, exp(fy-fx+qxy-qyx))

  if (runif(1)<rho){
    a.est[i] <- candidate[1]
    b.est[i] <- candidate[2]
  }else{
    a.est[i] <- a.est[i-1]
    b.est[i] <- b.est[i-1]
  }
}

plot(x, y, xlab="temp", ylab="oring", type="n")
for(i in (Niter-500):Niter){
  lines(x, 1/(1+exp(-(a.est[i]+b.est[i]*x))), lwd=2, col="grey")
}
curve(1/(1+exp(-(fit$coef[1]+fit$coef[2]*x))), from=53, to=81, add=T, col="orange", lwd=4)
points(x, y, pch=19)

c. 5000回のマルコフ連鎖によるサンプル列を生成し，サンプリングで得られた推定値を使ったロジスティック曲線を図示する

実際に上で書いたコードを動かして，事後分布からのサンプリングを行い，得られたサンプルのうち4500回から5000回までの反復で生成した当てはめ曲線を灰色の曲線で表している．

また，パラメータ $\alpha$ と $\beta$ の系列の変動，自己相関プロット，系列の分布を以下の図で示している．

par(mfrow=c(2,3))
plot(1:Niter, a.est, type="l", xlab="Iterations", main=expression(alpha))
acf(a.est)
plot(density(a.est))
plot(1:Niter, b.est, type="l", xlab="Iterations", main=expression(beta))
acf(b.est)
plot(density(b.est))

そして最後に，サンプリング推定値差の変動を示したのが以下の図である．glm()関数で推定した値を赤線で表示しており，そこからどれだけ離れているかを見ることができる．

plot(1:Niter, apply(cbind(a.est,b.est),1,function(x){loglik(x[1],x[2])+logq(x[1],x[2])}),
     type="l", xlab="Iteration", ylab="")
abline(h=(loglik(fit.param[1,1],fit.param[2,1])+logq(fit.param[1,1],fit.param[2,1])), col="red", lwd=2)

d. サンプリングで得られたパラメータの分布から，華氏60,50,40度のときの故障の確率と標準誤差を推定する

サンプリングで得られた値を用いて，ある温度における故障の確率と標準誤差を求める．なお，今回はBurn-inによる推定値の切り捨てを行なっておらず，サンプリング回数5000回で得られた推定値を全て使っている．

for(t in c(60,50,40)){
  p <- 1/(1+exp(-(a.est+b.est*t)))
  print(paste("temp:", t, "=>", "mean=", round(mean(p),digit=4), ", S.D.=", round(sd(p),digit=4) ))
}
[1] "temp: 60 => mean= 0.7476 , S.D.= 0.1282"
[1] "temp: 50 => mean= 0.9525 , S.D.= 0.0641"
[1] "temp: 40 => mean= 0.9895 , S.D.= 0.0248"

一般化線形モデルやロジスティックモデルなどの基本的な線形回帰を復習するために，久しぶりに本書を引っ張りだしてきてひと通り読み返した．自分としては理論や基本的な部分はある程度抑えているつもりだったけれども，検定であったり上手く回帰をするためのデータの条件であったりといった細かな部分に関しては，まだまだ不勉強なところが多い．とにかくデータに対して回帰をかけてこれで良しみたいな部分があったので，そういった部分を補うという意味でも良い復習になった．最近はMCMCにもっぱら専念しているが，基礎であったり古典の部分もしっかり抑えておかなければ…．

タイトルに「医療統計学」と難しそうな名前が付いている本書だが，原書のタイトルが”Introductory statistics with R”であるように，実質的には統計解析ソフトRの入門書である．例題で扱う題材として医療統計のデータが使われているだけで，内容自体に医療分野の専門性は無い．本書前半では，基本的なRの使い方とともに統計解析の基礎が説明され，後半になると，線形回帰や重回帰分析，ロジスティック回帰，生存分析などの具体的な解析手法が紹介される．全体を通して，サンプルデータを用いた豊富な解析実例とRのコード例がセットになっているため，必要ならば解析を再現することも可能である．そのため，Rを勉強するのも良し，解析手法を勉強するも良しと，初心者から中級者にかけての読者向けの非常にバランスの取れた構成となっている．また，理論よりも実戦向きの内容となっており，検定方法やモデル選択の比重が大きいのが本書の特徴でもある．そういえば，私がRを触り始めた頃は「The R Tips―データ解析環境Rの基本技・グラフィックス活用集」がまさにバイブルといった感じだったが，プログラミングとしてのRの解説がメインで統計解析の流れを説明した具体例は少なく，実用的な部分は本書をよく参照していた記憶がある．

ということで本書はR本の中でも結構オススメできる部類なのだが，どうやらAmazonや丸善の公式ページ，丸善&ジュンク堂のネットストアなどの主だったネット書店で品切れが相次いでいるようだ．ともなると，本書を手に入れるにはAmazonのマーケットプレイスで高い値段で買うか，大きな書店で売れ残っている本を見つけるかするしか無さそうだ．まあ今となってはRと統計に関する書籍は山ほど出版されているので，本書がこれから増刷などが掛かることは難しかもしれないが，このまま絶版になっていくには実に惜しい本だと思う．