Repeated measures ANOVA

同じ個体から何回もデータを取るなどの、同一対象からの複数回サンプリングは、研究の現場ではありふれているだろう。

にもかかわらず、このようなデザインは統計解析にとってかなりの曲者である。

なぜなら、同一対象からの複数回測定は、データが互いに独立であるという統計解析の大前提に抵触するからである。

Repeated measures ANOVA（反復測定分散分析）は、同一対象から複数回サンプリングされたデータという問題に対処できる手法の一つである。

※本稿では、同一対象から複数回サンプリングされたデータを解析するANOVAを、十把一絡げにrepeated measures ANOVAと呼びます。

※※Rの使い方と、ANOVAの理屈にある程度慣れてから挑戦することをオススメします。

要点

• 疑似反復（pseudoreplication）は統計解析にとって危険な存在

• 疑似反復の存在を見抜く基本キーワードは、「同じ個体」「同じ場所」「同じ時間」

• Repeated measures ANOVA（反復測定ANOVA）は、ランダム効果によって疑似反復に対処する

• Repeated measures ANOVAを実行する方法にはいろいろある（aov関数, nlmeパッケージ, lme4パッケージ+carパッケージまたはlmerTestパッケージなどなど）

• グループ間変数とグループ内変数を見分けよう

• Repeated measures ANOVAは、LMM（線形混合モデル）の入門編

例題１：疑似反復の問題点とランダム効果による解決

操作実験１

農薬散布が作物の収量に与える影響を調べるため、4つのプロットに各々10株の作物を植え、4つのうち2つのプロットに農薬を散布した。その後、プロットごとの作物の収量を測定した。その結果、以下のようなデータが得られた。（C: 無処理区、P: 農薬処理区を指す）

プロットの存在を無視した解析：図示

プロットの存在を無視した解析：one-way ANOVA

まず、作物の株を繰り返し（replicate）として、one-way ANOVAを行ってみる。

P値が0.001より（かなり）低く、農薬散布の影響は有意だという結果となる。

…ということで、農薬散布は作物の収量を増加させるのだ、と結論してよいのか？というのが、ここでの問題である。結論から先に言えば、これが疑似反復の落とし穴であり、このデータから「農薬散布が作物の収量を増加させる」と結論づけるのは危険である。

プロットの存在を考慮した解析：図示

同一対象からの複数回サンプリングなどによって、本来の解析のユニットよりも下の階層を繰り返し（replication）として扱ってしまうことを、疑似反復（pseudoreplication）と呼ぶ。今回の例で言えば、「本来なら解析のユニットをプロットにすべきところを、植物個体にしてしまうこと」である。同一のプロットから得られた複数のデータは、別のプロットから得られたデータ同士に比べて、どうしても互いに値が似通ってしまうと懸念される（これを、”データが相関している”と呼ぶ）。あらかじめ似通っていると想定されるデータをいくつも取るのは、データ数を水増ししたことになってしまう。

適切な解析のユニットを選ぶコツの一つは、何に対して実験処理を行ったか？である。農薬散布処理は植物個体ではなくプロットごとに行われている。こういう場合には、農薬処理の効果を検証するための群内変動は、プロット間の変動であり、植物個体間の変動ではないのである。

プロットの存在を考慮した解析：repeated measures ANOVA

プロットを解析のユニットにする方法として、プロットの効果を”ランダム効果”としてANOVAのモデルに組み込むという方法がある。ランダム効果を含めると、「データが互いに独立である」という統計解析の大前提を少し緩めることができる。平たく言うと、ランダム効果の投入により、同一プロットのデータはどれも似ている、という影響を取り除くことができる。本稿では、ランダム効果を組み込むことによって、同一対象からの複数回サンプリングに対処するANOVAを全て"repeated measures ANOVA"と呼ぶ。このように何でもかんでもrepeated measures ANOVAと呼ぶのは、たぶん教科書やWeb上の資料ではあまりやらないやり方なので、名前にこだわる人は色々調べてみてください。

