Hatena::ブログ(Diary)

驚異のアニヲタ社会復帰への道

Prima Project

2017-12-13

2017-12-12

予後因子 prognostic と予測因子 predictive は何が違うのか

興味深い話を教えてもらったので読んだ。

Understanding Prognostic versus Predictive Biomarkers

 

よく予後因子や予測因子となるバイオマーカーを探そう、という研究があるが、そもそも予後因子/予測因子とはなんぞや、という話。

現実世界においては、あるバイオマーカーがあったときに、厳密に予後因子であるか予測因子であるかは分けられないが、概念的には両者は異なる。

本文中では、

To identify a predictive biomarker, there generally should be a comparison of a treatment to a control in patients with and without the biomarker.

the biomarker illustrated in Figures 1A and 1B is prognostic but not predictive and will not be helpful in choosing between standard and experimental therapy.

とあるので、統計素人の理解としては、

臨床研究は新薬 vs 標準治療(もしくはプラセボ)の比較なので、比較した群において差があるかどうかを左右する因子であることが予測因子 predictive。比較しない群で(生存に)差があるかどうかが予後因子 prognostic。

因子のある/なしで、新薬を使うか標準治療を使うかの判断に影響するのが予測因子 predictive。

と思った。

 

例えばよく誤解されやすいのが、あるバイオマーカー(BM)のある/なしで、新薬 TRT+ (TRT- は標準治療でもプラセボでもなんでもいいが、とにかく比較される対照)についてこんなに差があったとする。(Fig 1A)

f:id:MikuHatsune:20171212152734p:image

このとき、BM のある/なしで生存に差があるから、BM を確認することで予後が予測できる、と思ってしまう。このとき、BM は予後因子ではあるが予測因子ではない。

 

なぜかというと、臨床研究はそもそも、新薬 vs 対照の比較研究だからである。BM のある/なしによる生存の差は検討課題ではない。実際の臨床研究の設定として、新薬 TRT+ と対照 TRT - の比較を、BM のある/なしで層別化したらこんな感じになる。

ここで、BM のある/なしは黒赤で、TRT は実線/破線で2x2 になっている。plot の重なりをずらすために値が微妙に変えてあるが、本質的には重なっていると思ってよい。(Fig 1B)

f:id:MikuHatsune:20171212152731p:image

さてここで、TRT + の集団に限ってみれば、BM のある/なしは生存に差があった。しかし、TRT - を見ても、BM のある/なしはTRT + の生存の差と同じ程度に生存に差が出るようになっている。

つまり、TRT + にしようがTRT - にしようが、BM によって生存に差が(ほぼ同等に)出るようになる。なので、BM があろうがなかろうが、「TRT + かTRT - にすべきかどうか」は悩む必要はないのである。だが、BM を見れば、その人が長生きしそうか、しなさそうなのかの「予後」には影響する。なのでこの状況ではTRT +/- の比較はされていないし、予後因子になるのである。

 

4本も直線があるとわかりにくいので、ばらばらに描いてみる。ここで予後因子らしさを出しているのは黄色のマスである。いま、TRT + の集団において、BM のある/なしは生存に影響しており、また、TRT - の集団においても、M のある/なしは生存に影響している

しかし、色のついていないマスをみると、BM がある集団に対してTRT + とTRT - は差がない。BM がない集団も同様である。臨床研究で比較研究しているのはこの状況であるので、BM のステータスはTRT の選択に影響しないから、TRT の選択結果(未来)を予測しない。たぶん。

f:id:MikuHatsune:20171212152727p:image

 

予測因子の話になるが、予後因子の最初の図とは対照的に、一見して予後にも予測にも役に立たなさそうな場合はある。これは、TRT + の集団においてBM のある/なしが生存に影響していない。(Fig 2A)

f:id:MikuHatsune:20171212152720p:image

 

しかしながら、何度も言うように、臨床研究は新薬 vs 対照の比較研究なので、TRT - が登場しないと結論は導けない。ここで、こんな感じの図になると、これは予測因子となる。(Fig 2B)

f:id:MikuHatsune:20171212152715p:image

 

