bayesmax.sblo.jp Open in urlscan Pro
59.106.18.137  Public Scan

Form analysis
1 forms found in the DOM

Name: form17288845 — GET http://bayesmax.sblo.jp/pages/user/search/

<form method="get" action="http://bayesmax.sblo.jp/pages/user/search/" class="sASform" name="form17288845" id="form17288845">
  <table class="seesaaArticleSearch">
    <tbody>
      <tr>
        <td colspan="5"><input value="" name="keyword" type="text" size="20" class="inputkeyword">&nbsp;<input value="検索" type="submit" class="input-submit"></td>
      </tr>
    </tbody>
  </table>
  <input type="hidden" name="tid" value="seesaa_hotspot"><input type="hidden" name="hid" value="167"><input type="hidden" name="c" value="12"><input type="hidden" name="search" value="1"><input type="hidden" name="ic" value="shift-jis">
</form>

Text Content

スマートフォン専用ページを表示

ベイズマックスへようこそ。
統計に関するあれやこれやを自由気ままに投稿します。
全記事一覧はこちらから。
-  1 2 3 4 5..  次の5件>>


2021年12月25日


10分ぐらいで分かる（？），有意性検定でサンプルサイズを事前に決めなければならない理由

　この記事は，Open and Reproducible Science Advent Calendar
2021の15日目の記事です（空いていたのでタイムリープして書きました）。いろんな人と話をしていると，p値を使って有意性検定をするときになぜサンプルサイズを事前に決めなきゃいけないのかがちゃんと周知されてないような気がしたので，この点について超簡単にまとめたいと思います。端的に言えば，「データの情報を使ってサンプルサイズを決めるのはまずい
(ただし例外もなくはない)」「データを集め切る前に途中のデータを覗くことそれ自体は必ずしも問題にはならない（問題になる場合もある）」というのが結論です。詳しく知りたい人は統計学の教科書を読んでください（『心理学統計法
(放送大学教材 1638)』や『瀕死の統計学を救え！』などが分かりやすい）。

　ここでは簡単のために1標本のt検定 (one-sample t-test) の例を挙げます。母集団の平均 (母平均) が0であるという帰無仮説 (μ = 0)
を棄却したい場合を考えてみます。その場合の手続きは，ごく簡単に書けば以下のようになります。
1. 確率モデルを決める。(e.g., 母集団分布が正規分布に従う。)
2. 標本抽出分布 (sampling distribution; 標本分布と呼ばれることが多い) を導出し，棄却域を決める。(e.g.,
自由度29のt分布の左側2.5%と右側2.5%。)
3. データを観測し，検定統計量の実現値 (定数) を計算して意思決定する。(e.g., 計算して得られたt =
5.62は棄却域に入っているので帰無仮説を棄却する。)
　なお，ここで「データを取得」と言わずに「データを観測」と言っているのは既に収集済みのデータを使って検定を実施する場合もあり得るからです。分野によっては「データを取得」と読み替えても差し支えないと思います。分かりやすさのために，以下の説明では新規にデータを取得する場合を考えます。
　ここで重要なのは，1と2の手続きは研究計画を立てた時点で（つまりデータを観測する前に）実施できるという点です。言い換えれば，データの情報を使わなくても（＝データとは独立に）標本抽出分布と棄却域を決定することができます。頻度主義統計学では確率的に変動するのはデータだけなので，標本抽出分布それ自体と棄却域はデータとは無関係に1つに固定されていることが要求されます。
　標本抽出分布は通常はサンプルサイズの関数であり，また棄却域を求めるためには帰無仮説と対立仮説 (両側か片側かを含む)
と有意水準が定まっていないといけないので，2の手続きを実施するためにはサンプルサイズ・仮説・有意水準を決めておく必要があります。t検定の場合，サンプルサイズが決まらなければ自由度が決まりません。棄却域が仮説と有意水準に依存するのは統計学をちょっとでも学んだことのある人であれば理解は難しくないと思います。
　ここで例えば，「データを1つずつ増やしてその都度検定を実施し，有意になった時点でデータの収集をやめる」という方略をとったとすると，(最終的な)サンプルサイズはデータの関数になってしまいます。その結果，サンプルサイズの関数である標本抽出分布もデータに依存してしまう
(独立でなくなる) ことになってしまうので，普通は検定の前提が満たされなくなります。例えば，t検定の結果が有意になるまでデータを1つずつ増やし続けて最終的にN
= 17で有意になったとします。このN = 17のデータに対するt検定で想定されている標本抽出分布は自由度16 (= 17-1)
のt分布ですが，確率的な標本抽出を行う度に自由度も確率的に揺らいでしまうので (e.g.,
もう1回データを取ったら自由度が23とか14とかになり得る)，この場合の自由度は定数ではなく確率変数です。したがって，本当の標本抽出分布は自由度の分布に依存するよく分からない分布になります
(少なくともt分布にはなりません)。それにもかかわらず無理やりt分布を適用してしまうと，一般には有意水準αが想定よりも高くなり，不当に有意になりやすくなります。サンプルサイズの確率的な揺らぎを考慮して導出した標本抽出分布を使えばこの問題は理論上回避できますが，一般にはかなり難しいです。
　一般に，サンプルサイズがデータに依存するような状況では標本抽出分布は想定通りにならないことが多いです。これこそが，データの情報を使ってサンプルサイズを決めてはいけない主な理由です。ときどき，検定の多重性が理由として挙げられることもありますが
(i.e.,
データを増やす度に検定を実施すると何度も検定をすることになるので有意水準がインフレする；これは結局は標本抽出分布と棄却域の決定に関する問題と本質的に変わらない)，先の説明で言う「データの情報」というのは検定結果だけを指すのではない点には留意が必要かもしれません。例えばデータを1つずつ増やしていき，その度に平均値を計算したり可視化したりして，目視で「そろそろ有意になりそうだな」と思ったタイミングでデータの収集を打ち切った場合でも，当然ながら標本抽出分布は歪みます
(いつ打ち切るかがデータに依存するため)。途中で検定を実施するか否かにかかわらず，通常はサンプルサイズの決定にデータの情報を使ってはいけません。（ただし，先にも述べた通り，サンプルサイズの確率的な変動を考慮した標本抽出分布をきちんと導出することができるのであれば，サンプルサイズがデータに依存する状況でも適切に検定を行うことは可能です。例えば，表が5回出るまでコインを投げ続けて表が出る確率についての検定を行う場合には，コインを投げた回数（＝サンプルサイズ）の標本抽出分布は負の二項分布という単純な分布になることが分かっているので簡単に検定できます。他にも，データを増やす度に有意水準を調整することで有意水準のインフレを防ぐ方法もあったりしますが，あまり現実的ではないように思います。）ちなみに，データを見てから対立仮説
(両側か片側か，片側ならどちらの向きにするか) を決めた場合は棄却域がデータの関数になってしまうのでこの場合も似たような問題が起こります。
　では，データを集め切る前にデータを要約・可視化したり検定したりしてはいけないかというと，実はそうとも限りません。途中の結果がどうであっても最初に決めたサンプルサイズを遵守することを決めていれば，サンプルサイズはデータに依存しないはずなので，標本抽出分布が歪むこともありません。データを集め切る前に，データが問題なく収集できているかを確認したり，検定や可視化のためのプログラムをテストしたりすることは，それはそれで必要な場合も多々ありますので，「データを集め切るまで絶対に中身を見てはいけない！」と強迫的になる必要は必ずしもありません。ただし，途中でデータを見ることによって実験者や調査者の構えが変わり，仮説に合うような結果を誘導してしまうようなバイアス
(いわゆる実験者効果とか要求特性)
が生じ得る場合や無自覚にデータの情報を使ってしまい得る場合にはこの限りではありません。そのようなリスクがある場合には，データを集め切るまで絶対に覗かないと決めておくことも選択肢になり得るでしょう。
　結局のところ，標本抽出分布を用いた検定のロジックをきちんと理解していれば，何が誤用で何が誤用でないかはすんなり分かるんじゃないかと個人的には思っています。特に標本抽出分布は伝統的な推測統計を理解する上で最も重要な考え方だと思うので，統計学をきちんと勉強したい人は何よりもまず標本抽出分布の理解に全力を尽くすのがいいと思います
(それができれば，様々な検定や推定の方法を統一的に理解できるはず……)。

