JP6758632B2 - 正確検定計算装置、正確検定計算方法、およびプログラム - Google Patents

Description

この発明は、フィッシャーの正確確率検定（Fisher's exact test）を効率良く計算する技術に関する。

2×2の分割表として与えられた説明変数と目的変数との関連の有無に関して仮説検定を行う科学的検定法の一つとして、フィッシャーの正確確率検定が広く知られている。非特許文献１には、フィッシャーの正確確率検定の利用例として、ゲノムワイド関連解析（GWAS: Genome-Wide Association Study）が記載されている。

フィッシャーの正確確率検定について説明する。表１はN人の被験者を変異の有無と特定の疾患の発症の有無で分類および数え上げした2×2分割表の例である。a, b, c, dは度数を、X, W, Y, Zは小計を表し、X=a+b, W=c+d, Y=a+c, Z=b+dである。a, b, c, d, X, W, Y, Zはすべて非負整数である。

フィッシャーの正確確率検定は、非負整数iについて、式（１）に示す確率p_iを計算し、

式（２）に示すp値と、有意水準と呼ばれる所定の値αとの大小関係により、説明変数（表１の例では変異の有無）と目的変数（表１の例では疾患の有無）との関連の有無を調査する。

ここで、p_aは実際の集計値であるa, b, c, dを度数とする分割表について式（１）により算出した確率値である。

非特許文献２には、度数の総和N=a+b+c+dと有意水準αとを固定して、すべてのN=a+b+c+dを満たす非負整数a, b, c, dの組それぞれについてp値がα未満であるか否かを記録した表を作成し、その表を参照することで効率的にフィッシャーの正確確率検定を行う方法が提案されている。

Konrad Karczewski, "How to do a GWAS", Lecture note in GENE 210: Genomics and Personalized Medicine, 2015. 千田浩司，長谷川聡，濱田浩気，荻島創一，三澤計治，長崎正朗，"秘密計算フィッシャー正確検定（1） 〜標本数が少ない場合" ，電子情報通信学会技術研究報告，Vol. 116，No. 129，pp. 245-251，2016年7月

しかしながら、非特許文献２に記載の従来技術は、N=a+b+c+dを満たす非負整数a, b, c, dの組すべてに対して検定の結果を記録しておくために、Ω(N³)の大きさの表を保持しておく必要があった。

この発明は、上述のような点に鑑みて、フィッシャーの正確確率検定の結果を効率的に求めることを目的とする。

上記の課題を解決するために、この発明の正確検定計算装置は、αをフィッシャーの正確確率検定の有意水準とし、Sをフィッシャーの正確確率検定の入力となる度数の組の集合とし、p(s)を度数の組s∈Sに対するフィッシャーの正確確率検定のp値とし、Tを各要素が対である集合とし、fを集合Sから集合Tへの全単写とし、f^-1を全単写fの逆写像とし、Y(x):={y | (x, y)∈T, p(f^-1(x, y))≧α}とし、iを集合Tの要素の番号とし、集合Tの各要素(x_i, y_i)のうち異なるx_iそれぞれについて、Y(x_i)の最小値y₀(x_i)とY(x_i)の最大値y₁(x_i)とをx_iに関連付けた参照表を記憶する参照表記憶部と、度数の組s∈Sに対して全単写f(s)を計算し、集合Tの要素(x_s, y_s)を得る全単写部と、参照表からx_sに関連付けられたY(x_s)の最小値y₀(x_s)とY(x_s)の最大値y₁(x_s)とを得る写像部と、y_sが最小値y₀(x_s)以上かつ最大値y₁(x_s)以下であるか否かを判定する判定部と、を含む。

この発明によれば、Nを固定しない場合にはO(M³)の大きさの表（MはNの上界）を、Nを固定する場合にはO(N²)の大きさの表を、それぞれ用いてフィッシャーの正確確率検定を計算できる。従来技術では、Nを固定する場合にはΩ(N³)の大きさの表を用いる必要があった。したがって、この発明によれば、フィッシャーの正確確率検定の結果を効率的に求めることができる。

図１は、正確検定計算装置の機能構成を例示する図である。 図２は、正確検定計算方法の処理手続きを例示する図である。

この発明では、自由度xの2×2の分割表に対するフィッシャーの正確確率検定において、x-1個の値を固定する場合に、p値が有意水準α以上となる範囲がある区間で表現できることを利用する。x-1個の値の組それぞれについて、p値が有意水準以上となるような残りの１個の値の範囲の上限と下限とを記録した参照表を用いることで、従来技術よりも小さい表を使ってフィッシャーの正確確率検定を効率よく計算することができる。

以下、この発明の実施の形態について詳細に説明する。なお、図面中において同じ機能を有する構成部には同じ番号を付し、重複説明を省略する。

