Pythonでt検定　2クラスの試験成績の比較

# 概要

「対応なしt検定」に関するメモ、実施のためのPythonコードです。
「対応なしt検定」のほか、その準備として「シャピロ・ウィルク検定」、「ルビーン検定」、「区間推定」に関する内容が含まれます。

### 関連エントリ
- [F#でt検定 2クラスの試験成績の比較](https://qiita.com/code0327/items/960c20ba51f9cf5430e7)

# 設定

問題例として、次のような設定（状況）を考えていきたいと思います。

>同質の2つのクラス（成績が均等になるように学生が振り分けられた学級）に対して、それぞれ異なる学習指導法を実施した（例えば、2つのクラスをA組、B組として、A組にはアクティブラーニング方式、B組には従来方式で授業を行なったなど…）。
>学期末に、各クラスで同じ内容の試験を実施し、その成績データが得られている。このデータから学習指導法の違いによる影響・効果を判断したい。

（補足）上記の`このデータから学習指導法の違いによる影響・効果を判断したい。` とは、厳密には「A組・B組と同質の学生らに対して、これらの異なる学習指導法を実施したら、試験成績に影響がでるのかを判断したい」ということになります。

　さて、学習指導法が試験成績に影響を与えているか？を測るには、通常、各クラスの成績の平均点を求めて比較をします。このとき、平均点に10点なり、15点なりの差異が確認できれば、それを根拠にどちらかの指導法の有効性を主張することができます。しかし、それでは「**たまたま、偶然による成績のばらつきの範囲内で差がでているだけでは？**」のように言われてしまうかもしれません。

　このとき、**統計学的仮説検定**（単に「**検定**」とも呼ばれる）を実施していれば「この差が偶然による成績のばらつきに起因するとは**統計学的に考えずらい**」ということが言えます（あるいは逆に「**統計学的にみると**偶然によるばらつきの可能性は捨てられない」との結論になることもあります）。いずれにしても、検定をすれば統計学的な裏付けを持って「ある学習指導法が有効／無効であるという判断」に説得力を持たせることができるようになります。なお、統計学では偶然によらない「意味を持つ差」のことを**有意差**と言います。

　検定には、様々な検定があって、**その目的**や**前提とする条件**により適切に使い分ける必要があります。ここでは、**2群の標本の平均に有意差があるかどうか**の検討に用いる「**t検定**」を取り上げていきます。

# 対応なしt検定

　先に問題例として挙げたものは「同質の2つの対象に対して異なる条件を与えた」というタイプのものなので、t検定のなかでも「**対応なしt検定**」というものを使います。

　一方で、同一の対象に対しての変化を見るような場合は「**対応ありt検定**」を使用します。例えば、英語コミュニケーション能力について、短期留学の前と後で変化しているか？などを検討するときに利用します。

　以下では、「**対応なしt検定**」を使って2つの学習指導法に有意差があるかどうかを評価していく手順をまとめています。

# 1. 実質的な差の評価

　A組の試験平均点 $\bar{x}\_A$、B組の試験平均点 $\bar{x}\_B$ を計算します。ここで、両者に**実質的な差がなければ**「2つの学習指導法には実質的な差はない」と判断して、「統計学的に有意差があるかどうか」のt検定を進めることは**終了**してもよいと思います。実質的な差がほぼ無くとも（十分な大きさの標本を用いていれば）検定の結果として「統計学的に有意差がある」という結論がでてくることもありますが、それはあまり意味を持ちません。

　例えば、100点満点で $\bar{x}\_A=89.0点$、$\bar{x}\_B=87.7点$ となったようなケースです（平均点の差は $1.3$ 点です）。こんなケースでも検定の結果、統計学的には有意差が認められる場合がありますが、それだけをもって「（A組で導入した）アクティブラーニングは試験成績向上に有効なので皆さん導入しましょう！」のような主張するのは厳しいですよね。労力をかけて検定する必要はないと思います。ただし「本当は差があるのだが、偶然のばらつきによって、今回は不運にも差がない結果となっただけだ」という可能性もあります。学習指導法による違いはきっとある！という信念があれば、翌年度にもう一度、標本を取り直すといったことができると思います。

　この処理を行なうPythonコードを以下に示します。各クラスの試験成績はpandasのSeriesに格納されていることを想定しています。なお、サンプル掲載の成績は正規乱数を使って生成した仮想的なものです（生成方法は一番最後に記載しています）。

```Python:標本平均の計算と比較（Python）
import math
import pandas as pd

xa = pd.Series( [75, 87, 89, 80, 84, 81, 88, 83, 88, 88, 82, 72, 74, 93, 77, 67, 88, 84, 68, 84, 80, 78, 75, 71, 82, 74, 84, 77, 79, 76, 83, 75, 86, 76, 80, 76, 68, 72, 75, 85])
xb = pd.Series( [64, 77, 79, 73, 89, 82, 59, 85, 80, 75, 65, 79, 65, 74, 73, 72, 69, 83, 90, 73, 88, 59, 62, 80, 64, 74, 81, 70, 69, 67, 81, 67, 72, 71, 72, 78, 78, 82, 72, 71])

print(f'A組の平均点 = {xa.mean():.2f}')
print(f'B組の平均点 = {xb.mean():.2f}')
print(f'2群の平均の差 = {(abs(xa.mean()-xb.mean())):.2f}')
```

```:実行結果
A組の平均点 = 79.60
B組の平均点 = 74.10
2群の平均の差 = 5.50
```

　なお、この時点で母集団について、いくつかの推定をすることができます。まずは各標本の母集団の平均 $\mu_A$ と $\mu_B$ の推定値（点推定）である $\hat{\mu}\_A$ と $\hat{\mu}\_B$ です。標本Aの母集団とは「（指導を受ける前の）A組と同質の学生らに対して、A組と同じ学習指導法を適用した後に、今回と同じ試験を受けさせたときの試験成績の集合」になります。

```math
\text{標本Aの母平均推定値}\quad \hat{\mu}_{A} = \bar{x}_{A} = \frac{1}{n}\sum_{i=1}^{n} x_{A_i}\\
\text{標本Bの母平均推定値}\quad \hat{\mu}_{B} = \bar{x}_{B} = \frac{1}{n}\sum_{i=1}^{n} x_{B_i}
```

　また、母集団の分散 $\sigma^2$ の推定値（点推定）である 不偏分散 $u^2$ についても計算できます。なお、母集団ではなくて、標本そのものの分散を計算する場合は、下記の計算式の分母の $(n-1)$ を $n$ に置き換えます。

```math
\text{標本Aの不偏分散}\quad u_A^2 = \frac{1}{n-1}\sum_{i=1}^{n} (x_{A_i}-\bar{x}_A)^2\\
\text{標本Bの不偏分散}\quad u_B^2 = \frac{1}{n-1}\sum_{i=1}^{n} (x_{B_i}-\bar{x}_B)^2
```

　これらを求めるPythonコードを以下に示します。

```Python:母集団の推定（Python）
import math
import pandas as pd

xa = pd.Series( [75, 87, 89, 80, 84, 81, 88, 83, 88, 88, 82, 72, 74, 93, 77, 67, 88, 84, 68, 84, 80, 78, 75, 71, 82, 74, 84, 77, 79, 76, 83, 75, 86, 76, 80, 76, 68, 72, 75, 85])
xb = pd.Series( [64, 77, 79, 73, 89, 82, 59, 85, 80, 75, 65, 79, 65, 74, 73, 72, 69, 83, 90, 73, 88, 59, 62, 80, 64, 74, 81, 70, 69, 67, 81, 67, 72, 71, 72, 78, 78, 82, 72, 71])

print(f'標本Aの母集団の平均の推定値 = {xa.mean():.2f}')
print(f'標本Bの母集団の平均の推定値 = {xb.mean():.2f}')
print(f'標本Aの母集団の標準偏差の推定値（不偏標準偏差）= {xa.std(ddof=1):.2f}')
print(f'標本Bの母集団の標準偏差の推定値（不偏標準偏差）= {xb.std(ddof=1):.2f}')
```

pandasでは標準偏差を求める `std()` がデフォルトで `std(ddof=1)` になります。一方、numpyではデフォルトでは `std(ddof=0)` になる点に注意します。この ddof (Delta Degrees of Freedom、**自由度**) は、標準偏差の計算の際、以下のように反映されます。

```math
\sqrt{\frac{1}{n-\text{ddof}}\sum_{i=1}^n (x_i-\bar{x})^2}
```

```:実行結果
標本Aの母集団の平均の推定値 = 79.60
標本Bの母集団の平均の推定値 = 74.10
標本Aの母集団の標準偏差の推定値（不偏標準偏差）= 6.42
標本Bの母集団の標準偏差の推定値（不偏標準偏差）= 7.89
```

# 2. 対応なしt検定の前提条件のチェック

　平均点に実質的な差（例えば、$\bar{x}\_A=79.60点$、$\bar{x}\_B=74.10点$）があっても、それは「たまたま」「偶然」「まぐれ」による誤差である可能性があります。この平均点の差が「偶然によるとは到底言えないこと」なのか「偶然によるものと言えてしまうこと」なのかを「対応なしt検定」によって定量的に見ることができます。

　ただし「対応なしt検定」を適用するためには、標本（＝ここでは各クラスの成績データ）について、以下の条件を満たす必要があります。

1. 2つの母集団から標本が無作為に抽出（ランダムサンプリング）されていること。
2. 比較する2群のデータ（=各標本）の母集団が正規分布、またはそれに近い分布となっていること（正規性）。
3. 比較する2群のデータ（=各標本）の分散が等しいこと（等分散性）。

1番目については「2つのクラスが同質」と言えれば、「A組については上位10名の成績を標本に、B組については下位10名を標本に使った」ようなことをしなければ問題ありません。2と3番目についてチェックする方法は、以下で詳細を示します。

### 補足
いまさらですが、t検定を行なう場合、データの種類が連続性のある「間隔尺度」または「比例尺度」である必要があります。「名義尺度」「順序尺度」の場合は適用できません。

## 2-1 正規性のチェック

　正規母集団（＝正規分布に従う母集団）から無作為に多数の標本を抽出すると、その標本の分布もまた正規分布になる、ということが知られています。よって、正規性は、各標本（各群のデータに）について**ヒストグラム**を作成して、その形状から判断することができます。ただし、標本数が少ないと（30未満ぐらい？）と形状から判断することは難しいです。ちなみに、ここで標本として使用している成績データは、各クラスとも40個となっています。

```Python:正規性チェックのためのヒストグラム描画（Python）
import pandas as pd
import matplotlib.pyplot as plt

xa = pd.Series( [75, 87, 89, 80, 84, 81, 88, 83, 88, 88, 82, 72, 74, 93, 77, 67, 88, 84, 68, 84, 80, 78, 75, 71, 82, 74, 84, 77, 79, 76, 83, 75, 86, 76, 80, 76, 68, 72, 75, 85])
xb = pd.Series( [64, 77, 79, 73, 89, 82, 59, 85, 80, 75, 65, 79, 65, 74, 73, 72, 69, 83, 90, 73, 88, 59, 62, 80, 64, 74, 81, 70, 69, 67, 81, 67, 72, 71, 72, 78, 78, 82, 72, 71])

fig = plt.figure(dpi=96)
ax1 = fig.add_subplot(211) 
ax1.set_xlim(55, 100)
ax1.set_ylim(0, 15)
ax1.hist(xa,bins=9,range=(55,100),color='r')

ax2 = fig.add_subplot(212) 
ax2.set_xlim(55, 100)
ax2.set_ylim(0, 15)
ax2.hist(xb,bins=9,range=(55,100),color='b')

plt.show()
```

実行結果は次のようになります。見た感じでは、どうやらいずれも正規分布に近い分布であり、それぞれの母集団も正規分布に従っていると言えそうです。

![無題.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/152805/22b3133e-9371-5ed3-96c6-4fa43b2e84c4.png)

　また、正規性は検定によりチェックすることもできます。正規性の検定には、いくつかの方法があり、ここでは**シャピロ・ウィルク検定**（Shapiro-Wilk検定）を使ってみます（標本数が少ないと分布に正規性があるとみなされる傾向が強く、判断には注意が必要となります）。

　シャピロ・ウィルク検定は「標本は、正規分布に従う母集団から無作為抽出されたものである」を帰無仮説としています。今回は**有意水準を $5\%$（$=0.05$）**とします。そして、検定により求まる**有意確率 $p$ が $5\%$ よりも大きいならば**「母集団は正規分布に従う」という仮説を採択します。ここでの有意確率 $p$ の具体的な意味は「正規母集団から無作為抽出したときに、今回の標本分布（＋今回の標本よりさらに正規分布と言えない標本分布）が観測される確率」を意味します。

　最初に「正規母集団から無作為に多数の標本を抽出すると、その標本の分布もまた正規分布になる」と書きましたが、標本数が限られていると（偶然のばらつきにより）正規分布とは言えないような分布が観測される場合もあります。次の図は、標準正規分布から40個を無作為抽出することを1セットとして、18セット実施したものです。これらの標本セットは概ね正規分布のようになっていますが、なかには#10のように正規性があるとは言えないような標本も得られています。

![無題.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/152805/8cc7e7fc-3202-d026-9302-9c82bc168a3e.png)

　今回の例において、A組の標本（標本数40）に対する検定の結果、仮に $p=50\%$ が得られれば、それは「正規母集団から40個の標本を無作為抽出する行為を1セットとして、それを何度も何度も行なうと『A組のような分布の標本』と『それ以上に乱れた分布の標本』が観測される確率は50%である」という意味となります。つまり、正規母集団から「A組のような分布の標本」が得られることは十分に起こり得ることであり、A組の標本の母集団は正規分布に従うとと考えても差し支えない、という結論になります。
　一方で、$p=0.01\%$ であれば、正規母集団から「A組のような分布の標本」が得られることは極めて稀であり偶然に起こるとは考えずらい、よってA組の標本の母集団は正規分布とは言えない、という結論になります。
　これらの結論の境界を定めたものが、**有意水準**であり分野によって$5\%$や、$1\%$、$0.1\%$などが慣例的に決まっているようです。

　まとめると以下のようになります。シャピロ・ウィルク検定により、$p$ 値（有意確率）を計算して有意水準（$5\%$）よりも・・・

- 大きければ（ $p\ge5\%$ なら）**正規性がある**と判断（帰無仮説を**採択**）
- 小さければ（ $p<5\%$ なら）**正規性がない**と判断（帰無仮説を**棄却**）

・・・します。正規性がない場合は「マン・ホイットニーのU検定」などの別の検定を利用することを検討します（ただ、参考資料[3]によれば、t検定は正規性に対して頑健であり、大きく外れていなければ、そのまま利用可能らしいです）。

 シャピロ・ウィルク検定を行なうPythonコードを以下に示します。

```Python:シャピロ・ウィルク検定による正規性のチェック（Python）
import pandas as pd
import scipy.stats as st

xa = pd.Series( [75, 87, 89, 80, 84, 81, 88, 83, 88, 88, 82, 72, 74, 93, 77, 67, 88, 84, 68, 84, 80, 78, 75, 71, 82, 74, 84, 77, 79, 76, 83, 75, 86, 76, 80, 76, 68, 72, 75, 85])
xb = pd.Series( [64, 77, 79, 73, 89, 82, 59, 85, 80, 75, 65, 79, 65, 74, 73, 72, 69, 83, 90, 73, 88, 59, 62, 80, 64, 74, 81, 70, 69, 67, 81, 67, 72, 71, 72, 78, 78, 82, 72, 71])

_, pa = st.shapiro(xa)
_, pb = st.shapiro(xb)

print('シャピロ・ウィルク検定')
print(f'標本A p値 = {pa:.3f}')
print(f'標本B p値 = {pb:.3f}')
```

```:実行結果
シャピロ・ウィルク検定
標本A p値 = 0.558
標本B p値 = 0.747
```

　実行例の場合、$p$ 値はいずれも $p\ge0.05$ $(5\%)$  なので、**帰無仮説を採択**して「各標本の母集団には正規性がある」と判断します。

　なお、この母集団が正規分布に従うということが言えると、各標本の母集団について平均値の**区間推定**（母平均の区間推定）ができます。よく利用されるのが**95%信頼区間**というもので、これは「$95\%$ の確率で母平均が存在する範囲」を意味します。これを統計学的に求めることができます。母集団が正規分布に従うときに、母平均の $95\%$ 信頼区間を求めるPythonコードを以下に示します。

```Python:母平均の95%信頼区間を求める（Python）
import math
import pandas as pd
import scipy.stats as st

xa = pd.Series( [75, 87, 89, 80, 84, 81, 88, 83, 88, 88, 82, 72, 74, 93, 77, 67, 88, 84, 68, 84, 80, 78, 75, 71, 82, 74, 84, 77, 79, 76, 83, 75, 86, 76, 80, 76, 68, 72, 75, 85])

u2 = xa.var(ddof=1)  # 母集団の分散推定値（不偏分散）
m = xa.mean()        # 標本平均
DF = len(xa)-1       # 自由度
SE = math.sqrt(u2/len(xa)) # 標準誤差

CI1,CI2 = st.t.interval( alpha=0.95, loc=m, scale=SE, df=DF )
print(f'標本Aの母平均の95%信頼区間CI = [{CI1:.2f} , {CI2:.2f}]')
```

```:実行結果
標本Aの母平均の95%信頼区間CI = [77.55 , 81.65]
```
 
　実行例の場合、母平均 $\mu$ が$95\%$の確率で、$77.55$点 から $81.65$点 の範囲にあることが分かりました。母集団の分布が分からないと母平均は点推定（$\bar{x}\_A=79.60$点）しかできませんが、分布がわかると区間推定までできるようになります。

## 2-2 等分散性のチェック

　等分散性は**ルビーン検定**（Levene検定）により判断することができます。この検定は、2つの母分散が等しいことを帰無仮説とします。ここでも、有意水準 $5\%$ で検定を行ないます。

 　ルービン検定を行なうPythonコードを以下に示します。

```Python:ルビーン検定による等分散性のチェック（Python）
import pandas as pd
import scipy.stats as st

xa = pd.Series( [75, 87, 89, 80, 84, 81, 88, 83, 88, 88, 82, 72, 74, 93, 77, 67, 88, 84, 68, 84, 80, 78, 75, 71, 82, 74, 84, 77, 79, 76, 83, 75, 86, 76, 80, 76, 68, 72, 75, 85])
xb = pd.Series( [64, 77, 79, 73, 89, 82, 59, 85, 80, 75, 65, 79, 65, 74, 73, 72, 69, 83, 90, 73, 88, 59, 62, 80, 64, 74, 81, 70, 69, 67, 81, 67, 72, 71, 72, 78, 78, 82, 72, 71])

_, p = st.levene(xa,xb,center='mean')

print('ルビーン検定')
print(f'p値 = {p:.3f}')
```

```:実行結果
ルビーン検定
p値 = 0.272
```

　実行例の場合、$p$ 値（有意確率）が $p\ge0.05$ $(5\%)$ なので、帰無仮説を採択して、2つの標本には等分散性があると判断します。

等分散性があるといえる場合は、**スチューデントのt検定** `st.ttest_ind(..., equal_var=True)` を使用します。また、等分散性がない場合は、**Welchのt検定**（Welch's method、**ウェルチの方法**）`ttest_ind(..., equal_var=False)` を用います。

## **<font color='red'>注意：検定の多重性</font>**

ここまで、正規性のチェックのために「シャピロ・ウィルク検定」、等分散性のチェックのために「ルービン検定」を適用し、さらに、このあと「対応なしt検定」を重ねて適用しようとしているわけですが、ここで**検定の多重性**という問題がでてきます。

- 引用:「[28-4. Welchのt検定](https://bellcurve.jp/statistics/course/9936.html)」＠BellCurve 統計WEB より・・・

> 統計学では検定を繰り返し行うと「多重性の問題」が生じるため、最近では2標本のt検定を行う場合には等分散性の検定は行わず、等分散かどうかを考慮する必要のない「Welchのt検定」を行ったほうが良いという考え方も一般的になりつつあります。 


# 3. 対応なしt検定

　対応なしt検定（対応なしのスチューデントのt検定）は、**2つの標本の平均値に有意差がない**こと帰無仮説とします。現在の問題設定では、これが棄却されること（＝学習指導法によって有意差がでること）を期待しています。

 　対応なしt検定を行なうPythonコードを以下に示します。

```Python:対応なしt検定（Python）
import math
import pandas as pd
import scipy.stats as st

xa = pd.Series( [75, 87, 89, 80, 84, 81, 88, 83, 88, 88, 82, 72, 74, 93, 77, 67, 88, 84, 68, 84, 80, 78, 75, 71, 82, 74, 84, 77, 79, 76, 83, 75, 86, 76, 80, 76, 68, 72, 75, 85])
xb = pd.Series( [64, 77, 79, 73, 89, 82, 59, 85, 80, 75, 65, 79, 65, 74, 73, 72, 69, 83, 90, 73, 88, 59, 62, 80, 64, 74, 81, 70, 69, 67, 81, 67, 72, 71, 72, 78, 78, 82, 72, 71])

t, p = st.ttest_ind(xa, xb, equal_var=True)
MU = abs(xa.mean()-xb.mean())
SE =  MU/t
DF = len(xa)+len(xb)-2
CI = st.t.interval( alpha=0.95, loc=MU, scale=SE, df=DF )

print('対応なしt検定')
print(f'p値 = {p:.3f}')
print(f't値 = {t:.2f}')
print(f'平均値の差   = {MU:.2f}')
print(f'差の標準誤差 = {SE:.2f}')
print(f'平均値の差の95%信頼区間CI = [{CI[0]:.2f} , {CI[1]:.2f}]')
```

```:実行結果
対応なしt検定
p値 = 0.001
t値 = 3.42
平均値の差   = 5.50
差の標準誤差 = 1.61
平均値の差の95%信頼区間CI = [2.30 , 8.70]
```

　まず、実行例の場合、$p$ 値（有意確率）が $p<0.05$ $(5\%)$ なので「2つの標本の平均値に有意差がない」としていた帰無仮説が**棄却**されます。つまり、t検定の結果として「2つの標本の平均差は（統計学的に見て）有意差があり、アクティブラーニング方式で授業を受けたA組のほうが有意に試験成績が良い」ということが言えます。

　ここでの有意確率 $p$ の具体的な意味は「**2つの標本がt検定の前提条件（無作為抽出・正規性・等分散性）を完璧に満たし、さらにそれら母集団の平均値の差がゼロであるとし、それら母集団から無作為に2群の標本を抽出して平均の差をとったとき、今回のような差、あるいは今回より大きな差が観測される確率**」を意味します。

　また、平均値の差の$95\%$信頼区間とは、$95\%$ の確率で母平均 $\mu$ を含んでいる範囲のことになります。

# おまけ

　平均80、標準偏差5の正規分布に従う40個の乱数（整数）を生成するコードです。

```Python:練習用データの生成（Python）
import numpy as np
import pandas as pd

a = (np.random.normal(80,5,size=40)).astype(np.int64).tolist()
print(a)

xa = pd.Series(a).hist()
```


# 参考資料

- [1] 小島寛之 著「完全独習 統計学入門」
- [2] 大村平 著「今日から使える統計解析 普及版」
- [3] 平井明代 編著「教育・心理系研究のためのデータ分析入門」
- [4] 向後千春 著「統計学がわかる ハンバーガーショップでむり無く学ぶ、やさしく楽しい統計学」