【統計の最新記事】

 * Stanでdrift diffusion..
 * 対数ロジスティック分布をStanに実装す..
 * p値で有意と言えない効果もベイズなら効果..
 * 一般化ガンマ分布をStanに実装する
 * StanでWiener分布からの乱数を得..

posted by mutopsy at 12:21 | 統計


2020年12月23日


STANでDRIFT DIFFUSION MODEL (DDM) を扱うときの実践上のTIPS (1)──まずは直観的に理解する──

　この記事は，Stan Advent Calendar 2020およびベイズ塾Advent Calendar
2020の23日目の記事です。この記事では，認知心理学で有名かつ頑健な現象の1つであるフランカー干渉を例に，drift diffusion model
(DDM) の解説を行うシリーズの第1弾です。数理的な説明というよりは，実践上知っておいたほうが良いと思われる点を解説するつもりです。

※2022/7/16追記：DDMについては，以前「基礎心理学研究」に寄稿したチュートリアル論文でもより詳しく解説しているので，そちらもご参照ください。

　本当はこの記事にもっといろいろな情報を盛り込みたかったのですが，最初の説明だけでそれなりの分量になってしまったので分割することにしました。第1弾である本記事ではDDMを直観的に理解するための解説を行います
(同様の解説記事は既にいろいろありますが，僕なりの視点で解説します)。第2弾ではランダムウォークからWiener過程を作る方法について，第3弾ではDDMのパラメータの単位やスケーリングの影響について，第4弾ではStanでDDMを実装する方法と注意点について解説する予定です。

　ちなみに，DDMについては過去にややマニアックな記事を2本書いていますが (「Drift Diffusion
Modelのパラメータから選択確率を計算する関数：RとStanによる実装」，「StanでWiener分布からの乱数を得る方法（＋累積分布関数ほか）」)，このシリーズは過去の記事とは独立に読めるようにしています。

目次
1. フランカー干渉の説明
2. DDMを直観的に理解する



1. フランカー干渉の説明

　このシリーズではフランカー干渉を例に挙げてDDMの説明をしたいので，まずはフランカー干渉の説明から。フランカー課題と呼ばれる実験課題にはいくつかのバリエーションがありますが，この記事では分かりやすさのために，矢印刺激を用いる課題を取り上げます。この課題の各試行では，図1のように左右どちらかを向いた矢印が複数提示されます。実験参加者の課題は中央にある矢印
(標的) が左右どちらを向いているのかをキー押しで回答することです。周辺にある矢印は無視すべき妨害刺激 (フランカー)
ですが，フランカーの矢印の向きが標的と逆向きの場合 (不一致条件) には，同じ向きの場合 (一致条件)
よりも反応時間が長くなり，誤答も増えるのが一般的です。このような現象をフランカー干渉と呼びます。


図1. フランカー課題の説明。

　この記事の文脈において重要なのは，(1) 実験参加者は提示された刺激を見て「左」か「右」かを回答するということと，(2)
条件の組み合わせは標的の向きが2通り，フランカーの向きが2通りの計4通りであることです。フランカーと標的の向きが一致しているか否かに基づいて分類すれば一致/不一致の2通りとなるので一致・不一致それぞれの平均反応時間と正答率を求めれば良さそうにも見えますが，統計モデリングの文脈では横着せずにまずは丁寧に全ての条件の組み合わせを考えたほうが間違いが起こりにくいです。今回の場合，参加者の反応のレパートリーは「左」「右」の2通りであって，「正解」「不正解」を参加者が選択する訳ではない
(正解/不正解は参加者の反応の結果によって決まる) ので，一致/不一致でまとめて正解/不正解を2値化してモデリングしようとするとおかしなことになりかねません
(一致/不一致で分けてモデル化しても結果として4通りの組み合わせを考えたモデルと等価になる場合はありますが，今回のDDMの例ではそうなりません)。


2. DDMを直観的に理解する

　それではDDMの直観的な説明をします。DDMはよく次のような図で説明されます。


図2. DDMの説明。

この図は，ある条件 (例えば，標的が「右」かつフランカーが「左」)
において実験参加者が反応のために証拠を集積する過程を表しています。横軸は刺激が提示されてからの経過時間を表します。縦軸は，提示された刺激に対して参加者が「右」と反応すべきか，あるいは「左」と反応すべきかを判断するための証拠の量を表しています。今回の例では，縦軸の数値が大きいほど「右」と反応するのに十分な証拠，数値が0に近いほど「左」と反応するのに十分な証拠があると考えてください。

　この例において，刺激が提示された時点では証拠の量の初期値zは真ん中の値 (灰色の破線)
よりも若干高い値ですが，これは，この参加者が始めから「右」と答えやすいバイアスを持っていることを表しています。

　刺激提示から少しの間は証拠の量に変化はなく，一定です。この変化がない期間は非決定時間と呼ばれ，刺激の符号化や反応のための運動 (e.g., キー押し)
といった，左右の判断とは無関係な処理に要する時間を表します。この時間の長さをτと表します。図では非決定時間は刺激提示の直後に描かれているのに，判断が終わった後に行われるはずのキー押しの時間が非決定時間に含まれるということに違和感を覚える方もいるかもしれませんが，実は非決定時間はどの位置にあっても数学的には全く同じです。例えば刺激提示直後と反応の直前の2か所に2つの非決定時間が別々に存在すると考えても問題ありません
(その場合，両者の和がτで表されます；その内訳は通常は識別不可能です)。説明の便宜のために，ここでは刺激提示直後のみに非決定時間が挿入されると考えることにします。

