はじめに
生物系のデータ解析というと、大抵ANOVAとt検定です。1 しかしANOVAは前提条件が多く存在し、データの正規性や等分散性が満たされていない場合は正しい結果が得られないことがあります。 正規性が満たされていない場合に選択肢として上がるのが、ノンパラメトリックな方法です。2郡比較の場合にはMann-WhitneyのU検定が便利に使えます。またANOVAのノンパラメトリックな代替としては、Kruskal-Wallis検定が選択されることが多いです。 ただ、これは教科書でもよく書いてありますが、二元配置デザインの場合にはノンパラメトリックな方法として広く受け入れられる手法がないです。
一方で二元配置デザイン (例えば2ストレインで2条件) でカウントデータをとることは結構ある気がします。 このデータに対して、いくつかの解析方法を試してみました。不勉強なところも多いので、間違いがあれば指摘していただけると幸いです。
データ生成
仮想データはStrain A, B に対して Condition A, Bを用意し、ポアソン分布からサンプルサイズ Nでデータを生成します。 Strain BのCondition Bのλだけ違う値にして、いくつかの解析を行ってみます。
library(tidyverse)
library(MASS)
cmap = c("#333333", "#FF4B00", "#005AFF", "#03AF7A",
"#4DC4FF", "#F6AA00", "#FFF100", "#B5EAD7")
sample_size <- 30
alpha <- 0.05
n_sim <- 1000
# make dataset
set.seed(1)
# number of samples for each group
# あとで再利用できるように関数にしておく
make_dataset <- function(N) {
StrainAConditionA <- rpois(N, 0.2)
StrainAConditionB <- rpois(N, 0.3)
StrainBConditionA <- rpois(N, 0.3)
StrainBConditionB <- rpois(N, 1)
# make tibble for dataset
data <- tibble::tibble(
Strain = factor(rep(c("A", "A", "B", "B"), each = N)),
Condition = factor(rep(c("A", "B", "A", "B"), each = N)),
Value = c(StrainAConditionA, StrainAConditionB, StrainBConditionA, StrainBConditionB)
)
return(data)
}
# box plot using ggplot2
# 簡単に図示してみる。
data <- make_dataset(sample_size)
ggplot2::ggplot(data, aes(x = Strain, y = Value, fill = Condition)) +
geom_boxplot() +
scale_fill_manual(values =cmap)+
theme_classic()
なんとなくBのデータだけ違う印象はあります。 ストレインによって、Conditionに対する応答が異なっているので、Interactionがありそうです。 実際に以下の解析方法でInteractionをどの程度有意として出してくるのかみてみます。
各手法での解析
Two-way ANOVA
まずはTwo-way ANOVAを行いました。ANOVAの場合、データには正規性が仮定されています。 カウントデータの場合、正規性の仮定は満たされないことが多い 2 気がしますが、実際こういう解析は論文でもよくみます。
ちょっと横道にそれますが、どんな統計手法を使うかは正しさの他に周りの環境にもよる気がします。単著ならなんでもいいと思うのですが、、、 個人的には、正しい統計モデルというのが存在しないのと同様に、正しい解析手法も存在しないかとおもいます。どうせ統計の専門家に見せたら生物の論文で使われている手法なんて間違いだらけなんでしょうし。 大切なのは勉強をし続けることと、その時点で自分の考えるベストな手法を選択することだと思っております。 私はまだまだ分からないことばかりですが。
閑話休題、Two-way ANOVAで解析してみます。 データを生成してTwo way ANOVA, interactionのp値を格納 これを1000回繰り返してプロットしてみます。
# Two-way ANOVA
# 1000回繰り返してp値を格納
pvalues <- replicate(n_sim, {
data <- make_dataset(sample_size)
aov_res <- aov(Value ~ Strain * Condition, data = data)
summary(aov_res)[[1]][["Pr(>F)"]][3]
})
# plot p-values
ggplot(tibble::tibble(pvalues = pvalues), aes(x = pvalues)) +
geom_histogram(bins = 30) +
theme_classic()
# count how many p-values are less than 0.05
sum(pvalues < alpha)
conf_interval1 <- prop.test(sum(pvalues < alpha), n_sim)$conf.int
P値の分布はこのようになり、
0.05を下回る率は65%程度でした。
glm with Poisson distribution
今度は一般化線形モデルで解析してみます。今回のようなカウントデータの場合にはデータが非負になります。カウントデータに対して正規分布を仮定して解析を行うと、推定値として負の値がでるおかしな解析をしていることになります。 この辺りはこの本が大変わかりやすくてよかったです。
Rではglmで一般化線形モデルを作成できるので、こちらで二元配置デザインのデータ解析を行います。 まずは分布としてはポアソン分布を使用します。
# glm with Poisson distribution
# 1000回繰り返してp値を格納
pvalues <- replicate(n_sim, {
data <- make_dataset(sample_size)
glm_res <- glm(Value ~ Strain * Condition, data = data, family = poisson(link="log"))
summary(glm_res)$coefficients[, 4][4]
})
# plot p-values
ggplot(tibble::tibble(pvalues = pvalues), aes(x = pvalues)) +
geom_histogram(bins = 30) +
theme_classic()
# count how many p-values are less than 0.05
sum(pvalues < alpha)
conf_interval2 <- prop.test(sum(pvalues < alpha), n_sim)$conf.int
glm with negative binomial distribution
先ほどの解析ではポアソン分布を用いました。ポアソン分布は平均と分散が等しくλとなるような分布ですが、実際に現場で得られるデータは分散が平均よりも大きくなることがしばしばあります。過分散と呼ばれる現象です。過分散に対する手法としてはいくつか知られていますが、ここでは負の二項分布を用いたモデリングを試してみます。 (ただし今回試しているのはポアソン分布からのデータ生成なので、負の二項分布を用いてよりよいモデリングになるかどうかは微妙なところ） MASSパッケージを利用して解析を行います。
# glm with negative binomial distribution
# 1000回繰り返してp値を格納
pvalues <- replicate(n_sim, {
data <- make_dataset(sample_size)
glm_res <- glm.nb(Value ~ Strain * Condition, data = data)
summary(glm_res)$coefficients[, 4][4]
})
# plot p-values
ggplot(tibble::tibble(pvalues = pvalues), aes(x = pvalues)) +
geom_histogram(bins = 30) +
theme_classic()
# count how many p-values are less than 0.05
sum(pvalues < alpha)
conf_interval3 <- prop.test(sum(pvalues < alpha), n_sim)$conf.int
Aligned rank transformation (ART)
最後に、Aligned rank transformation (ART) を使ってみます。
雑な理解ですが、ARTはデータをランクに変換してから解析する手法です。
ランクに変換する前に、align 3 という操作をしているのですが、これは各データからmain effectやinteractionを引いているようです。
その後にデータをランクに変換して、通常のANOVAをかける というのが手順になっているようです。
この手法が結構検出力が高いというので、試してみました。 RにはARToolというARTを実行するためのパッケージがあります。 かなり使い勝手は良く、anova()関数内にARTを指定するだけで良いので、簡単に使えます。 本来はart()で一度実行して、きちんと変換がうまくいっているかをチェックするべきですが、 今回は飛ばします。
# install.packages("ARTool")
library(ARTool)
# 1000回繰り返してp値を格納
pvalues <- replicate(n_sim, {
data <- make_dataset(sample_size)
art_res <- anova(art(Value ~ Strain * Condition, data = data))
art_res$`Pr(>F)`[3]
})
# plot p-values
ggplot(tibble::tibble(pvalues = pvalues), aes(x = pvalues)) +
geom_histogram(bins = 30) +
theme_classic()
# count how many p-values are less than 0.05
sum(pvalues < alpha)
conf_interval4 <- prop.test(sum(pvalues < alpha), n_sim)$conf.int
信頼区間の比較
それぞれの手法で得られた検出力の信頼区間を比較してみます。
# use conf_interval1, conf_interval2, conf_interval3 to compare
# make dataframe from conf_interval1, conf_interval2, conf_interval3
data <- data.frame(
Method = rep(c("Two-way ANOVA",
"glm with poisson",
"glm with negative binomial",
"ART"), each = 2),
Power = c(conf_interval1,
conf_interval2,
conf_interval3,
conf_interval4)
)
# 平均値の計算
means <- aggregate(Power ~ Method, data, mean)
# エラーバーの計算（最大値と最小値の差）
errors <- aggregate(Power ~ Method, data, function(x) max(x) - min(x))
# 平均値とエラーバーを結合
means$Error <- errors$Power / 2
# Means のデータ並び順を変更
# two-way ANOVA, glm with poisson, ARTの順に
means$Method <- factor(means$Method, levels = c("ART",
"glm with negative binomial",
"glm with poisson",
"Two-way ANOVA" ))
p <- ggplot(means, aes(x = Method, y = Power)) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = Power - Error, ymax = Power + Error), width = 0.2) +
labs(title = "Mean statistical powers with Error Bars", x = "Methods", y = "Powers") +
theme_classic()+
theme(text = element_text(size = 18)) +
coord_flip()
p
各手法で1000回中何回interactionが有意になったかと、その回数の95%信頼区間をグラフにしました。 この回数というのは検出力と考えられると思います。 (H0が正しくないときの棄却回数なので と書こうとしてこの辺りはかなり混迷を極めているのを思い出しました。ちょっと自分は理解しきれていませんが、以下の本に詳しかったと思います。)
こうしてみると、検出力からだけでいえば、ARTとTwo way Anova が高い検出力を持っているようです。 他の手法はのちに言及するとして、ARTは検出力が高いようで、使うのも簡単ですしかなり良い印象を持ちました。 ARTのデメリットとしては、モデリングの際に使用されているのがランク変換済みのデータのため、推定されたパラメータの意味は理解しにくいという面があります。例えばポアソン分布でglmを用いた場合には、パラメータ (λ) がどう異なるかなどを議論することができると思いますが、ARTの場合にはその部分の説明が難しそうです。 よく生物で見かける、「有意かどうかを知りたい」4 という状況下であればARTがかなり良い気がしますね。
異なるデータでの再検証
ARTは予想通りだったとして、Two way Anovaが高い検出力を持っている理由は何か、もう少し調べてみたいと思います。 カウントデータがTwo way ANOVAでの解析に向いていない点としては、非負の整数であることが挙げられますが、今回のデータは0が少ないため、正規分布で近似可能かもしれません。 また、サンプルサイズが小さい場合にはglmの推定がうまくいかない可能性もあります。 そこでもっとλを小さくして、0が多く含まれるデータにして、サンプルサイズを増やしてみます。
データはこんな感じでした。
信頼区間はこんな感じになりました。
どうも同じような結果になりました. 色々とパラメータを変えてみたりしましたが、どうも常にTwo way Anovaがglmよりも高い検出力を持っているようです。 もう少しリアルデータっぽい汚い分布にしたらどうなるのかというのも気になりますが、それは今度時間があるときに試してみようかと思います。
まとめ
今回は、カウントデータを使って、いくつかの解析手法の検出力を比較してみました。 試した手法の中では、ARTによるANOVAが結構いい線いっている気がしました。統計モデリングとしての解析ではないのであれば、かなり良い選択肢になると思います。コードを書くのも大変ではないですし、まだ試していませんが対比較も可能なようです。ドキュメントをしっかりよみつつ、今後使ってみようと思います。
もう一つわかったこととして、意外とポアソン分布ではTwo way anovaでも良く解析できるのかもしれません。 (単なる計算ミスや、コードの書き方のミスの可能性は大いにあるとはおもいますが) 他の分布の場合にはどうなるのかというのは興味のあるところですが、また時間のある時に試してみます。