少しややこしいので、先に2x2 に分割する。予測因子で注目するのは緑のマスである。いま、ひだり上のTRT + の集団においては、BM のある/なしは生存に影響しなかったが、みぎ下(予後因子では黄色にしていたマス)のTRT - の集団においてはBM のある/なしで生存に差があるので、これは実は予後因子でもある。しかし、TRT + とTRT - の比較はされていない。

注目すべきはひだり下の緑マスである。BM がない集団においては、TRT + とTRT - の生存に差がある。

みぎ上の緑マスは、BM がある集団において、TRT + とTRT - は差がない。

つまり、BM がない、と判定された集団においては、TRT の選択が生存に影響する。これが予測因子である。BM がある集団は、TRT の選択に関係がないが、TRT がない集団はTRT を選択しなければ生存が伸びない。というわけで選択を迫られる因子が予測因子である。たぶん。

f:id:MikuHatsune:20171212155957p:image

 

BM の付け方を逆にしてしまったので本文のbiomarker positive/negative は逆に読み替えてください。

# 生存曲線を適当に
x <- seq(0, 5, length=3000)
surv <- function(m){
  exp(-log(2)/m * x)
}

# plot 用
xl <- range(x)
yl <- c(0, 1)
m <- c(8, 3, 8.0-0.5, 3-0.3)
cols <- c("black", "red")
lty <- c(1, 3)
PN <- c("+", "-")
i <- c(1, 2, 1, 2)
j <- c(1, 1, 2, 2)


# prognostic
par(mar=c(4, 5, 2, 2), las=1, cex.axis=1.5, cex.lab=1.5)
plot(0, type="l", xlim=xl, ylim=yl, lwd=4, xlab="Time", ylab="Survival")
lines(x, surv(m[1]), lwd=3, lty=lty[1], col=cols[1])
lines(x, surv(m[2]), lwd=3, lty=lty[1], col=cols[2])
legend("bottomleft", legend=sprintf("BM %s;\tTRT +", PN), lty=lty[1], lwd=3, col=cols, text.font=2)

par(mar=c(4, 5, 2, 2), las=1, cex.axis=1.5, cex.lab=1.5)
plot(0, type="l", xlim=xl, ylim=yl, lwd=4, xlab="Time", ylab="Survival")
lines(x, surv(m[1]), lwd=3, lty=lty[1], col=cols[1])
lines(x, surv(m[2]), lwd=3, lty=lty[1], col=cols[2])
lines(x, surv(m[3]), lwd=3, lty=lty[2], col=cols[1])
lines(x, surv(m[4]), lwd=3, lty=lty[2], col=cols[2])
legend("bottomleft", legend=sprintf("BM %s;\tTRT %s", PN[j], PN[i]), lty=lty[i], lwd=3, col=cols[j], text.font=2)


# 2x2 comparison
par(mfrow=c(2, 2))
par(mar=c(4, 5, 2, 2), las=1, cex.axis=1.5, cex.lab=1.5)
plot(0, type="l", xlim=xl, ylim=yl, lwd=4, xlab="Time", ylab="Survival")
polygon(par()$usr[c(1,2,2,1)], par()$usr[c(3,3,4,4)], col="yellow")
lines(x, surv(m[1]), lwd=3, lty=lty[1], col=cols[1])
lines(x, surv(m[2]), lwd=3, lty=lty[1], col=cols[2])
#lines(x, surv(m[3]), lwd=3, lty=lty[2], col=cols[1])
#lines(x, surv(m[4]), lwd=3, lty=lty[2], col=cols[2])
legend("bottomleft", legend=sprintf("BM %s;\tTRT %s", PN[j], PN[i]), lty=lty[i], lwd=3, col=cols[j], text.font=2)

par(mar=c(4, 5, 2, 2), las=1, cex.axis=1.5, cex.lab=1.5)
plot(0, type="l", xlim=xl, ylim=yl, lwd=4, xlab="Time", ylab="Survival")
lines(x, surv(m[1]), lwd=3, lty=lty[1], col=cols[1])
#lines(x, surv(m[2]), lwd=3, lty=lty[1], col=cols[2])
lines(x, surv(m[3]), lwd=3, lty=lty[2], col=cols[1])
#lines(x, surv(m[4]), lwd=3, lty=lty[2], col=cols[2])
legend("bottomleft", legend=sprintf("BM %s;\tTRT %s", PN[j], PN[i]), lty=lty[i], lwd=3, col=cols[j], text.font=2)

