最終更新日 ： 2007.8.24
実験の計画‐簡潔に行おう
p77 BOX 5.1 完全に無作為化して計画した豆の収穫量の解析
p77最終行からp78にかけて「Box 5.1では、ブロックを無視して解析された。」とあるが、分かりやすくいえば、この解析における帰無仮説は「豆の種類（各群）によって収穫量（の平均値）に違いはない」という、いわゆる1要因の分散分析をしているわけである。
> YIELD
[1] 9.0 14.6 18.3 14.1 21.9 22.4 14.2 14.1 17.4 25.6 23.9 19.2 12.7 15.8 11.5
[16] 21.1 23.7 6.4 12.1 16.1 15.9 19.6 12.3 18.3
> BEAN
[1] 6 4 3 5 2 1 5 6 3 2 1 4 3 4 5 1 2 6 5 4 3 1 6 2
Levels: 1 2 3 4 5 6
> box5.1 <- glm(YIELD ~ BEAN)
> summary(aov(box5.1))
Df Sum Sq Mean Sq F value Pr(>F)
BEAN 5 444.43 88.89 14.586 8.579e-06 ***
Residuals 18 109.69 6.09
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p78 BOX 5.2 無作為化ブロック計画を用いた豆の収穫量の解析
これはいわゆる乱塊法を線形モデルとして回帰分析した結果である。MBLOCK（ブロック）という要因とMBEAN（豆種）という要因を考慮した分析で、この場合の帰無仮説は「ブロックによる収穫量の平均値に差はない」というものと「豆種によって収穫量の平均値に差はない」ということになる。なお、2要因の分散分析でも、この場合は繰り返し数が1回なので交互作用効果は検討できない（これについては7章で詳しく述べられている）。
> MYIELD
[1] 9.0 14.6 18.3 14.1 21.9 22.4 14.2 14.1 17.4 25.6 23.9 19.2 12.7 15.8 11.5
[16] 21.1 23.7 6.4 12.1 16.1 15.9 19.6
> MBLOCK
[1] 1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 3 3 4 4 4 4
Levels: 1 2 3 4
> MBEAN
[1] 6 4 3 5 2 1 5 6 3 2 1 4 3 4 5 1 2 6 5 4 3 1
Levels: 1 2 3 4 5 6
> box5.2 <- glm(MYIELD ~ MBLOCK + MBEAN)
> summary(aov(box5.2))
Df Sum Sq Mean Sq F value Pr(>F)
MBLOCK 3 49.29 16.43 6.415 0.00668 **
MBEAN 5 449.41 89.88 35.093 4.18e-07 ***
Residuals 13 33.30 2.56
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p80 BOX 5.3 ラテン方格法による解析
> COLUMN
[1] 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4
Levels: 1 2 3 4
> ROW
[1] 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4
Levels: 1 2 3 4
> TREATMT
[1] 1 4 3 2 2 1 4 3 4 3 2 1 3 2 1 4
Levels: 1 2 3 4
> SEEDS
[1] 120 88 131 126 138 119 87 134 93 107 98 107 128 90 87 89
> box5.3 <- glm(SEEDS ~ COLUMN + ROW + TREATMT)
> summary(aov(box5.3))
Df Sum Sq Mean Sq F value Pr(>F)
COLUMN 3 1332.25 444.08 8.7937 0.012910 *
ROW 3 1090.25 363.42 7.1964 0.020577 *
TREATMT 3 2650.25 883.42 17.4934 0.002270 **
Residuals 6 303.00 50.50
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p82 BOX 5.4 無作為化ブロック計画のANOVA表と係数表
YIELDに対する分散分析表はBOX 5.2で示したものと同じなので、ここでは係数表のみを提示する。
> summary(box5.2)
Call:
glm(formula = MYIELD ~ MBLOCK + MBEAN)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.87431 -1.11597 0.01597 0.81458 2.20903
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 21.7660 0.9891 22.007 1.13e-11 ***
MBLOCK2 2.3500 0.9240 2.543 0.024500 *
MBLOCK3 -1.5167 0.9240 -1.641 0.124664
MBLOCK4 -0.8972 1.0669 -0.841 0.415581
MBEAN2 1.6896 1.2440 1.358 0.197493
MBEAN3 -5.6750 1.1316 -5.015 0.000237 ***
MBEAN4 -5.3250 1.1316 -4.706 0.000411 ***
MBEAN5 -8.7750 1.1316 -7.754 3.14e-06 ***
MBEAN6 -12.2104 1.2440 -9.816 2.23e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 2.561245)
Null deviance: 532.000 on 21 degrees of freedom
Residual deviance: 33.296 on 13 degrees of freedom
AIC: 91.55
Number of Fisher Scoring iterations: 2
この結果を分かりやすく式でまとめると次のようになる。また、既に述べたことであるが、教科書ではBLOCKについては4番目の変数（ブロック4）を削除し、BEANについては6番目の変数（豆種6）を削除しているが、Rでは先頭の列（変数）を削除している（結果の表示は異なるが、実質的な意味は同じである）。
例えば、今回の解析結果であれば、ブロック1を基準としているので、ブロック2はブロック1よりも約2倍ほどの収穫量があるといえる。豆種についても豆種2は豆種1よりも1.7倍ほど収穫量が多いといえる。逆にマイナスの値をとっている場合は、基準となる変数（MBLOCK1およびMBEAN1）よりも収穫量が少ないということを意味している。
また、この式から条件ごとの収穫量を予測することも可能である。例えば、ブロック3と豆種5の組み合わせの場合においてはどの程度の収穫量が得られるか知りたいとする。その場合、該当する変数（つまりMBLOCK3とMBEAN5）に1を代入し、それ以外の変数にはすべて0を代入する。すなわち、
MYIELD = 21.76 - 1.51 * 1 - 8.77 * 1 = 11.48
と計算されるわけである。これを全ての組み合わせにおいて計算をして、収穫量の予測を立てれば、どのブロックでどの豆を育てればもっとも多い収穫が得られるかを見積もることが可能となるわけである。
ちなみに、係数表において、標準偏回帰係数のp値が高い場合（すなわち有意でない場合）はその係数は0であるとみなす。なぜなら、そこで立てている帰無仮説は「標準偏回帰係数は0である」というものであるので、この仮説を棄却できない場合は、その変数の係数は0であるということなのである。ただし、従来の検定のように有意水準を0.05、あるいはそれ以下に設定するほどの必要はないこともある（場合によってはp=0.20くらいでよいかもしれない）。
p86 BOX 5.5 無作為化ブロック計画において直行性を少し失う
特に補足的な説明しない。もし調整平方和と逐次平方和を比べたいのなら、次のようにすればよい。
> library(car) # carパッケージの呼び出し
> Anova(box5.2, test="F", type="III") # Type3の平方和を計算
Anova Table (Type III tests)
Response: MYIELD
SS Df F Pr(>F)
MBLOCK 49.69 3 6.4669 0.006482 **
MBEAN 449.41 5 35.0933 4.18e-07 ***
Residuals 33.30 13
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
また、既に述べたかもしれないが、carパッケージは最初から実装されているわけではないので、ダウンロードする必要がある。その場合は次のようにコマンドすればよい。これだけで自動的にダウンロード、インストールされるので、その後に上のようにライブラリからパッケージを呼び出せばよい。
> install.packages("car")
p88 練習問題
カーネーションの生長
(1)の問題はタイプ1の平方和とタイプ2の平方和を比べてみればよい。Rでタイプ2の平方和を出すためにはcarパッケージに含まれているAnova()を利用すればよい。なお、教科書ではタイプ3の平方和を採用しているが、私はタイプ2の平方和で良いと思うのでそのようにした。もし、タイプ3の平方和が出したければ、引数にtype="III"と指定すればよい。また、test="F"とするのも忘れないように！！
# データフレームの作成
> data5.6 <- data.frame(
+ SQBLOOMS = dat[,1], BED = as.factor(dat[,2]),
+ WATER = as.factor(dat[,3]), SHADE = as.factor(dat[,4]))
> data5.6
SQBLOOMS BED WATER SHADE
1 4.359 1 1 1
2 3.317 1 1 2
3 3.606 1 1 3
4 4.123 1 1 4
（途中省略）
> attach(data5.6) # 作成したデータフレームにパスを通す
> box5.6 <- glm(SQBLOOMS ~ BED + WATER + SHADE)
> summary(aov(box5.6))
Df Sum Sq Mean Sq F value Pr(>F)
BED 2 4.1323 2.0661 9.4570 0.0007277 ***
WATER 2 3.7153 1.8577 8.5029 0.0013016 **
SHADE 3 1.6465 0.5488 2.5120 0.0789451 .
Residuals 28 6.1173 0.2185
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> library(car) # パッケージcarの呼び出し
> Anova(box5.6, test="F") # タイプ2の平方和
Anova Table (Type II tests)
Response: SQBLOOMS
SS Df F Pr(>F)
BED 4.1323 2 9.4570 0.0007277 ***
WATER 3.7153 2 8.5029 0.0013016 **
SHADE 1.6465 3 2.5120 0.0789451 .
Residuals 6.1173 28
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> detach(data5.6) # パスの解除
(2)については教科書の解答通り、(3)はbox5.6とbox5.8のタイプ2の平方和を用いた分散分析表を見比べてみれば分かるはずである。box5.6においてはタイプ1とタイプ2とで結果は一致していた（つまり直行性が保たれているということ）が、box5.8においてはどうだろうか？
> dat <- excel.w(4)
Read 132 items
> data5.8 <- data.frame(
+ SQ2 = dat[,1], B2 = as.factor(dat[,2]),
+ W2 = as.factor(dat[,3]), S2 = as.factor(dat[,4]))
> data5.8
SQ2 B2 W2 S2
1 4.359 1 1 1
2 3.317 1 1 2
3 3.606 1 1 3
4 4.123 1 1 4
（途中省略）
> attach(data5.8)
> box5.8 <- glm(SQ2 ~ B2 + W2 + S2)
> summary(aov(box5.8))
Df Sum Sq Mean Sq F value Pr(>F)
B2 2 2.7626 1.3813 8.3790 0.001641 **
W2 2 5.0793 2.5397 15.4057 4.367e-05 ***
S2 3 0.8072 0.2691 1.6321 0.207156
Residuals 25 4.1213 0.1649
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> Anova(box5.8, test="F")
Anova Table (Type II tests)
Response: SQ2
SS Df F Pr(>F)
B2 2.6490 2 8.0346 0.002020 **
W2 4.6764 2 14.1836 7.644e-05 ***
S2 0.8072 3 1.6321 0.207156
Residuals 4.1213 25
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> detach(data5.8)
(4)は教科書とは違う図の描き方をしてみよう。以下のようにして描かれた図を見ても分かるだろうが、散水回数も遮光度もそれほど大きく影響しているようには見えない。例えば、週に1回の散水を基準（=0）として考えると、週2回の散水も週3回の散水もそれほど顕著に花数を増やすとはいえない。たとえ週2回の散水が最も良いだろうと予測したとしても、標準化係数が1にも満たないということは大した効果はない（つまり1つも花をつけない）という解釈もできる。私は花についてはからっきしの素人なので、この結果がどのようなものであるかはよく分からない。。。
素人考えでは遮光なしがもっとも花数をつけるような気がするが、そうでもないようである。いずれの条件においても特に差はないようである。ということは遮光量によって花数は影響しないのかもしれないが、この分野の専門家はこの結果をどう読み解くのだろうか？
> summary(box5.6) # BOX5.6の係数表
Call:
glm(formula = SQBLOOMS ~ BED + WATER + SHADE)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.23992 -0.18979 -0.01840 0.17471 0.74686
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.7765 0.2203 17.139 2.23e-16 ***
BED2 0.3185 0.1908 1.669 0.106242
BED3 -0.5044 0.1908 -2.643 0.013293 *
WATER2 0.7842 0.1908 4.109 0.000313 ***
WATER3 0.4489 0.1908 2.353 0.025905 *
SHADE2 0.1969 0.2203 0.894 0.379172
SHADE3 -0.2157 0.2203 -0.979 0.336068
SHADE4 -0.3673 0.2203 -1.667 0.106641
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 0.2184764)
Null deviance: 15.6114 on 35 degrees of freedom
Residual deviance: 6.1173 on 28 degrees of freedom
AIC: 56.357
Number of Fisher Scoring iterations: 2
# WATERの標準化係数のプロット
> WaterCoef <- c(0, 0.7842, 0.4489) # WATERの標準化係数を代入
> names(WaterCoef) <- c("Wat1", "Wat2", "Wat3") # x軸のラベル
> barplot(WaterCoef, ylim=c(0,1))
> abline(h=0)
# SHADEの標準化係数のプロット
> ShadeCoef <- c(0, 0.1969, -0.2157, -0.3673)
> windows() # 新しいグラフウィンドウを開く
> names(ShadeCoef) <- c("Sha1", "Sha2", "Sha3", "Sha4")
> barplot(ShadeCoef, ylim=c(-0.5, 0.5))
> abline(h=0)
雄のスムースイモリの背びれ
この問題はパス。