［第一実施形態］
第一実施形態の正確検定計算装置１は、図１に例示するように、参照表記憶部１０、参照表生成部１１、入力部１２、全単写部１３、写像部１４、判定部１５、および出力部１６を含む。この正確検定計算装置１が、図２に例示する各ステップの処理を行うことにより第一実施形態の正確検定計算方法が実現される。

正確検定計算装置１は、例えば、中央演算処理装置（CPU: Central Processing Unit）、主記憶装置（RAM: Random Access Memory）などを有する公知又は専用のコンピュータに特別なプログラムが読み込まれて構成された特別な装置である。正確検定計算装置１は、例えば、中央演算処理装置の制御のもとで各処理を実行する。正確検定計算装置１に入力されたデータや各処理で得られたデータは、例えば、主記憶装置に格納され、主記憶装置に格納されたデータは必要に応じて中央演算処理装置へ読み出されて他の処理に利用される。正確検定計算装置１の各処理部は、少なくとも一部が集積回路等のハードウェアによって構成されていてもよい。正確検定計算装置１が備える各記憶部は、例えば、RAM（Random Access Memory）などの主記憶装置、ハードディスクや光ディスクもしくはフラッシュメモリ（Flash Memory）のような半導体メモリ素子により構成される補助記憶装置、またはリレーショナルデータベースやキーバリューストアなどのミドルウェアにより構成することができる。

以下、図２を参照して、第一実施形態の正確検定計算装置１が実行する正確検定計算方法について説明する。

以下の説明では、αをフィッシャーの正確確率検定の有意水準とする。Sをフィッシャーの正確確率検定の入力となる度数の組の集合とする。Zは整数全体の集合を表す。正確検定計算装置１は、Sの要素である度数の組sを入力とし、度数の組sのp値が有意水準α未満であるか否かを示す検定結果を出力する。

第一実施形態では、記号を以下のように定義する。Mを0以上の整数とする。入力の集合Sを、S={(a, b, c, d) | a∈Z, b∈Z, c∈Z, d∈Z, a≧0, b≧0, c≧0, d≧0, a+b+c+d≦M}とする。p(s)を、入力s∈Sに対するフィッシャーの正確確率検定のp値とする。fを、f:(a, b, c, d)→((N, X, Y), a)である全単写とする。なお、N=a+b+c+d, X=a+b, Y=a+cである。全単写fの値域をTとする。すなわち、集合Tは、各要素が、ある度数と、その度数を含む度数の小計との対である集合である。f^-1を全単写fの逆写像とする。

ステップＳ１において、参照表生成部１１は、Y(x):={y | (x, y)∈T, p(f^-1(x, y))≧α}、Y(x)の最小の要素をy₀(x)、Y(x)の最大の要素をy₁(x)、写像gをg: x→(y₀(x), y₁(x))とし、写像gを表す参照表Gを生成する。すなわち、全単写fの値域Tの各要素((N_i, X_i, Y_i), a_i)から、異なる(N_i, X_i, Y_i)それぞれについて、p(f^-1((N_i, X_i, Y_i), a_i))≧αを満たすa_iの最大値y₀((N_i, X_i, Y_i))と最小値y₁((N_i, X_i, Y_i))とを求め、(N_i, X_i, Y_i)をキーとして最大値y₀((N_i, X_i, Y_i))と最小値y₁((N_i, X_i, Y_i))とを値とする参照表Gを生成する。参照表生成部１１は、生成した参照表Gを参照表記憶部１０へ記憶する。

ステップＳ２において、入力部１２は、度数の組s=(a, b, c, d)∈Sを入力として受け取る。入力部１２は、入力された度数の組sを全単写部１３へ入力する。

ステップＳ３において、全単写部１３は、入力部１２から度数の組sを受け取り、sに対する全単写fを計算する。すなわち、((N, X, Y), a)←f(s)を計算する。全単写部１３は、計算結果((N, X, Y), a)を写像部１４へ入力する。

ステップＳ４において、写像部１４は、全単写部１３から((N, X, Y), a)を受け取り、(N, X, Y)に対する写像gを計算する。すなわち、(y₀((N, X, Y)), y₁((N, X, Y)))←g((N, X, Y))を計算する。写像gの計算は、参照表記憶部１０に記憶された参照表Gを参照することで行う。具体的には、参照表記憶部１０に記憶された参照表Gから(N, X, Y)をキーとして最小値y₀((N, X, Y))と最大値y₁((N, X, Y))とを取得する。写像部１４は、取得した最小値y₀((N, X, Y))と最大値y₁((N, X, Y))とを判定部１５へ入力する。

