第六回（11月20日） Task Check and Weekly Assignment
To Do
□ t検定をやってみる
□ t検定の検定力をチェックしておく
□ 分散分析をやってみる(aov関数)
□ 分散分析をやってみる(anovakun関数)
Assignment
□ 国語の点数に男女差があるといってよいか，t検定してレポートしなさい
□ 算数の点数にクラス差があるといってよいか。分散分析してレポートしなさい
□ 理科の点数に性差とクラス差があるといってよいか。分散分析してレポートしなさい。
いつものファイルの読み込みと下処理。済んでいる場合はやらなくてよいです。
※Workspaceにsampleが入っている場合＝済んでいる場合です。
sample <- read.csv("sample(mac).csv", na.strings = "*")
sample$sex <- factor(sample$sex, labels = c("male", "female"))
t-testの例。まずはデータを二分割（別にしなくても良いけど，わかりやすくするため）。
height.male <- subset(sample$height, sample$sex == "male")
height.female <- subset(sample$height, sample$sex == "female")
正規性の検定。コルモゴロフスミルノフ検定。
# ks.test
ks.test(height.male, mean(height.male), sd(height.male))
##
## Two-sample Kolmogorov-Smirnov test
##
## data: height.male and mean(height.male)
## D = 0.56, p-value = 0.902
## alternative hypothesis: two-sided
ks.test(height.female, mean(height.female), sd(height.female))
##
## Two-sample Kolmogorov-Smirnov test
##
## data: height.female and mean(height.female)
## D = 0.52, p-value = 0.9804
## alternative hypothesis: two-sided
# each data set over 0.05 -> H0:this data seems to be normal curve
var.test(height.male, height.female)
##
## F test to compare two variances
##
## data: height.male and height.female
## F = 2.361, num df = 49, denom df = 49, p-value = 0.003206
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 1.34 4.16
## sample estimates:
## ratio of variances
## 2.361
問題なければt検定に進める。
result.ttest <- t.test(height ~ sex, data = sample, var.equal = F)
result.ttest
##
## Welch Two Sample t-test
##
## data: height by sex
## t = 10.58, df = 84.2, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 10.69 15.63
## sample estimates:
## mean in group male mean in group female
## 158.0 144.8
検定力分析。pwr関数を使う。
ちなみにt検定において検定力が大きいってどれぐらいの数字？
library(pwr)
cohen.ES(test = "t", size = "large")
##
## Conventional effect size from Cohen (1982)
##
## test = t
## size = large
## effect.size = 0.8
なるほど，では今回はどうだったか。
効果量esを計算してから関数に代入する。
es <- result.ttest$statistic * sqrt((50 + 50)/(50 * 50))
power.t.test(n = 50, d = es)
##
## Two-sample t test power calculation
##
## n = 50
## delta = 2.116
## sd = 1
## sig.level = 0.05
## power = 1
## alternative = two.sided
##
## NOTE: n is number in *each* group
次に分散分析。aov関数。
result.aov <- aov(height ~ class, data = sample)
summary(result.aov)
## Df Sum Sq Mean Sq F value Pr(>F)
## class 2 116 58.2 0.71 0.5
## Residuals 97 8000 82.5
複数の要因が入っている場合は独立変数を*で追加
result.aov2 <- aov(height ~ class * sex, data = sample)
summary(result.aov2)
## Df Sum Sq Mean Sq F value Pr(>F)
## class 2 116 58 1.51 0.23
## sex 1 4354 4354 112.75 <2e-16 ***
## class:sex 2 15 8 0.20 0.82
## Residuals 94 3630 39
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anovakunというと便利な関数でやってみよう。
まずはanovakunソースを読み込む
source("anovakun_431.txt")
画面上は何の変化もないが，workspaceのfuunctionsになにか追加されているはず。
実行にあたっては，まず分析専用のファイルを作るところから始める。
たとえば，性別とクラスを独立変数，国語の点数を従属変数とした分析をする。
sampleK <- subset(sample, select = c("sex", "class", "kokugo"))
anovakun(sampleK, "ABs", 2, 3, peta = TRUE)
##
## [ ABs-Type Design ]
##
## This output was generated via anovakun 4.3.1 at R version 2.15.1.
## It was executed on Tue Nov 20 14:26:14 2012.
##
##
## << DESCRIPTIVE STATISTICS >>
##
## ----------------------------------------------------
## A B N Mean S.D.
## ----------------------------------------------------
## a1 b2 17 68.5294 9.4281
## a1 b1 17 64.2941 14.5635
## a1 b3 15 62.5333 11.8072
## a2 b2 16 66.8750 13.0735
## a2 b1 17 63.9412 14.0822
## a2 b3 17 60.6471 14.5771
## ----------------------------------------------------
##
##
## << ANOVA TABLE >>
##
## == This data is UNBALANCED!! ==
## == Type III SS is applied. ==
## == The number of removed case is 1. ==
##
## ----------------------------------------------------------------------------------------------------
## Source SS df MS F-ratio p-value p.eta^2
## ----------------------------------------------------------------------------------------------------
## A 41.5966 1 41.5966 0.2432 0.6230 ns 0.0026
## B 612.9299 2 306.4650 1.7921 0.1723 ns 0.0371
## AxB 11.3897 2 5.6948 0.0333 0.9673 ns 0.0007
## Error 15904.0716 93 171.0115
## ----------------------------------------------------------------------------------------------------
## Total 16586.7273 98 +p < .10, *p < .05, **p < .01, ***p < .001
##
##
## output is over --------------------///
今回の例では下位検定まで行われないが，効果があれば自動的に算出してくれます。
最後に分散分析の検定力分析。例によって，どの程度の数字があれば普通の検定力があったのかを調べてから。
効果量f2はpetaから算出する。
cohen.ES(test = "f2", size = "medium")
##
## Conventional effect size from Cohen (1982)
##
## test = f2
## size = medium
## effect.size = 0.15
pwr.f2.test(u = 1, v = 93, f2 = 0.0026/(1 - 0.0026))
##
## Multiple regression power calculation
##
## u = 1
## v = 93
## f2 = 0.002607
## sig.level = 0.05
## power = 0.07822