par(mar=c(4, 5, 2, 2), las=1, cex.axis=1.5, cex.lab=1.5)
plot(0, type="l", xlim=xl, ylim=yl, lwd=4, xlab="Time", ylab="Survival")
#lines(x, surv(m[1]), lwd=3, lty=lty[1], col=cols[1])
lines(x, surv(m[2]), lwd=3, lty=lty[1], col=cols[2])
#lines(x, surv(m[3]), lwd=3, lty=lty[2], col=cols[1])
lines(x, surv(m[4]), lwd=3, lty=lty[2], col=cols[2])
legend("bottomleft", legend=sprintf("BM %s;\tTRT %s", PN[j], PN[i]), lty=lty[i], lwd=3, col=cols[j], text.font=2)

par(mar=c(4, 5, 2, 2), las=1, cex.axis=1.5, cex.lab=1.5)
plot(0, type="l", xlim=xl, ylim=yl, lwd=4, xlab="Time", ylab="Survival")
polygon(par()$usr[c(1,2,2,1)], par()$usr[c(3,3,4,4)], col="yellow")
#lines(x, surv(m[1]), lwd=3, lty=lty[1], col=cols[1])
#lines(x, surv(m[2]), lwd=3, lty=lty[1], col=cols[2])
lines(x, surv(m[3]), lwd=3, lty=lty[2], col=cols[1])
lines(x, surv(m[4]), lwd=3, lty=lty[2], col=cols[2])
legend("bottomleft", legend=sprintf("BM %s;\tTRT %s", PN[j], PN[i]), lty=lty[i], lwd=3, col=cols[j], text.font=2)


# predictive
par(mar=c(4, 5, 2, 2), las=1, cex.axis=1.5, cex.lab=1.5)
plot(0, type="l", xlim=xl, ylim=yl, lwd=4, xlab="Time", ylab="Survival")
lines(x, surv(m[1]), lwd=3, lty=lty[1], col=cols[1])
lines(x, surv(m[3]), lwd=3, lty=lty[1], col=cols[2])
legend("bottomleft", legend=sprintf("BM %s;\tTRT +", PN), lty=lty[1], lwd=3, col=cols, text.font=2)

m[2] <- m[1] - 0.3
par(mar=c(4, 5, 2, 2), las=1, cex.axis=1.5, cex.lab=1.5)
plot(0, type="l", xlim=xl, ylim=yl, lwd=4, xlab="Time", ylab="Survival")
lines(x, surv(m[1]), lwd=3, lty=lty[1], col=cols[1])
lines(x, surv(m[2]), lwd=3, lty=lty[1], col=cols[2])
lines(x, surv(m[3]), lwd=3, lty=lty[2], col=cols[1])
lines(x, surv(m[4]), lwd=3, lty=lty[2], col=cols[2])
legend("bottomleft", legend=sprintf("BM %s;\tTRT %s", PN[j], PN[i]), lty=lty[i], lwd=3, col=cols[j], text.font=2)

# 2x2 comparison
par(mfrow=c(2, 2))
par(mar=c(4, 5, 2, 2), las=1, cex.axis=1.5, cex.lab=1.5)
plot(0, type="l", xlim=xl, ylim=yl, lwd=4, xlab="Time", ylab="Survival")
lines(x, surv(m[1]), lwd=3, lty=lty[1], col=cols[1])
lines(x, surv(m[2]), lwd=3, lty=lty[1], col=cols[2])
#lines(x, surv(m[3]), lwd=3, lty=lty[2], col=cols[1])
#lines(x, surv(m[4]), lwd=3, lty=lty[2], col=cols[2])
legend("bottomleft", legend=sprintf("BM %s;\tTRT %s", PN[j], PN[i]), lty=lty[i], lwd=3, col=cols[j], text.font=2)

