PearsonもSpearmanもそのP値に、Fisher z変換からのP値を算出できないか？

[yoshi416]

現在JMPにおいては、PearsonもSpearmanもそのP値は、（z変換不要の）t分布近似の値である。しかし、Fisher z変換からのP値算出を切に熱望する。特に、最近の生物系の投稿論文では、相関係数PearsonもSpearmanも、よく見かける。その場合、Fisher z変換からのP値算出が圧倒的に多い。よろしくお願いしたい。また、信頼区間構成も可能な限りFisher z変換による95%信頼区間構成を望むものである。（なお、本件に関して正確法までは不要とする）

[yusuke_ono]

ご要望を寄せていただき、誠にありがとうございます。

また、弊グループ内でのやりとりに時間がかかり、返信が遅れましてすみません。

SAS Institute Japan株式会社 JMPジャパン事業部でテスターをしています小野と申します。

PearsonおよびSpearmanの相関係数に対するp値において、FisherのZ変換近似の方法もオプションとして追加してほしいとの要望を承りました。同方法を、要望リストに追加しました。さらに私たちのほうで検討を行いたいと思います。

現在の「多変量の相関」プラットフォームでは、PearsonおよびSpearmanの相関係数に対するp値は、Studentのt分布から計算されています。

Pearsonの相関係数に対しては、このStudentのt分布に基づくp値は、 次のような条件のもとで、正確なp値です。

1. 各観測ベクトルの分布が、多変量正規分布からの独立同一分布である。
2. 欠測値がない。
3. 帰無仮説がH0: ρij = 0である（ここでρij は、第i変数と第j変数のあいだのPearsonの母相関係数）

JMPの現在の「多変量の相関」プラットフォームでは、Pearsonの相関係数に対する信頼区間は、FisherのZ変換近似（バイアス項を含まないもの）から計算されています。そのため、信頼区間が、ある定数c（-1.0<c<1.0）を含まないのであれば、帰無仮説「H0: ρij = c」が棄却されることを意味します（デフォルトの信頼係数は95%、つまり、有意水準5%となっています）。

標本サイズがある程度、大きければ、帰無仮説「H0: ρij = 0」に対する検定として、Studentのt分布に基づく検定も、FisherのZ変換に基づく検定も似た結果となります。ただし、特に小標本の場合には、両者の結果は違ったものです。このため、現在の「多変量の相関」プラットフォームでは、検定と信頼区間が食い違った結果になる可能性もあります（統計的検定では「H0: ρij = 0」が棄却されるのに、信頼区間が0を含むこともありうる。もしくは、その逆の関係になることもありうる。）これは、ユーザーの方々には、やや混乱をもたらす動作かもしれません。しかし、上記したStudentのt分布に基づく方法は、上記した条件のもとで正確な方法であり、かつ、計算も簡単です。一方、「正確」な信頼区間は計算がやや難しいです。そのような理由から、p値としてはStudentのt分布を用いて、信頼区間としてはFisherのZ変換近似を用いるのは合理的だと考えています。

現在の「多変量の相関」プラットフォームでは、Spearmanの相関係数に対する信頼区間は計算されません。

「正確」な検定や信頼区間については、私たちの知る限り、2つの方法があります。１つは、Pearsonの相関係数に対して、多変量正規分布を前提として計算される信頼区間です。この方法はρij≠0の場合には、数値積分が必要です。計算が複雑なこともあり、こちらの考えに沿った、ρij≠0に対する正確な検定、および、ρijに対する正確な信頼区間はそれほど普及していないと思われます。もう１つの正確検定は、並び替え検定に基づく方法です。並び替え検定は、条件付き「正確」な方法です。ρij≠0のときに、どのように並び替えを行えばいいのかが難しいです。いずれの正確な方法も、JMPでは実装されていません。唯一、「H0: ρij=0」のときの、Studentのt分布に基づく（多変量正規分布であるときの）正確な方法が、Pearsonの相関係数に対して提供されているだけです。

まとめると、現在のJMPでは次のような作りになっています。

