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 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過程についてもっと詳しく説明します。

対数ロジスティック分布を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!

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];
}

一般化ガンマ分布をStanに実装する

　最近は生存時間解析 (survival analysis; 継続時間解析とも) を使った研究に勤しんでいるのですが，パラメトリックモデルについて勉強していたら，ガンマ分布・ワイブル分布・指数分布を包含する一般化ガンマ分布 (generalized gamma distribution) なる面白い分布が登場したので，これをStanに実装してみました。一般化ガンマ分布は非負の値を返す確率分布で，確率密度関数は以下の式で表されます。

この式ではrateパラメータをλ (逆数の1/λがscaleパラメータ)，2つの形状パラメータをg，wと置いています。gはガンマ分布の形状パラメータ，wはワイブル分布の形状パラメータに対応しています。g=1のときワイブル分布，w=1のときガンマ分布，g=w=1のとき指数分布と一致することが特徴です。他のパラメータ化の方法もありますが，他の分布との関係が分かりやすいのでこのパラメータ化を採用しています。分布間の関係は次のように図示できます。

一般化ガンマ分布を使えば，データがガンマ分布・ワイブル分布・指数分布のどれに従うのかが分からない場合でも分布形に当たりを付けることができたりします。もちろん，情報量規準などを用いて4つの分布のどれが適切かを判断することもできます。便利ですね。

　一般化ガンマ分布を使ったモデリングができるようにStanコードを書いてみました。

ggamma.stan：一般化ガンマ分布モデル

functions{
  //for文で1つずつ計算するときはこっち
  real ggamma_lpdf(real x, real rate, real g, real w){
    real lpdf;
    lpdf = log(w) + g*w*log(rate) + (g*w-1)*log(x) - (x*rate)^w - lgamma(g);
    return(lpdf);
  }
  
  //データをベクトルで渡すときにはこっち
  real sum_ggamma_lpdf(vector x, real rate, real g, real w){
    real sum_lpdf;
    sum_lpdf = sum(log(w) + g*w*log(rate) + (g*w-1)*log(x) - exp(w*log(x*rate)) - lgamma(g));
    return(sum_lpdf);
  }
}

data{
  int N;
  vector[N] Y;
}

parameters{
  real rate;
  real g;
  real w;
}

model{
  //事前分布（シミュレーションしやすいように狭めにしています）
  rate ~ uniform(0,10);
  g ~ uniform(0,10);
  w ~ uniform(0,10);
  
  //モデル
  target += sum_ggamma_lpdf(Y | rate, g, w);
  
  ////以下のような書き方でも同じ
  //for(i in 1:N) Y[i] ~ ggamma(rate, g, w);
}

　上記コードではggamma_lpdf()とsum_ggamma_lpdf()という2種類の関数を定義していますが，サンプルを要素ごとにreal型としてfor文で渡すときにはggamma_lpdf()を，サンプルをvectorとしてまとめて渡すときにはsum_ggamma_lpdf()を使います。後者では対数尤度をまとめて計算してその和を求めているという訳です。推定結果はどちらを使っても同じですが，データの数が多い時にはおそらくsum_ggamma_lpdf()のほうが計算が速いです（ちゃんと検証した訳ではありませんが）。では，このコードでうまくパラメータが推定できるか確認してみましょう。今回は，一般化ガンマ分布（λ = 0.5，g = 3，w = 5），ガンマ分布（λ = 0.5，g = 3），ワイブル分布（λ = 0.5，w = 5），指数分布（λ = 0.5）からの乱数を2,000個ずつ発生させて4つのデータを作り，それらに一般化ガンマ分布を当てはめてパラメータを推定してみます。Rコードと推定結果は以下の通りです。なお，一般化ガンマ分布をRで扱うにはggammaというパッケージが利用できます（gamlssというパッケージもありますが，パラメータ化の方法が異なります）。

ggamma.R：一般化ガンマ分布モデルの確認

#R 3.6.3を使用
library(tidyverse) #version 1.2.1
library(rstan) #version 2.19.2
library(ggamma) #version 1.0.1

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

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

N <- 2000
rate <- 0.5
shape_g <- 3
shape_w <- 5

set.seed(8823)
Y_ggamma <- rggamma(N, a = 1/rate, b = shape_w, k = shape_g)
Y_gamma <- rgamma(N, rate = rate, shape = shape_g)
Y_weibull <- rweibull(N, scale = 1/rate, shape = shape_w)
Y_exp <- rexp(N, rate = rate)

