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 at 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.