氏名：林　亜矢子

所属研究科：医療データ人材

学年：科目履修生

学生番号：0699320805

分散分析について講義で習ったことを演習しレポートとして提出しましょう

この Rmd ファイルは穴開きになっています．必要事項を自分で記入して完成させましょう

記入毎に Ctrl + s でこまめに保存すること

Rmd ファイルの出力は Knit を押して確認すると良い

提出方法

作業を完了させたこの Rmd ファイルと Knit してつくった html ファイルを PandA から提出すること

提出要領は PandA の「課題」を確認すること

注意

Rのコードチャンク以外にも記入すべきところがあるので，細部まできちんと確認すること

レポート問題

（参考：実践生物統計学，東京大学生物測定学研究室編，2004）

圃場試験で栽培条件の影響を調べたい．水稲の９品種をそれぞれ６区画の水田で栽培したときのアール当たりの玄米重量は以下のとおりであった．

品種 A1 A2 A3 B1 B2 C D1 D2 D3
34.3 33.5 35.2 36.0 37.3 38.3 35.1 36.2 34.6
35.9 35.2 36.5 37.5 38.5 38.0 37.0 38.0 35.5
36.8 36.1 37.9 38.1 39.4 36.8 38.6 39.9 36.4
37.4 34.8 38.6 39.5 40.5 37.4 36.5 37.9 37.8
36.0 35.3 37.1 38.6 40.1 35.9 38.4 39.1 35.9
37.0 36.5 38.1 36.9 38.8 36.6 35.7 37.5 35.6
  1. 収穫量の上のデータに対して，品種によって差があるかどうかの分散分析を有意水準 5% で行なえ．
  2. 多重比較を行い，有意な差がみとめられる品種の組を求めよ．

ただし，以下の設問のとおりにおこなうこと．

検定の設定

帰無仮説 \(H_0\) 品種による収穫量に差はない

有意水準 \(\alpha=0.05\) で一元配置分散分析をおこなう

データ格納を行なえ

R で標本をデータフレームとして格納せよ

# groupA1, groupA2, groupA3, ... , groupD3 までのデータを格納せよ
A1 <- c(34.3,35.9,36.8,37.4,36.0,37.0)
A2 <- c(33.5,35.2,36.1,34.8,35.3,36.5)
A3 <- c(35.2,36.5,37.9,38.6,37.1,38.1) 
B1 <- c(36.0,37.5,38.1,39.5,38.6,36.9) 
B2 <- c(37.3,38.5,39.4,40.5,40.1,38.8) 
C  <- c(38.3,38.0,36.8,37.4,35.9,36.6) 
D1 <- c(35.1,37.0,38.6,36.5,38.4,35.7) 
D2 <- c(36.2,38.0,39.9,37.9,39.1,37.5) 
D3 <- c(34.6,35.5,36.4,37.8,35.9,35.6) 



# data frame を構成し，DTライブラリを用いて描画せよ
datf = data.frame(group1=A1, group2=A2, group3=A3, group4=B1, group5=B2, group6=C, group7=D1, group8=D2, group9=D3)



# plot
DT::datatable(datf, extensions='FixedColumns', options = list(dom='t', scrollX=TRUE, scrollCollapse=FALSE))

各パラメータを代入せよ

一元配置分散分析で用いるパラメータの設定をおこなえ

# 全データを格納する配列 data をつくれ
data = c(A1,A2,A3,B1,B2,C,D1,D2,D3)
data
##  [1] 34.3 35.9 36.8 37.4 36.0 37.0 33.5 35.2 36.1 34.8 35.3 36.5 35.2 36.5 37.9
## [16] 38.6 37.1 38.1 36.0 37.5 38.1 39.5 38.6 36.9 37.3 38.5 39.4 40.5 40.1 38.8
## [31] 38.3 38.0 36.8 37.4 35.9 36.6 35.1 37.0 38.6 36.5 38.4 35.7 36.2 38.0 39.9
## [46] 37.9 39.1 37.5 34.6 35.5 36.4 37.8 35.9 35.6
# 総標本サイズ N を計算せよ
N = length(data)
N
## [1] 54
# 各群の標本サイズを計算せよ
n1 = length(A1); n1
## [1] 6
n2 = length(A2); n2
## [1] 6
n3 = length(A3); n3
## [1] 6
n4 = length(B1); n4
## [1] 6
n5 = length(B2); n5
## [1] 6
n6 = length(C); n6
## [1] 6
n7 = length(D1); n7
## [1] 6
n8 = length(D2); n8
## [1] 6
n9 = length(D3); n9
## [1] 6
# 群の数をmに代入せよ
m = 6
m
## [1] 6
m-1
## [1] 5
N-m
## [1] 48
N-1
## [1] 53

