多重比較検定と一元配置分散分析法(ANOVA)について
独立した群が３群以上あるとき、どの群とどの群の平均値に有意差があるか、評価をしたくなる場合が、研究の現場ではよくある話だと思います。
しかし、一般的に、3群以上ある場合にそれぞれの群同士を多重比較を用いて、2群での検定を何度も繰り返す多重比較検定は良くない手段とされています。なぜなら、有意水準が仮定の設定値以上に上昇してしまうからです。（多重性の問題） (なぜ、有意水準が上昇するかについては、私の実力では説明が難しい...泣)そのため、一元配置分散分析法(ANOVA)によって、群間同士に差があるかを評価するのが一般的だとされています。しかし、ANOVAでは、どの因子とどの因子同士で差がみられるのかを考察することができません。
多重性の問題を解決し、多重比較検定を行うための方法として、事前に有意水準を狭め、のちに全体の有意水準を5%に調整する手法がある。代表的なのがTukey‒Kramer法です。Tukey‒Kramer法では、各群のデータ数(n)が一致していなくても問題ないため、使い勝手がいい特徴があります。
それでは実際に検定を行っていきます
必要なライブラリーと、使用するデータを読み込む
#multcompのインストール
#install.packages("multcomp", repos="http://cran.ism.ac.jp/")
library(multcomp)
setwd("パスの指定")
data<-read.table("boxplot_SG_CL.txt",header = T)
head(data)
# sample_no sample_class heatmap_class SG CL
# 1 6 10a 3 0.5318429 0.8482768
# 2 7 10b 3 0.5245027 0.8566048
# 3 8 10c 3 0.5535782 0.8916458
# 4 9 10d 3 0.5012658 0.8107002
# 5 10 10e 1 0.5120243 0.6768474
# 6 11 10f 3 0.5340269 0.8464090
●sample_no : サンプルナンバー
●sample_class : 導入した遺伝子（伏せ記号で表す）
●heatmap_class : ある遺伝子的特徴によって、group1~5に分類分け(コントロールサンプルとして、WT:wild type, Vc: Vector controlを採用)
●SG : 化学分析の結果(S-lignin/G-lignin比)
●CL : 化学分析の結果(Carbohydrate/ total lignin比)
今回使用したデータは、遺伝子改変されたシロイヌナズナを用いたデータになります。遺伝子改変した植物が、どのように化学成分が変化するか、を測定した結果になっております。
グループ毎のSG比を比較
SG_data<-data.frame(group=factor(data$heatmap_class),SG=data$SG)
plot(SG~group,d=SG_data) #boxplotを作成
ANOVA
SG_res1<-aov(SG~group,d=SG_data)
summary(SG_res1)
# Df Sum Sq Mean Sq F value Pr(>F)
# group 6 0.02518 0.004197 3.72 0.00188 **
# Residuals 133 0.15005 0.001128
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
ANOVAにより、1%有意の結果を示す。しかし、box plotを見る限り、大きな違いが見られない。
多重比較検定（Tukey‒Kramer法）
SG_res2<-glht(SG_res1,linfct=mcp(group="Tukey"))
summary(SG_res2)
#
# Simultaneous Tests for General Linear Hypotheses
#
# Multiple Comparisons of Means: Tukey Contrasts
#
#
# Fit: aov(formula = SG ~ group, data = SG_data)
#
# Linear Hypotheses:
# Estimate Std. Error t value Pr(>|t|)
# 2 - 1 == 0 0.005649 0.012937 0.437 0.99943
# 3 - 1 == 0 0.009094 0.011589 0.785 0.98547
# 4 - 1 == 0 0.004076 0.012482 0.327 0.99989
# 5 - 1 == 0 -0.024931 0.010430 -2.390 0.20517
# VC - 1 == 0 -0.008298 0.011135 -0.745 0.98891
# WT - 1 == 0 0.007968 0.013533 0.589 0.99690
# 3 - 2 == 0 0.003444 0.011854 0.291 0.99995
# 4 - 2 == 0 -0.001574 0.012728 -0.124 1.00000
# 5 - 2 == 0 -0.030580 0.010723 -2.852 0.07043 .
# VC - 2 == 0 -0.013947 0.011410 -1.222 0.88036
# WT - 2 == 0 0.002319 0.013761 0.169 1.00000
# 4 - 3 == 0 -0.005018 0.011355 -0.442 0.99939
# 5 - 3 == 0 -0.034025 0.009052 -3.759 0.00426 **
# VC - 3 == 0 -0.017392 0.009855 -1.765 0.56564
# WT - 3 == 0 -0.001126 0.012502 -0.090 1.00000
# 5 - 4 == 0 -0.029007 0.010170 -2.852 0.07038 .
# VC - 4 == 0 -0.012374 0.010891 -1.136 0.91272
# WT - 4 == 0 0.003893 0.013333 0.292 0.99995
# VC - 5 == 0 0.016633 0.008462 1.966 0.43311
# WT - 5 == 0 0.032899 0.011436 2.877 0.06611 .
# WT - VC == 0 0.016266 0.012081 1.346 0.82341
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# (Adjusted p values reported -- single-step method)
cld(SG_res2)
# 1 2 3 4 5 VC WT
# "ab" "ab" "b" "ab" "a" "ab" "ab"
結果の見方としては、値が低い順にa, b, c...と文字をつけており、別の文字が書かれている試料間(例えば、group3 「b」 とgroup5 「a」)で0.5%有意差があることを示す。
グループ毎のCL比を比較
CL_data<-data.frame(group=factor(data$heatmap_class),CL=data$CL)
plot(CL~group,d=CL_data)
CL_res1<-aov(CL~group,d=CL_data)
summary(CL_res1)
# Df Sum Sq Mean Sq F value Pr(>F)
# group 6 1.804 0.30069 18.46 1.56e-15 ***
# Residuals 133 2.166 0.01628
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
0.1%有意の結果を示す。先ほどのS/Gデータよりは変動が大きい傾向にあることが読み取れる。
# Simultaneous Tests for General Linear Hypotheses
#
# Multiple Comparisons of Means: Tukey Contrasts
#
#
# Fit: aov(formula = CL ~ group, data = CL_data)
#
# Linear Hypotheses:
# Estimate Std. Error t value Pr(>|t|)
# 2 - 1 == 0 -1.896e-05 4.915e-02 0.000 1.00000
# 3 - 1 == 0 2.542e-01 4.403e-02 5.773 < 0.001 ***
# 4 - 1 == 0 4.047e-01 4.742e-02 8.535 < 0.001 ***
# 5 - 1 == 0 1.367e-01 3.963e-02 3.450 0.01247 *
# VC - 1 == 0 1.527e-01 4.230e-02 3.610 0.00735 **
# WT - 1 == 0 1.386e-01 5.142e-02 2.695 0.10438
# 3 - 2 == 0 2.542e-01 4.503e-02 5.645 < 0.001 ***
# 4 - 2 == 0 4.048e-01 4.836e-02 8.371 < 0.001 ***
# 5 - 2 == 0 1.367e-01 4.074e-02 3.356 0.01704 *
# VC - 2 == 0 1.527e-01 4.335e-02 3.523 0.00986 **
# WT - 2 == 0 1.386e-01 5.228e-02 2.651 0.11545
# 4 - 3 == 0 1.506e-01 4.314e-02 3.490 0.01088 *
# 5 - 3 == 0 -1.175e-01 3.439e-02 -3.416 0.01388 *
# VC - 3 == 0 -1.015e-01 3.744e-02 -2.710 0.10031
# WT - 3 == 0 -1.156e-01 4.750e-02 -2.434 0.18782
# 5 - 4 == 0 -2.680e-01 3.864e-02 -6.937 < 0.001 ***
# VC - 4 == 0 -2.520e-01 4.138e-02 -6.091 < 0.001 ***
# WT - 4 == 0 -2.662e-01 5.066e-02 -5.254 < 0.001 ***
# VC - 5 == 0 1.599e-02 3.215e-02 0.498 0.99879
# WT - 5 == 0 1.868e-03 4.345e-02 0.043 1.00000
# WT - VC == 0 -1.413e-02 4.590e-02 -0.308 0.99993
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# (Adjusted p values reported -- single-step method)
cld(CL_res2)
# 1 2 3 4 5 VC WT
# "a" "a" "c" "d" "b" "bc" "abc"
まとめ
*作図はRstudioではなく、Excel上で作っています。
*WT:wild typeのデータは省略
・S/G比については、どのグループにおいても大きな違いではないが、グループ3, 5間で5%有意差あり。
・炭化水素とトータルリグニンの比については、若干差が見られました。
また、コントロールとGroup3, 4で上昇がみられ、Group1,2で減少がみられた。
振り返り
ANOVA・多重比較検定用に設計した実験ではないため、結果の無理やり感がありますが...Rstudioによって、両者の検定を実際に行うことができました。
特に、多重比較検定の結果をa, b, c...で表す方法は、遺伝子発現量を比較した論文でよく見かけます。
最後にRstudioを使わずExcelで作図した理由は、使い勝手のいいグラフができなかったからです。案外Excelの方が、きれいなグラフを作れる時があります。