Agrestiの本で示されている例を実行してみたい。
コード集のみ、サプリメントとしてネット公開されている。
http://users.stat.ufl.edu/~aa/ordinal/R_examples.pdf
サンプルデータはこちらから作成できる。
http://users.stat.ufl.edu/~aa/cda/data.html
この中の、23. Insomnia data set of Table 12.3を利用する。なお、コードの一部を書き換えている。
insomnia<-read.csv("insomnia.csv", header=TRUE)
insomnia<-as.data.frame(insomnia)
head(insomnia)
データは以下のようなもの。
case treat occasion outcome count 1 1 1 0 1 7 2 1 1 1 1 7 3 2 1 0 1 7 4 2 1 1 1 7 5 3 1 0 1 7 6 3 1 1 1 7
repolrパッケージの説明はこちら
library(repolr)
model <- outcome ~ treat + occasion + treat * occasion
fit <- repolr(model, subjects="case", data=insomnia, times=c(1,2),
categories=4)
summary(fit)
結果。
repolr: 2016-02-26 version 3.4
Call:
repolr(formula = model, subjects = "case", data = insomnia,
times = c(1, 2), categories = 4)
Coefficients:
coeff se.robust z.robust p.value
cuts1|2 -2.2671 0.2091 -10.8422 0.0000
cuts2|3 -0.9515 0.1769 -5.3787 0.0000
cuts3|4 0.3517 0.1751 2.0086 0.0446
treat 0.0336 0.2369 0.1418 0.8872
occasion 1.0381 0.1564 6.6375 0.0000
treat:occasion 0.7077 0.2458 2.8792 0.0040
Correlation Structure: independence
Fixed Correlation: 0
Rで順序ロジスティック回帰を行うのは意外に難しいようだ。パッケージによって計算結果が違うと指摘するペーパーもあるようだ。どう違うかも把握しておかないといけないし、そもそもロバスト推定をしないとbrant検定でなんとも言えない状況に追い込まれることも判った。そちらの話はまた後ほど。そういう形で時代は進んでいるようなので、順序ロジスティック回帰はMASSパッケージで計算ができる、という説明は時代錯誤だそうだ。いろいろやってみたところ、計算に最適なパッケージはRでもStataでもなく、Mplusのようだ。今のところの印象だが。
追記：このパッケージは、推定値の符号が逆になるようだ。原因は不明。オプションも試してみたが、今のところわからない。またR version4系ではフリーズを起こす。2016年からアップデートがされていないこともわかった。推定結果をMplusで検算したエントリを入れている。
現在の個人的な見解として、このパッケージは使わない方がよいと思う。
