 |
はじめるまえに
本教科書で用いられているデータは著者のwebサイト（http://www.oup.com/grafenhails/）よりダウンロードできると書かれていますが、英語が苦手な人のために少しだけ補足しておきましょう。
まずリンク先へと行くと、Resourcesという項目の下にDatasets and supplementsという項目があるのでその部分をクリックしてください。するとリンク先のページにExcel, Minitab, SAS, SPSSという項目がありますから、Excelの部分をクリックしてください。そうすると、英語の注意書きのようなものが何行か書かれていて、Datasets one sheetとDatasets on separate sheetsという2つの項目があります。
[和訳] ここでは本書で用いたデータセットをExcelのワークシート形式でダウンロードすることができますが、全てのデータが入力されているワークシートはどちらのページからでも入手できます。"one sheet"のほうは1つのシートに各章のデータが全て入力されているので、MNITABやSAS、SPSSを使って解析する人はこちらをダウンロードした方が良いでしょう。一方で"separate sheets"のほうは各章の例データがそれぞれ別のシートに入力されているので、例えば、データをコピー&ペーストして使用する場合はこちらをダウンロードしたほうが便利でしょう。
というわけなので、結局はどちらのページからダウンロードしても構いませんが、separate sheetsのページからダウンロードしたものの方が便利だと思います。
第1章 分散分析への招待
p11 「BOX 1.1 1つの説明変数をもつ分散分析」の解析例
p11の表1.2における収穫量のデータを変数YIELDに代入し、肥料のデータを変数FERTILに代入します。このとき肥料のデータはカテゴリカルデータとして代入することを忘れないようにしてください。なお、Excelワークシートに表1.2と同じように入力してから（というか、webサイトからダウンロードしたワークシートから）、関数excel.w()を使用して読み込んでも構いません。
もしR上で直接やるなら、
> YIELD <- c(6.27, 5.36, 6.39, 4.85, 5.99, 7.14, 5.08, 4.07, 4.35, 4.95, + 3.07, 3.29, 4.04, 4.19, 3.41, 3.75, 4.87, 3.94, 6.28, 3.15, 4.04, + 3.79, 4.56, 4.55, 4.55, 4.53, 3.53, 3.71, 7.00, 4.61) > FERTIL <- as.factor(rep(c(1,2,3), c(10,10,10)))
もしExcelワークシートから読み込むなら、
> dat <- excel.w(2) > YIELD <- c(dat[,2]) > FERTIL <- as.factor(dat[,1])
のようにしてください。これで2つの変数に必要なデータを代入することができました。実際にBOX 1.1と同じ分散分析表を出力させるためには、関数lm()とanova()を使用します。
> result <- lm(YIELD ~ FERTIL)
> anova(result)
Analysis of Variance Table
Response: YIELD
Df Sum Sq Mean Sq F value Pr(>F)
FERTIL 2 10.8227 5.4114 5.7024 0.008594 **
Residuals 27 25.6221 0.9490
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ちなみに、表1.4のような各群の平均値の信頼区間は、Rで次のようにすれば求めることができます。標準誤差とは標本平均値の標準偏差のことで、これは s / sqrt(n) でもとめられるのでした（sは標本標準偏差でnはサンプルサイズ）。また、自由度27のt分布における片側面積が0.025（両側だと0.05）となるパーセント点は関数qt()を用いれば計算できます。それから、abs()は絶対値を返す関数です。
#肥料1の信頼区間を計算する > SE1 <- 0.976 / sqrt(10) #肥料1の標準誤差 > t.crit <- qt(0.025, 27) #qt(0.975, 27)としても同じ > 5.445 - abs(t.crit) * SE1 [1] 4.811727 > 5.445 + abs(t.crit) * SE1 [1] 6.078273
p18〜19 練習問題
メロンの解析例
p19の「表1.7 YIELDMに対する1元配置分散分析表」を出力するためには次のようにコマンドします。
> Melons <- excel.w(2)
> YIELDM <- c(Melons[,1])
> VARIETY <- as.factor(Melons[,2])
> Melons.result <- lm(YIELDM ~ VARIETY)
> anova(Melons.result)
Analysis of Variance Table
Response: YIELDM
Df Sum Sq Mean Sq F value Pr(>F)
VARIETY 3 1115.28 371.76 23.798 1.735e-06 ***
Residuals 18 281.19 15.62
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
問題(4)の平均の標準偏差というのは、標準誤差のことであるから、各群（各品種）のデータの標準偏差を求めて、それをsqrt(n)で割ってあげればよい。標準誤差と信頼区間の求め方は既に述べましたから省略します。
注意） 信頼区間の下側限界と上側限界の値が異なるのは、私は各群の不偏標準偏差を用いたのに対して、教科書では標準偏差をsqrt(EMS)としているからです。
> N1 <- length(VARIETY[VARIETY==1]) #品種1のデータ数 > N2 <- length(VARIETY[VARIETY==2]) > N3 <- length(VARIETY[VARIETY==3]) > N4 <- length(VARIETY[VARIETY==4]) > all.N <- c(N1, N2, N3, N4) > all.N [1] 6 6 4 6 > M1 <- mean(myData[VARIETY==1,][1]) #品種1の平均値 > M2 <- mean(myData[VARIETY==2,][1]) > M3 <- mean(myData[VARIETY==3,][1]) > M4 <- mean(myData[VARIETY==4,][1]) > all.M <- c(M1, M2, M3, M4) > all.M YIELDM YIELDM YIELDM YIELDM 20.49000 37.40333 20.46250 29.89667 > SD1 <- sd(myData[VARIETY==1,][1]) #品種1の標準偏差 > SD2 <- sd(myData[VARIETY==2,][1]) > SD3 <- sd(myData[VARIETY==3,][1]) > SD4 <- sd(myData[VARIETY==4,][1]) > all.SD <- c(SD1, SD2, SD3, SD4) > all.SD YIELDM YIELDM YIELDM YIELDM 4.694422 3.950497 4.764692 2.229894 > t.crit <- qt(0.975, 18) #自由度18のt分布における片側0.025となるパーセント点 > t.crit [1] 2.100922 > all.M - all.SD / sqrt(all.N) #95%信頼区間の下側限界 YIELDM YIELDM YIELDM YIELDM 18.57351 35.79055 18.08015 28.98632 > all.M + all.SD / sqrt(all.N) #95%信頼区間の上側限界 YIELDM YIELDM YIELDM YIELDM 22.40649 39.01612 22.84485 30.80702
メロンの解析例の補足
p18の図1.15のようなものを描きたい場合は以下のようにコマンドすればよい。なお、x軸のデータにカテゴリカルデータ（ VARIETY <- as.factor(Melons[,2]) ）を指定すると、自動的に箱ひげ図が描かれる。
> plot(Melons[,2], Melons[,1], ylab="収穫量", xlab="メロンの品種", ylim=c(0,45), + main="メロンの収穫量", yaxt="n", xaxt="n", bty="l", pch=18) > axis(1, at=c(1,2,3,4)) > axis(2, at=c(seq(0,45,5)))
|
またp19の表1.7のように各群の基本統計量を求めたい場合は、データフレームを作成することによって簡単に計算することができる。以下のコマンド例において、もし品種2の平均を求めたければ、VARIETY==1の部分をVARIETY==2とすればよい。
#データフレームの作成 関数data.fram()を使う > myData <- data.frame(YIELDM = YIELDM, VARIETY = VARIETY) > myData YIELDM VARIETY 1 25.12 1 2 17.25 1 3 26.42 1 4 16.08 1 5 22.15 1 6 15.92 1 7 40.25 2 8 35.25 2 9 31.98 2 10 36.52 2 11 43.32 2 12 37.10 2 13 18.30 3 14 22.60 3 15 25.90 3 16 15.05 3 17 28.05 4 18 28.55 4 19 33.20 4 20 31.68 4 21 30.32 4 22 27.58 4 > myData[VARIETY==1,][1] #myDataにおいて、VARIETYが1であるデータのみを抽出 YIELDM 1 25.12 2 17.25 3 26.42 4 16.08 5 22.15 6 15.92 > mean(myData[VARIETY==1,][1]) #品種1の収穫量の平均を求める YIELDM 20.49 > mean(myData[VARIETY==2,][1]) #品種2の収穫量の平均を求める YIELDM 37.40333 > mean(myData[VARIETY==3,][1]) #品種3の収穫量の平均を求める YIELDM 20.4625 > mean(myData[VARIETY==4,][1]) #品種4の収穫量の平均を求める YIELDM 29.89667
さらに、各群のサンプルサイズ（データ数）を求めたい場合、方法はいくつかあるが、最も単純な方法として次のようなものがある。
> length(VARIETY[VARIETY==1]) #品種1のデータ数 [1] 6 > length(VARIETY[VARIETY==2]) #品種2のデータ数 [1] 6 > length(VARIETY[VARIETY==3]) #品種3のデータ数 [1] 4 > length(VARIETY[VARIETY==4]) #品種4のデータ数 [1] 6
雌雄異株の樹木の解析例
> dio.trees <- excel.w(3)
Read 150 items
> SEX <- as.factor(dio.trees[,1])
> SEX
[1] 2 1 2 1 2 1 2 2 1 1 2 1 2 2 1 1 2 2 2 1 2 2 1 1 2 1 1 2 1 2 2 1 2 1 2 2 1 2
[39] 2 2 2 2 2 1 2 2 1 2 1 2
Levels: 1 2
> FLOWERS <- c(dio.trees[,3])
> FLOWERS
[1] 234 110 512 133 1107 624 8 333 637 225 512 47 164 1 36
[16] 524 198 225 548 272 164 5 133 133 419 442 3 1482 333 1138
[31] 561 488 282 536 11 148 512 70 1107 941 89 181 8 173 1
[46] 1091 561 3 312 1398
> plot(SEX, FLOWERS) #p274の図15.1
> anova(lm(FLOWERS ~ SEX)) #同じくBOX 15.1
Analysis of Variance Table
Response: FLOWERS
Df Sum Sq Mean Sq F value Pr(>F)
SEX 1 171841 171841 1.1754 0.2837
Residuals 48 7017255 146193