　非決定時間が過ぎると証拠の集積が始まります。図に描かれている4本のうねうねは各試行に対応しています
(1試行あたり1うねうね)。図では証拠の量がとても荒ぶっていますが，このうねうねを生み出すのがDDMの肝であるWiener過程です
(次回の第2弾でもう少し詳しく説明します)。この例では証拠の量が見かけ上不規則に変動しているように見えますが，平均的にはだんだん値が大きくなる傾向があるように作っています。この平均的な傾向すなわち傾きをドリフト率と呼び，δで表します。今回の例ではδが正なので，うねうねしながらも「右」と反応する方向に証拠が集まりやすい傾向があることが分かります。逆に，δが負であれば「左」と反応しやすくなります。つまり，ドリフト率δは証拠の集積の速さ
(絶対値の大きさ) と方向 (符号) を表すパラメータとして解釈できます。

　実験参加者は証拠を蓄積するだけではなくて，蓄積された証拠に基づいて最終的に左か右のキーを押さなければなりません。そこで，証拠の量がある閾値に達したときに反応が生成されることを仮定し，証拠の量がα以上になったら「右」，0以下になったら「左」と反応すると考えます。下側は0で固定なので，上側のαだけを決めれば閾値が確定します。αの値が小さければ証拠の量が不十分でも反応が起こり，逆にαの値が大きければ十分な証拠が集まるまで反応が起こらないことを意味するので，心理学的にはαは慎重さを表すパラメータとして解釈することができます。なお，先ほど説明した初期値zに関しては，下側の閾値
(0) と上側の閾値 (α) の間の相対的な位置として表現したほうが解釈をする上で便利なので，(0,1)の範囲をとるパラメータβを使ってz =
αβと表します。β = .5のとき反応バイアスなし，β > .5のとき上側 (この場合は「右」反応) へのバイアスあり，β < .5のとき下側
(この場合は「左」反応) へのバイアスありと解釈できます。

　このようなプロセスに基づいて反応が生成されると仮定すると，反応時間は図2の上側と下側に描かれているような確率分布に従います。上側は「右」と反応したときの反応時間分布，下側は「左」と反応した時の反応時間分布で，これら2つの分布の面積の和は1となります。もし正答が「右」であったなら，上側の確率分布の面積は正答率を，下側の確率分布の面積は誤答率を表します。このようにして，選択された反応
(「左」か「右」か) と反応時間の同時分布を表現できるのがDDMのいいところです。ちなみに，この同時分布 (図2の上側と下側に描かれている2つの分布の組)
はWiener過程にちなんでWiener分布と呼ばれます。このモデルはWiener過程が最初に閾値に達するときの時間の分布を表すのでWiener First
Passage Time分布と呼ぶ方がより正確かもしれません (Stan Functions Referenceではそう表記されています)。

　この記事は以上です。次回はWiener過程についてもっと詳しく説明します。

タグ：実験心理学 Drift Diffusion Model Wiener分布

posted by mutopsy at 14:03 | 統計


2020年12月21日


対数ロジスティック分布をSTANに実装する

　この記事は，Stan Advent Calendar 2020の21日目の記事です。この記事では，生存時間解析 (survival analysis;
継続時間解析とも)
でよく使われる対数ロジスティック分布を使ったモデリングをStanで行う方法について解説します。ちなみに，一般化ガンマ分布の記事も以前書きました。この記事を読んで得られる学びとしては，(1)ロジスティック分布と対数ロジスティック分布の関係が分かる，(2)自作関数をStanで作る方法が分かる，(3)対数〇〇分布の代わりに対数変換したデータ＋〇〇分布を使ってモデル化するとパラメータの推定はうまくいくがWAICの計算を間違えてしまうことがある，あたりが挙げられるかもしれません。

目次
1. ロジスティック分布と対数ロジスティック分布
2. 対数ロジスティック分布をStanに実装
3. 架空のデータを使って走らせてみる



1. ロジスティック分布と対数ロジスティック分布

　対数ロジスティック分布は，ある確率変数Xの対数ln(X)がロジスティック分布に従うときに，元の確率変数Xが従う分布です。逆にいえば，ロジスティック分布に従う変数を指数変換した変数は対数ロジスティック分布に従います。正規分布と対数正規分布の関係と同じですね。ロジスティック分布の台
(実現値が取り得る範囲)は全ての実数ですが，対数ロジスティック分布の台は正の実数です。それゆえに，対数ロジスティック分布は非負の値となる生存時間や継続時間の確率分布の候補になる訳ですね。ロジスティック分布と対数ロジスティック分布の見た目はこんな感じです。



という訳で，まずはロジスティック分布の説明から。ロジスティック分布の確率密度関数は



と表されます。μは位置パラメータ，σ (> 0) はスケールパラメータです。ちなみに，累積分布関数は



と書けます。この式に見覚えのある方も多いかと思いますが，μ = 0，σ = 1のとき (標準ロジスティック分布)，累積分布関数は



となりロジスティック関数と一致します。これがロジスティック分布という名前の由来ですね。標準ロジスティック分布とロジスティック関数の関係は，標準正規分布とその累積分布関数の関係と同じです。ロジスティック分布はStanではlogistic_lpdf()などの関数で実装されています
(サンプリングステートメントを使って「Y ~ logistic(mu, sigma);」と書ける)。

　ある確率変数Xがロジスティック分布に従うとき，exp(X)は対数ロジスティック分布に従います。スケールパラメータα(> 0)と形状パラメータβ(>
0)を使って，対数ロジスティック分布の確率密度関数は



累積分布関数は



と表されます。α = exp(μ)，β =
1/σを代入すればμとσを使った表現もできますが，αとβを使った方がシンプルな式で表せます。どちらのパラメータ化を使うかは目的によって変えると良いでしょう。


2. 対数ロジスティック分布をSTANに実装

　対数ロジスティック分布の関数はStanにはデフォルトで用意されていないので，(1) αとβでパラメータ化した対数ロジスティック分布を自分で定義する方法，(2)
μとσでパラメータ化した対数ロジスティック分布を自分で定義する方法，(3)
ロジスティック分布を利用する方法，の3通りの方法で実装してみます。3番目のロジスティック分布を使ったやり方は一番簡単ですが，後ほど示すように，WAICの計算などのために尤度を使う必要があるときには正しい結果を返さないので注意が必要です。結果を比較できるように，これらを1つのStanコードにまとめて書きます。WAICも計算します。Stanコードは以下の通り。

Stanコード: loglogistic.stan

functions{
  //1. αとβでパラメータ化
  real loglogistic1_lpdf(real x, real scale, real shape){
    real lpdf;
    lpdf = log(shape) - log(scale) + (shape-1)*log(x/scale) - 2*log(1+(x/scale)^shape);
    return(lpdf);
  }
  
  //2. μとσでパラメータ化
  real loglogistic2_lpdf(real x, real mu, real sigma){
    real lpdf;
    lpdf = -log(sigma) - mu + (1/sigma-1)*(log(x) - mu) - 2*log(1+(x/exp(mu))^(1/sigma));
    return(lpdf);
  }
}

data{
  int N;
  vector[N] Y1;
  vector[N] Y2;
  vector[N] Y3;
}

parameters{
  vector[3] mu;
  vector<lower=0>[3] sigma;
}

transformed parameters{
  vector<lower=0>[3] alpha;
  vector<lower=0>[3] beta;
  
  for(j in 1:3){
    alpha[j] = exp(mu[j]);
    beta[j] = 1/sigma[j];
  }
}

