第四回（10月22日） Task Check and Weekly Assignment
To Do □ 統計的仮説検定の手続きについて理解する。 □ χ二乗検定をやってみる。 □ t検定をやってみる □ t検定の検定力をチェックしておく □ 分散分析をやってみる(aov関数) □ 分散分析をやってみる(anovakun関数)
Assignment □算数のテストを四段階評価（優，良，可，不可）にわけたときクラス別の評価に有意な差があるといってよいか検定し，答えなさい。 □ 国語の点数に男女差があるといってよいか，検定し答えなさい。 □ 算数の点数にクラス差があるといってよいか，検定し答えなさい。 □ 理科の点数に性差とクラス差があるといってよいか，検定し答えなさい。
以下のような分割表があったとする
data <- matrix(c(20, 30, 10, 40, 30, 20), byrow = T, nrow = 3)
data
## [,1] [,2]
## [1,] 20 30
## [2,] 10 40
## [3,] 30 20
χ二乗検定をするにはchisq.test関数
chisq.test(data)
##
## Pearson's Chi-squared test
##
## data: data
## X-squared = 16.67, df = 2, p-value = 0.0002404
ファイルから読み込む場合。いつもの下処理をしてから。
sample <- read.csv("sample_win_.csv", na.strings = "*")
sample$sex <- factor(sample$sex, labels = c("male", "female"))
国語の点数を，優良可不可に区分してみる。
sample$kokugo.class <- ifelse(sample$kokugo > 80, 1, ifelse(sample$kokugo >
70, 2, ifelse(sample$kokugo > 60, 3, 4)))
sample$kokugo.class <- factor(sample$kokugo.class, labels = c("A", "B", "C",
"D"))
集計と検定，一通り。
tabs <- xtabs(~kokugo.class + sex, data = sample)
chisq.test(tabs)
##
## Pearson's Chi-squared test
##
## data: tabs
## X-squared = 2.345, df = 3, p-value = 0.5041
res.chisq <- chisq.test(tabs)
残差の分析。1.96が目安。
res.chisq$stdres
## sex
## kokugo.class male female
## A -0.5786 0.5786
## B 1.5043 -1.5043
## C -0.3968 0.3968
## D -0.5361 0.5361
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_433.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.3 at R version 3.0.2.
## It was executed on Tue Oct 22 11:48:26 2013.
##
##
## << 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