One-way ANOVA

（カテゴリ変数が連続変数に与える影響その１）

ANOVAとはanalysis of varianceの略で、日本語では分散分析と呼ばれる。

3群以上の平均値を比較するのに使われる。

3群以上のグループをまとめてカテゴリ変数と考えれば、これはカテゴリ変数が連続変数に与える影響の検証と同じである。

One-way ANOVAとは、カテゴリ変数が1つだけのANOVAのことを指す。

要点

例題1

実験

ある植食性昆虫の幼虫の成長が、摂食した植物の種によって異なるか調べるため、幼虫を植物A、植物B、植物Cを与えるグループに分けて飼育し、その後の体重を測定したところ、以下のようなデータが得られた（単位はmg）。

データファイル（csv）

> larval.data

      Weight Plant.Species

1   737.3546        PlantA

2   818.3643        PlantA

3   716.4371        PlantA

4   959.5281        PlantA

5   832.9508        PlantA

6   717.9532        PlantA

7   848.7429        PlantA

8   873.8325        PlantA

9   857.5781        PlantA

10  769.4612        PlantA

11 1502.3562        PlantB

12 1277.9686        PlantB

13 1075.7519        PlantB

14  757.0600        PlantB

15 1424.9862        PlantB

16 1191.0133        PlantB

17 1196.7619        PlantB

18 1388.7672        PlantB

19 1364.2442        PlantB

20 1318.7803        PlantB

21  782.7080        PlantC

22  770.3923        PlantC

23  706.7108        PlantC

24  520.9583        PlantC

25  755.7843        PlantC

26  694.9484        PlantC

27  685.9784        PlantC

28  567.6323        PlantC

29  656.9665        PlantC

30  737.6147        PlantC

> 

作図によってデータの傾向をつかむ（等分散性のチェックと変数変換）

ANOVAでは、誤差に”正規性”と”等分散性”を仮定している。ただ、これらの仮定が満たされていることをどうチェックしたらよいのかは、けっこう微妙な問題である。

実践的には、

のように、見た目での検証で十分のように思うし、そこまで気にする必要はないようにも思う。

なお、ANOVAは前者からの逸脱には頑健だが（正規性が満たされなくても結果を誤りにくい）、後者からの逸脱にはやや脆弱といわれている。そのため、等分散性のチェックは一応行っておいた方が良いかもしれない（見た目で）。特に、重さや体積などの比尺度の場合、平均値が大きくなるほど分散が大きくなることがよくある。このようなデータには、対数をかけて分散をグループ間で同じくらいにするという、対数変換がよく行われる。

対数変換の効果をチェックしたい場合は、以下のようにモデル式の左辺をlogで囲んでしまえばよい（logは自然対数（底がeの対数）を計算する関数）。

par(mfrow=c(1, 2), cex=1.5, las=1, lwd=2, family="sans")

plot(Weight ~ Plant.Species, larval.data)

plot(log(Weight) ~ Plant.Species, larval.data)

左のパネルが生値のデータ、右のパネルが対数変換後のデータである。まず左の生値のデータに注目すると、 B種の植物を摂食した方がほかのAやCに比べて体重が重くなっているようである。また、B種グループではAやCに比べて箱の高さが長いようにみえる。次に右の対数変換後のデータをみると、A、B、Cグループの箱の高さは、左のパネルよりは類似しているようである。

よって、対数変換後のデータを使って解析を進めることにする（必ず対数変換しなければならない、というわけではないことに注意、今回のデータならば生値でも変換しても結果はほとんど変わらなそう）。

One-way ANOVAの実践

ANOVAは線形モデルの枠組みによって扱うことができる。

Rでも、線形モデルを作成してから、そのモデルを使ってANOVAを実行できる。

まず線形モデルを保存し、それをanova関数に入れると、分散分析表が出力される。今回は対数変換後のデータを使って解析するためにモデル式の左辺をlogで囲んでいるが、生値で解析したい場合は単にlm(Weight ~ Plant.Species, larval.data)のように打ち込めばよい。

m <- lm(log(Weight) ~ Plant.Species, larval.data) #線形モデルの作成