model{
  //prior
  mu ~ uniform(-10,10);
  sigma ~ uniform(0,10);
  
  //model
  for(i in 1:N){
    //1. αとβでパラメータ化
    Y1[i] ~ loglogistic1(alpha[1], beta[1]);
    
    //2. μとσでパラメータ化
    Y2[i] ~ loglogistic2(mu[2],sigma[2]);
    
    //3. ロジスティック分布を利用
    log(Y3[i]) ~ logistic(mu[3],sigma[3]);
  }
}

generated quantities{
  vector[N] log_lik_1;
  vector[N] log_lik_2;
  vector[N] log_lik_3;
  
  //WAICの計算のために対数尤度を計算する
  
  for(i in 1:N){
    //1. αとβでパラメータ化
    log_lik_1[i] = loglogistic1_lpdf(Y1[i] | alpha[1], beta[1]);
    
    //2. μとσでパラメータ化
    log_lik_2[i] = loglogistic2_lpdf(Y2[i] | mu[2],sigma[2]);
    
    //3. ロジスティック分布を利用 (間違った結果になります)
    log_lik_3[i] = logistic_lpdf(log(Y3[i]) | exp(mu[3]), 1/sigma[3]);
  }
}



　データとして与えたY1，Y2，Y3に対し，それぞれ(1)～(3)の方法で対数ロジスティック分布を当てはめてパラメータを推定するコードです。4つの3次元ベクトルmu,
sigma, alpha,
betaの各要素にそれぞれの方法で推定されたパラメータの値が格納されます。functionsブロックでは，2通りのパラメータ化の方法で対数ロジスティック分布の対数確率密度関数
(lpdf) を定義しています (上述した確率密度関数の対数をとって整理したもの)。WAICを計算するために，generated
quantitiesブロックで対数尤度の計算も行っています。ここではあえてmodelブロックと対応した書き方をしています
(つまり，Y1～Y3を使って3通りの方法で対数尤度を計算している)。Y1～Y3を全て同じデータにしておけば，それぞれの方法の結果を比較することができます。


3. 架空のデータを使って走らせてみる

　それでは，Rで架空のデータを作って実際にパラメータを推定してみましょう。真値がμ=0.5 (α = 1.65)，σ = 2.0 (β = 0.5)
である対数ロジスティック分布からの乱数を5,000個発生させてみます。対数ロジスティック分布からの乱数を得るには，rlogis()で得たロジスティック分布からの乱数をexp()で指数変換するか，actuar::rllogis()などのパッケージに含まれる関数を使えばokです。Rコードと推定結果は以下の通り。

Rコード: loglogistic.R

library(tidyverse) #version 1.2.1
library(rstan) #version 2.19.2
library(loo) #2.3.1

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

## 架空のデータを作る ----

N <- 5000
mu <- 0.5
sigma <- 2.0

set.seed(8823)

Y <- rlogis(N, location = mu, scale = sigma) %>% exp()

#もしくは
#library(actuar)
#Y <- rllogis(N, shape = 1/sigma, scale = exp(mu))

## Stanに渡す ----
sm <- stan_model("loglogistic.stan")

warmup <- 1000
iter <- 2500 + warmup
t0<-proc.time()
fit_loglogis <- sampling(sm,data=list(N=N,Y1=Y,Y2=Y,Y3=Y)
                ,seed=8823,chains=4,iter=iter,warmup=warmup
)

#saveRDS(fit_loglogis,"fit_loglogis.rds")

## 推定結果 ----
s_loglogis <- summary(fit_loglogis)$summary
round(s_loglogis,3)[1:12,]

#           mean se_mean    sd  2.5%   25%   50%   75% 97.5%    n_eff Rhat
# mu[1]    0.534   0.000 0.048 0.440 0.503 0.535 0.566 0.628 13459.11    1
# mu[2]    0.534   0.000 0.049 0.437 0.500 0.534 0.567 0.629 13529.81    1
# mu[3]    0.534   0.000 0.049 0.438 0.501 0.533 0.567 0.630 12189.30    1
# sigma[1] 1.991   0.000 0.023 1.946 1.975 1.990 2.006 2.036 13718.75    1
# sigma[2] 1.991   0.000 0.024 1.945 1.975 1.991 2.007 2.038 12643.81    1
# sigma[3] 1.991   0.000 0.023 1.945 1.975 1.990 2.006 2.037 14590.76    1
# alpha[1] 1.708   0.001 0.082 1.553 1.653 1.707 1.762 1.874 13440.38    1
# alpha[2] 1.707   0.001 0.084 1.548 1.649 1.706 1.763 1.877 13436.23    1
# alpha[3] 1.708   0.001 0.084 1.549 1.650 1.705 1.763 1.878 12220.63    1
# beta[1]  0.502   0.000 0.006 0.491 0.498 0.502 0.506 0.514 13704.58    1
# beta[2]  0.502   0.000 0.006 0.491 0.498 0.502 0.506 0.514 12601.31    1
# beta[3]  0.502   0.000 0.006 0.491 0.498 0.502 0.506 0.514 14593.66    1

## WAICの計算
e_logliks <- rstan::extract(fit_loglogis, pars = c("log_lik_1", "log_lik_2", "log_lik_3"))
waic(e_logliks$log_lik_1) # WAIC = 32194.3 (SE = 527.7)
waic(e_logliks$log_lik_2) # WAIC = 32194.5 (SE = 527.8)
waic(e_logliks$log_lik_3) # WAIC = 53456.1 (SE = 662.6)



　パラメータの推定結果はどの方法でもほぼ同じ値でありかつ真値に近いので，うまく推定できていると見てよさそうです。一方，WAICの結果を見ると，方法(1)と方法(2)ではほとんど値が同じですが，ロジスティック分布を利用した方法(3)だけ結果が大きく異なっています。同じデータ・同じモデルなのにWAICが一致しないのはおかしいですね。これは，対数ロジスティック分布を使ったモデルではYの確率密度を使って尤度が計算されるのに対し，ロジスティック分布を使ったモデルではlog(Y)をデータとして与えたときの確率密度を使って尤度が計算されてしまうことに由来します。与えているデータが違うのだから，結果が変わるのは当たり前ですね
(もっとシンプルな例として，データのスケールを変えると確率密度や尤度が変わる，ということを考えると分かりやすいかもしれません)。サンプリングするときは方法3のようにlog(Y[i])
~ logistic(mu, sigma);
という書き方をしても問題ありませんし，おそらくその方が自作関数を使うよりも計算が高速で収束も安定すると思われますが，尤度を正確に計算する必要があるときには書き方に気を付けましょう。これは対数正規分布を使うときも同様で，サンプリング時にはlog(Y)
~ normal(mu, sigma); と書いてもYが対数正規分布に従うことを表せますが，尤度を計算するときにはlog_lik[i] =
lognormal(Y[i] | mu, sigma); という書き方をしなければなりません。Be careful and enjoy!

タグ：R STAN 生存時間解析 ロジスティック分布 対数ロジスティック分布

posted by mutopsy at 00:00 | 統計


2020年12月20日


心理学の再現性と科学性と「ベイズ」に関する自分の立ち位置を整理する──JPA2020再現可能性シンポジウムのスライドを添えて──

