お勧め本！
クイック・データアナリシス―
10 秒でできる実践データ解析法
田中 敏 (著), 中野 博幸 (著) 出版社: 新曜社 ; ISBN: 4788509199 ; (October 2004)
この本いいですね～。絶対おすすめします！
比率の検定
１ページの「今回の授業は今までの授業よりよかったですか」の結果を分析してみましょう。
生徒数が２０名で、ハイが１５名、イイエ５名です。
二項検定
binom .test
これは、正確な比率の検定です。
> binom.test(15,20)
Exact binomial test data: 15 and 20
number of successes = 15, number of trials = 20, p-value = 0.04139 alternative hypothesis: true probability of success is not equal to 0.5 95 percent confidence interval:
0.5089541 0.9134285 sample estimates:
probability of success 0.75
両側検定（デフォールト）で、5%水準で有意に「ハイ」が多いですね。
もし片側だと、その半分の水準になります。
> binom.test(15, 20, alternative="greater") Exact binomial test
data: 15 and 20
number of successes = 15, number of trials = 20, p-value = 0.02069 alternative hypothesis: true probability of success is greater than 0.5 95 percent confidence interval:
0.5444176 1.0000000 sample estimates:
probability of success 0.75
さて６人と１４人ではどうなるでしょうか？＞宿題 なお、オプションで母比率は変えることもできます。
P34の例を分析してみましょう。
女性 男性
４０代管理職の人数 １２人 １０人
英語教師の母比率 ７割 ３割
> binom.test(12,22, p=0.7,alternative="less") Exact binomial test
data: 12 and 22
number of successes = 12, number of trials = 22, p-value = 0.0916 alternative hypothesis: true probability of success is less than 0.7 95 percent confidence interval:
0.0000000 0.7286868 sample estimates:
probability of success 0.5454545 片側で、9.16%ですね。
この分析の実践は、上記の本にいろいろ参考になること（符号検定）が載っています。
加えて、説明のある自動集計検定は、とても便利そうです。
以下のリンクのJAVA-SCRIPT STAR を使ってみてください。
http://www.kisnet.or.jp/nappa/software/star/
~~おまけ
普通に行われているΧ二乗検定は、
chisq.test()
です。以下のリンクをみてください。
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/66.html
でも、近似値計算ではない、フィッシャーの検定をできるだけ行いましょう。
あと比率の検定には、
prop.test() というのもあります。
http://www.okada.jp.org/RWiki/index.php?cmd=read&page=Ķ%CC%F5%A1%A7%C8%E6Ψ%A4κ%B9%A4θ%A1%C 4%EA&word=prop.test
～～
クロス集計の際のフィッシャーの直接確率検定
次は、クロス集計表の分析です。
フィッシャーの直接確率検定 でやっていきましょう。
>fisher.test() は、引数にmatrix（行列）形式のデータをとります。
そこで、まずはmatrixを作成する方法を学びましょう。
1,2,3,4,5,6の数字を２行３列の行列にすることで学びます。
まず、ひとつのベクトルにします。これはc(1,2,3,4,5,6)です。それにmatrix関数を使います。
>matrix(c(1,2,3,4,5,6), nrow=2, ncol=3, byrow=T) 結果は以下のとおり。
[,1] [,2] [,3]
[1,] 1 2 3 [2,] 4 5 6
nrowは、number of row ということで行数を指定し、
ncolは、number of colum ということで列数を指定し、
byrow=T は、行ごとに、言い換えると、横ならびで入れていきます。省略すれば、縦ならびになります。（１から６ までの数字の並びを、1:6 とも指定できます。)
> matrix(c(1,2,3,4,5,6), nrow=2, ncol=3) [,1] [,2] [,3]
[1,] 1 3 5 [2,] 2 4 6 あるいは、
> matrix(1:6, nrow=2, ncol=3) [,1] [,2] [,3]
[1,] 1 3 5 [2,] 2 4 6
でも同じです。見やすいのは、byrow=Tのオプションをつけたほうでしょう。
Nrowとncolはどちらかを省略できます。
> matrix(1:6, nrow=2, byrow=T) [,1] [,2] [,3]
[1,] 1 2 3 [2,] 4 5 6
> matrix(1:6, ncol=3, byrow=T) [,1] [,2] [,3]
[1,] 1 2 3 [2,] 4 5 6
>
どちらも同じ行列が作られました。
さて、フィッシャーの直接確立検定をします。
次は、私の指導生が卒業論文でとったデータです。対象の内的想起の４タイプと、携帯電話依存の４タイプの 行列です。
> x <- matrix(c(11,13,6,8,12,7,5,6,7,8,4,1,4,4,2,2),ncol=4,byrow=T)
> x
[,1] [,2] [,3] [,4]
[1,] 11 13 6 8 [2,] 12 7 5 6 [3,] 7 8 4 1 [4,] 4 4 2 2
列に内的対象想起のタイプ、行に携帯電話依存のタイプとなっています。
> fisher.test(x)
以下にエラーfisher.test(x) : FEXACT error 6.
LDKEY is too small for this problem.
Try increasing the size of the workspace.
>
エラーが出ました。計算するのに必要なメモリーを確保していなかったようです。
workspace=3000000 というように、メモリー量を指定して行ってみましょう。
> fisher.test(x,workspace=3000000) Fisher's Exact Test for Count Data data: x
p-value = 0.8859
alternative hypothesis: two.sided
今度は、計算されました・・・・が、関係は認められなかったですね。
マクネマーの有意変化の検定
「テレビゲームの暴力表現は青少年に悪影響を及ぼすか？」を１００名の男子高校生に質問をしたあと、ゲー ムで覚えた技を使いたくて男性を殺害した事件が発生した。その事件の前と後で、同じ１００人に、それぞれ質 問を繰り返した結果が以下の表となったそうです。同じ人に質問を繰り返したので、４つの回答パターンにわか れます。その人数です。
例題から分かる心理統計学 この本もお勧めです。
影響あり（事件前） 影響なし（事件前） 計
影響あり（事件後） 49 16 65
影響なし（事件後） 7 28 35
計 56 44 100
例題からわかる心理統計学 p46より
> x <- matrix(c(49, 16, 7, 28), ncol=2, byrow=T)
> x
[,1] [,2]
[1,] 49 16 [2,] 7 28
> mcnemar.test(x, correct=F) McNemar's Chi-squared test data: x
McNemar's chi-squared = 3.5217, df = 1, p-value = 0.06057
>
6%ですので、殺人事件の報道をきいて意見が変わったとは言えないという結果ですね。
（correctのオプションの意味は自分で調べてください。)
カテゴリデータの検定を参考にしてください。
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/66.html
モザイク図を描きたいときには、下のページが参考になります。
http://www.okada.jp.org/RWiki/index.php?cmd=read&page=%A5%B0%A5%E9%A5 %A3%A5å%AF%A5%B9%BB%B2%ե B9ͼ%C2%CE
㽸 %A1%A7%A5⥶ %A5%A4%A5%AF%A5 ץ %ED%A5å%C8&word=%A5⥶ %A5%A4%A5%AF
心理学の実験や観察で、でてきる評点間一致について
κ (カッパ) 係数を計算する際は，パッケージ vcd を呼び出した後で関数 Kappa() を実行.
信頼性係数 ICC (Intraclass correlation coefficient) を求める場合は，パッケージ irr 中の関数 icc() を用いる．
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/66.html
＜おまけです＞
心理統計をする上での重要なポイント
http://www.okayama-u.ac.jp/user/le/psycho/member/hase/articles/1994/9407Hasegawa/9407Hasegawa.html 心理学研究法再考
(1)基礎的統計解析の誤用をなくすための３０のチェック項目
分散分析について
リンク
基本の説明があります。
http://www.shiga-med.ac.jp/~koyama/stat/com-anova.html
２要因以上になるとややこしくなる分散分析ですが、以下のリンクを参考にしてください。
http://www4.ocn.ne.jp/~murakou/anova.htm 上は、分散分析をするときの注意点です。
「心理学のためのデータ解析テクニカルブック」の例題を R で書いてみたというものが下にあります。
http://mat.isc.chubu.ac.jp/R/tech.html
雑誌に載ったテキストの分散分析の部分です。
http://www1.doshisha.ac.jp/~mjin/R/12.pdf 心理学者用のページです。
繰り返しのある場合の分散分析の例も載っています。
http://personality-project.org/r/r.guide.html#anova
特にそのなかの「R and Analysis of Variance」は以下のリンク。
http://personality-project.org/r/r.anova.html
#Statistical Computing Seminars Repeated Measures Analysis using SAS
#http://www.ats.ucla.edu/stat/sas/seminars/sas_repeatedmeasures/default.htm それを R で書いたコードが以下にあるます。
http://cat.zero.ad.jp/~zak52549/R/anova/repeated.txt
「R でラクラク分散分析」という研修テキストです。
http://cse.niaes.affrc.go.jp/minaka/R/Ranova-nlbc04.pdf
交互作用が出たときの理解のためのヒントです。
http://www.juen.ac.jp/psych/nakayama/anova.html Practical Statistics with R ・・・参考になるテキストです。
http://www.brandonu.ca/zoology/rutherford/RBookWebSite/RBook_pdf.html
以上のリンク先を参考にしましょう。
～～～～～～
さて、基本からはじめましょう。
一要因の分散分析
http://psy.isc.chubu.ac.jp/~oshiolab/teaching_folder/datakaiseki_folder/04_folder/da04_02.html にあるデータを使います。（被験者間データです。）
というデータです。ダブルクリックして、Calcで開いて（OpenOfficeを使ってる場合)、クリップボードにコピーしま しょう。
> x.df <- read.delim("clipboard")
> x.df
番号 条件 結果
1 1 4
2 1 1
3 1 3
4 1 2
5 1 2
6 1 4
7 1 3
8 2 6
9 2 8
10 2 5
11 2 9
12 2 8
13 2 7
14 2 7
15 3 4
16 3 3
17 3 4
18 3 6
19 3 5
20 3 5
21 3 5
番号 条件 結果 1 1 1 4 2 2 1 1 3 3 1 3 4 4 1 2 5 5 1 2 6 6 1 4 7 7 1 3 8 8 2 6 9 9 2 8 10 10 2 5 11 11 2 9 12 12 2 8 13 13 2 7 14 14 2 7 15 15 3 4 16 16 3 3 17 17 3 4 18 18 3 6 19 19 3 5 20 20 3 5 21 21 3 5
>
まず最初にやること、条件の変数を、因子タイプに変えること！
> x.df$条件 <- factor(x.df$条件)
> class(x.df$条件) [1] "factor"
分散分析は英語では、analysis of variance ですね。頭文字をとると、aov。その名前のオブジェクトを使います。
> result <- aov(結果~条件, data=x.df)
分析結果を、resultにいれておきました。どんな種類の結果がresultの中にあるかをnames()で調べます。
Aovの( )の中の式は、~ チルダを使います。
> names(result)
[1] "coefficients" "residuals" "effects" "rank"
[5] "fitted.values" "assign" "qr" "df.residual"
[9] "xlevels" "call" "terms" "model"
通常、必要な結果は、summary( )で、以下の出てきます。
> summary(result)
Df Sum Sq Mean Sq F value Pr(>F) 条件 2 69.238 34.619 25.964 4.961e-06 ***
Residuals 18 24.000 1.333
---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
繰り返し強調しますが、要因となるのは、
factor( ) で、因子タイプに変えてから分析してください。でないと、おか
しなことになります。
要因ごとの、基本統計量を describeとbyを使って確認しておきましょう。
> by(x.df$結果, x.df$条件, describe) INDICES: 1
V n mean sd median min max range skew se 1 1 7 2.71 1.11 3 1 4 3 -0.15 0.42
--- INDICES: 2
V n mean sd median min max range skew se 1 1 7 7.14 1.35 7 5 9 4 -0.22 0.51
--- INDICES: 3
V n mean sd median min max range skew se 1 1 7 4.57 0.98 5 3 6 3 -0.17 0.37
>
有意な結果になっていますので、Tukeyの下位検定をしてみます。分散分析の結果を入れたresultを
TukeyHSD( )の中にいれます。
> TukeyHSD(result)
Tukey multiple comparisons of means 95% family-wise confidence level Fit: aov(formula = 結果 ~ 条件, data = x.df)
$条件
diff lwr upr p adj
2-1 4.428571 2.8533421 6.0038007 0.0000032 3-1 1.857143 0.2819136 3.4323722 0.0196279 3-2 -2.571429 -4.1466579 -0.9961993 0.0015998
>
どの条件間においても差が認められますね。
フィッシャーのLSDを使いたいならば、
> pairwise.t.test(x.df$結果, x.df$条件, p.adjust.method=”none”) Bonferroniの手法ならば、
> pairwise.t.test(x.df$結果, x.df$条件, p.adjust.method=”bonferroni”) Holmの手法ならば
> pairwise.t.test(x.df$結果, x.df$条件, p.adjust.method=”holm”)
グラフにしたいのならば、まず条件ごとの平均値を、y にいれましょう。
> y <- by(x.df$結果, x.df$条件, mean)
> y
INDICES: 1 [1] 2.714286
--- INDICES: 2
[1] 7.142857
--- INDICES: 3
[1] 4.571429
次に、plot( )を使います。
> plot(y)
折れ線にしたければ、line の”l” をtypeで指定します。
> plot(y , type="l")
～～ちょっと細かい話～～
分散分析は、各群の分散が等質性を前提としています。その前提をチェックするには、
bartlett.test()
を用いることができます。
> bartlett.test(x.df$結果, x.df$条件)
Bartlett test of homogeneity of variances data: x.df$結果 and x.df$条件
Bartlett's K-squared = 0.5877, df = 2, p-value = 0.7454
上の結果は等分散ということですが、もし等分散でなかったならば、oneway( )で検定できます。オプションで、
var=F と明示するか、書かずに分析してみましょう。
なお、正規性の検定は、shapiro.test() で Shapiro-Wilk 検定 を, ks.test() で Kolmogorov-Smirnov 検定 を行え ます。分散分析は、正規性については頑強だといわれてます。
～～～～
さて、次は被験者内要因の場合をみていきましょう。
同じく小塩のデータを使います。
「要因の分散分析（被験者内計画）
•
６名の被験者にミュラー・リヤーの錯視実験を行った。羽の角度を 30 度，60 度，90 度，120 度とした４つの条件を設け，6 名の被験者は繰り返し 4 つの条件の試行を行った。実験で は，２つの図形の線分の長さが等しいと判断した際の，２つの線分の長さの差を測定した。
•
角度によって錯視量に差があるかどうか，差があるかとすればどこにあるかを検定したい。
」
という例です。
> x.df <- read.delim("clipboard")
> x.df
subject r30 r60 r90 r120 1 1 42 39 36 34 2 2 22 17 15 8 3 3 35 32 25 25 4 4 34 30 22 20 5 5 40 33 28 23 6 6 34 30 23 28
>
この形で R は分析しませんので、データの形を変えます。
> subject <- c(rep(1,4), rep(2,4), rep(3,4), rep(4,4), rep(5,4), rep(6,4))
> subject
[1] 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 5 5 5 5 6 6 6 6
>
> condition <- rep(c(1:4), 6)
> condition
[1] 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4
>
> sakusi <- c(42,39,36,34, + 22,17,15, 8, + 35,32,25,25, + 34,30,22,20, + 40,33,28,23, + 34,30,23,28)
> sakusi
[1] 42 39 36 34 22 17 15 8 35 32 25 25 34 30 22 20 40 33 28 23 34 30 23 28
>
> subject <- factor(subject)
> condition <- factor(condition) 要因は、factorタイプに変えます。
> x <- data.frame(subject, condition, sakusi)
> x
subject condition sakusi 1 1 1 42 2 1 2 39 3 1 3 36 4 1 4 34 5 2 1 22
6 2 2 17 7 2 3 15 8 2 4 8 9 3 1 35 10 3 2 32 11 3 3 25 12 3 4 25 13 4 1 34 14 4 2 30 15 4 3 22 16 4 4 20 17 5 1 40 18 5 2 33 19 5 3 28 20 5 4 23 21 6 1 34 22 6 2 30 23 6 3 23 24 6 4 28
>
被験者内要因が入るときには、+Error( / )を付け加えます。
～～～ちょっと詳しく～～
Mauchlyの球面性検定のやり方は謎？？なので、仮定を満たしたものとしてやってみます。
mauchly.test(stats) Mauchly's Test of Sphericity があり、SSD(stats)を使って前処理したあとで使うよ うですが・・ちょっと使い方を分かりかねています。
また、http://www.psych.upenn.edu/~baron/rpsych/rpsych.html#SECTION000710000000000000000 6.10 Sphericity にthe Greenhouse-Geisser and the Huynh-Feldt corrections のイプシロン算出の方法の記載があ ります。Greenhous-Geisser と Huygh-Feldt の修正は、自由度にイプシロンをかけるということです。それでＦ 値が決まり、pf( ）関数を使って、p値を求めます。
http://www.psych.upenn.edu/~baron/rpsych/rpsych.html にErrorの使い方を含めて分散分析のやり方が詳し くあります。
～～～
> result <- aov(sakusi~condition+Error(subject/condition))
> summary(result) Error: subject
Df Sum Sq Mean Sq F value Pr(>F) Residuals 5 1058.37 211.67 Error: subject:condition
Df Sum Sq Mean Sq F value Pr(>F) condition 3 491.46 163.82 32.855 7.766e-07 ***
Residuals 15 74.79 4.99
---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>