[R] predict.lrm ( Design package)

Chris Mcowen chrismcowen at me.com
Mon Sep 20 16:13:41 CEST 2010 

A few comments; sorry I don't have time for any more.

- Combining categories is almost always a bad idea
- It can be harder to discriminate more categories but that's only because the task is more difficult
- Split-sample validation is not reliable unless you have say 10,000 samples to begin with; otherwise 100 fold repeats of 10-fold cross-validation, or (faster) 200 or so bootstraps is better.  Be sure to re-do all modeling decisions from scratch for each re-sample.  If you did no statistical variable selection this may not be an issue.
- You don't need ".lrm" after predict
- Not sure why your variable names are all upper case; harder to read this way

Good luck
Frank

Frank E Harrell Jr   Professor and Chairman        School of Medicine
                    Department of Biostatistics   Vanderbilt University

On Mon, 20 Sep 2010, Chris Mcowen wrote:

> Dear Professor Harell
> I am familier with binary models, however i am now trying to get predictions from a ordinal model and have a
> question. 
> I have a data set made up of 12 categorical predictors, the response variable is classed as 1,2,3,4,5,6, this
> relates to threat level of the species ( on the IUCN rating). 
> Previously i have combined levels 1 and 2 to form = non threatened and then combined 3-6 to form threatened, and run
> a binary model. I have tested the performance of this based on the brier score (0.198) and the AUC or C value
> (0.75). I also partitioned the data set into training and test data and used the predict function to get a predicted
> probability for the newdata. When visualising the results with a cutoff value calculated with epi, roughly 75% of
> the time the prediction was correct. 
> Now i am interested in predicting the threat level of a species not purely as threatened or not but to specific IUCN
> levels. I have used the predict.lrm function (predict.lrm(model1, type="fitted.ind"))  to generate probabilities for
> each level,  and also (predict(model1, traist, type="fitted")) see below. 
> When i call the model the Brier score is  0.201  and C value 0.677. However  when i inspect the output and relate it
> to the corresponding species ( for which i know the true IUCN rating) the model performs very badly, only getting
> 43% correct. Interestingly i have noticed the probabilities are always highest for levels 1 and 6, on no occasion do
> levels 2,3,4 or 5 have high probabilities?
> I am unsure if this is just because the model can not discriminate between the various levels due to insufficient
> data? Or if i am doing something wrong, if this is the case i would greatly appreciate any advice or suggestion of a
> different method. 
> Thanks in advance,
> Chris 
> model1 <- lrm(EXTINCTION~BS*FR+WO+LIF+REG+ALT+BIO+HAB+PD+SEA, x=TRUE,y=TRUE)
> predict.lrm(model1, type="fitted.ind")
>     EXTINCTION=1 EXTINCTION=2 EXTINCTION=3 EXTINCTION=4 EXTINCTION=5
> 1     0.19748393   0.05895033   0.12811186  0.086140778  0.068137104
> 2     0.27790178   0.07247496   0.14384976  0.087487677  0.064584865
> 3     0.24605628   0.06777939   0.13931242  0.087996215  0.066625303
> 4     0.24605628   0.06777939   0.13931242  0.087996215  0.066625303
> 5     0.24605628   0.06777939   0.13931242  0.087996215  0.066625303
> 6     0.24605628   0.06777939   0.13931242  0.087996215  0.066625303
> 7     0.13928899   0.04558050   0.10636220  0.077770389  0.065500459
> 8     0.24605628   0.06777939   0.13931242  0.087996215  0.066625303
> 9     0.24605628   0.06777939   0.13931242  0.087996215  0.066625303
> 10    0.33865077   0.07915126   0.14744522  0.083923247  0.059387585
>     EXTINCTION=6
> 1     0.46117600
> 2     0.35370096
> 3     0.39223038
> 4     0.39223038
> 5     0.39223038
> 6     0.39223038
> 7     0.56549746
> 8     0.39223038
> 9     0.39223038
> predict(model1, traist, type="fitted")
>           y>=2       y>=3       y>=4       y>=5       y>=6
> 1   0.80251607 0.74356575 0.61545388 0.52931311 0.46117600
> 2   0.72209822 0.64962327 0.50577351 0.41828583 0.35370096
> 3   0.75394372 0.68616432 0.54685190 0.45885569 0.39223038
> 4   0.75394372 0.68616432 0.54685190 0.45885569 0.39223038
> 5   0.75394372 0.68616432 0.54685190 0.45885569 0.39223038
> 6   0.75394372 0.68616432 0.54685190 0.45885569 0.39223038
> 7   0.86071101 0.81513051 0.70876831 0.63099792 0.56549746
> 8   0.75394372 0.68616432 0.54685190 0.45885569 0.39223038
> 9   0.75394372 0.68616432 0.54685190 0.45885569 0.39223038
> 10  0.66134923 0.58219797 0.43475276 0.35082951 0.29144192
On 20 Sep 2010, at 14:50, Frank Harrell wrote:


You sent a private note about this which I just took the time to answer. 
Please send only one note, and please post my reply to you to r-help.
Frank


-----
Frank Harrell
Department of Biostatistics, Vanderbilt University
-- 
View this message in context: http://r.789695.n4.nabble.com/predict-lrm-Design-package-tp2546894p2547011.html
Sent from the R help mailing list archive at Nabble.com.

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

More information about the R-help mailing list