　この記事は，Open and Reproducible Science Advent Calendar 2020およびベイズ塾 Advent Calendar
2020の20日目の記事です。2020年9月9日に行われた，日本心理学会第84回大会 (JPA2020)
の若手の会企画シンポジウム「若手が聞きたい再現可能性問題の現状とこれから」で発表したスライドの紹介と裏話 (経緯の話とか)
をしつつ，再現性に関する現時点での自分の立ち位置について整理をしたいと思います。最後に，自分がなぜ「ベイズ」を推すのかについても少し。当初は「スライド紹介＋おまけ」という感じの短い記事にするつもりだったのですが，書き始めたら止まらなくて，最終的にはスライドがおまけになってしまいました。なお，シンポジウムの動画は若手の会のWebサイトにまとめられています。また，同シンポジウムで登壇された山田先生のスライド
(紹介記事)，平石先生・中村先生のスライドもSlideShareに公開されています。
※12/20 12:40ごろ記事タイトルを若干修正しました。

目次
1. スライドの紹介と裏話
2. 再現性に関する自分の立ち位置
3. 再現性・科学性と「ベイズ」
4. まとめ



1. スライドの紹介と裏話

まずはオープンデータから！──高めよう信用性，広めよう二次分析── by @mutopsy

　心理学においてオープンデータを実践することの意義について，(1) 科学的知見の信用性 (credibility) の向上と (2)
データの二次分析の促進という2つの側面からお話をしました。オープンデータにしないことの弊害というよりは，オープンデータを実践することで得られるメリットについて，具体的かつ説得的かつポジティブに説明することを目指しました。

　僕が再現性に関するトークをしたのは今回が初めてでした。グリーンバックなのにバーチャル背景を使わない某Y先生にお声かけ頂いたことがきっかけです。再現性に関する話題にはもともと関心があり，いろいろな情報を聞きかじってはいたものの，こういう場で話せるほど引き出しがある訳でもないし，もっと他に適任がいるだろうと思っていましたが，これを機にちゃんと勉強するのも悪くないかなーということで勢いで殻を破って引き受けることにした次第です。それで何の話をしようかと思ったときに，自分がまともに実践しているのはオープンデータと統計学の勉強とプレプリントの公開ぐらいしかなかったので，せっかくだしオープンデータの意義についてちゃんと考えるかーと思ってこんなトークになりました。この手の話題に関連する論文やプレプリントはTwitterのタイムラインでよく流れてくるので，文献を探すのは比較的楽でした。いい時代ですね。といってもトーク内容はほとんど自分の経験ベースの考察でしたが，企画の趣旨からして若手研究者による実践例みたいなのが求められていたように思うのでこれで良かったかもしれません。この発表のお陰で，自分がもともと実践していたオープンデータの営みにちゃんと意味があることを再確認できたし，今までのやり方では不十分であったことなども分かり
(例えば，データだけでなく分析スクリプトも共有したほうがいいこととか)，とても良い学びになりました。せっかく殻を破ったので，これからも再現性についての情報発信をしていきたいなーと思っています。


2. 再現性に関する自分の立ち位置

　このシンポジウムは再現性に関する自分の立ち位置を整理するのにもとても役立ちました。自分のこれまでを振り返ってみると，実は「再現性問題」が明るみになる以前から僕は再現性に関心があったんだなーと気づかされました。学部時代から僕は心理学の科学的方法論に関心があって，学部入学当初はむしろ
(無知ゆえに)
心理学という学問を疑ってかかっていたからこそ必然的に，心理学の科学性の拠り所である研究法や統計学についてそれなりにちゃんと勉強できたという経緯があります。当時の指導教員であった千野直仁先生の数理的・論理的な厳密さへのこだわりや科学観にもすごく影響を受けています。千野先生は当時からオープンデータの重要性を訴えられており，また統計学の誤用に関する啓発論文も執筆されていたので，心理学の科学性や再現性は僕にとって当時から身近なトピックでした。千野先生の先見の明は流石としか言いようがないです。

　僕は心理学には科学であってほしいと思っています
(全ての心理学がそうあるべきと思っている訳ではありませんが)。なので，心理学者が自ら「心理学は科学ではない」という旨の発言をするのを聞くと悲しい気持ちになります。ただ，僕は心理学における「再現性問題」をどうにかしたい，というような使命感はあまり持っていなくて，どちらかといえば，自分と自分の周り
(近接分野) の研究の科学らしさ (再現性を含む)
さえ高められればいいかなと今は思っています。心理学は広大であり良くも悪くも一枚岩ではないので，自分の狭い見識だけに基づいて全体についてどうこう言うのは無理筋だし，別の分野で起こった問題について外側から首を突っ込むのも正直面倒だし
(自分の分野を棚に上げている訳じゃなく，いちいち他の分野のことまで気にしていられない)，そもそも自分はそういうことができる立場にもいないので
(制度そのものを整えるとか)，とりあえず目の前の自分の仕事がちゃんとできればいいんじゃないかと思っています。実際，僕は科学の俎上に乗りやすそうなテーマ
(心的回転とか) ばかりを意図的に選んで研究しているので
(自分のテーマも再現性問題と無縁ではありませんが)，厳密さを妥協せざるを得ないこともありうる他の研究テーマに対してまで厳密な方法を強要するようなことはとてもできません。そういう意味では，僕は「再現性問題」という心理学全体や他の分野にも波及した時事ネタ自体に興味があるというよりは，科学性・再現性を高める方法論とその実践にもともと関心があったと言うべきかもしれません
(もともと疑ってましたし)。もちろん，再現性問題をきっかけに明るみになったこともたくさんありますし，再現性問題それ自体を軽視している訳では決してないのですが，再現性問題という言い方をしてしまうとQRPsのようなネガティブな面が強調されやすくなって，これまで誠実にやってきた研究者からの反感をいたずらに買ってしまうリスクがあるんじゃないかと危惧しています
(危機感を煽りたい気持ちも分かりますが，やりすぎは逆効果だと思います)。そういう視点も大事ですが
(現にそれがなければこんなに盛り上がらなかったでしょう)，再現性を高めるという意識を広く共有したいなら，もうちょっとポジティブな視点で語る人が増えてもいいんじゃないかなあと思っています。再現性は高ければ高いほど良いに決まっていますし，そのことに反対する研究者はいないでしょうから。


3. 再現性・科学性と「ベイズ」

　この記事はベイズ塾Advent
Calendarの記事でもあるし，そもそもこれは統計ブログなので，ベイズの話もしましょうか。再現性の文脈では「ベイズ」の話題もよく上がります。p値を使った帰無仮説検定の代替案としてベイズ的な仮説評価が有用である，といった内容であることが多いように思います。これについてもたくさん議論がありますが，僕自身はベイズと再現性を直接的に結びつけるようなことは普段はあまり考えていません
(後述するように，結果的に再現性の向上に繋がることは十分ありうると思っていますが)。再現性を高めたくてベイズに傾倒しているという訳ではないんですよね。p値をやめてベイズにすれば再現性が上がる，という単純な話でもなくて，使い方を間違えればむしろ再現性を損なうリスクもあります
(前回の記事「p値で有意と言えない効果もベイズなら効果があると言える？──事後分布に基づく仮説評価について──」も参照)。p値だってちゃんと使えばとても役立ちますし，なんなら僕は修論でも博論でも古典的な分析しかしていません。どちらかといえば，僕がベイズを推す理由は再現性のためというより科学としてのスコープを広げるためと言ったほうが良いかもしれません。僕は物理学や科学史が好きで，物理学
(特に古典力学)
のような数理モデルを心理学にも適用したいということを学部入学当初から考えていたのですが，大学院に入ってベイズ統計モデリングと出会ったことでまさにそれに近いことができるようになりました。幸せです。