Repeated measures ANOVAを行うためには、ランダム効果が含まれる線形モデルを作らないといけないが、この作り方には大きく分けて最小二乗法（ランダム効果なしの線形モデルで普通使われる方法）と最尤法の二つがある。ランダム効果入りの線形モデルはランダム効果なしのものに比べて作るのが難しいので、最尤法を使った方が結果が安定するようである。

最小二乗法アプローチ（関数aov）

最小二乗法でrepeated measures ANOVAを行うためには、関数aovを使う。関数aovでは、モデル式中でError()という記法を用いることができる。この中で指定した変数がランダム効果となる。プロットをランダム効果に含めるためには、以下のようにすればよい。

##Plotをランダム効果に含めることで、疑似反復に対応したANOVA

summary(aov(Yield~Treatment+Error(Plot), yieldData))

###Error: PlotにおけるResidualsのDfに注意

これまで見てきたANOVAの結果と異なり、Errorという”階層”を示す結果が表示されている。Treatmentの行はError: Plotの階層だけに現れており、これはF=5.78, P=0.14となっている。よって、農薬処理が収量に与える効果は有意ではないと結論づけられる。

注意すべきはError: Plotの階層におけるResidualsの自由度である。プロットの存在を無視したANOVAでは38（=40(植物個体の数)-1(処理の自由度)-1）だったものが、2と、大幅に減少している。これは、農薬処理を検証するための群内変動として、プロット間変動が選ばれたことを示している。プロットは合計4つだけなので、Residualsの自由度は4(プロットの数)-1(処理の自由度)-1=2となるのである。

最尤法アプローチ（パッケージlme4）

次に、より現代的な方法である最尤法を用いて解析する。最尤法を用いてランダム効果入りモデルを作るパッケージには、主にnlmeとlme4があるが、ここでは後者のlme4パッケージを用いた方法を紹介する。

lme4パッケージには、ランダム効果入りの線形モデル（正規分布仮定）を作るlmer関数が用意されており、これを使ってrepeated measures ANOVAを実行することができる。lmer関数の使い方はlm関数に準じたものだが、モデル式のランダム効果の記法がlmer関数に独自のものである。今回の実験の場合、lmer関数でプロットをランダム効果に含めるためには、モデル式の中に"(1|Plot)"という項を加える（なんでこんなにややこしい書き方をするのかと思うかもしれないが、まずは丸覚えで良い）。

また、ポイントは、carパッケージのAnova関数の中で、test="F"と指定することである。これを指定すると、pbkrtestというパッケージを経由して、適切なF検定を自動で行ってくれる。

##Plotをランダム効果に含めることで、疑似反復に対応したANOVA

library(lme4) #lmer関数を含むパッケージ

library(car) #Anova関数を含むパッケージ

lmer.m <- lmer(Yield~Treatment+(1|Plot), yieldData)

Anova(lmer.m, test="F")

aov関数のアウトプットと比較すると、F値、P値ともに同一である。Df.resは、TreatmentのResidualsの自由度を表しており、これも2であるから、aov関数の結果と同一である。

このくらいのシンプルなデザイン、データであれば、最小二乗法と最尤法で結果が異なることはほとんどない。ただし、データによっては両者には違いが出てくる。また、aov関数は処理間でデータ数が異なる（unbalanced）場合は使えない。今日では基本的に最尤法を用いる方が推奨される。

Least square meansを用いた、repeated measures ANOVAの結果の図示

ANCOVAの図示で活躍したleast square meansは、repeated measures ANOVAでも大いに有用である。

Repeated measures ANOVAのモデルにはランダム効果が含まれているので、単に植物個体ベースでSEを計算しても、正しいSEにならない。lsmeans関数では、モデルにランダム効果が含まれていることを考慮のうえ、平均値だけでなく正しいSEも算出が可能である。

関数の使い方は、ANOVAやANCOVAのときとまったく同じである。

library(lsmeans)

