inst/tests/orm-example.r

In rms: Regression Modeling Strategies

## See http://stats.stackexchange.com/questions/65548/which-model-should-i-use-to-fit-my-data-ordinal-and-non-ordinal-not-normal-an
pred_1 = rep(c(10,20,50,100),30)
pred_2 = rep(c('a','b','c'),40)

resp = c(0.08666667, 0.04000000, 0.13333333, 0.04666667, 0.50000000, 0.04000000, 0.02666667, 0.24666667, 0.15333333, 0.04000000, 0.06666667, 0.06666667, 0.03333333,
    0.04000000, 0.26000000, 0.04000000, 0.04000000, 1.00000000, 0.28666667, 0.03333333, 0.06666667, 0.15333333, 0.06666667, 0.28000000, 0.35333333, 0.06000000,
    0.06000000, 0.05333333, 0.96666667, 0.06666667, 0.03333333, 0.22000000, 0.04666667, 0.04666667, 0.05333333, 0.05333333, 0.05333333, 0.08000000, 0.48666667,
    0.08666667, 0.02666667, 0.21333333, 0.45333333, 0.04666667, 0.36000000, 0.06666667, 0.04000000, 0.06000000, 0.07333333, 0.06000000, 0.04000000, 0.04666667,
    0.30000000, 0.08666667, 0.07333333, 0.06666667, 0.29333333, 0.36000000, 0.17333333, 0.04000000, 0.09333333, 0.11333333, 0.03333333, 0.08000000, 0.27333333,
    0.08666667, 0.03333333, 0.04000000, 0.02666667, 0.07333333, 0.07333333, 0.02000000, 0.02666667, 0.08000000, 0.07333333, 0.02666667, 0.06666667, 0.07333333,
    0.95333333, 0.05333333, 0.04000000, 0.11333333, 0.04000000, 0.07333333, 0.06666667, 0.05333333, 0.04000000, 0.04000000, 0.06000000, 0.12666667, 0.04666667,
    0.04000000, 0.21333333, 0.05333333, 0.97333333, 0.11333333, 0.02666667, 0.04000000, 0.03333333, 0.37333333, 0.25333333, 0.06000000, 0.06000000, 0.06000000,
    0.04666667, 0.26666667, 0.98000000, 0.02000000, 0.26000000, 0.06000000, 0.05333333, 0.28000000, 0.99333333, 0.04666667, 0.02666667, 0.04000000, 0.12666667,
    0.04666667, 0.18000000, 0.03333333)

require(rms)
row <- 0
png('/tmp/lookdist.png')
for(gvar in list(pred_1, pred_2)) {
  row <- row + 1; col <- 0
  for(fun in list(qlogis, qnorm, function(y) -log(-log(y)))) {
    col <- col + 1
    cat(row, col, '\n')
    print(Ecdf(~resp, groups=gvar, fun=fun,
               main=paste(c('pred_1','pred_2')[row],
                 c('logit','probit','log-log')[col])),
          split=c(col,row,3,2), more=row < 2 | col < 3)
  }
}
dev.off()
f <- orm(resp ~ pred_1 + pred_2)
f
anova(f)
dd <- datadist(pred_1, pred_2); options(datadist='dd')
bar <- Mean(f)
png('/tmp/Predict.png')
plot(Predict(f, fun=bar), ylab='Predicted Mean')
dev.off()

png('/tmp/or.png')
plot(summary(f), log=TRUE)
dev.off()