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

[Adam Peer]

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