## Stanに渡す ----
sm <- stan_model("stanmodel/ggamma.stan")

warmup <- 1000
iter <- 1000 + warmup
t0<-proc.time()
fit_ggamma <- sampling(sm,data=list(N=N,Y=Y_ggamma)
                ,seed=8823,chains=4,iter=iter,warmup=warmup
)

fit_gamma <- sampling(sm,data=list(N=N,Y=Y_gamma)
                       ,seed=8823,chains=4,iter=iter,warmup=warmup
)

fit_weibull <- sampling(sm,data=list(N=N,Y=Y_weibull)
                       ,seed=8823,chains=4,iter=iter,warmup=warmup
)

fit_exp <- sampling(sm,data=list(N=N,Y=Y_exp)
                       ,seed=8823,chains=4,iter=iter,warmup=warmup
)

## 一般化ガンマ分布 (rate = 0.5, g = 3, w = 5) からの乱数の結果 ----
s_ggamma <- summary(fit_ggamma)$summary
round(s_ggamma,3)

#          mean se_mean    sd     2.5%      25%      50%      75%    97.5%   n_eff  Rhat
# rate    0.511   0.002 0.037    0.458    0.484    0.504    0.532    0.598 488.966 1.008
# g       3.271   0.030 0.667    2.276    2.777    3.161    3.663    4.825 492.162 1.008
# w       4.992   0.023 0.511    3.982    4.628    4.996    5.364    5.974 493.785 1.009
# lp__ -332.568   0.040 1.214 -335.789 -333.102 -332.264 -331.696 -331.168 920.305 1.001

## ガンマ分布 (rate = 0.5, g = 3, w = 1) からの乱数の結果 ----
s_gamma <- summary(fit_gamma)$summary
round(s_gamma,3)

#           mean se_mean    sd      2.5%       25%       50%       75%     97.5%   n_eff  Rhat
# rate     0.873   0.031 0.595     0.372     0.569     0.731     0.979     2.186 370.082 1.017
# g        3.695   0.039 0.830     2.426     3.142     3.575     4.111     5.644 464.461 1.012
# w        0.878   0.004 0.099     0.686     0.813     0.877     0.941     1.081 496.439 1.010
# lp__ -5125.935   0.038 1.201 -5129.165 -5126.483 -5125.623 -5125.061 -5124.577 982.829 1.004

## ワイブル分布 (rate = 0.5, g = 1, w = 5) からの乱数の結果 ----
s_weibull <- summary(fit_weibull)$summary
round(s_weibull,3)

#         mean se_mean    sd      2.5%       25%       50%       75%     97.5%    n_eff  Rhat
#rate     0.502   0.001 0.017     0.475     0.491     0.501     0.512     0.540  515.125 1.007
#g        1.036   0.006 0.131     0.814     0.942     1.027     1.110     1.327  515.630 1.008
#w        4.999   0.017 0.388     4.234     4.738     4.987     5.254     5.772  529.559 1.007
#lp__ -1076.843   0.040 1.281 -1080.196 -1077.396 -1076.519 -1075.901 -1075.407 1020.983 1.003

## 指数分布 (rate = 0.5, g = 1, w = 1) からの乱数の結果 ----
s_exp <- summary(fit_exp)$summary
round(s_exp,3)

#           mean se_mean    sd      2.5%       25%       50%       75%     97.5%   n_eff  Rhat
# rate     0.574   0.004 0.103     0.425     0.506     0.560     0.624     0.814 553.282 1.005
# g        1.066   0.006 0.136     0.838     0.975     1.059     1.144     1.362 581.750 1.005
# w        0.982   0.003 0.076     0.837     0.932     0.980     1.030     1.132 598.027 1.004
# lp__ -3293.357   0.042 1.288 -3296.725 -3293.945 -3293.026 -3292.422 -3291.911 932.671 1.002

　サンプルサイズや乱数の関係でEAP (mean) がやや真値から離れているように見える箇所もありますが，いずれの場合でも95%確信区間に真値が含まれていることが確認できます。特に，ガンマ分布ではw≃1，ワイブル分布ではg≃1，指数分布ではg≃w≃1と正しく推定できていることが分かります。というわけで，ぜひみなさんも一般化ガンマ分布で遊んでみてください。※この記事で紹介したコードを研究等に使用される場合は自己責任でお願いします。