被験者内（参加者内）の分散分析を行うときには、球面性の仮定が満たされている必要があります。また、満たされていない場合は、ε（イプシロン）統計量による自由度の補正を行って、p値を少し大きくすることによって有意になりづらくなるようにしてあげます。
SPSSだと反復測定分散分析を行うと自動的に球面性の仮定、自由度補正を行った共分散分析の結果を行ってくれます。
Rで共分散分析を行うときに、同じように球面性の仮定や自由度補正を行ってくれる方法はないかなぁと思って調べてみたところ、見つけたので忘れないうちにメモをしたいと思います。
理論的なことはさておき、RでSPSSと同じような結果を出力したいというときにこちらの記事が参考になれば幸いです。
Rで球面性の仮定と共分散分析を行えるパッケージと関数
afexパッケージに含まれる、aov_car()という関数を使うことで、球面性の仮定を確認しながら共分散分析を行うことができます。
分散分析だったり、重回帰分析の多重共線性で確認するときによく使われるcar パッケージではないので、注意が必要です。ちなみに、carの意味はCompanion to Applied Regressionらしいので、車(car)とは関係ないっぽいです。
aov_car()で被験者内要因のある共分散分析を行うには、以下の引数を入れていきます。
aov_car(従属変数 ~ 独立変数 + Error(idに関する変数/被験者内要因),data=分析対象のデータフレーム名,
factorize = FALSE,type = 3,return = “afex_aov”)
基本的には、通常のaov関数のように入れていく形です。
後半にあるfacorize=FALSEを入れることで、連続変数の共変量の変数が因子化されるのを防ぎます。
return=”afex_aov”は、afexパッケージの関数と互換性のある分散分析の結果を含むオブジェクトを返してくれます。ほかにも、aov、lm、Anova.mlmなど、いくつかの戻り値が選べますが、それぞれの特徴にあった結果オブジェクトを返してくれるようです。
戻り値として、lmやaovのオブジェクトにすれば、gtsummaryなどで可視化する際に一気に出力できるのかもしれません（まだ試してないですが）。
また、afex以外にも、他には多重比較とか単純主効果検定をやる場合は、emmeansパッケージも必要となります。
実際にaov_car関数を使ってみる
実際にやってみたいとおもいます。
イメージしやすいように、健康得点（10点満点）に影響を与える謎の新薬の効果を確認する研究として、2つの群(drug / placebo）に5名ずつ、測定時期（befor（投与前）/after1（投与直後）/after2（投与から1日後）の３回測定した結果があるとします（2×3の混合計画）。さらに、共変量として、年齢（age）があります。
データとしては、以下のような形です。
Rでは、被験者内要因がある場合は上記のような横長のデータではなく、縦長のlong型データにする必要があります。
縦長データと横長データについてはこちらを確認すると良いかもです。
今回はExcel上のデータを読み込むのではなく、Rのスクリプト上で入れていきたいと思います。
Rを立ち上げて、以下のスクリプトを入れることで、上記のデータのlong型データを作ることができます。
id<-rep(c(1:10),3)
btwn<-rep(c(“drug”,”placebo”),each=5,time=3) #投薬条件（被験者間要因）
wtin<-rep(c(“before”,”after1″,”after2″),each=10)
age<-rep(c(45,44,40,48,49,41,46,46,43,38),time=3)
score<-c(3,3,2,4,4,3,4,2,3,3,
5,7,5,6,5,3,5,3,4,3,
4,5,4,4,3,5,3,2,3,4)
dat<-data.frame(id,btwn,wtin,age,score)
読み込むとこんな感じになります。
このデータに対して、2×3の混合計画の年齢を共変量として、分散分析を行うときは以下のスクリプトを実行します。
library(afex)
library(emmeans)
options(contrasts = c(“contr.sum”, “contr.poly”)) #type3平方和の計算の前に入れておく呪文
#共分散分析の実行
mixAncova<-aov_car(score~btwn*wtin + age + Error(id/wtin),data=dat,
factorize = FALSE,type = 3,return = “afex_aov”)
mixAncova
共変量となる変数は、要因と同じように計算式に組み込む形でOKです。
mixAncovaの出力を見ると以下のようになります。
Sphericity correction method: GG と書かれているように、Greenhouse-Geisserの方法で自由度の調整をした結果が出力されています（SPSSともすべてのF値、P値が一致していました）。
今回の都合の良いダミーデータの場合は交互作用（btwn:wtin）が有意だったので、単純主効果検定を行いたい場合は、以下のように行います。
mfi_Ancova<-emmeans(mixAncova,~wtin|btwn, cov.reduce=FALSE)
# 各参加者内水準における参加者間の比較
contrast(mfi_Ancova, method = “pairwise”,by=”wtin”)
#各参加者間水準ごとに、参加者内水準の比較（bonferroni補正）
pairs(mfi_Ancova, adjust = “bonferroni”)
※mfi_Ancovaの変数を作成するときに、cov.reduce=FALSEとすると、spssの結果と一致しました。
１つめのconrstnの部分で、各測定時期間で被験者間要因の比較をして、2つめのpairsで、投与群ごとに被験者内要因の水準間比較を行っています。
まとめ
参加者内要因のある共分散分析を行うときに、球面性の仮定についても考慮しながら行うときには、afexパッケージに含まれるaov_car関数を使えばうまくできそうです。
なお、どういう場合に共分散分析を行ってもよいのか、球面性の仮定で使うべき補正方法などについては触れていませんので、専門書籍などを参考にすると良いかもしれません。
SPSSの結果をRで再現したいなぁというときなどに参考になれば幸いです。