(lsm.lmer.m <- lsmeans(lmer.m, specs="Treatment"))

このようにして計算されたleast square meanとSEを用いて、棒グラフを作成することができる。

まず、summary関数をかませてlsmeansをデータフレーム型に変換し、平均値とエラーバーの区間を計算する。続いて、計算した平均値とエラーバーの区間を用いて、棒グラフとエラーバーを描く。

例題２：repeated measures ANOVAによる経時的測定データの解析

操作実験２

ある昆虫の成長が、食べる餌の質によってどのように変化するか調べるため、8個体を4個体ずつ餌Aと餌Bを与えるグループに分けて5日間飼育し、1日目（給餌前）、3日目、5日目の体重を測定した。ただし、実験の都合上、いくつかの欠損値が生じた。その結果、以下のようなデータが得られた。

このデザインは、一見典型的な2×3の二要因デザインに見える（各餌と各測定日の全ての組み合わせが存在する）。しかし、ここでも同一対象、つまり昆虫の各個体に対して、複数回の測定が行われている。したがって、同じ個体からのデータはどれも似ている、という影響を、餌→体重、測定日→体重の影響から取り除く必要がある。

繰り返し数の確認

まず、要因の各組み合わせにおける繰り返し数を確認しておこう。

要因FoodとDayの各組み合わせにおいて、繰り返し数が異なっている、unbalancedなデザインであることが分かる。

普通の箱ひげ図

次に、とりあえず個体はムシして、要因DayとFoodの組み合わせによって体重がどう変化しているか、普通の箱ひげ図でチェックしてみる。ついでに、ln変換した方が良いかも一応みておく。

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

boxplot(Mass ~ Day + Food, growthData, ylab="Mass")

boxplot(log(Mass) ~ Day + Food, growthData, ylab="log(Mass)")

普通にtwo-way ANOVAしてみる

比較のため、まずは個体を考慮することなく、普通にtwo-way ANOVAを行ってみる。農薬散布の例でも説明したように、同一個体からのデータは”疑似反復”のはずなので、本来はおかしい方法である。

なお、変数Dayは数値変数になっているはずなので、解析の前にカテゴリカル変数に変換しておく。

すべての項が有意になっている。注意すべきは、分母自由度（ResidualsのDfのこと）がすべての項で15になっていることである。分母自由度の問題はあとでもう一度解説する。

各個体のデータを対応させた図の作成

さて、ここまでは各個体の違いは考慮せずに解析を進めてきた。しかし、実際はこの実験は同一の幼虫を経時的に複数回測定している。そのため、「ある個体は別の個体より常に重い」のようなことが起こっているかもしれない。

以下のようなスクリプトで、幼虫個体別の体重増加の過程を図示することができる（詳細ははぶきます）。

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

library(dplyr)

library(tidyr)

growthData %>%

  select(-Food) %>%

  mutate(Mass=log(Mass)) %>%

  spread(key=ID, value=Mass) %>%

  select(-Day) %>%

  as.matrix(.) %>%

  matplot(., type="b", col=rep(c(1, 2), each=4), lty=1, xlab="Day", ylab="ln(Mass)")

Repeated measures ANOVAの実践

というわけで、次は個体を考慮した（正しい）解析であるrepeated measures ANOVAを行う。今回、unbalancedなデザインなのでaov関数による最小二乗法は使えない。よって、lme4::lmer関数を用いて解析を行う。グループ間変数とグループ内変数が混在していても、指定の仕方は農薬実験のときと全く同じであり、グループ変数（ID）について、"+(1|ID)"と指定することで解析が行われる。

※注意！！：この手のデザインの場合、特に時間に関する説明変数が数値変数になっていることが多いが、これはカテゴリカル変数にしておかないと結果がおかしくなる。事前に確認が必要。

library(lme4)

library(car)

growth.lmer.m <- lmer(log(Mass) ~ Food*Day + (1|ID), growthData)

Anova(growth.lmer.m, test="F")

