AICなどの情報量規準を用いた多重比較の方法(Dayton 1998)を実装するRコードの紹介です。以前、共著論文で使用する機会があったのですが、久しぶりに掘り返して、だいぶスッキリしたコードへと改定できたので公開することにしました。手作業だと面倒くさかったのでだいぶ手間が省けるはずです。
非常によく見かける解析方法に、ANOVAの後で多重比較(post-hoc tests)というコンボがあリます。いわゆる多重比較の方法は検定(頻度論)に基づいており、Type I errorの回避に特化しているのが気になっていました。一方で情報量規準に基づくモデル選択による解析を行うケースも多いですが、旧来の多重比較(頻度論)も併用するとダブルスタンダードになってしまいます。実際に査読で突っ込まれたこともあります。また多重比較は正規分布を仮定できるデータにしか適用できないことも制約になります。この点についてはHolmの方法など原理的に分布型に制約されない方法もありますが、組み合わせ数が多くなると検出力が低下する問題もあります。
こういった問題をすっきりと解決できる、情報量規準に基づいたモデル選択による多重比較の方法があります。以前、北大・久保先生のスライドで知った手法ですが( https://kuboweb.github.io/-kubo/ce/2004/index.html )、下記のDayton(1998)が同様の手法(PCIC)を発表しており、根拠論文として参照可能です。特に手法名が提唱されているとは知らずに使っている例も多いかもしれませんが。また、名前が知られていないからこそ、統計ソフトにもRパッケージでも実装されてこなかったのかもしれません。
Dayton CM (1998) Information criteria for the paired-comparisons problem. The American Statistician, 52(2):144-151. DOI: 10.1080/00031305.1998.10480554(ResearchGateからも取得可能)
下記の10数行のコードで実装できました。lm, glm, glm.nbなどのモデル式を下記のpcic関数に食わせれば実行できるはずです。
require(partitions)
# pcic関数の定義、ここから:
pcic <- function(Model, Data=Model$model, Expl=names(Model$model)[2]) {
groupPTs <- apply(setparts(nLevels), 2, function(x) as.numeric(factor(x, levels=unique(x))))
# グルーピングの全パターン、setpartsの標準ではb,a,aなどが出るので正順に並べ替え
AICs <- c(AIC(update(Model,
Data$groupLv <- factor(letters[1:nLevels][groupPTs[,i]][Data[,Expl]], levels=letters[1:nLevels])
AIC(update(Model, reformulate(c(setdiff(names(Data)[-1],Expl),"groupLv"),
comps <- data.frame(AIC=AICs, deltaAIC=AICs-min(AICs), Grouping=groupPTA)
print(paste(c(paste(levels(Data$Species), collapse=", "), "are labeled in this order.",
Expl, "sharing the same alphabet do not differ in the", names(Data)[1]), collapse=" "))
# 関数定義ここまで(これを一度R上で実行するとカスタム関数として利用可能になる)
・lm, glm, glm.nbなどのモデルに適用可能です。それ以外でもModel$modelを実行してモデルに使用しているデータ(y, x)が返ってくるならば使用可能、offset項にも対応。
実行例:
data(iris) # Sepal.Lengthが種によって異なるかについて多重比較を実行
Model <- glm(Sepal.Length ~ Species, data=iris)
require(partitions)