情報量規準による多重比較のRコード

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) 

# 組み合わせを表現する関数setpartsを使うのに必要なパッケージ、要インストール

# pcic関数の定義、ここから：
pcic <- function(Model,  Data=Model$model, Expl=names(Model$model)[2]) { 

    nLevels <- length(levels(Data[,Expl])) # グループのリスト
    groupPTs <- apply(setparts(nLevels), 2, function(x) as.numeric(factor(x, levels=unique(x))))
    # グルーピングの全パターン、setpartsの標準ではb,a,aなどが出るので正順に並べ替え

    groupPTA <- apply(matrix(letters[1:nLevels][groupPTs], nrow=nLevels), 2, paste, collapse=" ") 

    # アルファベット変換&結合
    AICs <- c(AIC(update(Model, 

        reformulate(c(setdiff(names(Data)[-1],Expl),1), response = names(Data)[1]))),     # 帰無モデル

    sapply(2:ncol(groupPTs), function(i) { # 帰無モデル以外
    Data$groupLv <- factor(letters[1:nLevels][groupPTs[,i]][Data[,Expl]], levels=letters[1:nLevels]) 

    # グループ分けの変数をDataに追加
    AIC(update(Model, reformulate(c(setdiff(names(Data)[-1],Expl),"groupLv"), 

        response = names(Data)[1]), data=Data)) } ) ) # 各モデル
    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=" "))

    return(comps[order(comps$AIC),]) }
# 関数定義ここまで（これを一度R上で実行するとカスタム関数として利用可能になる）

・lm, glm, glm.nbなどのモデルに適用可能です。それ以外でもModel$modelを実行してモデルに使用しているデータ（y, x）が返ってくるならば使用可能、offset項にも対応。

・注意点として、グルーピングのための説明変数は "~"のすぐ後ろの変数と見なす。異なる時は文字で指定する（この例の場合は"Species"）

・Data: Model$model # モデルに使用しているデータ：y, x。これで大丈夫なはずだが、自前で指定も可能

実行例：
data(iris) # Sepal.Lengthが種によって異なるかについて多重比較を実行
Model <- glm(Sepal.Length ~ Species, data=iris)
require(partitions)

pcic(Model)

[1] "setosa, versicolor, virginica are labeled in this order. Species sharing the same alphabet do not differ in the Sepal.Length"

       AIC  deltaAIC Grouping

5 231.4520   0.00000    a b c

4 265.6359  34.18393    a b b

3 295.6778  64.22581    a a b

1 372.0795 140.62758    a a a

2 373.1310 141.67904    a b a

一番上がベストモデルです。setosa: a, versicolor: b, virginica: c、のように異なるアルファベットで識別されており、種によってSepal.Lengthが異なることが示されました。このアルファベットによる識別は一般的な多重比較でよく見かける表現なので、多くの人が馴染みがあるかと期待します。

結果を表形式で取り出すには、こんな感じでOKです。

write.csv(pcic(Model), "保存場所/ファイル名.csv")

最後に、post-hoc testの用途について。実験系など個々のグループの意味が明確な場合に適した手法だと思いますが、その場合は頻度論の方が適していそうです。ここまで書いておいて何ですが、個人的にはpcicが必要なケースってそんなに多く無いと思っています。生態学系でよく見かける使用例の多くは回帰分析の方が適しているケースだし（説明変数が連続変数になっている等）、2way ANOVAなどの方が適しているケースも多いように感じます。何の仮定も置かずにpost-hoc testをやってみて、その結果を見てから演繹的に考察をするような使い方をしようとしていたら、少し立ち止まって考え直してみませんか。仮説を立ててそれを検証するための解析を行うという設計さえしていれば、より適切な解析方法が見つかるはずです。

cf. lme4は今のところ非対応です。Data=Model@frameでデータを引っ張って来れそうですが、モデルのアップデートのところの改定がすぐにはできそうもなく。