まずはホームディレクトリを設定します。
home <- [配布したフォルダをコピーしたパス]
setwd(home)
分散分析について説明します。t検定は分散分析の特殊ケースで、t検定では二群の平均値を比較する問題でしたが、分散分析では、より一般的に二つ以上の群間で平均値を比較することができます。例えば、以下のように、3つの条件(例えば、薬剤の投与条件など)における各遺伝子の発現量を3サンプルずつ測定す、条件間で平均値の異なるような遺伝子群を同定したいときに、分散分析を用いることができます。
set.seed(100)
n <- 10
p <- 5
cond1 <- matrix(rnbinom(n*p,2,mu = 100),ncol=p,dimnames=list(paste0("Gene",1:n),paste0("Condition",rep(1,p))))
cond2 <- matrix(rnbinom(n*p,2,mu = 100),ncol=p,dimnames=list(paste0("Gene",1:n),paste0("Condition",rep(2,p))))
cond3 <- matrix(rnbinom(n*p,2,mu = 300),ncol=p,dimnames=list(paste0("Gene",1:n),paste0("Condition",rep(3,p))))
data_example <- cbind(cond1,cond2,cond3)
Condition1 | Condition1 | Condition1 | Condition1 | Condition1 | Condition2 | Condition2 | Condition2 | Condition2 | Condition2 | Condition3 | Condition3 | Condition3 | Condition3 | Condition3 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Gene1 | 36 | 28 | 155 | 49 | 186 | 87 | 58 | 87 | 46 | 33 | 366 | 67 | 321 | 253 | 542 |
Gene2 | 140 | 78 | 66 | 47 | 149 | 163 | 246 | 54 | 96 | 135 | 392 | 50 | 863 | 450 | 437 |
Gene3 | 90 | 44 | 20 | 49 | 16 | 102 | 3 | 233 | 58 | 49 | 543 | 206 | 543 | 510 | 1015 |
Gene4 | 100 | 65 | 87 | 68 | 71 | 127 | 64 | 128 | 52 | 284 | 696 | 343 | 293 | 175 | 370 |
Gene5 | 81 | 104 | 58 | 192 | 25 | 345 | 62 | 189 | 39 | 45 | 696 | 361 | 137 | 151 | 282 |
Gene6 | 72 | 113 | 29 | 156 | 61 | 89 | 37 | 294 | 125 | 169 | 1162 | 102 | 547 | 137 | 204 |
Gene7 | 40 | 144 | 124 | 91 | 28 | 24 | 141 | 70 | 29 | 166 | 186 | 324 | 195 | 64 | 272 |
Gene8 | 172 | 263 | 76 | 36 | 85 | 37 | 96 | 228 | 7 | 195 | 394 | 186 | 1046 | 459 | 209 |
Gene9 | 25 | 47 | 5 | 100 | 92 | 86 | 75 | 46 | 329 | 121 | 540 | 47 | 561 | 67 | 288 |
Gene10 | 99 | 52 | 73 | 85 | 210 | 93 | 43 | 89 | 62 | 119 | 64 | 126 | 586 | 357 | 345 |
一つ目の遺伝子(Gene1)に対して
(groups <- colnames(data_example))
[1] "Condition1" "Condition1" "Condition1" "Condition1" "Condition1" "Condition2" "Condition2" "Condition2" "Condition2" "Condition2" "Condition3" "Condition3"
[13] "Condition3" "Condition3" "Condition3"
(x <- data_example[1,])
Condition1 Condition1 Condition1 Condition1 Condition1 Condition2 Condition2 Condition2 Condition2 Condition2 Condition3 Condition3 Condition3 Condition3 Condition3
36 28 155 49 186 87 58 87 46 33 366 67 321 253 542
(x2 <- data.frame(expr=x,groups=groups))
(lm_fit <- lm(expr~0 + factor(groups),x2))
Call:
lm(formula = expr ~ 0 + factor(groups), data = x2)
Coefficients:
factor(groups)Condition1 factor(groups)Condition2 factor(groups)Condition3
90.8 62.2 309.8
lm_fit$contrasts
$`factor(groups)`
[1] "contr.treatment"
(result <- anova(lm_fit))
Analysis of Variance Table
Response: expr
Df Sum Sq Mean Sq F value Pr(>F)
factor(groups) 3 540448 180149 15.052 0.0002274 ***
Residuals 12 143620 11968
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
result$`Pr(>F)`[1]
[1] 0.0002274161
いま、Gene1に対して、
groups <- colnames(data_example)
pvals <- NULL
for(i in 1:n){
x <- data_example[i,]
x2 <- data.frame(expr=x,groups=groups)
lm_fit <- lm(expr~0 + factor(groups),x2)
result <- anova(lm_fit)
pvals[i] <- result$`Pr(>F)`[1]
}
(qvals <- p.adjust(pvals,method="BH"))
[1] 0.0005685403 0.0011965310 0.0003635731 0.0003635731 0.0029166489 0.0170762568 0.0003635731 0.0032017675 0.0053416020 0.0015868761
もし、二群間の比較であれば、t-検定と同じ結果を与えます。
data_example2 <- cbind(cond1,cond3)
groups <- colnames(data_example2)
pvals_anova <- NULL
for(i in 1:n){
x <- data_example2[i,]
x2 <- data.frame(expr=x,groups=groups)
lm_fit <- lm(expr~factor(groups),x2)
result <- anova(lm_fit)
pvals_anova[i] <- result$`Pr(>F)`[1]
}
(qvals_anova <- p.adjust(pvals_anova,method="BH"))
[1] 0.07687011 0.07687011 0.03995054 0.04594354 0.08388471 0.12430750 0.07687011 0.09120758 0.08388471 0.09269908
pvals_t <- NULL
for(i in 1:n){
x <- data_example2[i,1:5]
y <- data_example2[i,6:10]
pvals_t[i] <- t.test(x,y,var.equal = TRUE)$p.value
}
(qvals_t <- p.adjust(pvals_t,method="BH"))
[1] 0.07687011 0.07687011 0.03995054 0.04594354 0.08388471 0.12430750 0.07687011 0.09120758 0.08388471 0.09269908
round(pvals_anova,3)==round(pvals_t,3)
[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
round(qvals_anova,3)==round(qvals_t,3)
[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
つまりt検定は分散分析の特殊ケースです。また、分散分析では複数条件間の差異を一度に検定できるだけでなく、条件間の交互作用なども自在にモデリングできるため、実験デザインに合わせてモデル化することが可能である点も重要です。