平均の計算をせよ

各群の平均を計算せよ

# 群内平均
mean1 = mean(A1)
mean1
## [1] 36.23333
mean2 = mean(A2)
mean2
## [1] 35.23333
mean3 = mean(A3)
mean3
## [1] 37.23333
mean4 = mean(B1)
mean4
## [1] 37.76667
mean5 = mean(B2)
mean5
## [1] 39.1
mean6 = mean(C)
mean6
## [1] 37.16667
mean7 = mean(D1)
mean7
## [1] 36.88333
mean8 = mean(D2)
mean8
## [1] 38.1
mean9 = mean(D3)
mean9
## [1] 35.96667

総平均を計算せよ

# 総平均
meanT = mean(data)
meanT
## [1] 37.07593

平方和の計算をせよ

群内平方和と群間平方和をそれぞれ計算せよ

# 各群の平方和を計算せよ
SS1 = sum((A1-mean1)^2)
SS1
## [1] 6.173333
SS2 = sum((A2-mean2)^2)
SS2
## [1] 5.553333
SS3 = sum((A3-mean3)^2)
SS3
## [1] 7.753333
SS4 = sum((B1-mean4)^2)
SS4
## [1] 7.753333
SS5 = sum((B2-mean5)^2)
SS5
## [1] 6.74
SS6 = sum((C-mean6)^2)
SS6
## [1] 4.093333
SS7 = sum((D1-mean7)^2)
SS7
## [1] 9.988333
SS8 = sum((D2-mean8)^2)
SS8
## [1] 8.26
SS9 = sum((D3-mean9)^2)
SS9
## [1] 5.773333
# 群内平方和 SSW
SSW = SS1+SS2+SS3+SS4+SS5+SS6+SS7+SS8+SS9
SSW
## [1] 62.08833
# 群間平方和 SSA
SSA= n1*(mean1-meanT)^2 + n2*(mean2-meanT)^2 +n3*(mean3-meanT)^2+n4*(mean4-meanT)^2 + n5*(mean5-meanT)^2 +n6*(mean6-meanT)^2+n7*(mean7-meanT)^2 + n8*(mean8-meanT)^2 +n9*(mean9-meanT)^2
SSA
## [1] 66.17037

総平方和を計算して確認せよ

# 総平方和 SS を計算せよ
SS = sum((data-meanT)^2)
SS
## [1] 128.2587
# 総平方和が群内平方和と群間平方和の和になっていることを確認せよ
SSA + SSW
## [1] 128.2587

平均平方の計算を計算せよ

# 群内平均平方 MSW を計算せよ
MSW = SSW / (N - m)
MSW
## [1] 1.293507
# 群間平均平方 MSA を計算せよ
MSA = SSA/(m-1)
MSA
## [1] 13.23407
# 総平均平方 MS を計算せよ
MS = SS / (N-1)
MS
## [1] 2.419976

\(F\) 統計量を計算せよ

# F統計量 Fstar を計算せよ
Fstar = MSA / MSW
Fstar
## [1] 10.23116

棄却限界を求めよ

# F(m-1, N-m) の上側5％点 F0 を計算せよ
F0 = qf(0.95, m - 1, N - m)
F0
## [1] 2.408514

検定の結果を述べよ

この部分に検定の結果を述べよ

\(F^*\) が \(F(m-1, N-m)\) の上側5％点 \(F_0\) より上回っているので，有意水準 \(\alpha=0.05\) の \(F\) 検定によって帰無仮説を棄却する
(Fstar = 10.04676, F0 = 2.408514)

分散分析表を作成せよ

上で計算した値を下の表に（手入力はダサいけど）埋めて，分散分析を完成させよ

変動要因 平方和 自由度 \(\mathrm{df}\) 平均平方 分散比 \(F\)
群間変動(処理) \(\mathrm{66.17}\) \(5\) \(\mathrm{13.23}\) \(\mathrm{10.23}\)
群内変動(誤差) \(\mathrm{62.09}\) \(48\) \(\mathrm{1.294}\)
全体 \(\mathrm{128.3}\) \(53\) \(\mathrm{2.420}\)

検算

R の aov 関数を用いて良いので上で行った計算を検算せよ

# 群のラベルを表す配列 group を作成せよ
group = c(rep(1,n1),rep(2,n2),rep(3,n3),rep(4,n4),rep(5,n5),rep(6,n6),rep(7,n7),rep(8,n8),rep(9,n9))