anova(m) #One-way ANOVA

> anova(m) #One-way ANOVA

Analysis of Variance Table


Response: log(Weight)

              Df  Sum Sq Mean Sq F value    Pr(>F)    

Plant.Species  2 1.83524 0.91762  41.964 5.192e-09 ***

Residuals     27 0.59041 0.02187                      

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> 

Plant.Speciesの行の一番右にある数字がP値であり、P=0.0000000005である。なので、幼虫の体重は摂食する植物の種によって有意に異なる、と結論できる。

Sum Sq（Sum of Squares: 平方和）とは、線形モデルの枠組みにおける”変動”のことである（変動の分解を思い出そう）。変動をDf（自由度）で割ると、表中のMean Sq（Mean Square: 平均平方）になっているはずである（1.835/2=0.9176,  0.5904/27=0.02187）。この値はそれぞれ、Plant.Speciesの分散と、Residuals（誤差）の分散の不偏推定量である。さらに、Plant.SpeciesのMean SqをResiduals（誤差）のMean Sqで割ると、表中のF valueの値になるはずだ（0.9176/0.02187=41.96）。

このように、ANOVAにおけるF値とは、説明変数の分散と誤差の分散の比のことである。そのため、F値には分子自由度（説明変数の自由度）と分母自由度（誤差の自由度）という2つの自由度が存在することに注意する必要がある。

ANOVA表を論文でレポートしないのであれば、最低限、F値と2つの自由度とP値を報告すればよいだろう。論文では、"F2,27 = 41.96, P < 0.001"のように、Fの下付き文字として2つの自由度を報告することが多い。

ここで、典型的なANOVA表（分散分析表）のサンプルを2つ載せておく。

表1は詳しめに報告するタイプ、表2はシンプルなタイプのANOVA表である。両者では、sum of squaresとmean squaresの有無が異なる。

前述の通り、sum of squaresとは、データの変動（バラツキ）を示している。よって、sum of squaresを見れば、植物の種と残差がそれぞれどの程度のバラツキを持っているかがわかる。例えば、1.84 / (1.84+0.59) = 0.76のように計算すれば、全体の変動のうち、76％を植物の種という要因が説明していることがわかる。また、前述の通り、mean squaresを見れば、植物の種と残差の分散（の不偏推定量）もわかる。

一方、各要因の有意性だけに注目してる場合には、少なくとも自由度・F値・P値があれば十分である。このような場合、自由度の列には分子自由度と分母自由度をコンマでつないで並べることが多い。

表1. 詳しめに報告するタイプのANOVA表。Residuals: 残差、df: 自由度、Sum of squares: 平方和、Mean squares: 平均平方、F: F値、P: P値。
表2. シンプルに報告するタイプのANOVA表。df: 自由度、F: F値、P: P値。

事後検定（post-hoc test）

植物の種によって幼虫の体重が異なることがわかると、次はどの植物で体重が重くなるのかということが知りたくなる。それを知るための検定が事後検定と呼ばれるもので、要はANOVAの後にグループ間の個々のペアで多重比較するというものである。しかし、そもそもやるべきなのかどうかから含めて議論の多い手法でもある。

統計的仮説検定は、P値（帰無仮説が真のときデータが得られる確率）が5％以下であれば帰無仮説を棄却する。なので、検定を100回行えば帰無仮説が真であってもそのうち5回くらいは帰無仮説が棄却される（本当は差がないのに差があると判定してしまう）と期待される。多重比較では検定を多く行うので、本当は差がないのに差があると判定するという誤りを犯しにくくなるように有意水準を調整するということが行われる。

しかしこの調整方法にTukey法だのBonferoni法だのHolm法だのと色々あって大変ややこしい。さらに、データや方法の組み合わせによってはANOVAの結果と事後検定の結果が食い違うこともある。

筆者自身はANOVAの（つまりF検定の）結果の方が本質的であって、事後検定の有無や方法にこだわる必要はないという立場で、あまり背後の理屈に深入りして理解しているわけでもない。検定を何度も行うと「本当は差がないのに差があると判定してしまう」危険性が高まることと、それを調整する方法が存在することを押さえておけば十分だろう。


