Data Blue

データの海で遊んでます。

多腕バンディット問題の解(Dynamic Programming)

多腕(k)バンディット問題をDirect programmingで解きます。
現在の状態sについて、将来の最大期待値であるV(s)は下記のように書け、
下記を満たすようなアームjを選択するのが、次の一手になります。

f:id:mako_shark:20170902025017p:plain

ここで、

V(s)は、全てのアームについてこれまで得られた情報sから考えられる最大期待値

Rj(sj)は今、アームjを選んだ時に直ちに得られる報酬

si (i = 1, 2,..k)は、k個のスロットそれぞれについてのこれまでの情報(成功、失敗)

s'はそれから一つアクションをおこした時の情報(jが選ばれたらs'j)

γは将来に対しての割引率です。

 

この式の大切なところは、V(s)が、これまで得られた情報sと、
一アクション進んだ時の最大期待値であるV(s')で直接計算できることです。

つまりT回引くと最初に決まっていれば、最後の時点での
V(last) =0と置いておけば、その一個手前のV(last-1)は、
どれか一つが選ばれる前の可能性として求められ、、、
と逆算逆算によって、V(s)が直接求められます。

問題は、Tと、kの数によって計算量が指数的に増えていくところです。

これに解決策を与えるために色々な方法があります。

Gittins Indexについて次は書いてみたいと思います。

臨床試験の最適化

医師はこれまでの経験・エビデンスから、
最も良い治療を患者さんに提供することを使命としています。
しかし、治療法AとBがどちらが良いかわからないときは
どうすればよいでしょうか*1

そのような状況に限り、患者さんに最適の医療を提供するという
本来の目的から外れ、治療法AとBを比較するという臨床試験
することが正当化されます。
臨床試験によって得られた情報(例:AがBより有効)は、
同じ疾患を持つ将来の患者さんの治療指針となります。

しかしながら、このメリットを被験者となる目の前の患者さんには
提供することができないというのをいつも申し訳なく思っています。

今回、試験中に得られる結果をシームレスに活用することで、
臨床試験の「比較」性能を担保しながら、「治療」的側面も最大化
するようなデザインについて報告している下記論文

A Bayesian adaptive design for clinical trials in rare diseases - ScienceDirect

を読みました。

 

具体的には、治療結果が分かるごとに
各治療の事後分布から将来的に得られる期待値を考慮して、
期待値が大きいほうに次の人を割り振るという考え方です。
ただし、これを*2完璧にやってしまうと初期の偶然によって
一方の治療法にほとんど人がいなくなってしまい、
そもそもの比較が成り立たなくなる(検出力不足)ので、
それを防ぐために

1.期待値が低い方にも少な目(0.1など)の可能性で割り振る

2.試験終了時までに最低限各群必要な症例数を設定する

という2つの制約をつけることを提案しています*3

割当決定プログラムも論文の参考としてつけられており*4、参考になります。

一番の問題はデザイン上、盲検化が出来ないことです。
その他にもいろいろと制限はありますが、臨床試験を受ける
患者さんにも最大限メリットを提供しようというベースの考え方
自体は大変共感できました。

 

使っているフレームワークは古くからある(多腕)バンデット問題です。
これは複数のスロットマシーンがあり、それぞれの当たり確率が
わからないときに、どういう戦略をとれば一番報酬を大きくできるか
という問題です。
治療法の選択に有効そうですので、これから勉強していきたいと思います。

*1:Aを新規治療、Bをプラセボと読み替えていただいても結構です

*2:Direct programmingで

*3:Constrained randomized direct programming

*4:一か所訂正が必要ですのでもし気になる方は、ご連絡ください

多重検定問題

第一種の過誤(Type I error)はアルファエラー(alpha error)とも呼ばれ、
実際に差が無いものを誤って差があるといってしまうことです。
アワテンボウのエラーと呼ばれる所以です。

よく行われるのが α (有意水準)=  0.05として、p値*1
これを下回ると、2群には差があると結論付けます。

「20回に1回以下しか見られないような現象であれば、偶然と考えるより、
帰無仮説が間違っていると考え、対立仮説を信じよう。」

という宣言が α=0.05 なのです。

ただ、これは逆に言うと、何も差がなくても20回に1回は有意と
なってしまうことを容認しています。
極端な話、独立な検定を100回行えば5回は偶然によって
p値が0.05より小さいものが出てしまうのです。

例えばその中で有意であったものだけ結果として報告されても
読んでいる人はわかりません。

私の受けた統計レクチャーの講師は、
「私はどんな試験でもp値が0.05より小さくなる検定を
見つけることができる!」
と豪語していました*2

このようなことを防ぐために、2つの方法があります。

1.事前に検定する項目を決めておく

恣意的に良い結果が得られたもののみを報告することを防げます。

2.試験全体のαを調整する。

複数の解析を行う際には、解析回数に合わせてα補正を行います。
有名なボンフェロー二の補正は4回解析を行う場合は
「α/4」を検定一回当たりの有意水準とします。

*1:p値:もし帰無仮説が正しい時に、観察された結果以上に極端な結果が得られる可能性

*2:冗談です

False Discovery Rateのメモ

FDR: False Discovery Rateは、q-valueとも呼ばれ、
SNPやMetabolomeなどで多重検定を行いシグナル検出を
するときに、とりあえず可能性のあるシグナルは多めに
拾っておきたいというときに使用されます。

 

例えば q-value < 0.05だとすると
検定で有意としたもののうち約5%は疑陽性であることを
容認するということになります。

 

p-value補正(ボンフェローニなど)では、行う全ての検定を対象
Type I errorの割合をalphaにコントロールするのに対して、

q-value補正では有意と結論付けた検定のみを対象として、
Type I error をalphaにコントロールといったような理解をしています。

 

具体的な計算方法*1は下記
p-valueを小さい順に並べて k 番目のものをp(k)とすると

  • q-value = p(k) * 全検定回数 /  k  < alpha *2

となる最大のkをみつけ、p(1) ~ p(k) までを有意とします。
例えば1000回の検定を行い、p(3) = 0.00012のとき(alpha = 0.05)

  • q-value = 0.00012 * 1000 / 3 = 0.04 < 0.05

となるのでp(3)は (p(1), p(2)が有意であれば)有意です*3

*1:Benjamini-Hochberg法

*2:q * = alpha / 1000 * kとしてp(k)とq*を比べるという方法もありますが結果は同じです

*3:ちなみにBonferonniの補正では、0.00012 * 1000 = 0.12 > 0.05ですので、これは、有意とはなりません

Permutation テスト

遺伝子やメタボロームなどのOmicsデータの解析でよく聞く
パーミュテーションテストについて概要をメモします。

 

結論からいうと、

”知りたい因子の相関係数正規分布が仮定できないときでも使える”

検定です。

 

例えばデータX(5列1000行)が、下記のような構造をしているとします。

  1. 興味のある変数 E
  2. アウトカム Y(連続変数)
  3. 交絡因子3つ V1, V2, V3

次にEとYの相関を調べるため下記モデルを当てはめます

  • Y = b0+ b1*E + b2*V1 + b3*V2+ b4*V3 + error

通常はEの相関係数"b1"の t-test において

p < 0.05(alpha)

であれば、有意にEとYは相関すると結論するわけです。

 

しかし、ここで問題。

t-testの前提として、b1が正規分布する

ことが必要です。

 

 

この前提が成り立たなければどうしましょう。

 

ここで、データXについて、

”Yのみランダムに入れ替えて同じモデルを作ること”

を繰返し行います。

例えば1000回行うとb1も1000個できますね。

この1000個のb1の中で、

”観察データでのb1より極端な値は10回”

ということであれば、p = 10/1000  =0.01となります。

 

ポイントはE, V1-3のデータ構造を保持してYのみ並べ替えるので、
E, V1-3の同士の相関が保たれている点だと思います。

勾配ブースティング(Gradient Boosting )について

前回は、CARTについてまとめてみましたが、
ブースティングの中では勾配ブースティング(gradient boosting)が
特によく用いられます。

【木1+残差→木1+木2+残差】

と残差について(小さな)木を次々に適用していく(Gradient)モデルです。

 

残差は連続数であり、これに回帰木をあてはめていくことから
Multiple additive regression trees: MARTと呼ばれることがあります*1
他にはGradient boosting machine やGeneralized boosting model、
TreeNetなどと記載してあることもあり、混乱させてくれます。

 

アルゴリズムの中でXGBoost *2は、使いやすくて、性能もよく、大変重宝します。

 

XGBoostのpythonインストールに困った方は

mako-shark.hatenablog.com

をご覧下さい。少しでも参考になれば幸いです。

*1:これはSalford Systemのソフトウェアの商標

*2:Extreme Gradient Boosting の短縮

CARTまとめ(decision tree method)

決定木:Decision treeは2種類

  • 分類木:Classification Tree(目的変数がカテゴリー変数)
  • 回帰木:Regresion Tree (目的変数が連続変数)

 

これら決定木を用いた分析が

Classification and Regression Tree: CART分析と呼ばれます。

色々な解析方法がありますが、混乱しやすいので特徴等まとめてみました。

 

ポイントは何をどうアンサンブル*1するかです。

 

バッギングツリー(bootstrap aggregating = bagging trees):

【データをブートストラップ→木作成】を繰り返してそれをアンサンブル。

 

ランダムフォレスト(random forest):

baggingに加え【変数をサンプル→木作成】を繰り返してアンサンブル。

outlierに強い。

 

ブースティング (boosting):

【データにウェイトをかけて木作成→残差でウェイトを更新】を繰り返し、

正確性にウェイトをつけてアンサンブル。

予測性能の悪い変数が多い時でも効率的に良いモデルができる。

 

CARTの発展:

一般的にはSingle tree → Bagging → Random Forest → Boostingと性能が良い。

*1:色々なClassificationを組み合わせることをアンサンブル(Ensemble)と呼びます