# group をカテゴリカル配列に変更せよ
group = factor (group)
group
##  [1] 1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 3 3 4 4 4 4 4 4 5 5 5 5 5 5 6 6 6 6 6 6 7 7
## [39] 7 7 7 7 8 8 8 8 8 8 9 9 9 9 9 9
## Levels: 1 2 3 4 5 6 7 8 9
# aov の結果を result に格納せよ
result = aov(data ~ group)

# aov の結果の要約を出力せよ
summary(result)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## group        8  66.17   8.271   5.995 3.09e-05 ***
## Residuals   45  62.09   1.380                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Tukey HSD 検定を行え

群の平均値に統計的有意な差があるかどうかを調べるための多重比較を Tukey HSD 法を用いて行え

有意水準を0.05とする

# 各群の箱ひげ図を描画し，各群の要約統計量を可視化せよ
boxplot(data ~ group)

# TukeyHSD 関数で多重比較を行い，結果を THSD に格納せよ
THSD = TukeyHSD(result)
THSD
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = data ~ group)
## 
## $group
##            diff        lwr          upr     p adj
## 2-1 -1.00000000 -3.2088892  1.208889184 0.8611598
## 3-1  1.00000000 -1.2088892  3.208889184 0.8611598
## 4-1  1.53333333 -0.6755559  3.742222518 0.3863696
## 5-1  2.86666667  0.6577775  5.075555851 0.0033493
## 6-1  0.93333333 -1.2755559  3.142222518 0.9008132
## 7-1  0.65000000 -1.5588892  2.858889184 0.9877618
## 8-1  1.86666667 -0.3422225  4.075555851 0.1589303
## 9-1 -0.26666667 -2.4755559  1.942222518 0.9999804
## 3-2  2.00000000 -0.2088892  4.208889184 0.1039413
## 4-2  2.53333333  0.3244441  4.742222518 0.0140468
## 5-2  3.86666667  1.6577775  6.075555851 0.0000288
## 6-2  1.93333333 -0.2755559  4.142222518 0.1291023
## 7-2  1.65000000 -0.5588892  3.858889184 0.2919204
## 8-2  2.86666667  0.6577775  5.075555851 0.0033493
## 9-2  0.73333333 -1.4755559  2.942222518 0.9740672
## 4-3  0.53333333 -1.6755559  2.742222518 0.9967255
## 5-3  1.86666667 -0.3422225  4.075555851 0.1589303
## 6-3 -0.06666667 -2.2755559  2.142222518 1.0000000
## 7-3 -0.35000000 -2.5588892  1.858889184 0.9998447
## 8-3  0.86666667 -1.3422225  3.075555851 0.9326029
## 9-3 -1.26666667 -3.4755559  0.942222518 0.6382878
## 5-4  1.33333333 -0.8755559  3.542222518 0.5737355
## 6-4 -0.60000000 -2.8088892  1.608889184 0.9927373
## 7-4 -0.88333333 -3.0922225  1.325555851 0.9253890
## 8-4  0.33333333 -1.8755559  2.542222518 0.9998924
## 9-4 -1.80000000 -4.0088892  0.408889184 0.1938180
## 6-5 -1.93333333 -4.1422225  0.275555851 0.1291023
## 7-5 -2.21666667 -4.4255559 -0.007777482 0.0485871
## 8-5 -1.00000000 -3.2088892  1.208889184 0.8611598
## 9-5 -3.13333333 -5.3422225 -0.924444149 0.0009926
## 7-6 -0.28333333 -2.4922225  1.925555851 0.9999688
## 8-6  0.93333333 -1.2755559  3.142222518 0.9008132
## 9-6 -1.20000000 -3.4088892  1.008889184 0.7010313
## 8-7  1.21666667 -0.9922225  3.425555851 0.6856156
## 9-7 -0.91666667 -3.1255559  1.292222518 0.9094981
## 9-8 -2.13333333 -4.3422225  0.075555851 0.0657200

帰無仮説の下でどの組が有意水準以下の出現確率となるか調べよ

# 有意水準以下となる有意確率を出す群の組を表示せよ
which(THSD$group[,4] < 0.05)
## 5-1 4-2 5-2 8-2 7-5 9-5 
##   4  10  11  14  28  30

Tukey HSD 検定の結果をここに述べよ

品種による収穫量の差で、
A1とB2の差（p = 0.0033493)
A2とB1の差（p = 0.0140468)
A2とB2の差（p = 0.0000288)
A2とD2の差（p = 0.0033493)
B2とD1の差（p = 0.0485871)
B2とD3の差（p = 0.0009926)
の6ペアは有意差が認められました（P < 0.05）。

以上