　統計モデリングの利点については，心理学ワールド86号の「この人をたずねて」に寄稿した清水先生へのインタビュー記事で清水先生がいみじくも説明されている内容にとても納得しており，また第3回犬山認知行動研究会議
(2020/1/11)
では当時の自分の考えを整理して発表したこともありますが，せっかくの機会なので改めて考えをまとめてみます。実験計画法＋有意性検定という伝統的な研究法では効果の有無や関連の有無といった質的な評価はできますが，関数形の数理表現のような量的な説明・予測は難しく，また，研究が進むにつれてどんどん要因が増えてどんどん小さな効果を検証する方向に進みがちなので，それが再現性の低下を招いてしまう面もあると思います
(小さな効果を検出するには大サンプルが必要なので追試のコストも上がる)。こういう事情を考えると，頑健な現象の背後にあるメカニズムを統計モデリングによって検証するアプローチを取り入れることで結果的に再現性の向上に繋がることはあるだろうと考えています。一方で統計モデリングも万能ではなくて，分析するデータに含まれていない情報はモデル化できないし
(あるパラメータを入れたくてもそれを識別できる情報がなければ推定できない)，モデルに含まれなかった重要な変数の存在を見落とすリスクもある上に (e.g.,
疑似相関)，モデル内のパラメータの妥当性 (実在性や解釈性を含む)
の確証はモデリングだけではできないことが多いです。だからこそ，それぞれの方法で手の届かないところを補うために，従来の実験計画法＋有意性検定と統計モデリングの両方を組み合わせたり，両者を交互に繰り返したり
(e.g., パラメータを操作する実験を行って条件間差を見る)
することが必要になるのではないかと思います。統計モデリングはベイズじゃなくてもできますが，ベイズには推定の容易さであったり優れた仮説評価指標を使えたりといった利点があるので個人的にはベイズ推しです。

　あと，実験計画法＋有意性検定だけだとデータが持っている情報のほとんどを捨ててしまうことになってもったいないなと思います。それだけで終わらずに，統計モデリングや探索的アプローチなどの二次分析をしたほうが，有限の資源を有効活用できて無駄がありません。まだ食べられるデータがあるならぜひオープンにしてください。骨までしゃぶりつくします。


4. まとめ

　存外長くなってしまいましたが，まとめると，僕の関心は再現性に特化しているというよりは，心理学における科学的方法論とその実践に興味があるということです。自分の研究が依拠している研究法について無知なのは怖いですし，「よい」科学的方法論に根差した研究であれば再現性も高くなるはずですしね。科学的方法を洗練させていく営みはとても健全かつ前向きな営為なので，無駄に悲観的になる必要はないと思います。考えるだけじゃなくてちゃんとその方法を自分で実践することで規範になれたらもっといいですね。僕はこれまでいろんな研究テーマや方法論
(ベイズとかGLMMとか多変量解析とか，実験とか調査とか疫学研究とか二次分析とか)
を扱ってきましたが，こうやって考えればこれまでの自分の研究のスタンスをだいたいうまくまとめられる，ということに気付いたのが今年の大きな収穫でした。

タグ：再現性 オープンデータ ベイズ スライド紹介

posted by mutopsy at 00:00 | 雑記


2020年12月10日


P値で有意と言えない効果もベイズなら効果があると言える？──事後分布に基づく仮説評価について──

　この記事は，Open and Reproducible Science Advent Calendar 2020およびベイズ塾 Advent Calendar
2020の10日目の記事です。この記事では，ベイズ統計学と再現性に関わる次のようなケースについて考えてみます。

　研究者のAさんは，ある実験操作をしたときとしなかったときで変数Vの平均値に差が生じるかどうかを調べるために，60名の実験参加者を30名ずつ無作為に統制群と実験群に割り当てて実験データを取得しました。変数Vの平均値を計算してみたところ，統制群は-0.028，実験群は0.348であったので，実験操作によって変数Vの平均値が数値的には増加しているように見えます(図1)。


図1. Aさんが収集したデータ。黒い横線は各群の平均値。

　ところが，この変数VについてWelchのt検定 (両側検定，α = .05) を行ったところp =
.163という結果となり，有意差は認められませんでした。諦めきれなかったAさんはふと，ベイズで分析し直せば実験操作の効果があると主張できるのではないかと思い立ち，Welchのt検定に相当するモデルを使って平均値差μdiff
(= μ実験群 - μ統制群) の事後分布を推定しました(図2)。


図2. 2群の平均値差 (μdiff) の事後分布。

　この事後分布を利用して，μdiffが0よりも大きい確率 (図2の赤い領域の面積) を計算したところ，p(μdiff > 0) =
91%という高い値が得られました。この結果をもってAさんは，「この実験操作は変数Vを高める効果があるということが高い確信を持って言える」と結論づけました。

　この例は筆者による創作でありフィクションですが，これに類する報告を学会や研究会，論文等でときどき見聞きします。このような推論には問題があるのですが，無自覚にこれを行ってしまっていると思われる実例としばしば遭遇してしまうので，この状況はまずいと思い，この推論の何が問題なのかについて自分なりに整理してまとめようと思いました。あくまでも私見なので，鵜呑みにするのではなく参考意見として受け止めて頂ければと思います。ちなみに，Aさんの実験データ
(架空) の生成と分析を行うR・Stanコードはこの記事の末尾に掲載しています。

目次
1. ベイズの枠組みにおける仮説評価の方法
2. 何が問題なのか
3. 他の指標（95%確信区間・ベイズファクター）との比較
4. まとめ
5. RとStanのコード



1. ベイズの枠組みにおける仮説評価の方法