par(mar=c(4, 5, 2, 2), las=1, cex.axis=1.5, cex.lab=1.5)
plot(0, type="l", xlim=xl, ylim=yl, lwd=4, xlab="Time", ylab="Survival")
polygon(par()$usr[c(1,2,2,1)], par()$usr[c(3,3,4,4)], col="lightgreen")
lines(x, surv(m[1]), lwd=3, lty=lty[1], col=cols[1])
#lines(x, surv(m[2]), lwd=3, lty=lty[1], col=cols[2])
lines(x, surv(m[3]), lwd=3, lty=lty[2], col=cols[1])
#lines(x, surv(m[4]), lwd=3, lty=lty[2], col=cols[2])
legend("bottomleft", legend=sprintf("BM %s;\tTRT %s", PN[j], PN[i]), lty=lty[i], lwd=3, col=cols[j], text.font=2)

par(mar=c(4, 5, 2, 2), las=1, cex.axis=1.5, cex.lab=1.5)
plot(0, type="l", xlim=xl, ylim=yl, lwd=4, xlab="Time", ylab="Survival")
polygon(par()$usr[c(1,2,2,1)], par()$usr[c(3,3,4,4)], col="lightgreen")
#lines(x, surv(m[1]), lwd=3, lty=lty[1], col=cols[1])
lines(x, surv(m[2]), lwd=3, lty=lty[1], col=cols[2])
#lines(x, surv(m[3]), lwd=3, lty=lty[2], col=cols[1])
lines(x, surv(m[4]), lwd=3, lty=lty[2], col=cols[2])
legend("bottomleft", legend=sprintf("BM %s;\tTRT %s", PN[j], PN[i]), lty=lty[i], lwd=3, col=cols[j], text.font=2)

par(mar=c(4, 5, 2, 2), las=1, cex.axis=1.5, cex.lab=1.5)
plot(0, type="l", xlim=xl, ylim=yl, lwd=4, xlab="Time", ylab="Survival")
#lines(x, surv(m[1]), lwd=3, lty=lty[1], col=cols[1])
#lines(x, surv(m[2]), lwd=3, lty=lty[1], col=cols[2])
lines(x, surv(m[3]), lwd=3, lty=lty[2], col=cols[1])
lines(x, surv(m[4]), lwd=3, lty=lty[2], col=cols[2])
legend("bottomleft", legend=sprintf("BM %s;\tTRT %s", PN[j], PN[i]), lty=lty[i], lwd=3, col=cols[j], text.font=2)

2017-12-09

Python っぽい感じでうまくクラスの概念を使ってR を書きたい

写経しながらR で書き換えつつやってみている。

素のR ユーザーはクラスとかそういう概念に触れない or 真面目に勉強しないままきているのでよくわからなくなってる。

ここらへんを参考にしてクラスっぽいものを作ってみる。

R6 class さっそく調査 #rstatsj - Qiita

[翻訳]R6 vignette: R6クラス入門 - Qiita

Rでオブジェクト指向っぽくクラスを使ってみる(R6メソッド) - Qiita

 

ことの発端としては、コンストラクタを作ってそれに対して関数を使ったとき、もとのコンストラクタが更新されて値を持っているという事態に遭遇したことである。

例として、上本のPython コードを出す。

# Python
class MulLayer:
  def __init__(self):
    self.x = None
    self.y = None
  
  def forward(self, x, y):
    self.x = x
    self.y = y                
    out = x * y
    return out
  
  def backward(self, dout):
    dx = dout * self.y
    dy = dout * self.x
    return dx, dy

apple = 100
apple_num = 2
tax = 1.1

# layer
mul_apple_layer = MulLayer()   ### 1
mul_tax_layer = MulLayer()

# forward
apple_price = mul_apple_layer.forward(apple, apple_num)   ### 2
price = mul_tax_layer.forward(apple_price, tax)

# backward
dprice = 1
dapple_price, dtax = mul_tax_layer.backward(dprice)
dapple, dapple_num = mul_apple_layer.backward(dapple_price)

MulLayer というクラスが全体に定義されて、その中で__init__, forward, backward という関数(?) が定義されている。

###1 の部分で、mul_apple_layer というオブジェクトがMulLayer() で作られるが、MulLayer そのものは関数ではなくクラスであって、MulLayer() でオブジェクトが作られる意味がわからない。

