The help files do not make this clear, but predict and Mean in rms do not support computation of predicted means for the continuation ratio model. The help file for cr.setup shows you how to compute the unconditional probabilities you would need for doing that. Frank
apeer wrote: > > Dear list, > > I'm currently trying to use the rms package to get predicted ordinal > responses from a conditional ratio model. As you will see below, my > model seems to fit well to the data, however, I'm having trouble > getting predicted mean (or fitted) ordinal response values using the > predict function. I have a feeling I'm missing something simple, > however I haven't been able to determine what that is. Thanks in > advance for your help. > > Adam > > > dd <- datadist(all.data2.stand) > options(datadist='dd') > bp.cat2 <- all.data2.stand$bp.cat2 > u <- cr.setup(bp.cat2) > u > > b.mean <-rep(all.data2.stand$b.mean, u$reps) > r.mean <-rep(all.data2.stand$r.mean, u$reps) > mean.ova.energy <- rep(all.data2.stand$mean.ova.energy, u$reps) > > y <- (u$y) # constructed binary response > cohort <- u$cohort > attach(all.data2.stand[u$subs,]) > dd <- datadist(dd, cohort) > > ord.cr <- lrm(y ~ cohort + mean.ova.energy + b.mean + r.mean, x=TRUE, > y=TRUE, na.action=na.delete) > summary(ord.cr) > > > > p.cr <- predict(ord.cr, all.data2.stand, type='mean', codes=TRUE) > pred.mean2 <- data.frame(p.cr) > pred.mean2 > >> ord.cr <- lrm(y ~ cohort + mean.ova.energy + b.mean + r.mean, x=TRUE, >> y=TRUE, na.action=na.delete) >> summary(ord.cr) > Effects Response : y > > Factor Low High Diff. Effect S.E. > mean.ova.energy 0.36902 1.00810 0.63906 -2.732000e+01 11.74 > Odds Ratio 0.36902 1.00810 0.63906 0.000000e+00 NA > b.mean -0.98219 0.18109 1.16330 -6.760000e+00 3.14 > Odds Ratio -0.98219 0.18109 1.16330 0.000000e+00 NA > r.mean -0.50416 0.89758 1.40170 1.175000e+01 4.84 > Odds Ratio -0.50416 0.89758 1.40170 1.270308e+05 NA > cohort - bp.cat2>=2:all 1.00000 2.00000 NA 4.307000e+01 18.37 > Odds Ratio 1.00000 2.00000 NA 5.055545e+18 NA > cohort - bp.cat2>=3:all 1.00000 3.00000 NA 5.538000e+01 23.52 > Odds Ratio 1.00000 3.00000 NA 1.130317e+24 NA > Lower 0.95 Upper 0.95 > -50.32 -4.310000e+00 > 0.00 1.000000e-02 > -12.92 -6.100000e-01 > 0.00 5.400000e-01 > 2.27 2.124000e+01 > 9.66 1.671337e+09 > 7.07 7.907000e+01 > 1171.10 2.182447e+34 > 9.29 1.014700e+02 > 10876.06 1.174706e+44 >> ord.cr > > Logistic Regression Model > > lrm(formula = y ~ cohort + mean.ova.energy + b.mean + r.mean, > na.action = na.delete, x = TRUE, y = TRUE) > > Model Likelihood Discrimination Rank > Discrim. > Ratio Test Indexes Indexes > > Obs 182 LR chi2 174.09 R2 0.953 > C 0.998 > 0 143 d.f. 5 g 33.065 > Dxy 0.996 > 1 39 Pr(> chi2) <0.0001 gr 2.290780e+14 > gamma 0.996 > max |deriv| 6e-07 gp 0.338 > tau-a 0.337 > > Brier 0.013 > > > Coef S.E. Wald Z Pr(>|Z|) > Intercept -20.6064 8.5979 -2.40 0.0165 > cohort=bp.cat2>=2 43.0670 18.3684 2.34 0.0190 > cohort=bp.cat2>=3 55.3845 23.5159 2.36 0.0185 > mean.ova.energy -42.7469 18.3663 -2.33 0.0199 > b.mean -5.8150 2.6984 -2.16 0.0312 > r.mean 8.3840 3.4523 2.43 0.0152 > >> p.cr <- predict(ord.cr, all.data2.stand, type='mean', codes=TRUE) > Error in model.frame.default(Terms, newdata, na.action = na.action, ...) : > variable lengths differ (found for 'mean.ova.energy') > In addition: Warning message: > 'newdata' had 72 rows but variable(s) found have 182 rows >> pred.mean2 <- data.frame(p.cr) > Error in data.frame(p.cr) : object 'p.cr' not found >> pred.mean2 > > ______________________________________________ > 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. > ----- Frank Harrell Department of Biostatistics, Vanderbilt University -- View this message in context: http://r.789695.n4.nabble.com/help-with-predict-for-cr-model-using-rms-package-tp3723475p3723763.html Sent from the R help mailing list archive at Nabble.com. ______________________________________________ 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.