　いわゆるベイズ的な仮説評価には複数の方法がありますが，研究者によって推奨する方法が異なっており，議論もたくさん行われています（この記事ではそれぞれの長短について深入りはしません）。仮説評価の方法は大雑把に，モデル比較に基づく方法とパラメータ推定に基づく方法の2種類に分けることができます。モデル比較に基づく方法にはベイズファクター
(周辺尤度・自由エネルギー) やWAICによるモデル比較が代表的です。先ほどの例に出てきた事後分布を使った仮説の評価はパラメータ推定に基づく方法に含まれます
(他にもROPEを使った方法や95%確信区間が0をまたぐかどうかを基準にする方法などがあります；ついでに，ROPEを利用したベイズファクターの計算方法も存在します)。Makowski
et al. (2019) は，事後分布においてパラメータが中央値と同じ符号になる確率 (つまり，0を下回る確率と0を上回る確率のうち大きいほう)
のことをprobability of direction (pd)
と呼んでいますが，先ほどの例はこのpdを使った仮説評価を行っていることと等しいです。pdは[.5, 1] の範囲を取り得る指標で，値が大きいほど効果の存在
(effect existence)
の強い証拠とみなすことができます。一方，値が小さかったとしても効果がないことの証拠とはみなせないという非対称性があります。Makowski et al.
(2019) も言及しているように，(少なくとも単純なモデルにおいては) pdはp値と一対一に対応します (効果の事前分布に(-∞, +∞)
の範囲の非正則な一様分布を仮定したとき，pd = 1 –
p/2が成り立つ；ただし，この記事の分析では各群の平均の事前分布に一様分布を仮定しているためμdiffの事前分布は三角分布となり，この式とは一致しない)。よくよく考えてみたら，上述したpdの性質はp値の性質と何ら変わりませんよね。帰無仮説が正しい時
(真の効果の大きさが0のとき) にはサンプルの変動によってp値は (0, 1)の範囲の一様分布に従いますが，pdも同様に(.5, 1)の範囲の一様分布に従います
(事前分布にもよりますが)。そう考えると，p = .163で有意じゃないのにpd = .91で効果があるように見えるのは不思議な気がしませんか？


2. 何が問題なのか

　Aさんのpdの使い方には，少なくとも2つの問題点があるように思われます（もっとあるかもしれませんが，ここでは2つに絞ります）。1つ目は，pdが示す確率の値を素朴に解釈してしまっている点でしょう。例えば，天気予報で降水確率が80%であることを知れば，多くの人は傘を持って出かけるでしょう。降水確率が90%なら尚のこと，間違いなく雨が降るだろうと確信することでしょう。僕だったら降水確率が50%でも傘を持って出かけます（何なら毎日折り畳み傘を持参していますけど）。これと同じように，pd
= .91を根拠に「実験操作の効果がある」と判断するのは妥当でしょうか。例えば，100回同じ実験をしたら91回はμdiff > 0となる，といった具合に，
pdが実験の再現率を表すと解釈できるのであればこの判断は正当化できるかもしれませんが，残念ながらそうではありません
(実際，今回の設定においては100回同じ実験を繰り返した時にpd >
.90となるのはたった20回だけです)。そもそも，pdは必ず.5以上の値を取りますしね。パラメータの事後分布とそこから計算された指標は，あくまでも今回のサンプルで条件付けられた結果に過ぎません。言い換えれば，別の実験参加者60人に対して同じ実験を行った場合には異なる事後分布が得られるということです
(あるいは，今回の実験参加者のうち10名がたまたま何らかの事情で実験に参加できず，たまたま都合の付いた他の10名が実験に参加していた平行世界があったとしたら，異なる事後分布が得られていたでしょう)。たまたまサンプルが偏ってしまった場合には，真の効果の大きさが0であっても大きなpdが得られてしまうことは十分あり得ます
(偶然p < .05になり得るのと同じように)。このように，pdはそもそも直感的に解釈するのが難しい指標であるように思われます
(少なくとも僕にとっては)。少なくとも降水確率のように素朴に解釈してしまうのはまずい，というのが1点目の問題でした。

　2点目の問題は，基準を明確にせずに効果の有無について二値判断をしてしまっている点です（i.e., 「この実験操作は変数Vを高める」）。今回はpd =
.91だったので素朴に確率を解釈すれば効果がありそうな感じがしますが，pd = .84だったらどうでしょうか。じゃあpd =
.73だったら？二値判断をせずに，pdの値をそのまま量的に解釈するのであればこの点は問題にならないのですが，前述した通りpdはそれほど直感的に解釈しやすい指標とは言えず，量的な解釈は意外と難しいように思います。それに，効果の存在を示すことはなんやかんやで学術的にも社会的にも求められることが多く，二値判断せずに論文を書くことは現状かなり難しいと思います。たとえ留保付きだとしても，効果の有無についての言及はほとんどの場合避けられないし，意図せずとも暗黙に二値判断をしてしまうこともあるし，論文で明示しなかったとしても読者が勝手に効果の有無を判断してしまうこともあり得えます
(余談ですが，ベイズファクターによるモデル比較の場合にも，「モデルXがモデルYよりも支持された」といった具合に離散的・質的な判断が行われることがほとんどだと思います)。それならばせめて，事前に「pdが〇〇以上であれば効果があったと判断する」という基準を決めておかないとまずい訳です。残念ながら，そういった基準を明示せずにpdの結果を見て「効果があった」と主張している報告を複数回観測したことがあるのですが……
(しかもpdの値が.80ぐらいのものもありました)。帰無仮説検定の問題点として二値判断である点が批判されることは多いですが，基準を明確にせずに事実上二値判断をするのに比べたら，α
= .05と事前に決めた上でp値を使って二値判断をするほうが遥かに健全だと思います。pdの基準については特別な理由やこだわりがないのであれば両側検定のα =
.05に相当するpd > .975を採用するのが無難だと思います。※もちろん，二値判断せずにpdを量的に解釈するという方針を採用するのであれば
(それができるのであれば)，このような基準を決める必要はありません。ここで批判の対象はあくまでもpdを恣意的に解釈して二値判断する行為に対するものであり，pdそのものを批判している訳ではありません。


3. 他の指標（95%確信区間・ベイズファクター）との比較

　他の指標を使ってAさんの仮説を評価したらどうなるでしょうか。試しに確信区間を計算してみたところ，μdiffの95%確信区間は[-0.177,
0.925]であり，0をまたいでいるようです。95%確信区間が0をまたがないことはpd > .975と同値なので，pd =
.91のとき0をまたぐのは当たり前ですね。

　続いてベイズファクターではどうでしょうか。ベイズファクターは事前分布に依存しますが，今回は標準化効果量の事前分布をスケールパラメータr =
1/√2のコーシー分布としました (JASPのデフォルトの設定と同じ；Morey & Rouder, 2011およびRouder et al.,
2009も参照)。帰無仮説 (効果ゼロ; μdiff = 0) の周辺尤度を分母とするベイズファクターを計算したところBF =
0.605となり，どちらかといえば帰無仮説に軍配が上がるが今回のデータからはどちらともいえない，という結果が得られました
(Jeffreysの基準でいえばanecdotalな証拠)。なんにせよ，今回のAさんのデータから「効果がある」と結論づけるのは難しそうですね。
※ここでの分析方法の詳細が気になる方は後述のR & Stanコードを参照してください。


4. まとめ

　p値を使ったら有意にならないけどベイズなら効果があると主張できる，という幻想をぶち壊したくてこの記事を書きました。少なくとも僕は経験上，他の条件がほぼ同じであるにもかかわらずp値では有意にならないがベイズなら効果があると言えるケースに遭遇したことは一度もありません（もしそういう事例があったら教えてほしいです）。そういうことが起こった場合には，分析の設定か解釈を誤っていることを疑ったほうが良いと思います。例えば事前分布を不適切に設定した場合，ベイズファクターが誤ったモデルを選択してしまうといったことは起こり得ます
(清水先生のこの記事が参考になります)。逆に，p値では有意なのにベイズファクターで評価するとむしろ帰無仮説に軍配が上がるといったようなことはときどき起こります。ベイズの利点は存在しない効果をあたかも存在するかのように主張できることではなくて，統計モデリングとの相性の良さや優れた仮説評価の方法を提供してくれるところにあると僕は考えています。今回紹介したpdも，適切に使用・解釈する分には豊かな情報を与えてくれる指標の1つだと思われます
(議論はあるでしょうが)。