ステップＳ５において、判定部１５は、写像部１４から最小値y₀((N, X, Y))と最大値y₁((N, X, Y))とを受け取り、全単写部１３の計算結果として得たaが最小値y₀((N, X, Y))から最大値y₁((N, X, Y))までの範囲に入っているか否かを判定する。すなわち、y₀((N, X, Y))≦a≦y₁((N, X, Y))が成り立つか否かを判定する。判定部１５は、判定結果を出力部１６へ入力する。

ステップＳ６において、出力部１６は、判定部１５から受け取った判定結果をフィッシャーの正確確率検定の検定結果として正確検定計算装置１から出力する。

［第二実施形態］
第二実施形態では、記号を以下のように定義する。Nを0以上の整数とする。入力の集合Sを、S={(a, b, c, d) | a∈Z, b∈Z, c∈Z, d∈Z, a≧0, b≧0, c≧0, d≧0, a+b+c+d=N}とする。p(s)を、入力s∈Sに対するフィッシャーの正確確率検定のp値とする。fを、f(a, b, c, d)→((X, Y), a)である全単写とする。なお、X=a+b, Y=a+cである。全単写fの値域をTとする。すなわち、集合Tは、各要素が、ある度数と、その度数を含む度数の小計との対である集合である。f^-1を全単写fの逆写像とする。

ステップＳ１において、参照表生成部１１は、Y(x):={y | (x, y)∈T, p(f^-1(x, y))≧α}、Y(x)の最小の要素をy₀(x)、Y(x)の最大の要素をy₁(x)、写像gをg: x→(y₀(x), y₁(x))とし、写像gを表す参照表Gを生成する。すなわち、全単写fの値域Tの各要素((X_i, Y_i), a_i)から、異なる(X_i, Y_i)それぞれについて、p(f^-1((X_i, Y_i), a_i))≧αを満たすa_iの最大値y₀((X_i, Y_i))と最小値y₁((X_i, Y_i))とを求め、(X_i, Y_i)をキーとして最大値y₀((X_i, Y_i))と最小値y₁((X_i, Y_i))とを値とする参照表Gを生成する。参照表生成部１１は、生成した参照表Gを参照表記憶部１０へ記憶する。

ステップＳ２において、入力部１２は、度数の組s=(a, b, c, d)∈Sを入力として受け取る。入力部１２は、入力された度数の組sを全単写部１３へ入力する。

ステップＳ３において、全単写部１３は、入力部１２から度数の組sを受け取り、sに対する全単写fを計算する。すなわち、((X, Y), a)←f(s)を計算する。全単写部１３は、計算結果((X, Y), a)を写像部１４へ入力する。

ステップＳ４において、写像部１４は、全単写部１３から((X, Y), a)を受け取り、(X, Y)に対する写像gを計算する。すなわち、(y₀((X, Y)), y₁((X, Y)))←g((X, Y))を計算する。写像gの計算は、参照表記憶部１０に記憶された参照表Gを参照することで行う。具体的には、参照表記憶部１０に記憶された参照表Gから(X, Y)をキーとして最小値y₀((X, Y))と最大値y₁((X, Y))とを取得する。写像部１４は、取得した最小値y₀((X, Y))と最大値y₁((X, Y))とを判定部１５へ入力する。

ステップＳ５において、判定部１５は、写像部１４から最小値y₀((X, Y))と最大値y₁((X, Y))とを受け取り、全単写部１３の計算結果として得たaが最小値y₀((X, Y))から最大値y₁((X, Y))までの範囲に入っているか否かを判定する。すなわち、y₀((X, Y))≦a≦y₁((X, Y))が成り立つか否かを判定する。判定部１５は、判定結果を出力部１６へ入力する。

ステップＳ６において、出力部１６は、判定部１５から受け取った判定結果をフィッシャーの正確確率検定の検定結果として正確検定計算装置１から出力する。

上述のように構成することにより、この発明の正確検定計算技術によれば、Nを固定しない場合にはO(M³)の表（MはNの上界）を、Nを固定する場合にはO(N²)の表を、それぞれ用いてフィッシャーの正確確率検定を計算できる。従来の方法では、Nを固定する場合にΩ(N³)の表を用いる必要があった。

以上、この発明の実施の形態について説明したが、具体的な構成は、これらの実施の形態に限られるものではなく、この発明の趣旨を逸脱しない範囲で適宜設計の変更等があっても、この発明に含まれることはいうまでもない。実施の形態において説明した各種の処理は、記載の順に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。

［プログラム、記録媒体］
上記実施形態で説明した各装置における各種の処理機能をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記各装置における各種の処理機能がコンピュータ上で実現される。

この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。

また、このプログラムの流通は、例えば、そのプログラムを記録したDVD、CD-ROM等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。

このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。そして、処理の実行時、このコンピュータは、自己の記憶装置に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるASP（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。