• PearsonおよびSpearmanの相関係数に対するp値は、Studentのt分布に基づいて計算されます。
• Pearsonの相関係数に対する信頼区間は、FisherのZ変換（バイアス項なし）による近似で計算されます。
• Spearmanの相関係数に対する信頼区間は計算されません。

現在のJMP上にて、FisherのZ変換近似からp値を計算したい場合には、JMPスクリプト言語（JSL; JMP Scripting Language）にてプログラミングをする必要があります。コードは以下のようなものになります（以下のコードは、バイアス項なしのFisherのZ変換によるものです)。もし、以下のようなコードによるFisherのZ変換の計算に興味がある場合はご返信ください。コードについて解説いたします。また、以下は計算方法を分かりやすく提示するためのものですが、さらに汎用性があるコードにすることもできます。

r = 0.4;
n = 15;

z = 0.5 * Log((1+r)/(1-r));
up_z = z + Normal Quantile(0.975)*Sqrt(1/(n-3));
low_z = z - Normal Quantile(0.975)*Sqrt(1/(n-3));
Show((Exp(2*low_z)-1)/(Exp(2*low_z)+1));
Show((Exp(2*up_z)-1)/(Exp(2*up_z)+1));
p_z = 2*Normal Distribution(-Abs(z*Sqrt(n-3)));
Show(p_z);

t = r * Sqrt((n-2)/(1-r^2));
p_t = 2 * T Distribution(-Abs(t), n - 2);
Show(p_t);

長文、失礼いたしました。

ご要望を寄せて頂いたことを、重ねてお礼申し上げます。

上記の返信にて不十分・不明な点などございましたら、お手数ですがご返信ください。

[yoshi416]

小野さま

ご連絡ありがとうございました。

急ぎではないので、むしろ、丁寧に回答ありがとうございます。

つきましては、以下となります。

(1)Spearmanの相関係数に対する、Fisher Z変換によるバイアス補正ありの95%信頼区間構成は諦めます。

(2)PearsonおよびSpearmanの相関係数に対するp値に関しては、「Fisher Z変換」、「Studentのt分布（現搭載済み）」、「正確法」のうち、やはり、「Fisher Z変換」によるP値とうのを、論文に記載したいため、

PearsonおよびSpearmanの相関係数に対するp値に関しては、

1)FisherのZ変換の計算に興味がありコードについて解説希望

2)FisherのZ変換の計算に興味がありさらに汎用性があるコード解説希望

となります。

上の(2)で「正確法」は諦めます。

背景は、20-30例の臨床画像解析で用いたいため、小数例20-30例であるためです。

私は、現在、JMP15.2(win10の64bit版)を利用しています

引き続き、よろしくお願い致します。

[yusuke_ono]

以下のリンクのページに、解説のPDF文書と、プログラム例をアップロードしました。

プログラムの例（Fisherのz変換） 

PDF文書では、JMPスクリプト言語の概要をまず説明し、次にFisherのz変換近似によるp値や信頼区間を求める方法を説明しました。最後に、より汎用的に、それらを求めるプログラム例を提示しました。

拡張子がjslとなっているファイルは、プログラム例です。このファイルをまずダウンロードしておき、分析対象のデータテーブルをJMPで開いたあと、［ファイル］→［開く］で同ファイルを開いてください。右クリックで呼び出されるメニューにて［スクリプトの実行］を選択すると、同プログラムが実行されます。呼び出された起動ダイアログにて、分析対象の列を選択してください。

ご不明な点や不備などを見つましたら、返信して頂けるとありがたいです。

[yoshi416]

テクニカルグループ　小野さま

この度は、丁寧なPDFファイル御礼申し上げます。

非常に判り易く、見やすいです。

さて、PDF P.4 z変換p値の計算式で、分散を1/(n-3)としたのは正規分布により従うようにする為、と理解しましたが、

その認識で合っておりますでしょうか？

さて、その理解を示すため、以下のグラフ作成は、JMPのグラフビルダーで作成可能でしょうか？

