累積ロジットの汎用Rパッケージ {ordinal}

累積ロジット（cumulative logit model）を使う際に、今ひとつ使い勝手のいいRパッケージがないのが気になっていた。
（cf. 累積ロジット：よい、ふつう、わるい、のような段階的な現象についての推定に用いる。等間隔に近ければ二項分布のGLMで構わないだろうけれど、そうでない場合にはこれが適しているようだ）

例えば、vglmではランダム効果が入れられない。

mixcatでは、ランダム効果が入れられるが、対数尤度までは出力できても、AICなど情報量規準は出してくれない（手計算すれば良い話ではあるが…）。

最近、ordinalというパッケージを見つけた。１つのパッケージでGLM版、GLMM版の両方が含まれているし、使用法もglm()やglmer()と同じようにしてくれていて使いやすい。

ちなみに、GLM版にstepAICは使えたが、dredgeはダメだった。
GLMM版では、モデル選択関数はdrop1が使用できた。

require(ordinal)
example(ordinal) # パッケージで用意されている例を出力

# まずはGLM版
ordinl> ## A simple cumulative link model:
ordinl> fm1 <- clm(rating ~ contact + temp, data=wine)

ordinl> summary(fm1)
formula: rating ~ contact + temp
data:    wine

 link  threshold nobs logLik AIC    niter max.grad cond.H
 logit flexible  72   -86.49 184.98 6(0)  4.01e-12 2.7e+01

Coefficients:
           Estimate Std. Error z value Pr(>|z|)    # 回帰係数部分（この例ではカテゴリカルだが）
contactyes   1.5278     0.4766   3.205  0.00135 **
tempwarm     2.5031     0.5287   4.735 2.19e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Threshold coefficients:           # 各ランクの切片
    Estimate Std. Error z value
1|2  -1.3444     0.5171  -2.600
2|3   1.2508     0.4379   2.857
3|4   3.4669     0.5978   5.800
4|5   5.0064     0.7309   6.850


# 比較対照用にvglmでの結果
require(VGAM)
vm1 <- vglm(rating ~ contact + temp, family=cumulative(parallel=T), data=wine)

Coefficients:
Estimate Std. Error z value
(Intercept):1 -1.3444 0.50850 -2.6438
(Intercept):2 1.2508 0.43908 2.8487
(Intercept):3 3.4669 0.59711 5.8061
(Intercept):4 5.0064 0.72906 6.8669
contactyes -1.5278 0.47362 -3.2258 # 回帰係数部分の正負がclmと逆なことに注意
tempwarm -2.5031 0.53199 -4.7052

# ordinalに戻って、こちらはGLMM版
ordinl> ## A simple cumulative link mixed model:
ordinl> fmm1 <- clmm(rating ~ contact + temp + (1|judge), data=wine) # glmer()同様の構造

ordinl> summary(fmm1)
Cumulative Link Mixed Model fitted with the Laplace approximation

formula: rating ~ contact + temp + (1 | judge)
data:    wine

 link  threshold nobs logLik AIC    niter    max.grad cond.H # ちゃんとAICも算出する
 logit flexible  72   -81.57 177.13 331(996) 1.05e-05 2.8e+01

Random effects:                                     # ランダム効果
 Groups Name        Variance Std.Dev.
 judge  (Intercept) 1.279    1.131
Number of groups:  judge 9

Coefficients:
           Estimate Std. Error z value Pr(>|z|)  
contactyes   1.8349     0.5125   3.580 0.000344 ***
tempwarm     3.0630     0.5954   5.145 2.68e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Threshold coefficients:
    Estimate Std. Error z value
1|2  -1.6237     0.6824  -2.379
2|3   1.5134     0.6038   2.507
3|4   4.2285     0.8090   5.227
4|5   6.0888     0.9725   6.261

（次の記事に続く）