2015年生物測定基礎実験

統計解析３

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

この実験の目的

　統計解析ソフトＲを用いて，統計解析の理論と実践を学ぶ

分散分析（Analysis of Variance）

因子と水準

一元配置（one-way layout）

構造モデル

平方和分解

平方和の期待値

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

分散分析表と F 検定

多重比較（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 である．

多重比較法の比較

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

a <- 10; 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) 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) # チューキー HSD 　 which(th$x[,4] < 0.05)

pairwise.t.test(y, x, p.adj = "none") # 対比較補正なし
だと，本来差が無いはずなのに，差があるとしてしまいゴーストを拾ってしまう．従って， 多くの比較を行う場合，補正しないといけない．

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

品種データダウンロード

品種収量の分散分析

hinsyu <- read.csv("hinsyu.csv"); hinsyu # csv データ読み込み boxplot(hinsyu) x <- rep(names(hinsyu), each=n) y <- as.vector(as.matrix(hinsyu))

問題１

1. コメデータの分散分析を行え．
2. コメデータの品種効果を固定効果とみなして多重比較を行い，有意な差が認められる組を求めよ．

多元配置

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="Rice yield to fertility levels（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] は交互作用である．

複数年次にわたる解析

構造モデル

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

  yield <- rice$gy # gy 列を取り出す dens <- factor(rice$density) #密度水準ラベル化 fert <- factor(rice$fert) year <- factor(rice$year) rya.aov <- aov(yield ~ dens + fert + year + dens:fert + fert:year + dens:year + dens:fert:year) # 3 次交互作用まで計算 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")

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

形質間相関

　収量と収量構成要素間での形質間相関をみてみる．相関の高い形質の組みはみられなかった．

• コメ形質間相関
  

  cor(rice[,7:11]) # 形質間相関 pairs(rice[,7:11]) # 形質間対散布図

コメ形質間対散布図

問題２

1. 田無農場の４年次にわたるデータで，収量構成要素のどれか２つに対し，分散分析を行い， 収量の結果との共通点や相違点を考察せよ．(主効果の他，交互作用の出方の違いなど．)

回帰分析（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)

  

ブートストラップ法による信頼区間の構成

• 自動車の速度と制動距離の回帰係数における，ブートストラップによる信頼区間の構成

  n <- length(y) M <- 1000 coeff.vec <- NULL for(i in 1:M){ nx <- sample(1:n, replace=TRUE) bb <- (lm(y[nx] ~ x[nx]))$coefficients[2] coeff.vec <- c(coeff.vec, bb) } hist(coeff.vec) quantile(coeff.vec, c(0.025, 0.975)) # 95％信頼区間

  

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

データダウンロード

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

問題３

1. 帰無仮説　H₀：b = 1
   の検定を y を頭部，前胸部，交尾器としたときでそれぞれ行え．
2. 回帰係数の 95％信頼区間をそれぞれ求めよ．
3. 回帰係数の 95％信頼区間をブートストラップ法で求め，先ほどの信頼区間と比較せよ．
4. 頭部，前胸部，交尾器の回帰係数が異なっている生物学的な理由を考察せよ．

参考文献（古い順）

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