R演習_03_統計処理
2024_応生科実験(植物実験)
Hiroto Yamashita
Last update: 2024-10-29
Licenses
CC BY 4.0
はじめに：統計処理について
実験データの信頼性・再現性を確かめるために統計処理は必要不可欠である．実験データは科学的証拠であり、統計学により実験データを証拠に変えることができる． (ただ、統計学で生命現象をわかった気になるのは注意である)．生物実験では、必然的にデータはばらつきを伴うので (100%正確な測定、観察、標本は存在しない)、標本の反復データを統計解析することで母集団を推定し、科学的な考察を行う．Rは、統計解析向けに開発されたプログラミング言語であるため、統計処理のツールが充実している．ここでは、時間の都合上、細かい統計検定法については解説しないが、本実験結果の検定について、Rでの実装方法について解説する．
本稿で、ggplot2関数を用いるために、library関数を用いて、tidyverseを呼び出しておく．
Step01. 二群間比較 (t検定)
対照区とある処理区の比較などの、一番シンプルな検定である．
まず、デモデータとして、Cont区とTreatA区の実験データからなるデータフレームを作成する．
# データフレーム作成
d <- data.frame(Exp = c("Cont", "Cont", "Cont", "Cont", "Cont", "TreatA", "TreatA",
"TreatA", "TreatA", "TreatA"), Value = c(3.1, 5, 5.2, 4.3, 2.8, 5.9, 6, 6.5,
7.1, 7.3))
# データの中身を確認する
head(d)## Exp Value
## 1 Cont 3.1
## 2 Cont 5.0
## 3 Cont 5.2
## 4 Cont 4.3
## 5 Cont 2.8
## 6 TreatA 5.9データの分布を可視化して、確認する．
# ggplotによるグラフ作成
ggplot(d, aes(x = Exp, y = Value)) + geom_point(size = 3) + theme(text = element_text(size = 20))
この実験群間のデータを眺めてみると、両者の分布の差がありそうに見える．これを統計検定で確かめる．
2群間の検定の場合、t検定が使われる (パラメトリック、ノンパラメトリックの違いもあるが、今回のサンプルサイズではあまり問題とならないためパラメトリック検定として考える．詳細は勉強すべし．) ．
t検定は大別すると、対応のないt検定と対応のあるt検定に分けられ、さらに対応のないt検定はStudentの t検定とWelchのt検定に分けられる．Studentのt検定は等分散している場合、Welchのt検定は等分散していない場合 (不等分散)に適用される．
まず、この実験データの当分散性を決めるために、var.test関数を用いてF検定を行う．データフレームを参照して、Cont (d[1:5,2])とTreatA (d[6:10,2]) の実験データを指定している． F検定の帰無仮説は、対象の2群に差はないこと、つまり有意水準を0.05とすると、p<0.05 (p値：帰無仮説が成立する確率) ならば帰無仮説は棄却され、対象の2群は不等分散となる．
##
## F test to compare two variances
##
## data: d[1:5, 2] and d[6:10, 2]
## F = 2.9824, num df = 4, denom df = 4, p-value = 0.315
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.3105214 28.6446658
## sample estimates:
## ratio of variances
## 2.982412F検定の結果、p値 (p-value) は0.315なので、帰無仮説は棄却されず、2群間に等分散性があることが示唆された．なので、この実験データの検定には、Studentのt検定を適用する．t検定は、t.test関数で実行できる．t検定の帰無仮説は、2群間の平均値に差がないことである．
※var.equalオプション (データ間に等分散性が仮定できるかどうか) が、TRUEの場合Student’s t-test、FALSEの場合Welch’s t-test
##
## Two Sample t-test
##
## data: d[1:5, 2] and d[6:10, 2]
## t = -4.4048, df = 8, p-value = 0.002272
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -3.778344 -1.181656
## sample estimates:
## mean of x mean of y
## 4.08 6.56検定結果は上記のように出力される．p値が0.002272なので、有意水準p<0.05とすると、帰無仮説が棄却され、この2群間には差があることが示される．一般的に、有意水準は5%または1%以下 (p<0.05、p<0.001) とすることにされているがこの限りでない．
したがって、この実験データでは、Cont処理区と比較してTreatA処理区では有意に差があることが示される(厳密には有意水準5%以下で有意差あり)．
Step02. 多群間比較
続いて、多群間比較の場合の統計検定について説明する． 農学、生物実験のデザインでは、多群間比較に該当する場合が多い．
例えば、対照区 (Cont) に対して、様々なタイプの肥料 (FertA, FertB, FertC) を施用して比較する実験、品種間 (CultivarA, CultivarB, CultivarC) での栽培試験を比較する実験などが挙げられる．
前者の場合、Dunnett検定が適用できる．Dunnett検定は、1つの対照区と2つ以上の処理群がある場合、母平均について対照区と処理群の対比較のみを同時に検定するための多重比較法．
後者の場合、Tukey検定が適応でき．Tukey検定は、群間ですべての対比較を同時に検定するための多重比較法．
以下では、RでのDunnett検定とTukey検定の方法について説明する．
Dunnett検定：対照区と処理群の対比較のみを同時に検定するための多重比較法
Dunnett検定用のデモデータをデータフレームで作成する． 対照区 (Cont) に対して、処理A (TreatA)，処理B (TreatB) があり、それぞれ5反復 (n=5) の実験データである．
# データフレーム作成
d.Dun <- data.frame(Exp = c("Cont", "Cont", "Cont", "Cont", "Cont", "TreatA", "TreatA",
"TreatA", "TreatA", "TreatA", "TreatB", "TreatB", "TreatB", "TreatB", "TreatB"),
Value = c(3.1, 5, 5.2, 4.3, 2.8, 5.9, 6, 6.5, 7.1, 7.3, 3, 4.3, 2.2, 2.1, 3.9))
# factor型に変換
d.Dun <- dplyr::mutate(d.Dun, Exp = fct_inorder(Exp))
# d.Dunの中身を確認
head(d.Dun, n = 15)## Exp Value
## 1 Cont 3.1
## 2 Cont 5.0
## 3 Cont 5.2
## 4 Cont 4.3
## 5 Cont 2.8
## 6 TreatA 5.9
## 7 TreatA 6.0
## 8 TreatA 6.5
## 9 TreatA 7.1
## 10 TreatA 7.3
## 11 TreatB 3.0
## 12 TreatB 4.3
## 13 TreatB 2.2
## 14 TreatB 2.1
## 15 TreatB 3.9データの分布を可視化して、確認する．ContとTreatAは差がありそうだが、ContとTreatBは差がなさそう．これをDunnett法で統計検定する．
# ggplotによる作図
ggplot(d.Dun, aes(x = Exp, y = Value)) + geom_point(size = 3) + theme(text = element_text(size = 20))
多群間検定に便利なパッケージであるmulticompをインストールして、呼び出しておく． multcompのインストールが未実施の場合、install.packages関数により、まずはインストールを行う．
library関数を用いて、multcompを呼び出す．
まず、従属変数であるValue列と独立変数であるExp列を指定して、一元配置の分散分析 (one-way ANOVA) をaov関数で行う．帰無仮説はすべての群間において平均値に差がないことである．
以下の結果の通り、p値は0.000209であり、有意水準p<0.001としても、帰無仮説は棄却され、すべての群間の平均値が等しいとはいえない、つまり群間に差があることが示された．
## Df Sum Sq Mean Sq F value Pr(>F)
## Exp 2 31.80 15.902 18.64 0.000209 ***
## Residuals 12 10.24 0.853
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1続いて、どの群間に差があるのかをDunnett検定で確かめる． glht関数にANOVAの結果res1を入力し、linfct=mcp(Exp="Dunnett")で独立変数に対する検定法を指定する．ここでは、Dunnettを指定している．
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Dunnett Contrasts
##
##
## Fit: aov(formula = Value ~ Exp, data = d.Dun)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## TreatA - Cont == 0 2.4800 0.5842 4.245 0.00214 **
## TreatB - Cont == 0 -0.9800 0.5842 -1.677 0.20284
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)各群間比較の組合せごとに統計表が出力される．この検定結果では、ContとTreatAの比較でp値が0.00214の値であり、ContとTreatBの比較でp値が0.20284の値である． 有意水準p<0.01とすると、ContとTreatAの群間には有意差があり、ContとTreatBの群間には有意差が認められないことが示された．
なお、Dunnett検定やt検定での統計検定結果をグラフに反映させる際、グラフ内において、有意差が認められた処理区に「＊」を付与して表現するのが慣例的である．
Tukey検定：群間ですべての対比較を同時に検定するための多重比較法
Tukey検定用のデモデータをデータフレームで作成する． 品種A (CultivarA) 、品種B (CultivarB)，品種C (CultivarC) があり、それぞれ5反復 (n=5) の実験データである．
d.Tuk <- data.frame(Cultivar=c("CultivarA","CultivarA","CultivarA","CultivarA","CultivarA","CultivarB","CultivarB","CultivarB","CultivarB","CultivarB","CultivarC","CultivarC","CultivarC","CultivarC","CultivarC"),
Value=c(99,110,151,150, 80, 170, 149, 152, 173, 147, 82, 93, 90, 100, 80))
#factor型に変換
d.Tuk <- dplyr::mutate(d.Tuk, Cultivar=fct_inorder(Cultivar))
#データの中身を確認
head(d.Tuk, n=15)## Cultivar Value
## 1 CultivarA 99
## 2 CultivarA 110
## 3 CultivarA 151
## 4 CultivarA 150
## 5 CultivarA 80
## 6 CultivarB 170
## 7 CultivarB 149
## 8 CultivarB 152
## 9 CultivarB 173
## 10 CultivarB 147
## 11 CultivarC 82
## 12 CultivarC 93
## 13 CultivarC 90
## 14 CultivarC 100
## 15 CultivarC 80データの分布を可視化して、確認する．CultvarAとCultivarB、CultivarBとCultivarCは差がありそう．
# ggplot2による作図
ggplot(d.Tuk, aes(x = Cultivar, y = Value)) + geom_point(size = 3) + theme(text = element_text(size = 20))
先ほどと同様に、まず分散分析を行う．p値の値は低く、有意水準p<0.001においても有意差あり．
## Df Sum Sq Mean Sq F value Pr(>F)
## Cultivar 2 12076 6038 14.92 0.000557 ***
## Residuals 12 4857 405
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1では、どの群間の組合せに差があるのかをTukey検定で確認する．方法はDunnettの時とほぼ同じで、glht関数を用いて、オプションに'Tukey'を指定して、実行する．
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = Value ~ Cultivar, data = d.Tuk)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## CultivarB - CultivarA == 0 40.20 12.72 3.159 0.0210 *
## CultivarC - CultivarA == 0 -29.00 12.72 -2.279 0.0976 .
## CultivarC - CultivarB == 0 -69.20 12.72 -5.439 <0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)各群間比較の組合せごとに統計表が出力される．この検定結果では、CultivarAとCultivarBの比較でp値が0.0208、CultivarAとCultivarCの比較でp値が0.0976、CultivarBとCultivarCの比較でp値が<0.001の値である．つまり、有意水準p<0.05とすると、CultivarAとCultivarBの群間ないし、CultivarBとCultivarCの群間に有意差が認められ、CultivarAとCultivarCの群間には有意差が認められないことが示された．
Tukeyの統計検定結果をテーブルやグラフに反映させる際、アルファベットを用いて、有意差があるかないかを示す場合が多い．異なるアルファベットの場合、有意差があることを意味し、同一のアルファベットの場合、有意差がないことを表現する．
cld関数を用いることで、自動的にアルファベットのパターンを付与できる．
## CultivarA CultivarB CultivarC
## "b" "a" "b"簡易的に、このアルファベットのパターンと箱ひげ図のグラフを統合させて、可視化してみる．このように示すことで、どの群間が統計的に有意な差が認められているか表現することできる．
実習. データフレームdemoを対象に、グラフの作成および統計検定を実施せよ．PおよびK濃度のデータに対して行うこと．
ここでは、read.csv関数を用いて、csvファイルを読みこむ． デモデータとして、コマツナを様々な栄養欠乏条件で栽培した際の多元素濃度 (イオノーム) データがまとめられたcsvファイル demo.csv を読み込む．
Treat列を横軸として、データフレーム内での順番にするために、factorの水準を設定する．
ggplot2による作図
分散分析
Dunnett検定