ロジスティック回帰分析
ロジスティック回帰分析について、まず分析例の説明をする。その後、後半において、理論とプログラムの使用法について解説する。
シンプソンのパラドックスと調整変数の関係も扱っている。なお、本書「心理学データ分析と測定」の「第4章 比率(確率)の条件差を見る」は、本ウェブサイトの内容を基にしている。「心理学データ分析と測定」では、回帰分析の基本的なことも説明している。統計学の基本的なことは、拙著「データ分析のための統計学入門:統計学の考え方」で説明している。
ロジスティック回帰分析は、比率あるいは確率を従属変数とする回帰分析である。例えば、入学試験の合格確率と模擬試験の点数および体重の組に対して、次のロジスティック回帰モデル
を設定する。
表1.(合否、模試成績、体重)の人工データ。合否は合格を1、不合格を0で表している。 |
(1, 790, 67) (0, 626, 57) (1, 869, 53) (1, 688, 54) (0, 687, 53) (0, 614, 48) (1, 763, 55) (0, 744, 58) (0, 629, 59) (0, 660, 62) (0, 680, 62) (1, 651, 59) (0, 771, 59) (0, 787, 72) (0, 641, 60) (1, 706, 62) (1, 821, 60) (1, 840, 59) (0, 644, 60) (0, 595, 48) (1, 753, 49) (1, 697, 72) (1, 743, 56) (0, 593, 56) (1, 828, 72) (1, 726, 55) (0, 713, 67) (1, 728, 58) (1, 748, 51) (0, 521, 56) (1, 770, 71) (1, 751, 77) (1, 752, 59) (0, 615, 48) (1, 956, 65) (1, 698, 62) (0, 668, 50) (1, 781, 54) (1, 583, 51) (1, 768, 72) (0, 705, 66) (1, 769, 62) (1, 744, 67) (1, 826, 62) (1, 864, 57) (0, 646, 70) (1, 896, 75) (1, 879, 56) (0, 664, 46) (0, 594, 73) (1, 696, 68) (0, 630, 65) (0, 729, 60) (1, 930, 54) (1, 816, 74) (0, 645, 62) (1, 777, 68) (1, 809, 47) (1, 761, 51) (0, 740, 60) (1, 816, 50) (1, 660, 53) (0, 711, 64) (1, 681, 61) (1, 650, 71) (1, 913, 68) (1, 616, 63) (0, 743, 50) (0, 687, 53) (0, 564, 63) (1, 620, 59) (1, 724, 61) (0, 650, 67) (0, 631, 69) (0, 748, 60) (0, 612, 58) (0, 522, 48) (0, 699, 61) (0, 561, 68) (1, 879, 59) (1, 863, 60) (1, 761, 55) (1, 803, 52) (1, 791, 56) (1, 786, 63) (1, 738, 57) (0, 707, 63) (1, 680, 49) (1, 695, 48) (1, 873, 67) (1, 580, 53) (1, 799, 57) (0, 614, 46) (0, 607, 55) (1, 783, 54) (1, 728, 47) (1, 734, 55) (0, 667, 79) (0, 721, 57) (1, 734, 78) |
表1のデータ(DataTable411.txt)に対しての事後分布を求め、その平均値、95% HPD区間、の値を求めると、以下のようになる。
ただし、事前分布は十分に広い範囲で一様分布とした。本Websiteにおける分析では、事前分布を十分に広い範囲で一様分布であるとする。また、本Websiteにおける計算は、このWebsiteの後半に用意されているプログラムで行うことができる。
上の結果では、係数の値は負であり、95% HPD 区間は0を含んでいない。の値は0の左側の確率を示すものであるが、係数の値が負であるので0の右側の確率が、有意に係数の値が0と異なるかどうかの基準になる。
であるので、係数の値は5%水準で有意に0と異なるといえる。上の式において「1.000」あるいは「0.000」と書かれているのは、小数点以下第4位で四捨五入すると小数点以下第3位までの数値が0になるということを表す。
係数の場合は、値が正であることから、上のの場合と逆の方向に見ていけばよい。0の左側の領域の確率であるので、係数の値は5%水準で有意に0と異なるといえる。すなわち、模擬試験の点数が高いほど合格確率も高くなるといえる。
係数の場合は、95% HPD 区間が0を含んでいるので、両側検定の考え方では、5%水準で有意に0と異なるとはいえない。0の両側の領域の確率も
であり、確率的にはの値がそれほど0から離れているとはいえない。以上の結果より、表1のデータにおいては、合格確率と体重との有意な関係は認められないといえる。
表1のデータは、個人単位(ケース単位)のものであるが、データが条件ごとの頻度で与えられている場合も同様に分析できる。
表2 数学・算数の授業法と数学・算数の好き嫌いの人数. |
|||
|
従来の教授法 |
新しい教授法 |
|
小学3年生 |
好き |
20 |
30 |
嫌い |
20 |
10 |
|
好きの比率 |
0.5 |
0.75 |
|
中学2年生 |
好き |
5 |
60 |
嫌い |
35 |
200 |
|
好きの比率 |
0.125 |
0.231 |
例えば、表2のデータ(DataTable421.txt)は、数学(算数)の2種類の授業法(従来の方法と新しい方法)の効果を比較したもの(人工データ)であるが、従来の方法による授業を受けた場合、小学3年生では20人が数学(算数)が好きと答え、20人が嫌いと答えたが、新しい授業法を受けた生徒は、30人が好きと答え、10人が嫌いと答えている。中学2年生の場合は、従来の授業法を受けた場合は5人が数学が好きと答え、35人が嫌いと答えたが、新しい授業法では60人が好きと答え200人が嫌いと答えている。好きと答えた者の比率は、小学3年生では従来の方法の場合は0.5であるのに対して、新しい方法の場合は0.75であり、中学2年生では従来の方法の場合は0.125であり、新しい授業法では0.231である。小学3年生、中学2年生とも、同じ学年同士の比較では、新しい授業法を受けた生徒の方が従来の授業法を受けた生徒より、数学・算数が好きと答えた者の比率が高い。このデータをロジスティック回帰モデルで分析してみる。
好きと答える確率をとおき、ロジスティック回帰モデルを
とおく。ここで、
とおいた。データがどの条件のものであるかを示し、それらの効果を表すために導入された変数、、はダミー変数と呼ばれている(岡本(2009b)、pp. 130-133)。
式(1)の各条件における値を表3に示す。
表3 ロジスティック回帰モデル(1)の値。 |
||
|
従来の授業法 |
新しい授業法 |
小学3年生 |
|
|
中学2年生 |
|
|
係数は小学3年生における従来の授業法を受けた場合の「好き」の比率に対する式(1)の値を表し、この条件での値を基準にして授業法の効果および学年の効果を表している。は授業法の主効果、は学年の主効果、は交互作用である。
モデル(1)を表2のデータに適用して、の値の事後分布における平均値、95%HPD区間、0以下の値を取る確率()を求めると以下のようになる(これらの分析は、このWebsiteの後半に用意されているプログラムで行うことができる)。
授業法の主効果の値は正であり、95%HPD区間は0を含まず、事後分布の0以下の領域の確率は5%よりも小さい。したがって有意水準5%で有意に正である。すなわち、新しい授業法を受けた生徒の方が「好き」と答える比率は有意に高い。
学年の主効果の値は負であり、95% HPD区間は0を含まず、事後分布における0以上の値の領域の確率はと5%よりも小さい。したがって有意水準5%で有意に負の値である。すなわち、学年が上がると有意に「数学が好き」と答える比率が下がっている。
交互作用の値は負であるが、95% HPD区間は0を含み、事後分布における0以下の領域および0以上の領域はいずれも5%より大きい。したがって、交互作用は有意ではない。
調整変数(moderator variable)とシンプソンのパラドックス(Simpson’s Paradox)
独立変数と従属変数の関係に影響を与える変数は、調整変数と呼ばれている。表2のデータの場合、学年が調整変数である。学年を無視してデータをまとめると、表4(DataTable441.txt)のようになる。
表4 授業法における数学の好き嫌いの人数(学年を無視) |
||
|
従来通り |
新しい教授法 |
好き |
25 |
90 |
嫌い |
55 |
210 |
好きの比率 |
0.31 |
0.30 |
表2の分析では、新しい教授法の方が有意に数学・算数が好きになる比率が高いことが認められたが、学年を無視して集計した表4では新しい授業法の方が「好き」の比率が小さくなっている。
表4のデータをロジスティック回帰モデルで分析するために、好きと答える確率に対して
とおく。ここで、
である。従来通りの授業法の場合の好きと答える確率は
で、新しい授業法の場合に好きと答える確率は
で表される。新しい授業法の効果がによって表されている。
モデル(2)を表4のデータに適用して係数の事後分布を求め、それらの平均値、95%HPD区間、0以下の値を取る確率()を求めると以下のようになる。
授業法の効果の値は負であるが、95%HDP区間が0を含むこと、事後分布における0以下および0以上の領域の確率が5%より大きいことから、有意に負の値であるとはいえない。すなわち、調整変数である学年を考慮した表2の分析では授業法の効果が有意であることが検出できたが、調整変数を無視した表4の分析では授業法の効果は認められない。
上の比率の場合における調整変数の効果は、次のように説明できる。
データが1あるいは0の2値の1つを取るものとする。問題としている要因をAとし、2つの水準A1とA2の間でデータの値が1である比率に差があるかどうかを調べるものとする。調整変数の要因をBで表し、2つの水準B1とB2があるものとする。B1のときのA1のデータ数をで表し、値が1であるものの比率をで表す。同様に、B1のもとでのA2のデータ数をで、値が1であるものの比率をで表す。B2のもとでも同様に、A1(A2)の条件でのデータ数を()で表し、データの値が1であるものの比率を()で表す。このとき、調整変数Bを無視してデータをまとめると、条件A1でのデータ数はであり、データの値が1であるものの個数はである。したがって、1であるものの比率は
となる。同様に、A2のもとでの値が1であるデータの比率は
となる。
比率は、との加重平均である。したがって、サンプル数との比率を調整することによりとの間の任意の値をとりうる。も同様にサンプル数との比率を調整することによりとの間の任意の値をとりうる。このため、例えば
であっても、
と関係が逆になることもありうる(シンプソンのパラドックス)。
表5 数学・算数の授業法と数学・算数の好き嫌いの人数(仮想データ). |
|||
|
従来の教授法 |
新しい教授法 |
|
小学3年生 |
好き |
200 |
280 |
嫌い |
200 |
120 |
|
好きの比率 |
0.5 |
0.7 |
|
中学2年生 |
好き |
40 |
1200 |
嫌い |
160 |
2800 |
|
好きの比率 |
0.2 |
0.3 |
表5のデータ(DataTable442.txt)は、調整変数である学年を無視すると授業法の効果に有意な逆転が見られる例である。式(1)のモデルで分析すると以下の結果をえる。
授業法、学年の効果が有意で、新しい授業法を受けた方が「好き」の比率が高く、中学2年生の方が「好き」と答える比率が低い。交互作用は有意ではない。
表5のデータを学年の違いを無視してまとめると表6(DataTable443.txt)になる。
表6 授業法における数学の好き嫌いの人数(学年を無視) |
||
|
従来通り |
新しい教授法 |
好き |
240 |
1480 |
嫌い |
360 |
2920 |
好きの比率 |
0.4 |
0.34 |
表6では新しい授業法を受けた条件の方が「好き」と答える比率が低い。表6のデータをモデル(2)で分析すると
となる。の値は負であり、95% HPD区間は0を含まない。の事後分布における0以上の領域の確率は
である。したがって、は5%水準で有意に負の値であるといえる。すなわち、新しい授業法を受けた群の方が「好き」と答える生徒の比率が有意に低いといえる。調整変数である学年を考慮した場合と無視した場合とで、逆の有意な差が見られる(シンプソンのパラドックス)。
単回帰モデルあるいは重回帰モデルでは、従属変数は連続量の変数であって、あらゆる実数値をとりうるものである。従属変数が、「はい/いいえ」などのカテゴリ変数の場合、あるいは比率のように0から1までの値しかとらない場合は、従属変数が連続量であるモデルをそのまま適用するのではなく、従属変数に合わせたモデルが用いられる。
従属変数が連続量である場合の重回帰モデルは、次式で与えられる。
これに対して、従属変量が確率である場合は、その変換値
を重回帰モデル
で表す。ここで、変換関数は次式
で与えられるものであり、ロジット変換(logit transformation)と呼ばれている。ロジット変換による回帰モデル(4)は、ロジスティック回帰モデル(logistic regression model)と呼ばれている。確率のロジスティック変換による変換値は、−∞から+∞までの値をとりうる。ここで、
はオッズ(odds)と呼ばれているもので、ロジット変換とはオッズの対数値を与えるものである。の値が0から1まで変わるとき、オッズの値は0から+∞まで変わるので、このときロジット変換の値は−∞から+∞まで変わる。
式(4)の係数の大きさが確率に及ぼす影響は、以下のように評価できる。
式(4)において、の値を1増やしてとしたときのの値をとおく。
である。このとき、
となる。式(5)より
を得る。したがって、
となる。すなわち、係数は、独立変数の値が1増加したときの確率の変化量をオッズ比の対数で表したものに等しい。
確率の変換を、プロビット変換(probit transformation)と呼ばれている次の変換
によって行い、値をモデル(4) においての代わりに用いたものは、プロビット回帰モデルと呼ばれている。ここで、は標準正規累積分布関数の逆関数であり、確率の値が0から1まで動くとき、の値は−∞から+∞まで動く。
ここでは、ロジット変換を用いたロジスティック回帰分析について説明する。
モデル(4) に対しては、重回帰分析の仮定「」に基づいた分析が適用できないので、最尤法においてデータ数が十分に多いときに成り立つ漸近定理に基づいて分析が行われる(Kutner et al., 2005)。しかし、心理学におけるデータの場合、漸近定理に基づく分析が可能なだけの十分な数のデータが収集されているかどうか疑わしいことが多い。最尤法とは異なり、データ数に依存しない分析法としてベイズ的方法がある。このベイズ的方法による分析について説明する。なお、ベイズ的方法については、岡本(2009a、pp. 153-157)、岡本(2009b、pp. 160-167)、岡本(2006、pp. 252-262)、および本書「心理学データ分析と測定」の「付録C ベイズ的分析法」を参照されたい。
ベイズ的方法では、データから得られる情報は事後分布として表される。事後分布は、事前分布と尤度関数の積に比例する形で与えられる。データが2値をとる従属変数の場合、それらの値を0と1で表しても一般性は失われない。従属変数が次の確率で1あるいは0の値をとるとき、
尤度関数は次式で与えられる。
ここで、は式(4)により決まるもので
である。
このとき、事前分布を、事後分布をとおくと
を得る。式(7)で与えられた事後分布は、MCMCというシミュレーションによって推定することができる。MCMCについては、岡本(2008、pp. 110-113)、岡本(2006、pp. 257-261)、および本書「心理学データ分析と測定」の「付録C ベイズ的分析法」を参照されたい。
モデル(6)はデータが1つ1つ個別に与えられた場合のものであるが、表2のように条件別にそれぞれいくつのケースが値1あるいは0であるか個数の形式でデータが集められた場合は以下のようになる。
いま、条件のもとでデータが個集められ、そのうち個が値1のものであったとする。このとき、尤度関数は
で与えられる。
1つ1つのデータが1または0の値をとる従属変数の各値に対してデータが与えられている表1のような場合は、プログラムLogisticRegIndividual.exeで分析することができる。プログラムは、プログラム名LogisticRegIndividual.exeをクリックしてダウンロード・実行することができる。実行開時に、「*.dllファイルがインストールされていません」というようなエラーメッセージが表示されたときは、Visual Studio Express for Windows Desktopを
http://www.microsoft.com/visualstudio/jpn/downloads
からダウンロード・インストールすれば、必要なファイルがインストールされ、プログラムLogisticRegIndividual.exeがダウンロード・実行できるようになる。なお、Windows 8/8.1におけるダウンロード・実行の手順については、特にウェブサイトを用意して説明した。
ソース・コード・ファイルなどは圧縮ファイルLogisticRegIndividual.zipとしてまとめられているので、圧縮ファイル名をクリックしてダウンロードすることができる。ダウンロードした圧縮ファイルは、圧縮ファイル名LogisticRegIndividual.zipをマウスの右ボタンでクリックして表示されるメニュ「解凍」あるいは「展開」などを選ぶと解凍される。解凍されたプログラムファイルは、Visual Studio 2013で開くことができる。Visual C++プログラミングの入門書として<岡本安晴「大学生のための心理学VC++プログラミング入門」勁草書房>を用意した。
プログラムを実行すると、下図のフォームが表示される。
「GO」ボタンのクリックで、入力データファイル名の設定を求めるダイアログ・フォームが表示される。
入力データファイルは、表1の場合、次の形式で用意する。
==========================
/
3
Result Const Score Weight
1
1
790
67
0
1
626
57
1
1
869
53
1
1
688 54
0
1
687
53
・
・
・
1
1
783
54
1
1
728
47
1
1
734
55
0
1
667
79
0
1
721
57
1
1
734
78
/
/
1
1 0 0
/
1
0 1 0
/
-1
==========================
スラッシュ/で始まる最初の行までは、無視される(読み飛ばされる)。
このスラッシュ/で始まる最初の行の次の行に、独立変数の数pを書く。
その次の行には、変数のラベルを、従属変数、独立変数1、・・・、独立変数pの順に並べる。
このラベルを書いた行の次の行から、1行に1組のデータを、従属変数、独立変数1、・・・、独立変数pの値の順に書いていく。
すべての組のデータを書いた最後は、スラッシュ/で始まる行を置いて、データの終わりを表す。
次のスラッシュ/で始まる行の次の行に書かれている整数値が0以下の値であれば、プログラムの実行終了となる。
正の値が書かれているときは、その次の行に各独立変数に対するウェイトを並べる。例えば、
と書かれているときは、係数の重み付き和
の事後分布が求められ、その平均値、中央値、95% HPD区間、p値が算出される。
上の例の場合、ウェイト
1 0 0
に対しては
の事後分布が求められ、ウェイト
0 1 0
に対しては、
の事後分布が求められる。
入力ファイル名を設定して「開く」ボタンをクリックすると、次に出力ファイル名の設定を求めるダイアログ・フォームが表示される。
出力ファイルはテキストファイルとして作成されるので、ファイル拡張子は「.txt」としておく。出力ファイル名の設定後、「保存」ボタンをクリックすると計算が始まる。計算が終了すると、出力ファイルが保存された旨を通知するメッセージ・ボックスが表示される。
「Ok」ボタンのクリックで、プログラムの実行終了となる。実行終了後、出力ファイルはテキスト・エディタで開くことができる。表1のデータを上の入力ファイルの説明のように用意した場合、出力ファイルは以下のようになる。
=============================
Input File = F:\yasuharu\www\LaCoocan\booksetc\pm2010\LogisticRegression\LogisticRegIndividual\LogisticRegIndividual\Data.txt
Y -> Result
X[0] -> Const
X[1] -> Score
X[2] -> Weight
n = 100
Data...
Y
X[0]
X[1]
X[2]
b[0] -> -12.70369/mean -12.56033/median 95% HPD
Interval:[-18.85963, -6.71347]
P(v <= 0.0) = 1.00000
b[1] -> 0.02041/mean
0.02015/median
95% HPD Interval:[0.01258, 0.02850] P(v <= 0.0) = 0.00000
b[2] -> -0.02307/mean
-0.02219/median
95% HPD Interval:[-0.08876, 0.03659] P(v <= 0.0) = 0.76761
Weight...
1 0 0
-12.70369/mean
-12.56033/median
95% HPD Interval:[-18.85963, -6.71347] P(v <= 0.0) = 1.00000
Weight...
0 1 0
0.02041/mean
0.02015/median
95% HPD Interval:[0.01258, 0.02850] P(v <= 0.0) = 0.00000
=============================
上の出力結果の解釈については、本書「心理学データ分析と測定」の「第4章 比率(確率)の条件差を見る」を参照されたい。
表2のように、データが条件ごとの集計の形式で与えられているときは、プログラムLogisticRegGroup.exeで分析することができる。プログラムはプログラム名LogisticRegGroup.exeをクリックしてダウンロード・実行することができる。実行開始時に、「*.dllファイルがインストールされていません」というようなエラーメッセージが表示されたときは、Visual Studio Express for Windows Desktopを
http://www.microsoft.com/visualstudio/jpn/downloads
からダウンロード・インストールすれば、必要なファイルがインストールされ、プログラムLogisticRegGroup.exeがダウンロード・実行できるようになる。なお、Windows 8/8.1におけるダウンロード・実行の手順については、特にウェブサイトを用意して説明した。
ソース・コード・ファイルなどは圧縮ファイルLogisticRegGroup.zipとしてまとめられている。圧縮ファイル名をクリックしてダウンロードすることができる。ダウンロードした圧縮ファイルは、圧縮ファイル名LogisticRegGroup.zipをマウスの右ボタンでクリックして表示されるメニュ「解凍」あるいは「展開」などを選ぶと解凍される。解凍されたプログラムファイルは、Visual Studio 2013で開くことができる。
プログラムを実行すると下図のフォームが表示される。
「GO」ボタンのクリックで、入力データファイル名の設定を求めるダイアログ・フォームが表示される。
入力データファイルは、表2の場合、次の形式で用意する。
==========================
Y=1の数 Y=0の数 独立変数・・・
/
4
Math(+) Math(-)
Const Dummy1 Dummy2 Dummy3
20
20
1
0 0 0
30
10
1 1 0 0
5
35
1
0
1 0
60
200
1
1
1
1
/
/
1
0 1 0 1
/
-1
==========================
スラッシュ/で始まる最初の行までは無視される(読み飛ばされる)。
このスラッシュ/で始まる最初の行の次の行に独立変数の数pを書く。
その次の行には、変数のラベルを、従属変数Y(1のデータ数を表すもの)、従属変数N(0のデータ数を表すもの),独立変数1、・・・、独立変数pの順に並べる。
このラベルを書いた行の次の行から、1行に1組のデータを、従属変数Y(1のデータ数)、従属変数N(0のデータ数)、独立変数1、・・・、独立変数pの値の順に書いていく。
すべての組のデータを書いた最後は、スラッシュ/で始まる行を置いて、データの終わりを表す。
次のスラッシュ/で始まる行の次の行に書かれている整数値が0以下の値であれば、プログラムの実行終了となる。
正の値が書かれているときは、その次の行に各独立変数に対するウェイトを並べる。例えば、
と書かれているときは、係数の重み付き和
の事後分布が求められ、その平均値、中央値、95% HPD区間、p値が算出される。
上の入力データファイルの例における重み
0 1 0 1
に対応する重み付きパラメタの和は
を表す。これは、中学2年生における新しい授業法と従来の授業法の差
を表すものである(表3参照)。
入力ファイル名を設定して「開く」ボタンをクリックすると、次に出力ファイル名の設定を求めるダイアログフォームが表示される。
出力ファイルはテキストファイルとして作成されるので、ファイル拡張子は「.txt」としておく。出力ファイル名の設定後、「保存」ボタンをクリックすると計算が始まる。計算が終了すると、出力ファイルが保存された旨を通知するメッセージボックスが表示される。
「Ok」ボタンのクリックでプログラムの実行終了となる。実行終了後、出力ファイルはテキスト・エディタで開くことができる。表2のデータを上の入力ファイルの説明のように用意した場合、出力ファイルは以下のようになる。
=============================
Input File =
F:\yasuharu\www\LaCoocan\booksetc\pm2010\LogisticRegression\LogisticRegGroup\LogisticRegGroup\DataMdr2.txt
Y_P -> Math(+)
Y_N -> Math(-)
X[0] -> Const
X[1] -> Dummy1
X[2] -> Dummy2
X[3] -> Dummy3
n = 4
Data...
Y_P
Y_N
X[0]
X[1]
X[2]
X[3]
b[0] -> -0.00123/mean -0.00054/median 95% HPD Interval:[-0.62859,
0.63074] P(v <= 0.0) =
0.50072
b[1] -> 1.13625/mean 1.12987/median 95% HPD
Interval:[0.18136, 2.11460]
P(v <= 0.0) = 0.00919
b[2] -> -2.02913/mean -2.00535/median 95% HPD
Interval:[-3.20173, -0.87131]
P(v <= 0.0) = 0.99995
b[3] -> -0.31604/mean -0.33049/median 95% HPD
Interval:[-1.70589, 1.11333]
P(v <= 0.0) = 0.67752
Weight...
0 1 0 1
0.82022/mean 0.79114/median 95% HPD
Interval:[-0.17132, 1.86460]
P(v <= 0.0) = 0.04668
=============================
上の出力結果の解釈については、本書「心理学データ分析と測定」の「第4章 比率(確率)の条件差を見る」を参照されたい。
引用文献
Kutner, M. H., Nachtsheim, C. J., Neter, J. and Li, W. (2005). Applied linear statistical models, 5th ed. McGraw-Hill.
岡本安晴(2008).統計学を学ぶための数学入門[上]―算数から数学へ―.培風館.
岡本安晴(2009a).統計学を学ぶための数学入門[下]―データ分析に活かす―.培風館.
岡本安晴(2009b).データ分析のための統計学入門―統計学の考え方―.おうふう.
岡本安晴、2011.8, 2014.6