勉強 4日目
2つのグループの平均値に差があるかどうかを検定するために,t検定を行う。 社会科学では,2グループに対応関係があるケースはほとんど無く,また2グループの分散も通常等しくはない。 そこで,対応関係にない分散の異なる2グループの平均差の検定を行う際に用いるウェルチのt検定(Welch’s t-test)について説明する。
検定は次のステップで行う。
2グループの平均値に差は無い,という帰無仮説を考える。 ウェルチのt値は次のように定義される。
\[ t = \frac{\bar x_1 - \bar x_2}{\sqrt{\frac{s^2_1}{n_1} + \frac{s^2_2}{n_2}}} \]
平均の差を見るための作図パッケージgplots
を導入する。
#install.packages("gplots") # first time only
library(gplots)
ID,性別,結婚満足度尺度からなるデータセットを読み込む。
d6 <- read.csv("chap6.csv", header=T)
まず,結婚満足度の男女差を図示する。
plotmeans(d6$m_satis ~ d6$sex)
このままだと縦軸の範囲が狭く,グラフが見づらいので,縦軸などを少し調整する。 plotmeans
関数のオプションとして,ylim
でY軸の範囲を,connect
で平均をつなげる直線を引くかどうかを決める。
plotmeans(d6$m_satis ~ d6$sex, ylim=c(1,5), connect = F)
これより,平均にはあまり差が無いように見える。 次に等分散を仮定しないウェルチのt検定を行い,統計的に差を評価する。
res <- t.test(d6$m_satis~ d6$sex)
print(res)
Welch Two Sample t-test
data: d6$m_satis by d6$sex
t = 0.93704, df = 17.989, p-value = 0.3611
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.4968692 1.2968692
sample estimates:
mean in group 1 mean in group 2
4.4 4.0
分析の結果から,グループ1の平均は4.4,グループ2の平均は4.0であった。 t値は,0.94となり, このとき帰無仮説の下でこの差が出てくる確率は0.36となる。 よって帰無仮説の下でこの差がでてくることがあり得ないとはいえず,2群の平均値には差があるとはいえない。 したがって,結婚満足度の男性平均と女性平均に差があるかどうか分からない。
悩みの相談相手の人数が性別により異なるのかを調べる。 まず,ID,性別,相談相手の人数についてのデータセットを読み込む。
d6.hw <- read.csv("chap6.hw.csv", header=T)
まずはデータの特徴を掴むためにplotmeans()
を用いて図示する。
plotmeans(d6.hw$friend ~ d6.hw$sex, ylim=c(0,5), connect = F)
今度は女子の相談相手人数の平均は男性のを大きく上回っていることが図より明らかである。 さらにウェルチの平均の差の検定を行い,統計的にこの差が有意なものなのかどうか,を検証する。
res <- t.test(d6.hw$friend ~ d6.hw$sex)
print(res)
Welch Two Sample t-test
data: d6.hw$friend by d6.hw$sex
t = -2.3062, df = 16.439, p-value = 0.03443
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-2.8545203 -0.1232574
sample estimates:
mean in group 1 mean in group 2
1.400000 2.888889
t値は,-2.31となり, このとき帰無仮説の下でこの差が出てくる確率は0.03となる。 よって帰無仮説の下でこの差がでてくる確率は極めて低いことから,帰無仮説を棄却し,2群の平均値には差があるといえる。 したがって,悩みの相談相手の人数について,女性は男性よりも悩みを相談する相手が多い,といえる。
n1 <- length(d6.hw$friend[d6.hw$sex==1])
n2 <- length(d6.hw$friend[d6.hw$sex==2])
m1 <- mean(d6.hw$friend[d6.hw$sex==1])
m2 <- mean(d6.hw$friend[d6.hw$sex==2])
s1 <- var(d6.hw$friend[d6.hw$sex==1])
s2 <- var(d6.hw$friend[d6.hw$sex==2])
welch_t <- ( m1 - m2 ) / sqrt( (s1/n1) + (s2/n2) )
ウェルチのt値を計算すると,-2.31となる。 ウェルチのt検定で用いられる自由度は
\[ v \sim \frac{\left ( \frac{s_1^2}{n_1} + \frac{s_2^2}{n_2} \right )^2 }{\frac{s_1^4}{n_1^2 (n_1 - 1)} + \frac{s_2^4}{n_2^2 ( n_2 -1)}} \]
で計算される。
v <- ((s1/n1) + (s2/n2) )^2 /#
#---------------------------------------------
( (s1^2 / (n1^2*(n1-1))) + ( s2^2 / (n2^2*(n2-1))) )
自由度は,16.44となる。 t値が-2.31,自由度が16.44であるとき,t分布のもとで差が無いという結果が出る確率を表すp値は, 0.0341となり,有意水準5%のもとで帰無仮説が棄却される。
以上の計算を関数としてオブジェクトにしておく。
welch <- function(d1,d2){
n1 <- length(d1)
n2 <- length(d2)
m1 <- mean(d1)
m2 <- mean(d2)
s1 <- var(d1)
s2 <- var(d2)
w <- ( m1 - m2 ) / sqrt( (s1/n1) + (s2/n2) )
v <- ((s1/n1) + (s2/n2) )^2 / ( (s1^2 / (n1^2*(n1-1))) + ( s2^2 / (n2^2*(n2-1))) )
p <- dt(w,v)
if (p < 0.01) {
print("1%水準で有意!")
} else if (p < 0.05) {
print("5%水準で有意!")
} else if (p < 0.10) {
print("10%水準で有意!")
} else
print("なにもいえない。")
}
先ほどの結果を再現してみる。
d1 <- d6.hw$friend[d6.hw$sex==1]
d2 <- d6.hw$friend[d6.hw$sex==2]
welch(d1,d2)
[1] "5%水準で有意!"