Cross validation and ordinal logistic regression

Question

I am trying to understand cross-validation for ordinal logistic regression. The aim of the game is to validate the model used in an analysis...

I first construct a toy data set:

set.seed(1)
N <- 10000
# predictors
x1 <- runif(N)
x2 <- runif(N)
x3 <- runif(N)

# coeffs in the model
a <- c(-2,-1)
x <- -x1+2*x2+x3

# P( y ≤ i ) is given by logit^{-1} ( a[i]+x )
p <- outer(a,x, function(a,x) 1/(1+exp(-a-x)) )

# computing the probabilities of each category
q <- 1 - p[2,]
p[2,] <- p[2,] - p[1,];
p <- rbind(p,q);

# outcome
y <- ordered( apply( p, 2, function(p) which(rmultinom(1,1,p)>0) ) ) 

Now, I fit the model it using lrm in the package rms.

require("rms")
fit <- lrm(y~x1+x2+x3, x=TRUE,y=TRUE)


> fit

Logistic Regression Model

lrm(formula = y ~ x1 + x2 + x3, x = TRUE, y = TRUE)

                      Model Likelihood     Discrimination    Rank Discrim.
                         Ratio Test            Indexes          Indexes   
Obs         10000    LR chi2    1165.46    R2       0.126    C       0.664
 1           2837    d.f.             3    g        0.779    Dxy     0.328
 2           2126    Pr(> chi2) <0.0001    gr       2.178    gamma   0.329
 3           5037                          gp       0.147    tau-a   0.203
max |deriv| 4e-10                          Brier    0.187                 

     Coef    S.E.   Wald Z Pr(>|Z|)
y>=2  2.1048 0.0656  32.06 <0.0001
y>=3  1.0997 0.0630  17.45 <0.0001
x1    0.8157 0.0675  12.09 <0.0001
x2   -1.9790 0.0701 -28.21 <0.0001
x3   -1.0095 0.0687 -14.68 <0.0001

I understand the second part of the outcome: the coefficients I put in the model are here (it is almost perfect with N = 100000). The sign is reversed because in my model I used the coeffs to compute the odds of being $\le 1$, and $\le 2$, here it’s the other way, I think there’s not much issues there.

However I don’t understand the discrimination and rank discrimination indexes. Can you help me?! Some pointers?

Things are worse when we turn to cross validation...

> validate(fit, method="cross")
          index.orig training    test optimism index.corrected  n
Dxy           0.3278   0.3278  0.3290  -0.0012          0.3291 40
R2            0.1260   0.1260  0.1313  -0.0053          0.1313 40
Intercept     0.0000   0.0000 -0.0072   0.0072         -0.0072 40
Slope         1.0000   1.0000  1.0201  -0.0201          1.0201 40
Emax          0.0000   0.0000  0.0056   0.0056          0.0056 40
D             0.1164   0.1165  0.1186  -0.0021          0.1186 40
U            -0.0002  -0.0002 -0.8323   0.8321         -0.8323 40
Q             0.1166   0.1167  0.9509  -0.8342          0.9509 40
B             0.1865   0.1865  0.1867  -0.0001          0.1867 40
g             0.7786   0.7786  0.7928  -0.0142          0.7928 40
gp            0.1472   0.1472  0.1478  -0.0007          0.1478 40

Mmffff? What’s this? How do I interpret this? The man page gives few explaination, I don’t have access to this paper... and I feel overwhelmed by an ocean of complexity. Please help!

Answer 1

Concentrate on a few of the indexes right now. index.orig is the apparent predictive ability/accuracy score when you evaluate it on the data used to fit the model. index.corrected is the cross-validation-corrected version of the same index, i.e., corrected for overfitting (de-biased). Dxy is Somers' $D_{xy}$ rank correlation coefficient - a measure of pure discrimination. See original paper or nonparametric texts for details. $D_{xy} = 2(C - \frac{1}{2})$ where $C$ is the generalized ROC area (concordance probability). Intercept and Slope pertain to the calibration curve on the logit scale. Emax is the estimated maximum calibration error using that slope and intercept. B is the Brier accuracy score (combines discrimination and calibration).

Methods are described in my book or the course notes on the book's web site: http://biostat.mc.vanderbilt.edu/rms

Answer 2

To your first question: they are different model performance measurements. You want some of them to be big, others to be small. In fact they are somehow related so suggest you focus on one or two

To your second question: what you have in your first R output is your model performance on training sample. When you do validation by cross-validation, you get those measurements on multiple test samples and their average. This give you more realistic estimate of how your model performs.

HTH