これは、__init__ というものがself というものを勝手に作るようにできていて、MulLayer() によって .x と .y によって呼び出せるものが出来上がる、というようになっている。確かに、

mul_apple_layer.x
mul_apple_layer.y

は、エラーをはかない(かわりに何も持ってないので出力もない)が、

mul_apple_layer.z
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AttributeError: 'MulLayer' object has no attribute 'z'

と、持っていない .z についてはエラーをはく。

 

さてこれを何も考えずにR で書くとこうなる。

MulLayer <- function(){
  self <- list(x=NULL, y=NULL)
  return(self)
}

forward <- function(self, x, y){
  self$x <- x
  self$y <- y
  out <- x * y
  return(list(self=self, out=out))
}

backward <- function(self, dout){
  dx <- dout * self$y
  dy <- dout * self$x
  return(list(dx=dx, dy=dy))
}

apple <- 100
apple_num <- 2

# layer
mul_apple_layer <- MulLayer()

# forward
tmp <- forward(mul_apple_layer, apple, apple_num)
apple_price <- tmp$out
mul_apple_layer <- tmp$self   ### 更新作業が必要

MulLayer という関数を作っただけなので、引数なしで、出力としてself というものを返す。

ただし、forward 関数では引数としてのself は、関数内でごそごそいじっても入力としてのself は変更されない。たぶん <<- とか使うといいのかと思ったけどこれでも .x と .y にあたる $x, $y は変更されない。

というわけで、tmp という中間変数的なオブジェクトを作って、再度オブジェクトを作る、というような二度手間をやるはめになる。

 

そこで、R6 というパッケージを使って、クラスなるものを定義して書いてみる。

上のPython と同じように書くと、こうなる。

library(R6)

MulLayer = R6Class("MulLayer",
  public = list(
    x = NULL,
    y = NULL,
    initialize = function(){
      self$x = NA
      self$y = NA
    },
    forward = function(x, y){
      self$x <- x
      self$y <- y
      out <- self$x*self$y
      return(out)
    },
    backward = function(dout){
      dx <- dout * self$y
      dy <- dout * self$x
      return(list(dx=dx, dy=dy))
    }
  )
)

m <- MulLayer$new()
<MulLayer>
  Public:
    backward: function (dout) 
    clone: function (deep = FALSE) 
    forward: function (x, y) 
    initialize: function () 
    x: NA
    y: NA

いま、MulLayer というクラスができて、それで作られたオブジェクトには、$x (.x と同じ) と$y (.y と同じ)という変数が付属している。また、MulLayer に適応できる関数として、forward とbackward という関数がある。

initialize は、これだけを発動すると、$x と$y を持つオブジェクトができる。

m$x
[1] NA
m$y
[1] NA

である。

さてここで、Python っぽいことをしてみよう。MulLayer には.forward という関数が使えて、これの引数は(self, x, y) だった。ただし、self はself.forward とすることで、self の値が x, y に応じて変更されながら、return out するという機能がある。

R6 ではこれをm$forward() という書き方でできて、

m$forward(100, 2)
[1] 200

と、確かにx * y の値が出力される。このとき、m は

m
<MulLayer>
  Public:
    backward: function (dout) 
    clone: function (deep = FALSE) 
    forward: function (x, y) 
    initialize: function () 
    x: 100
    y: 2

となっていて、確かにm$x (m.x と同じ)、m$y (m.y と同じ)とすると、それぞれx, y の値が入って更新されたm が出力される。

というわけで、tmp を作っていちいち更新、という作業がなくなる、と思われる。

 

R6 の書き方を教えてくれたSK 氏に感謝。

2017-12-04

MikuHatsune2017-12-04

ラブライブ!サンシャイン!! の名前呼び合いグラフをかく 9話目

saint snow の姉妹の仲がギクシャクしたのでルビィら1年生組が函館に残ってsaint snow のために歌を作ったけど、結局千歌がセンターを奪った話。

ダイヤさんの姉みがやばい。DQN ネームじゃなかったら推しにできたのに…

 

果南ちゃんは誰も名前を呼ばず、また誰からも名前を呼ばれなかったので孤立しました。

f:id:MikuHatsune:20171204223916p:image

 

