こちらでの計算をMplusで検算する。
Mplus用のデータの書き出し
insomnia<-read.csv("insomnia.csv",header=TRUE)
insomnia<-as.data.frame(insomnia)
head(insomnia)
library(MplusAutomation) variable.names(insomnia) # 変数名を書き出し
prepareMplusData(insomnia, filename="insomnia.dat", overwrite=T)
repolrのコード
library(repolr)
model <- outcome ~ treat + occasion + treat * occasion
fit <- repolr(model, subjects="case", data=insomnia, times=c(1,2),
categories=4)
summary(fit)
Mplusのコード
TITLE: Your title goes here
DATA: FILE = "insomnia.dat";
VARIABLE:
NAMES = case treat occasion outcome count;
USEVARIABLES = case treat occasion outcome count Int;
CATEGORICAL = outcome;
IDVARIABLE = case;
MISSING=.;
DEFINE:
Int = treat * occasion;
ANALYSIS:
ESTIMATOR = MLR;
MODEL:
outcome on treat occasion Int
plot:
type = plot3;
もともとの式に積項があるので、DEFINEで定義しておく必要があるDEFINE:Int = treat * occasion;。ここで定義したIntはUSEVARIABLESに入れておく必要がある。これは忘れがち。
結果
repolrの結果。
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
Mplusの結果。
MODEL RESULTS
Two-Tailed
Estimate S.E. Est./S.E. P-Value
OUTCOME ON
TREAT -0.034 0.240 -0.140 0.889
OCCASION -1.038 0.252 -4.117 0.000
INT -0.708 0.333 -2.123 0.034
Means
COUNT 12.908 0.291 44.390 0.000
Thresholds
OUTCOME$1 -2.267 0.217 -10.449 0.000
OUTCOME$2 -0.951 0.188 -5.061 0.000
OUTCOME$3 0.352 0.180 1.953 0.051
Variances
COUNT 40.421 1.854 21.804 0.000
解釈
推定値の符号が逆である。ちなみに正解はMplusの方だ。repolrでオプションで変えられないか試してみたが、無理だった。
ロバスト推定をした標準誤差だが、閾値の方は許容範囲としても、変数の方がかなり違う。これはまずいな、という感じだ。
順序ロジッスティック回帰分析をした際には、少なくとも、計算にどのソフトウェア・パッケージを使用したかという記述は必要だし、使ってはいけないパッケージもあるようだ。
サプリメント
Mplusでのbrant検定の結果
BRANT WALD TEST FOR PROPORTIONAL ODDS
Degrees of
Chi-Square Freedom P-Value
OUTCOME
Overall test 7.827 6 0.251
TREAT 0.624 2 0.732
OCCASION 0.155 2 0.926
INT 2.346 2 0.309
