東京農工大学

2020.10.30

生物統計学３

東京大学大学院農学生命科学研究科　大森宏

分散分析（Analysis of Variance）

因子と水準

一元配置（one-way layout）

構造モデル

平方和分解

平方和の期待値

帰無仮説のもとでの平方和の比の分布

分散分析表と F 検定

対比（contarst）

　対比 C = Σ_ic_iμ_i は C^{^} = Σ_ic_iX^-_i. で推定されるが， 帰無仮説（α_i = 0）のもとでは， 対比推定量の平均と分散は，

直交対比

多重比較（multiple comparison）

• なにもしない
  推奨されない方法である．すなわち，補正なしの t 検定を行う．昔の LSD である．
  R では，|d_ij| の p 値が出力 される．
• ボンフェローニ（Bonferroni）補正
  いま，有意水準 α' のそれぞれ独立な検定を r 回行ったとすると，1 回の検定で正しい判断を行う確率が 1 - α' なので，r 回の検定で正しい判断を行う確率は，(1 - α')^r となる．よって，正しい判断 を行わない（第 1 種の過誤の）確率は，
  1 - (1 - α')^r ≒ 1 - (1 - rα') = rα'，ただし，α' ≒ 0
  となる．これが，r 回の検定全体での有意水準となる．よって，検定全体での有意水準を α にするには， 1 回の検定の有意水準を α' = α/r にすればよい．これがボンフェローニ補正である．しかし， 多重比較における検定は独立な検定ではないので，この補正は厳しすぎ（保守的）て， 有意な組み合わせが見つからない恐れがある．
  R の多重比較では，補正なしの p 値を r 倍した p 値を出力する．ただし，これが 1 を超えた場合は 1 とする．
• ホルム（Holm）補正
  ボンフェローニ補正を改良したものである．すべての比較組み合わせ（対比）の t 値を計算し， それを大きさの順に並べる．一番大きな t 値 t₍₁₎ の有意確率を α/r，次の大きさ の t₍₂₎ の有意確率を α/(r - 1)，というように有意確率を調整する．
  R ではホルム補正がデフォルトで， p 値を大きさの順に並べ最も小さな p 値を r 倍し，次に大きな p 値を r - 1 倍して出力するようである．
• チューキー（Tukey）の HSD（honestly significant difference）
  今までは，t 検定の有意確率を補正することにより，多重比較の問題に対処していたが， スチィーデント化された範囲の分布 （Studentized range distribution）という多重比較専用の分布 を用いて検定する．2 つの処理 i，j 間の比較 を行うときに用いる検定統計量は，先ほどの t_ij である．

多重比較法の比較

シミュレーションデータに対する分散分析と対比、多重比較

# start of script a <- 5; n <- 10 # 処理水準数 a，処理内標本数 n x <- rep(1:a, each=n) x <- factor(x) # ラベル化 y <- rnorm(n*a) # 標準正規乱数 cbind(y, x) # boxplot(y ~ x) av <- aov(y ~ x) # 分散分析 summary(av) av2 <- lm(y ~ x) # 線形モデル anova(av2) # 分散分析 # 対比 c1 <- c(1/3,1/3,1/3,-1/2,-1/2) # 対比係数ベクトル c2 <- c(0, 0, 0, 1, -1) c3 <- c(1, -1/2, -1/2, 0, 0) contrasts(x) <- cbind(c1, c2, c3) # 対比の定義 contrasts(x) # 定義した残りの直交対比は R が勝手につくる fc <- lm(y ~ x) # 対比つき線形モデル summary(fc) # 結果表示（最初の3つの対比のみに着目） # 多重比較 pt <- pairwise.t.test(y, x); pt # 対比較ホルム補正 which(pt$p.value < 0.05, arr.ind=T) pt <- pairwise.t.test(y, x, p.adj = "bonf"); pt # 対比較ボンフェローニ補正 which(pt$p.value < 0.05, arr.ind=T) pt <- pairwise.t.test(y, x, p.adj = "none"); pt # 対比較補正なし which(pt$p.value < 0.05, arr.ind=T) th <- TukeyHSD(av); th # チューキー HSD which(th$x[,4] < 0.05) # end of script

品種によるコメの収量の違い

データダウンロード

データの品種ごとのちらばりと，品種平均と標準誤差を視覚化する

yy <- read.csv("hinsyu.csv"); yy # csv データ読み込み n <- nrow(yy) a <- ncol(yy) nn <- n*a # 総データ数 xm <- 1:a ym <- mean(yy) # 品種ごとの平均 se <- sum(diag(var(yy)*(n-1))) # 誤差平方和 ve <- se/(nn - a); ve # 誤差平均平方（誤差分散） vm <- ve/n # 品種平均の分散 svm <- sqrt(vm) # 品種平均の標準誤差 stripchart(yy, xlab="Yield (kg/a)") # 品種ごとの個体分布表示 text(ym, xm,"|", col="red") # 品種平均 text(ym+svm, xm,"|", col="red") # 品種平均 ＋ 標準誤差 text(ym-svm, xm,"|", col="red") # 品種平均 － 標準誤差 segments(ym-svm, xm, ym+svm, xm, col="red")

多元配置

実験計画法

因子と水準

　作物において，収量や抵抗性などの調査研究したい特性を形質(character)という． 着目した形質に影響を与えると考えられる品種，温度，施肥量などを 要因または因子(factor)という．要因の影響を調べるためいくつかの品種を用いたり， 施肥量に段階を設けたりするが，それを水準(level)という．
　要因と水準を決定したら，その反応を調べるためには，圃場で実際に栽培してみる 必要がある場合が多い． しかしながら，圃場の環境条件は場所により異なり，均一の条件で栽培管理を行うことは 難しい．環境条件に多少のばらつきがあると想定される場合には，ばらつきの影響を無作為化 する必要がある．これを行うのが，実験計画法である．

Fisher の3原則

　実験計画の創始者である R.A. Fisher が提唱した原則で， 反復(replication)，無作為化(randomization)，局所管理(local control) の3つからなる．これらは， 処理の違いに基づく差が偶然の誤差を超えて有意なものかどうか統計的に判定するため， 実験を行うにあたり必要な管理原則である．この原則に基づいて圃場での実験配置が設計される．

• 反復
  1つの処理組み合わせに対して，少なくとも2回以上の繰り返しが必要である． 繰り返しにより，誤差の大きさが評価できる．
• 無作為化
  実験順序や，空間配置の無作為化．ある特定の処理条件を圃場のある特定の場所にかためて しまうと，その場所がたまたま水はけがよかったり肥沃だったりした場合，処理効果によるものか 場所の効果によるものかの判別ができなくなる． このようなある特定の場所による効果を系統誤差(systematic error)というが， 無作為に処理条件を配置させることにより，系統誤差を偶然誤差(random error)に転化させる ことができる．
• 局所管理
  実験が大規模で，実験全体を無作為化したら誤差が大きくなってしまうような場合， 実験をある程度細分化してブロックを構成し，ブロック内で処理条件を無作為化する． ブロック内のバックグラウンドはなるべく均一になるように管理できれば， 系統誤差の一部がブロック間変動として除去できるので，実験誤差が 小さくなる．

よく用いられる試験配置

• 乱塊法(randomized block design)
  Fisher の3原則を完全に満足する実験手法で，すべての処理組み合わせの実験を 1回ずつ集めたもので 1つのブロックを形成する．ブロック数が反復数になる．
• 分割区法(split-plot design)
  潅水区や非潅水区のような水管理は，比較的広い区画で行わないと労力，費用等が かかる．このような処理を用いた実験の場合，乱塊法ででブロック内ですべての組み合わせの試験区 を設定することは現実的でない．水管理のように 広い区画が必要な因子 A を一次因子，その区画内で水準を変化させる品種や施肥量などの 因子 B を二次因子と呼ぶ．

2000年度試験配置
データダウンロード

特定年次の解析

構造モデル

処理平方和

F 検定

分散分析表の見方

• コメ収量（2000年）の分散分析
  

  rice <- read.csv("ricecul.csv") # データ読み込み rice yield <- rice$gy dense <- factor(rice$density) # 水準のラベル化 fert <- factor(rice$fert) blk <- factor(rice$rep) tapply(yield[1:12], dense[1:12], mean) # 密度水準ごとの平均 tapply(yield[1:12], fert[1:12], mean) tapply(yield[1:12], blk[1:12], mean) cm <- tapply(yield[1:12], dense[1:12]:fert[1:12], mean); cm # 処理組み合わせごとの平均 # 施肥量水準に対する収量のグラフ plot(1:3, cm[1:3], type="b", lwd=2, cex=1.5, xaxt="n", xlab="fertility", ylab="yield", ylim=c(300, 600), pch=0, cex.lab=1.0, cex.axis=1.0, col="blue") axis(1, 1:3, labels=c("fert1","fert2","fert3"), cex.axis=0.8) points(1:3, cm[4:6], type="b", lwd=2, cex=1.5, pch=2, col="red") legend(1.2, 580, legend=c("density1","density2"), pch=c(0,2), col=c("blue","red"), cex=0.8) title(main="施肥量水準に対するコメ収量（2000年）") # 分散分析 ry.aov <- aov(yield[1:12] ~ blk[1:12] + dense[1:12] + fert[1:12] + dense[1:12]:fert[1:12]) summary(ry.aov) #分散分析表の表示

  注）blk[1:12] などは主効果であり，dens[1:12]:fert[1:12] は交互作用である．

複数年次にわたる解析

構造モデル

X_ijkl = μ + α_i + β_j + γ_k + (αβ)_ij + (βγ)_jk + (αγ)_ik + (αρ)_il + ε_ijkl

　分散分析表をみると，年次効果と施肥量効果は主効果と交互作用ともに有意 となったが，密度に関しては有意とはならなかった．どのような効果が存在するかをみるため， 年次と施肥量水準の 組み合わせごとの平均を図にプロットした．
　これをみると， 年次ごとに収量に変動があることがわかる．また， 無肥料区(fert1)と少肥料区(fert2)は 年次ごとの変動が似ていて，少肥料区の方が収量が高くなっているが， 多肥料区(fert3)は2000年のように高収量のときもあるが，1996年や1997年のようにそれほど高い収量が得られない場合 もあり，これが交互作用として検出されている． これは，多肥料区で繁茂しすぎて 倒伏が起こった年があったためと考えられる．

「日本晴」収量の年次効果と施肥水準によるプロット（1996 ～ 2000年度）

• コメ収量複数年次の分散分析
  

  yield <- rice$gy # gy 列を取り出す dens <- factor(rice$density) #密度水準ラベル化 fert <- factor(rice$fert) year <- factor(rice$year) blk <- factor(rice$rep) rya.aov <- aov(yield ~ dens + fert + year + dens:fert + fert:year + dens:year + year:blk) summary(rya.aov) # #「日本晴」収量の年次効果と施肥水準によるプロット x <- tapply(rice$gy, fert:year, mean) plot(1:4, x[1:4], type="b", lwd=2, cex=1.5, xaxt="n", xlab="year", ylab="yield", ylim=c(300,600), pch=0, lty=1,cex.lab=1.2, cex.axis=1.2, col="red") axis(1, 1:4, labels=levels(year), cex.axis=1.2) box(lwd=3) points(1:4, x[5:8], type="b", pch=2, lwd=2, lty=1, col="blue", cex=1.5) points(1:4, x[9:12], type="b", pch=3, lwd=2, lty=1, col="green", cex=1.5) legend(2, 600, legend=c("fert1","fert2","fert3"), pch=c(0,2,3), col=c("red","blue","green"), cex=1.2) title(main="Yield of 'Nipponbare' to year effects and fertility levels")

回帰分析（Regression Analysis）

直線回帰

モデル

最小２乗法

　データに最もよくあてはまる直線回帰式を得るには，データ点 （x_i ，y_i ）， と回帰による推定点， （x_i ，y^_i ）， y^_i ＝ a ＋ b x_i ， の間の距離の２乗和 S が最小になるような回帰係数 a ，b を求める．つまり，

　これは，S を a ，b で偏微分して 0 とおくことによって得られる．つまり，

回帰式の統計モデル

　推定された直線回帰式がどの程度現実のデータに適合しているかを調べるために， 回帰式が従う統計モデルを考える．標本の格データ点， （x_i ，y_i ）， が，

残差分散と回帰係数の標準誤差

　回帰で説明がつかない残差平方和 S_e は，

　一般に，Var(y_i ) ＝ σ² であるとき，その定数（c） 倍の分散は，

回帰係数の標準誤差による t 検定

目的変数 y が説明変数 x との回帰関係にないという 帰無仮説，

分散分析

平方和分解

　回帰式により， 従属変数 y のデータ y_i は，

決定係数（重相関係数の２乗）

　データが直線回帰式でよく説明できるのは，回帰平方和が大きく，残差平方和 が小さい場合である．総平方和のうち回帰平方和で説明される割合を決定係数，もしくは 重相関係数の２乗といい，

回帰分析の分散分析表

変動因	平方和	自由度	平均平方	F 値
回帰	SR	1	SR	F ＝ SR/se2
残差	Se	n－2	se2 ＝ Se/n－2
全体	ST	n－1

F 検定

従属変数 y が説明変数 x の回帰関係にないという 帰無仮説，

自動車の速度と制動距離
自動車の速度 x mile/h とブレーキをかけてからの静止するまでの距離 y ft のデータに対し， y が x の累乗で表せると考えられるときは， 両者の対数変換を行い，直線回帰の問題に帰着させることができる．すなわち，
y = ax^b －＞ logy = loga + b logx
の関係を利用する．これの回帰分析を行ってみる．

cars y <- log(cars$dist) x <- log(cars$speed) plot(x, y, xlab="log of velocity", ylab="log of braking distance") car.lm <- lm(y ~ x) abline(car.lm, col="red") summary(car.lm)

回帰係数に対する検定

回帰式の信頼区間

回帰係数の信頼区間

　回帰係数の標準誤差 s _b を用いて， 回帰係数 b の信頼区間がつくれる．すなわち， 自由度 n－2 の t 分布の 97.5％点を t₀ とすると， 回帰係数 b の 95％信頼区間の幅 d は，d ＝ t₀ s _b となるので， 95％信頼区間は，

回帰直線の信頼区間

回帰予測値の信頼幅

• 自動車の速度と制動距離の回帰係数に対する検定
  　物理法則によると，自動車の速度の２乗と制動距離が比例するので，求めた回帰係数は 2 となることが 期待される．このため，H₀：b ＝ 2，の検定を行ってみる．
  　これは，回帰係数の信頼区間を求め，これが 2 を含むかどうかでも調べることもできる．

  bh <- summary(car.lm)$coefficient[2,1] # 回帰係数 sb <- summary(car.lm)$coefficient[2,2] # 回帰係数の標準誤差 tv <- abs(bh - 2)/sb; tv edf <- summary(car.lm)$df[2] # 残差分散の自由度 1 - pt(tv, df=edf) # 有意確率 t0 <- qt(0.975, df=edf) d <- t0*sb bh - d # 95％信頼区間の下限 bh + d # 95％信頼区間の上限 new <- data.frame(x=seq(1.3, 3.5, by=0.02)) dc <- predict(car.lm, new, interval="confidence", level=0.95) # 回帰推定値（回帰直線）の 95 ％信頼幅 dp <- predict(car.lm, new, interval="prediction", level=0.95) # 回帰予測値の 95 ％信頼幅 matplot(new$x, cbind(dc,dp[,-1]), lty=c(1,2,2,3,3), type="l", col=c("blue","blue","blue","red","red"), xlab="log of velocity", ylab="log of braking distance") points(x, y)

  

ミヤマクワガタの相対成長解析

データダウンロード

成長段階の異なる 47 頭のミヤマクワガタのパーツ別の重量データ(g)．パーツは， 頭部（WHEAD），前胸部（WTHORAX），中胸～腹部（WABDOM），交尾器（WGENI）からなる．

参考文献（古い順）

1. Introduction to the Theory of Statistics, Mood, A. M., Graubill, F. A. & Boes, D. C., 1974, McGRAW-HILL
2. 「実験」生産環境生物学，東京大学大学院農学生命科学研究科生産・環境生物学専攻編，1999，朝倉書店
3. 工学のためのデータサイエンス入門－フリーな統計環境Rを用いたデータ解析－，間瀬茂ら，2004， 数理工学社
4. 実践生物統計学－分子から生態まで－（第 1 章，第 2 章）， 東京大学生物測定学研究室編（大森宏ら）， 2004，朝倉書店
5. The R Tips データ解析環境 R の基本技・グラフィックス活用集，船尾暢男，2005，九天社
6. R で学ぶデータマインニング I －データ解析の視点から－，熊谷悦生・船尾暢男，2007，九天社
7. R で学ぶデータマインニング II －シミュレーションの視点から－，熊谷悦生・船尾暢男，2007，九天社