[R] help with predict for cr model using rms package

[Frank Harrell]

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.