どの項も有意だという結果は変わらないのだが、統計量FやP値には違いがみられる。特に、変数Foodの有意性がだいぶ下がっている。注目すべきはDf.res（residualの自由度、分母自由度とも）で、なんと整数ではなくなっている。これは、Kenward-Roger法による分母自由度の近似値で、unbalancedなデザインの時には整数ではなくなる。

また、two-way ANOVAの出力ではFood, Day, Food:Dayの分母自由度はすべて15だったが、今回は順に5.76, 9.56, 9.63となっている。仮に、欠損値がなかったと仮定して計算してみると、幼虫の数は8、餌処理の水準数は2、日の水準数は3なので、データの数は（欠損がなければ）8×3=24になる。すると、グループ間変数の繰り返し数は8なので、8-1(餌処理の自由度)-1=6(グループ間階層の分母自由度)と計算されるはずである。さらに、グループ内階層の分母自由度は、24-1(餌処理の自由度)-6(グループ間階層の分母自由度)-2(日の自由度)-2(餌×日の自由度)-1=12となる。実際には欠損があったのでこの自由度よりやや小さくなっているが、グループ間階層＜グループ内階層という大小関係は変わらない。

このように、グループ間階層の分母自由度と、グループ内階層の分母自由度は、それぞれグループ間階層の繰り返し数とグループ内階層の繰り返し数を反映しているものだと考えることができる。よって、repeated measures ANOVAを行ったときには、分母自由度をチェックして、自分が考えている繰り返し数を自由度が反映していそうか確かめるのが有用だろう。

Least square meansを用いた事後検定

さて、今回のモデルの場合でも、lsmeansパッケージを用いた更なる解析が可能であるが、図示は農薬実験のやつと同じようにやればよいので、ここでは経時的な測定データにありがちな事後検定を行ってみる。

普通、この手の実験では、「同じ日の中で、処理間で体重が異なるか」は知りたいが、「異なる日の間で体重が異なるか」についてはそこまで注意が払われないことが多いだろう。幼虫がある程度成長するのは当たり前だからである。そこで、上記のlmerで作ったモデルを使って、「同じ日の中で、処理間で体重が異なるか」を検証してみよう。

lsmeansは引数specsに変数を取ることで、その変数の組み合わせについて平均値の推定値を計算してくれる。

これを、関数pairsで囲むと、すべての組み合わせで事後検定を行ってくれる。しかし、この方法の場合、興味のない比較もすべて行われてしまう。

そこで、Foodに注目するために、Dayで場合分けしたい。その場合は、引数byが使える。

これをpairsで囲むと、byで指定した水準ごとに、繰り返し事後検定をやってくれる。

Day1では餌間で差がないが、Day3、Day5では餌間で差があることがわかる。

ただし、この方法では、byで指定した水準内ではP値の調整が行われるが、水準間では調整が行われていない。この3つの比較全体でP値の調整を行いたい場合には、上をさらにrbindで囲む。

ちょっとだけP値が大きくなっていることがわかるはずである。

Repeated measures ANOVAはmixed modelingの入門編

さて、ここまで、repeated measures ANOVAによる反復測定デザインの解析方法を解説してきた。Repeated measures ANOVAは、ランダム効果を導入することで、pseudoreplicationの問題を回避する解析法であった。このように、ランダム効果が導入されたモデルを”mixed model”と呼ぶ。（なぜmixedなのかと言えば、線形回帰やANOVAなどに含まれる普通の説明変数を”固定効果”と呼ぶが、mixed modelは固定効果とランダム効果が混在しているからである。）

Mixed modelingはなかなか理解するのが難しいが、上記の2つは非常に典型的なmixed modelingの例であり、入門編として有用だと思う。グループ変数・グループ間変数・グループ内変数を認識できるようになってくれば、もう少し複雑なmixed modelingも使いこなせるようになるだろう。

ついでに、ランダム効果とcorrelation structureの記事も読むと、より理解が深まると思うので、ぜひ続けて読んでみてください。

'21 1/26