Dear list,
I am relatively new to ordinal models and have been working through the example
given by Frank Harrell in the predict.lrm {Design} help
All of this makes sense to me, except for the responses, i,e how do i interpret
them? i would be extremely grateful if someone could explain the results?
First i establish the date and model -
> y <- factor(sample(1:3, 400, TRUE), 1:3, c('good','better','best'))
> x1 <- runif(400)
> x2 <- runif(400)
> f <- lrm(y ~ rcs(x1,4)*x2, x=TRUE)
Get 0.95 confidence limits for Prob[better or best
# How do i interpret this on the y scale i.e good,better,best?
>
> L <- predict(f, se.fit=TRUE) #omitted kint= so use 1st intercept
> plogis(with(L, linear.predictors + 1.96*cbind(-se.fit,se.fit)))
> se.fit
> 1 0.6430994 0.8305201
> 2 0.5812662 0.7919122
> 3 0.5692593 0.7976906
> 4 0.5600308 0.7278637
> 5 0.6845250 0.8819143
> 6 0.5518848 0.7228657
> 7 0.5876031 0.7717215
> 8 0.6291766 0.8354423
> 9 0.5839353 0.8333790
> 10 0.5631326 0.8314051
Get Prob(better) than all others -
# Does this mean that for data point 1, y= best as it has the higher
probability?
> predict(f, type="fitted.ind")[1:10,]
y=good y=better y=best
1 0.2517915 0.3469692 0.4012392
2 0.3031733 0.3554471 0.3413796
3 0.3046236 0.3555365 0.3398398
4 0.3514780 0.3546880 0.2938340
5 0.1989827 0.3251784 0.4758390
6 0.3581265 0.3540297 0.2878438
7 0.3130150 0.3559091 0.3310759
8 0.2541324 0.3476007 0.3982669
9 0.2740127 0.3519713 0.3740160
10 0.2839907 0.3535331 0.3624763
Establish data frame to use as newdata
> d <- data.frame(x1=c(.1,.5),x2=c(.5,.15))
Predict newdata - Prob(Y>=j) for new observation
> predict(f, d, type="fitted")
# Does this mean that for data point 1, y= better as it has the higher
probability?
y>=better y>=best
1 0.6800290 0.3239935
2 0.5846743 0.2409657
# Prob(Y=j)
# Again - Does this mean that for data point 1, y= better as it has the higher
probability?
predict(f, d, type="fitted.ind")
y=good y=better y=best
1 0.3199710 0.3560355 0.3239935
2 0.4153257 0.3437086 0.2409657
predict mean(y) using codes 1,2,3
# How do i interpret this on the y scale i.e good,better,best?
> predict(f, d, type='mean', codes=TRUE)
1 2
2.004022 1.825640
Thanks for any advice it is greatly appreciated
John
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
