男女で心理学テストの平均値に差があるかどうか
統計の好き嫌いで統計テストの平均値に差があるか---統計が好きである人の統計テストの点数の分布と、統計が嫌いな人の統計テストの点数の分布の比較
この内容は山田、杉澤、村井(2008)「R」によるやさしい統計学を参考にしています。
この講義では、「２つの平均値を比較する」をとりあげます。 具体的には、２つの平均値を比較する方法の学習、つまり、 (1)独立な2群の平均値差の検定と、(2)対応のある2群の平均値差の検定、 について学びます。
学習項目です:
物語: ある心理学の教授は考えました。「今までの講義では教育効果があがらない。 そこで指導方法を変えてみよう（例えば、アクティブラーニングを取り入れる、 毎回宿題を出して学習を促進する）」。
一方、試験では今までと同程度の問題を出して、指導方法を変えたことによる効果があったか、検証することにしました。
本当に効果があったか調べるには、今までの指導法による試験のデータと、
指導方法を変えた後の試験のデータ(「指導法データ」と呼びます)
を比較すればよいのだが、
具体的にはどうすればよいのだろう???
2つの平均値を比較することが必要とされる場合:
2群が独立である具体的な例: 男女それぞれの心理学テストの成績、統計が好きな人と嫌いな人それぞれの統計学テストの成績、というように2群がそれぞれ別々の標本から得られたデータが対象--- t検定を実施するための条件については、あとで述べるt検定の前提条件を参照のこと
平均がそれぞれ$\mu_1, \mu_2$、かつ分散$\sigma^2$が等しい、正規分布に従う母集団から無作為抽出された標本$X_1, X_2$を考える：式で書くと $X_1 \sim N(\mu_1, \sigma^2)$ ,$X_2 \sim N(\mu_2, \sigma^2)$
その平均値の差(ここで$n_1$と$n_2$はそれぞれの群のサンプルサイズとする)$\bar{X_1}-\bar{X_2}$
はやはり正規分布に従う：
これを次のように標準化する: $$ \frac{(\bar{X_1}-\bar{X_2}) -(\mu_1 - \mu_2)}{\sigma \sqrt{\frac{1}{n_1}+\frac{1}{n_2}} }$$ これは$N(0,1)$に従う(なぜだろうか？)。
ここで母分散$\sigma^2$が未知であるなら標本から推定する:(ここで、$\hat{\sigma}_1^2$と$\hat{\sigma}_2^2$はそれぞれの群の不偏分散) $$ \hat{\sigma}^2_{pooled} = \frac{(n_1-1)\hat{\sigma}_1^2+(n_2-1)\hat{\sigma}^2_2}{n_1+n_2 - 2} $$
この値によって次の検定統計量を導く(注：これは$\mu_1=\mu_2$という帰無仮説のもとでの検定統計量だから$\mu_1-\mu_2(=0)$が省略されている)：
$$t = \frac{\bar{X_1}-\bar{X_2}}{\sqrt{\frac{(n_1-1)\hat{\sigma}_1^2+(n_2-1)\hat{\sigma}^2_2}{n_1+n_2 - 2}(\frac{1}{n_1}+\frac{1}{n_2})}}$$
こうして導かれた検定統計量$t$の標本分布は、帰無仮説$H_0：\mu_1 = \mu_2$のもとで、
自由度$ n_1 + n_2 - 2$の$t$分布にしたがう(この理由はCramer(1946)参照)。この検定統計量を用いて、２つの平均値の差に関する検定を行うことができるようになる(問題: 対立仮説はどうなる？):
例題:「統計テスト１」(StatTest1)の得点の平均値に男女で有意な差があるかどうかを有意水準5%、両側検定で検定せよ
import numpy as np
StatTest1Male = np.array([6,10,6,10,5,3,5,9,3,3])
StatTest1FeMale = np.array([11,6,11,9,7,5,8,7,7,9])次のステップで行う:
print(np.mean(StatTest1Male))
⇒ 6.0
print(np.mean(StatTest1FeMale))
⇒ 8.0
print(np.var(StatTest1Male, ddof=1)
⇒ 7.777778
print(np.var(StatTest1FeMale, ddof = 1))
⇒ 4.0次に「プール標準偏差」(PooledStd)$\sigma_{pooled}$を求める。
n1 = len(StatTest1Male)
n2 = len(StatTest1FeMale)
PooledStd = (((n1-1)*np.var(StatTest1Male,ddof=1)+(n2-1)*np.var(StatTest1FeMale,ddof=1)) / (n1+n2-2))**0.5
print(PooledStd)
⇒ 2.426703最後に検定統計量の実現値を計算する:
tDivider =PooledStd * (1.0/n1+1.0/n2)**0.5
tDividend = np.mean(StatTest1Male)-np.mean(StatTest1Female)
t = tDividend / tDivider # 検定統計量
⇒ -1.842885検定統計量の実現値がt = -1.84と求まった
import scipy.stats as st
st.t.ppf(0.025,18) # 検定統計量 df=18、下側確率0.05/2 = 0.025となるtの値を求める
⇒ -2.100922
# 下側確率であるから、この値よりもt値が小さければ棄却される
st.t.ppf(0.975,18) # df=18、上側確率 1 - 0.05/2 = 0.975となるtの値を求める
⇒ 2.100922
# 上側確率であるから、この値よりもt値が大きければ棄却されるこの結果、棄却域は $t < -2.10$ または $t > 2.10$ となるので、tの値は棄却域に入らないため、帰無仮説は棄却されない。ゆえに、検定の結果は「有意水準5%で有意差が見られなかった」となる
なおcdf関数 を用いて、直接p値を求めることもできる:
st.t.cdf(-1.842885,18) # 下側確率
⇒ 0.04093903
# 下側確率とすれば、p値は0.04という値( < 0.05)
2*st.t.cdf(-1.842885,18) # 両側検定なので2倍する
⇒ 0.08187807
# 両側検定であるから２倍したp値は0.08という値(> 0.05)なお、scipy.statsモジュールのttest_ind関数によって簡単に求めることができる:(注意: 分散が等しい場合のみ。そうでなければ、後述するようにオプションにequal_var=Falseを指定する)
st.ttest_ind(StatTest1Make,StatTest1Female)
⇒ Ttest_indResult(statistic=-1.8428853505018534, pvalue=0.081878014624106793) # 検定統計量とp値# 以上の実行
import numpy as np
StatTest1Male = np.array([6,10,6,10,5,3,5,9,3,3])
StatTest1Female = np.array([11,6,11,9,7,5,8,7,7,9])
print(np.mean(StatTest1Male))
print(np.mean(StatTest1Female))
print(np.var(StatTest1Male, ddof=1))
print(np.var(StatTest1Female, ddof = 1))
n1 = len(StatTest1Male)
n2 = len(StatTest1Female)
PooledStd = (((n1-1)*np.var(StatTest1Male,ddof=1)+(n2-1)*np.var(StatTest1Female,ddof=1)) / (n1+n2-2))**0.5
print(PooledStd)
tDivider =PooledStd*(1.0/n1+1.0/n2)**0.5
tDividend = np.mean(StatTest1Male)-np.mean(StatTest1Female)
t = tDividend / tDivider # 検定統計量
print(t)
import scipy.stats as st
print(st.t.ppf(0.025,18) ) # 検定統計量 df=18、下側確率0.05/2 = 0.025となるtの値を求める
# 下側確率であるから、この値よりもt値が小さければ棄却される
print(st.t.ppf(0.975,18)) # df=18、上側確率 1 - 0.05/2 = 0.975となるtの値を求める
# 上側確率であるから、この値よりもt値が大きければ棄却される
print(st.t.cdf(-1.842885,18)) # 下側確率
# 下側確率とすれば、p値は0.04という値( < 0.05)
print(2*st.t.cdf(-1.842885,18)) # 両側検定なので2倍する
# 両側検定であるから２倍したp値は0.08という値(> 0.05)
st.ttest_indを使って一気に実行する場合（注意：これもst.ttest_1samp関数と同じく標本データがないと使えない)
st.ttest_ind(StatTest1Male,StatTest1Female)
重要
結論を書くこと。この例の場合、$t= -1.84$の値は棄却域に入らないため、帰無仮説は棄却されない($p=0.08 > 0.05$)。ゆえに、検定の結果は「有意水準5%で有意差が見られなかった」
独立な２群のｔ検定では、「それぞれ$N(\mu_1, \sigma^2)$と$N(\mu_2, \sigma^2)$に従う２つの標本$X_1$と$X_2$に対し、その平均値の差が$N(\mu_1-\mu_2, \sigma^2(\frac{1}{n_1}+\frac{1}{n_2}))$に従う」ということが前提となっていた（ただしここで、$n_1, n_2$はそれぞれの標本のサイズ)。このことを、今まで学んだことから示せ。
[ヒント] 正規分布に従う標本平均の標本分布がどのような分布に従うかは前項の「標本分布を求める」節で述べられている。また、それぞれ正規分布にしたがう２つの標本データの和(差)がどのような分布に従うかは「Pythonで確率を学ぶ」の主要な確率分布の項で述べられている。 これらを組み合わせえて考えよ。
解答欄
t検定を実行するには次の３つの条件が必要
2群の分散が等しいことを「この2群の分散は等質である」といいます。t検定では2群の分散が等質であることが前提です。それには「分散の等質性の検定」を行う必要があります。これは検定統計量をFとしているためF検定とも呼ばれています。
まずそれぞれの分散($\sigma_1^2, \sigma_2^2$とする)を求め、その比($F=\frac{\sigma_1^2}{\sigma_2^2}$)を検定統計量とします。ここで、分散の大きい方を分母にとる必要があります(今の例では$\sigma_2 > \sigma_1$と仮定した）。この$F$値は、2群のサイズをそれぞれ$n_1, n_2$とすると、分子の自由度$df_1＝n_1-1$、分母の自由度$df_2＝n_2-1$のF分布に従うので、scipy.stats.fモジュールのcdf関数によってp値を求めます(このとき、帰無仮説は「等質である」、対立仮説は「等質ではない」となることに注意)。
以下では例を用いて示します:
import scipy.stats as st
import numpy as np
ClassA = np.array([54,55,52,48,50,38,41,40,53,52])
ClassB = np.array([67,63,50,60,61,69,43,58,36,29])
print("variance of ClassA = %f, ClassB = %f" % (np.var(ClassA), np.var(ClassB)))
# ClassBの分散が大きいので、それを分母にとる
F = np.var(ClassA)/np.var(ClassB)
print(st.f.cdf(F, len(ClassA)-1, len(ClassB)-1)*2) # 両側検定のため2倍する
この例では、p値が 0.03であるため、帰無仮説が棄却されます。 つまりこの例では2つの変数の分散が等質であるという仮定は成り立ちません。そのため、この2つの変数に対してはt検定ができない、ということになります。その場合でも2つの群の平均値の比較の検定は可能です。その一つの方法が、次で紹介するWelch(ウェルチ)の検定です。
[なぜ分散が等質でないと判断されたか] 今の説明で「p値が 0.03であるため、帰無仮説が棄却された(つまりこの例では2つの変数の分散が等質であるという仮定は成り立たない)」という点に疑問を持った人のために、解説します。 ここで行ったのは「2つの群が等質である」ことを示すための検定ではないのです。行ったのは「等質でない」(これが『対立仮説』、帰無仮説は「等質である」）ことの検定でした。 そして有意水準を5%とすると、$p<0.05$ですから、帰無仮説が棄却され、対立仮説が採択されたのです。つまり、 「この例では2つの変数の分散が等質であるという仮定は成り立たない」という結論になります。
?st.f.cdf
母分散が等質でないときは、t検定は使えないのでWelch検定を使う---st.ttest_ind関数で、equal_var=Falseというオプションをつけて実行します:
st.ttest_ind(ClassA, ClassB, equal_var=False)
この結果はp値が0.28なので、帰無仮説が棄却されない、つまりクラスAとクラスBの平均に有意な差はない、という結論になります。
差得点$D \sim N(\mu_D, \sigma^2_D)$ と仮定ｓ，データ数を$n$とすると、その標本平均$\bar{D} \sim N(\mu_D, \frac{\sigma^2_D}{n})$となります。
この標本分布を標準化すると $$Z = \frac{\bar{D}-\mu_D}{\sigma_D / \sqrt{n}} \sim N(0,1)$$
ここで、$\sigma_D$が未知なので、これを標本から求めた標準偏差$\hat{\sigma}_D$で代用すると、 $t = \frac{\bar{D}-\mu_D}{\hat{\sigma}_D / \sqrt{n}}$ は自由度 $df = n - 1$ のt分布にしたがいます。
例題:指導の前後に行った「統計テスト１」(StatTest1)と「統計テスト２」(StatTest2)で有意な差があるかどうかを 有意水準5%、両側検定で検定せよ
import numpy as np
StatTest1 = np.array([6,10,6,10,5,3,5,9,3,3,11,6,11,9,7,5,8,7,7,9])
StatTest2 = np.array([10,13,8,15,8,6,9,10,7,3,18,14,18,11,12,5,7,12,7,7])次のステップで行う:
Difference = StatTest2 - StatTest1
print(np.std(Difference, ddof=1))
⇒ 2.772041
tDivider = np.std(Difference, ddof=1)/(len(Difference))**0.5
tDividend = np.mean(Difference)
t = tDividend / tDivider
print(t)
⇒ 4.839903この結果、検定統計量の実現値は t= 4.84
st.t.ppf(0.025,19) # 検定統計量 df=19、下側確率0.05/2 = 0.025となるtの値を求める
⇒ -2.093024
# 下側確率であるから、この値よりもt値が小さければ棄却される
st.t.ppf(0.975,19) # df=19、上側確率 1 - 0.05/2 = 0.975となるtの値を求める
⇒ 2.093024
# 上側確率であるから、この値よりもt値が大きければ棄却されるこの結果、棄却域は $t < -2.09$ または $t> 2.09$ となるので、tの値は棄却域に入るため、帰無仮説は棄却される。ゆえに、検定の結果は「5%水準で有意差が見られた」となる
なお、scipy.statsモジュールのttest_rel関数によっても対応のあるt検定を求めることができる:
st.ttest_rel(StatTest1, StatTest2)
⇒ Ttest_relResult(statistic=-4.8399026368560945, pvalue=0.00011378795008861286) # 検定統計量とp値この結果、p値が0.0001という小さな値であるので、帰無仮説が棄却され、検定の結果は「有意水準5%で、統計テスト1と統計テスト2では有意な差があった」となる。
import scipy.stats as st
import numpy as np
StatTest1 = np.array([6,10,6,10,5,3,5,9,3,3,11,6,11,9,7,5,8,7,7,9])
StatTest2 = np.array([10,13,8,15,8,6,9,10,7,3,18,14,18,11,12,5,7,12,7,7])
Difference = StatTest2 - StatTest1
print(np.std(Difference, ddof=1))
tDivider = np.std(Difference, ddof=1)/(len(Difference))**0.5
tDividend = np.mean(Difference)
t = tDividend / tDivider
print(t)
print(st.t.ppf(0.025,19)) # 検定統計量 df=19、下側確率0.05/2 = 0.025となるtの値を求める
# 下側確率であるから、この値よりもt値が小さければ棄却される
print(st.t.ppf(0.975,19) ) # df=19、上側確率 1 - 0.05/2 = 0.975となるtの値を求める
# 上側確率であるから、この値よりもt値が大きければ棄却される
標本データがある場合、一つの関数を用いて実行することもできます。ただ計算値を表示するだけではだめで、結論を「言葉で」書くことを忘れないようにしよう
st.ttest_rel(StatTest1, StatTest2)
今の例題のデータに対して「対応なしの場合の検定」を行え。そして対応ありの場合と対応なしの場合では、どのような違いがあるか、具体的に述べよ。
解答欄
注: numpyをnp, numpy.randomを random、 matplotlib.pyplotをplt、pandasをpd、scipy.statsをstと略記する
| 目 的 | 関数名とモジュール | 使い方 |
|---|---|---|
| 分散の等質性の検定(F検定)におけるp値 | st.f.cdf(F,分子の自由度,分母の自由度) | st.f.cdf(np.var(ClassA)/np.var(ClassB),len(ClassA)-1,len(ClassB)-1) |
| # ClassAとClassBの分散の等質性の検定のためのp値(np.var(B) > np.var(A)と仮定) | ||
| 独立な２群の検定 | st.ttest_ind(データ1,データ2) | st.ttest_ind(StatTest1Male,StatTest1Female) |
| Welchの検定(分散が等質でない場合) | st.ttest_ind(データ1,データ2, equal_var=False) | st.ttest_ind(ClassA, ClassB, equal_var=False) |
| 対応のあるt検定 | st.ttest_rel(データ1,データ2) | st.ttest_rel(StatTest1, StatTest2) |
以下は、統計学が好きな人の得点と嫌いな人の得点のデータである。統計学が好きか嫌いかという2群について、平均値に有意な差があるかどうかを、有意水準5%の両側検定で調べよ。
StatLovers = np.array([6, 10, 6, 10, 11, 6, 11, 7])
StatHaters = np.array([5, 3, 5, 9, 3, 3, 9, 5, 8, 7, 7, 9])[ヒント]「統計学好き」(StatLovers)と「統計学嫌い」(StatHaters) の2変数は、分散が等質かどうか不明である。そこで、まず(1)分散の等質性の検定を行い、その結果(2a)分散が等質であればt検定を行い、(2b)分散が等質でなければWelchの検定を行う。
import numpy as np
import scipy.stats as st
StatLovers = np.array([6, 10, 6, 10, 11, 6, 11, 7])
StatHaters = np.array([5, 3, 5, 9, 3, 3, 9, 5, 8, 7, 7, 9])
以下は、心理学テストにおける、男と女のデータである。男女で平均値に有意な差があるかどうかを、有意水準5%の両側検定で調べよ。
PsycTestMale= np.array([13, 14, 7, 12, 10, 6, 8, 15, 4, 14])
PsycTestFemale = np.array([9, 6, 10, 12, 5, 12, 8, 8, 12, 15])[ヒント] 「心理学テスト男」(PsycTestMale)と「心理学テスト女」(PsycTestFemale) の2変数は、分散が等質かどうか不明である。そこで、まず(1)分散の等質性の検定を行い、その結果(2a)分散が等質であればt検定を行い、(2b)分散が等質でなければWelchの検定を行う。
# 計算用
import numpy as np
import scipy.stats as st
PsycTestMale= np.array([13, 14, 7, 12, 10, 6, 8, 15, 4, 14])
PsycTestFemale = np.array([9, 6, 10, 12, 5, 12, 8, 8, 12, 15])
以下はビデオ視聴によるダイエットプログラムに参加した10名の、プログラム参加前後の体重のデータである(単位はkg)。 このデータから、ダイエットプログラムは効果があったと言えるかどうかを判定せよ。有意水準を5%、両側検定とする。
Before = np.array([61, 50, 41, 55, 51, 48, 46, 55, 65, 70])
After = np.array([59, 48, 33, 54, 47, 52, 38, 50, 64, 63])[ヒント] 同じ人についてのプログラム参加前(Before)と後(After)のデータなので、「対応のあるt検定」を行う。
# 計算用
import numpy as np
import scipy.stats as st
Before = np.array([61, 50, 41, 55, 51, 48, 46, 55, 65, 70])
After = np.array([59, 48, 33, 54, 47, 52, 38, 50, 64, 63])