　ちなみに，今回使用したAさんのデータは，帰無仮説が正しい (真の効果がゼロ) という前提のもとで生成したデータです
(ちょうどいいp値とpdになるようにseedの調整はしましたが)。つまり，本当は効果がないのにもかかわらず素朴にpdを解釈すると効果があるように見えてしまうことが起こりうるという例示でもありました。p値との対応を考えると，例えばpd
> .90で効果ありとみなすという基準を設けたとしたら，単純なモデルにおいてこの基準はα =
.20と等価になってしまいます。つまり，帰無仮説が正しかったとしても20%の確率でpd >
90となり「効果あり」と誤って判断してしまうことを意味します。分析の目的によってはこれで良いこともあるかもしれませんが，科学的な研究においてはこの二値判断はかなりリスキーですね。ベイズ的なアプローチには再現性の向上が期待されている側面もありますが，よく考えずになんとなーくベイズを使ってしまうと，p値を使う以上に再現性を低めてしまう危険性があります
(査読で弾ければ良いのですが，査読者が必ずしも理解しているとは限らない)。ベイズに限らず，自分が依拠する研究法についてきちんと理解しておくことは非常に大切なことだと思います。僕も勉強し続けます。


5. RとSTANのコード

　架空データの生成・分析に使用したRコードとStanコードは以下の通りです。

Rコード: run.R

library(tidyverse)
library(BayesFactor)
library(rstan)

set.seed(21)

# Making data -------------------------------------------------
N <- 60
Y1 <- rnorm(N/2,0,1)
Y2 <- rnorm(N/2,0,1)

# Visualization of data -------------------------------------------------
dg <- data.frame(Y1 = Y1, Y2 = Y2) %>%
  gather(key = "group", value = "value") %>%
  mutate(group = ifelse(group=="Y1", "統制群", "実験群")) %>%
  mutate(group = factor(group, levels = c("統制群", "実験群")))

g <- ggplot(data = dg, aes(x = group, y = value)) +
  theme_bw() +
  geom_jitter(width = 0.2, color = "99BB99", alpha = 0.8, size = 2) +
  stat_summary(fun = "mean", fun.min = function(x)mean(x), fun.max = function(x)mean(x), geom = "crossbar", color = "black", width = 0.3)
plot(g)
ggsave("dataplot.png", g, width = 3, height = 3, dpi = 600)

# Result of t-test -------------------------------------------------
t.test(Y1,Y2) # p = .1625

# MCMC sampling by Stan -------------------------------------------------
sm_t <- stan_model(file = "ttest.stan")

standata <- list(
  N1 = length(Y1),
  N2 = length(Y2),
  Y1 = Y1,
  Y2 = Y2
)

fit <- sampling(sm_t, data = standata, seed = 8823,
                chains = 4, iter = 25500, warmup = 500, thin = 1)

s <- summary(fit)$summary
e <- rstan::extract(fit) %>% as.data.frame()

# Summary of the posterior of mu_diff -------------------------------------------------
s["mu_diff",] %>% round(3)

# Result of probability of direction (pd) -------------------------------------------------
max(mean(e$mu_diff > 0), 1-mean(e$mu_diff > 0)) # pd = .91

# Visualization of posterior of mu_diff -------------------------------------------------
dg <- data.frame(
  mu_diff = density(e$mu_diff)$x,
  density = density(e$mu_diff)$y
)

g <- ggplot(dg, aes(x = mu_diff, y = density)) +
  theme_bw() +
  geom_area(fill = "#888888") +
  geom_area(data = dg[dg$mu_diff>0,], fill = "#DD8888") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  ylab("事後確率密度")

plot(g)

ggsave("posterior.png", g, width = 3, height = 3, dpi = 600)

# Result of Bayes factor -------------------------------------------------
ttestBF(x = Y1, y = Y2, rscale = 1/sqrt(2)) # BF = 0.605, in favor of the null hypothesis



Stanコード: ttest.stan

data{
  int N1;
  int N2;
  real Y1[N1];
  real Y2[N2];
}

parameters{
  vector[2] mu;
  vector<lower = 0>[2] sigma;
}

model{
  Y1 ~ normal(mu[1], sigma[1]);
  Y2 ~ normal(mu[2], sigma[2]);
}

generated quantities{
  real mu_diff = mu[2] - mu[1];
}



タグ：R 再現性 ベイズ STAN 事後分布

posted by mutopsy at 00:00 | 統計
-  1 2 3 4 5..  次の5件>>

検索ボックス

 

プロフィール

管理人
mutopsy (武藤拓之)

自己紹介
大学で研究をしている小さな生き物です。心理学の科学的方法（数理&統計モデリング・実験法・心理測定論・仮説検定・ベイズ統計学・再現性と信用性の向上・科学哲学）とその実践（特に知覚・認知・数理心理学）に関心があります。

研究関連のWebサイト
mutopsy.net

Twitter
@mutopsy

よく読まれている記事
1. MCMCとともだちになろう [スライド紹介]
2. p値で有意と言えない効果もベイズなら効果があると言える？──事後分布に基づく仮説評価について──
3. 実験心理学者のためのベイジアンモデリング──正答率データをどうやって分析するか──
4. 金のエンゼルが当たる確率のベイズ推定
5. 認知心理学への実践：データ生成メカニズムのベイズモデリング[スライド紹介とWS感想]
6. 銀のエンゼルが当たる確率のベイズ推定：金のエンゼルの何倍当たりやすいか？
7. 回帰分析の悩みどころ (「アヒル本」7.1-7.5) [スライド紹介]
8. StanとRで折れ線回帰──空間的視点取得課題の反応時間データを説明する階層ベイズモデルを例に──[スライド紹介]
※全記事一覧はこちらにあります
最近の記事
(12/25)10分ぐらいで分かる（？），有意性検定でサンプルサイズを事前に決めなければならない理由
(12/23)Stanでdrift diffusion model (DDM) を扱うときの実践上のtips (1)──まずは直観的に理解する──
(12/21)対数ロジスティック分布をStanに実装する
(12/20)心理学の再現性と科学性と「ベイズ」に関する自分の立ち位置を整理する──JPA2020再現可能性シンポジウムのスライドを添えて──
(12/10)p値で有意と言えない効果もベイズなら効果があると言える？──事後分布に基づく仮説評価について──

カテゴリ
統計(20)
Rに関するTips(3)
雑記(3)

過去ログ
2021年12月(1)
2020年12月(4)
2020年10月(1)
2020年09月(3)
2020年03月(1)
2019年12月(2)
2019年01月(3)
2018年12月(2)
2018年07月(1)
2017年12月(8)

著書の紹介

Twitter