また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、本装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。

上述のゲノムワイド関連解析への応用結果から明らかなように、この発明は、フィッシャーの正確確率検定を利用した解析、例えば、ゲノムワイド関連解析、ゲノム分析、臨床研究、社会調査、学術研究、実験結果の分析、マーケティング調査、統計計算、医療情報分析、顧客情報分析、売り上げ分析において、フィッシャーの正確確率検定を効率良く行うことに適用することが可能である。

１ 正確検定計算装置
１０ 参照表記憶部
１１ 参照表生成部
１２ 入力部
１３ 全単写部
１４ 写像部
１５ 判定部
１６ 出力部

Claims (5)

1. αをフィッシャーの正確確率検定の有意水準とし、Sをフィッシャーの正確確率検定の入力となる度数の組の集合とし、p(s)を度数の組s∈Sに対するフィッシャーの正確確率検定のp値とし、Tを各要素が対である集合とし、fを集合Sから集合Tへの全単写とし、f^-1を上記全単写fの逆写像とし、Y(x):={y | (x, y)∈T, p(f^-1(x, y))≧α}とし、
   iを集合Tの要素の番号とし、集合Tの各要素(x_i, y_i)のうち異なるx_iそれぞれについて、Y(x_i)の最小値y₀(x_i)とY(x_i)の最大値y₁(x_i)とを上記x_iに関連付けた参照表を記憶する参照表記憶部と、
   度数の組s∈Sに対して上記全単写f(s)を計算し、集合Tの要素(x_s, y_s)を得る全単写部と、
   上記参照表から上記x_sに関連付けられたY(x_s)の最小値y₀(x_s)とY(x_s)の最大値y₁(x_s)とを得る写像部と、
   上記y_sが上記最小値y₀(x_s)以上かつ上記最大値y₁(x_s)以下であるか否かを判定する判定部と、
   を含む正確検定計算装置。
2. 請求項１に記載の正確検定計算装置であって、
   上記集合Sは、Mを0以上の整数とし、Zを整数全体の集合とし、a, b, c, dをフィッシャーの正確確率検定の入力となる度数とし、S={(a, b, c, d) | a∈Z, b∈Z, c∈Z, d∈Z, a≧0, b≧0, c≧0, d≧0, a+b+c+d≦M}を満たす度数の組(a, b, c, d)の集合であり、
   上記集合Tの要素xは、N=a+b+c+dとし、X=a+bとし、Y=a+cとし、N, X, Yの組(N, X, Y)であり、
   上記集合Tの要素yは、上記度数aである、
   正確検定計算装置。
3. 請求項１に記載の正確検定計算装置であって、
   上記集合Sは、Nを0以上の整数とし、Zを整数全体の集合とし、a, b, c, dをフィッシャーの正確確率検定の入力となる度数とし、S={(a, b, c, d) | a∈Z, b∈Z, c∈Z, d∈Z, a≧0, b≧0, c≧0, d≧0, a+b+c+d=N}を満たす度数の組(a, b, c, d)の集合であり、
   上記集合Tの要素xは、X=a+bとし、Y=a+cとし、X, Yの組(X, Y)であり、
   上記集合Tの要素yは、上記度数aである、
   正確検定計算装置。
4. 参照表記憶部と全単写部と写像部と判定部とを備える正確検定計算装置が実行する正確検定計算方法であって、
   αをフィッシャーの正確確率検定の有意水準とし、Sをフィッシャーの正確確率検定の入力となる度数の組の集合とし、p(s)を度数の組s∈Sに対するフィッシャーの正確確率検定のp値とし、Tを各要素が対である集合とし、fを集合Sから集合Tへの全単写とし、f^-1を上記全単写fの逆写像とし、Y(x):={y | (x, y)∈T, p(f^-1(x, y))≧α}とし、
   上記参照表記憶部に、iを集合Tの要素の番号とし、集合Tの各要素(x_i, y_i)のうち異なるx_iそれぞれについて、Y(x_i)の最小値y₀(x_i)とY(x_i)の最大値y₁(x_i)とを上記x_iに関連付けた参照表が記憶されており、
   上記全単写部が、度数の組s∈Sに対して上記全単写f(s)を計算し、集合Tの要素(x_s, y_s)を得、
   上記写像部が、上記参照表から上記x_sに関連付けられたY(x_s)の最小値y₀(x_s)とY(x_s)の最大値y₁(x_s)とを得、
   上記判定部が、上記y_sが上記最小値y₀(x_s)以上かつ上記最大値y₁(x_s)以下であるか否かを判定する、
   正確検定計算方法。
5. 請求項１から３のいずれかに記載の正確検定計算装置としてコンピュータを機能させるためのプログラム。

Images