Rによる多重比較の方法にはいくつかあるが、のちのち応用が効く方法としてここではパッケージ"lsmeans"を用いたものを紹介する。（lsmeansはその名前が誤解を呼びやすいとして開発が終了していて、作者が新しいパッケージである"emmeans"への乗り換えを勧めている。lsmeansは今後はサポートされなくなる可能性大）

そもそもlsmeansは、説明変数の水準ごとの平均を（線形モデル等の）モデルを使って推定するためのパッケージである。

lsmeans関数の第1引数にモデルを入れ、第2引数"specs"に推定したい説明変数の名前を指定する。

library(lsmeans)

lsmeans(m, specs="Plant.Species")

> lsmeans(m, specs="Plant.Species")

 Plant.Species   lsmean         SE df lower.CL upper.CL

 PlantA        6.696894 0.04676206 27 6.600946 6.792841

 PlantB        7.114888 0.04676206 27 7.018940 7.210836

 PlantC        6.526092 0.04676206 27 6.430144 6.622040


Results are given on the log (not the response) scale. 

Confidence level used: 0.95 

> 

すると、各植物ごとにlsmean, SE, df, lower.CL, upper.CLが計算される（SEがすべて同じ値であることを不審に思うかもしれないが、線形モデルでは等分散が仮定されているので、SEがすべて同じなのはむしろ当然）。このうちのlsmeanが、推定された平均値である。

今、この3つの値のうち、どことどこに差があるのかを知りたいわけだ。

これらの値の大小関係を総当たりで検定するには、このlsmeansをpairsで囲めばよい。また、引数adjustによって有意水準の調整法が指定できる。

lsmeans.Plant <- lsmeans(m, specs="Plant.Species")

pairs(lsmeans.Plant, adjust="tukey")

> pairs(lsmeans.Plant, adjust="tukey")

 contrast          estimate         SE df t.ratio p.value

 PlantA - PlantB -0.4179942 0.06613154 27  -6.321  <.0001

 PlantA - PlantC  0.1708016 0.06613154 27   2.583  0.0399

 PlantB - PlantC  0.5887957 0.06613154 27   8.903  <.0001


Results are given on the log (not the response) scale. 

P value adjustment: tukey method for comparing a family of 3 estimates 

> 

以上のような結果が出力されるはずである。上から順に、A種-B種、A種-C種、B種-C種の差が有意かどうかを検定している。今回はTukey法による有意水準の調整を行っているが、adjustの指定を変えれば別の方法で調整が可能である。

また、以下の2つの方法でも同じ結果を出力する。

contrast(lsmeans.Plant, adjust="tukey")


lsmeans(m, specs=pairwise~Plant.Species, adjust="tukey")

さらに、multcomp::cld関数を使えば、有意差に基づく水準のグルーピングも可能である。

cld(lsmeans.Plant, adjust="tukey", sort=FALSE, Letters=letters)

> cld(lsmeans.Plant, adjust="tukey", sort=FALSE, Letters=letters)

 Plant.Species   lsmean         SE df lower.CL upper.CL .group

 PlantA        6.696894 0.04676206 27 6.577881 6.815906  a    

 PlantB        7.114888 0.04676206 27 6.995875 7.233900   b   

 PlantC        6.526092 0.04676206 27 6.407080 6.645104    c  


Results are given on the log (not the response) scale. 

Confidence level used: 0.95 

Conf-level adjustment: sidak method for 3 estimates 

P value adjustment: tukey method for comparing a family of 3 estimates 

significance level used: alpha = 0.05 

> 

一番右のアルファベットが異なるペアでは有意差があることを示す。したがって、植物A、植物B、植物Cの間ではどのペアでも有意に体重が異なると結論できる。

結果の図示（back transform）

t検定のときと同じように生値でグループごとに平均値とSEを計算してグラフを作成しても良いのだが、今回はデータを対数変換して解析しているため、生値の平均±SEのグラフは解析で比較した値とは異なるものになってしまう。