ラブライブ!サンシャイン!! の名前呼び合いグラフをかく 1話目と2話目 - 驚異のアニヲタ社会復帰への道

ラブライブ!サンシャイン!! の名前呼び合いグラフをかく 3話目 - 驚異のアニヲタ社会復帰への道

ラブライブ!サンシャイン!! の名前呼び合いグラフをかく 4話目 - 驚異のアニヲタ社会復帰への道

ラブライブ!サンシャイン!! の名前呼び合いグラフをかく 5話目 - 驚異のアニヲタ社会復帰への道

ラブライブ!サンシャイン!! の名前呼び合いグラフをかく 6話目 - 驚異のアニヲタ社会復帰への道

ラブライブ!サンシャイン!! の名前呼び合いグラフをかく 7話目 - 驚異のアニヲタ社会復帰への道

ラブライブ!サンシャイン!! の名前呼び合いグラフをかく 8話目 - 驚異のアニヲタ社会復帰への道

2017-12-03

MikuHatsune2017-12-03

きららフェスタ2017 に出演していた声優たちの集客力をRstan で推定する

この記事は

RStudio Advent Calendar 2017 - Qiita

まんがタイムきらら Advent Calendar 2017

Stan Advent Calendar 2017 - Qiita

R Advent Calendar 2017 - Qiita

の3日目の配当記事です。

 

声優統計第9号で、きららフェスタ2016に出演した声優の集客力について検討したが、2017年度もキララフェスタに参加したのでまったく同じことをやってみる。

 

2017年では4作品のアニメが出て、東山奈央大久保瑠美が総合司会ということで出演していた。この二人は4作品には出演していなかった。

声優アニメ出演パターンは以下の通りである。茅野愛衣が3作品に出ている。水瀬いのりも意外にサクラクエストに出演していた。

f:id:MikuHatsune:20171203001821p:image

 

声優の集客力の推定結果は以下の通りである。サクラクエスト声優陣の推定がグダグダである。

2016年度ほどに、いろいろなアニメにほどよく散らばって出演している声優が少なくて推定結果はグダグダになったようである。

f:id:MikuHatsune:20171203002118p:image

 

さで、Rstudio アドベントカレンダーでもあるので、Rstan の推定結果をRstudio からRpubs へ出すときの注意点を述べる。Rstudio でknitr すると、library の読み込み時に関連メッセージがでてきて、これがhtml に紛れ込んでしまう。

library のオプションでquietly=FALSE にしても、チャンクオプションでinclude=FALSE にしても出てきてしまうのでイラッ☆とする。

ここで、library で読み込まずに、namespace::function として呼び出すとこの関連メッセージを呼び起こすことなく関数を実行できるっぽいので、こうする。

さて、rstan の実行結果を使ってhtml に計算結果や図を出したいが、rstan のsampling はいちいち実行メッセージを出してしまうので、これもすべてhtml に紛れ込んでしまう。なので、sampling した結果をsave 関数によって保存して、knitr するときにはload することによってknitr でコンパイルするときにはrstan による演算はしないことで回避できるっぽい。

ということで、Rstudio でknitr して、Rpubs にpublish ボタンでアップロードすれば、html やrevealjs がいろいろ作成できる。

JapanR#5 お疲れ様でした。こんな感じでrevealjs のスライドができます。

RPubs - 20171202japanr#5

anime,bd,安済知佳,上田麗奈,大久保瑠美,茅野愛衣,久保ユリカ,小松未可子,高田憂希,竹尾歩美,田中ちえ美,東山奈央,徳井青空,戸田めぐみ,七瀬彩夏,原田彩楓,日笠陽子,本渡楓,水瀬いのり,村川梨衣,山口愛,佳村はるか
ご注文はうさぎですか,11038,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,1,1,0,0
NEWGAME,6452,0,0,0,1,0,0,1,1,0,0,0,1,0,0,1,0,0,0,1,0
サクラクエスト,1165,1,1,0,0,0,1,0,0,1,0,0,0,1,0,1,0,1,0,0,0
うらら迷路貼,694,0,0,0,1,1,0,0,0,0,0,0,0,0,1,0,1,0,0,0,1