つきましては、更なるオプションですが、以下を希望します。

（ここまでは大変でしたら、無理にとは申しません）

FisherのZ変換：Z=1/2log(1+r/1-r) ~N(0, 1/n-3)

を、以下のロジット変換

ロジっト変換：Z=log (p / 1-p)

と対比させて、以下のX軸Y軸からなる、グラフ作成機能があると理解が助かります。

FisherのZ変換：X軸は r(-∞<r<+∞)、Y軸はZ

ロジット変換：X軸は p(p>0のみ)、Y軸はZ

FisherのZ変換とロジット変換のグラフ機能は、さらなるオプション、という事で構いません。

よろしくお願い致します。

[yusuke_ono]

ご返信ありがとうございます。

「分散を1/(n-3)としたのは正規分布により従うようにするため」ではないと思います。「標本相関係数をz変換したものの分散は1/(n-1)よりも、1/(n-3)のほうが近い。だから、1/(n-1)ではなく1/(n-3)を分散として用いたほうがよい」ということだと思います。

つまり、正規分布に対する近似を良くするためではなく、正規分布の分散に対する評価を良くするために、（1/(n-1)ではなくて）1/(n-3)を用いているのだと思います。

z変換近似の分散を1/(n - 3)とするのは、Fisher(1925, p.162)や増山（1943, p.71）といった古い教科書（ハウツー本）にも登場しますので、ある程度は普及している近似法ではないかと思います。なお、H0: ρ=0の検定について、Fisher(1925)と増山(1943）は両方ともt分布による正確法を最初に紹介しています。

Fisher自身がどこかに1/(n-3)にした理由を書いているとは思いますが、最近の書籍では竹村（1991,pp.259-261）に理由が書かれています。私自身は数式や「Hotelling（1953）」といった参考文献を追えていませんが、竹村（1991,pp.259-261）にて、Var[Z]の分散をn^{-2}の項まで求め、それを少し変形すると1/(n-3)にしたほうが、n^{-2}の項の絶対値が小さくなっているのが分かる、と書かれています。

後半部分の話ですが、お手軽な方法として、乱数シミュレーションによってz変換近似で分散を1/(n-3)としたほうが（1/(n-1)としたときよりも）近似が良くなることを確認してみたいということでしょうか？ JMPスクリプト言語では、反復処理や乱数生成の関数も用意されていますので、次のぐらいのコード量でプログラミングできます。（各関数については、［ヘルプ］→［スクリプトの索引］でご確認ください。）

Names Default To Here(1);
Random Reset(111111111);
n = 20;
nsim = 1000000;
rho = 0.4;
c = Sqrt(1 - rho^2);
result = J(nsim, 1, .);
For(i = 1, i <= nsim, i++,
  x1 = J(n, 1, Random Normal());
  x2 = rho * x1 + c * J(n, 1, Random Normal());
  r = Correlation(x1 || x2)[1, 2];
  result[i] = 0.5 * Log((1+r)/(1-r));
);

dt = As Table(result);
dt << Distribution( Y([1]));
Show(Stddev(result), Sqrt(1/(n -3)), Sqrt(1/(n - 1)));

上記のプログラムの実行結果は、次の通りです（Random Reset関数に指定している乱数シード値や、変数nsimで指定しているシミュレーション回数により、結果は変わってきます）。この例では、sqrt(1/(n-3))のほうがsqrt(1/(n-1))よりも真の標準誤差に近くなっています。

Std Dev(result) = 0.241861636994535;
Sqrt(1 / (n - 3)) = 0.242535625036333;
Sqrt(1 / (n - 1)) = 0.229415733870562;

正規密度関数を重ね書きしたものは次図のようになります。

シミュレーション結果.png

参考文献
Fisher, R.A. (1925) Statistical Methods for Research Workers, Oliver and Boyd

増山元三郎（1943）『少數例の纏め方と實驗計畫の立て方』河出書房

竹村彰通（1991）『多変量推測統計の基礎』共立出版