生値の平均値を対数変換しても、対数変換後のデータの平均値とは一致しないからである。これは、実際に対数変換後のデータの平均値を式変形してみればわかる。

  (log(a)+log(b)+log(c))/3 = log((a×b×c)1/3) ≠ log((a+b+c) / 3)

ここで、下線を引いた部分は”(a×b×c)1/3”であり、(a+b+c) / 3とは違うことがわかる。

この”(a×b×c)1/3”は”a, b, cの幾何平均”と呼ばれるものだ（それに対して、(a+b+c) / 3は算術平均と呼ぶ）。

このことは、対数変換を施したANOVAは、生値の算術平均を比較しているのではなく、生値の幾何平均を比較していることを意味している。したがって、生値の算術平均を図示すると、解析では幾何平均を比較しているのに図では算術平均を比較していることになり、図と解析が一致しなくなってしまう。


解析とグラフを一致させるための方法として、lsmeansで推定した平均±SEを指数関数で”back transformする”方法を紹介する。これはモデルが複雑になればなるほど威力を発揮する方法なので、覚えておくととても応用が効く。

まず、全段で保存したlsmeans.Plantをsummaryで囲み、さらに保存する。

summary.lsmeans.Plant <- summary(lsmeans.Plant) 

summary.lsmeans.Plant #中身の確認

> summary.lsmeans.Plant <- summary(lsmeans.Plant)

> summary.lsmeans.Plant

 Plant.Species   lsmean         SE df lower.CL upper.CL

 PlantA        6.696894 0.04676206 27 6.600946 6.792841

 PlantB        7.114888 0.04676206 27 7.018940 7.210836

 PlantC        6.526092 0.04676206 27 6.430144 6.622040


Results are given on the log (not the response) scale. 

Confidence level used: 0.95 

> 

lsmeans関数でつくったオブジェクトは、summaryをかませるとデータフレーム形式に変換される。そのため、上の表の各列を$列名で呼び出せるようになる。

これを使い、上表のlsmean列とSE列、すなわち対数スケールの平均値と対数スケールのSEを呼び出して保存する。

log.m <- summary.lsmeans.Plant$lsmean #対数スケールの平均

log.e <- summary.lsmeans.Plant$SE #対数スケールのSE

これが対数スケールの平均値であることは、以下のように手動で対数値を平均すれば確かめられる。

> with(larval.data, tapply(log(Weight), Plant.Species, mean))

  PlantA   PlantB   PlantC 

6.696894 7.114888 6.526092 

> log.m

[1] 6.696894 7.114888 6.526092

あとはこれらを使って棒グラフを描く。

平均値、平均値-SE、平均値+SEの3つの値をそれぞれ対数の逆関数つまり指数関数（exp()）で実数スケールに直し、barplot2関数に渡してやればよい。

library(gplots) #ライブラリの呼び出し


par(cex=1.5, las=1, lwd=2, family="sans")

barplot2(exp(log.m), plot.ci=TRUE, 

         ci.l=exp(log.m-log.e), ci.u=exp(log.m+log.e), 

         names.arg=c("PlantA", "PlantB", "PlantC"), 

         xlab="Treatment", ylab="Larval weight (mg)")

このような方法を使うことにより、解析と図を完ぺきに一致させることが可能である。論文では、back-transformed least square means ± SE を図示した、という表現がよく使われる。

なお、有意差を示すアルファベットを描きたければ、以下のようにすればよい。

par(cex=1.5, las=1, lwd=2, family="sans")

x <- barplot2(exp(log.m), plot.ci=TRUE, 

              ci.l=exp(log.m-log.e), ci.u=exp(log.m+log.e), 

              names.arg=c("PlantA", "PlantB", "PlantC"), 

              xlab="Treatment", ylab="Larval weight (mg)", ylim=c(0, 1500))

text(x, exp(log.m+log.e)+100, labels=c("a", "b", "c"))

text関数は図に文字を描き加える関数で、第1引数に文字のx座標、第2引数に文字のy座標、labelsに文字を指定する。ここではx座標に棒の位置、y座標にエラーバー上限値+100を指定している。

上図のようなグラフを描くことができる。

'19 2/27