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.
