# Proportional odds assumption in ordinal logistic regression in R with the packages VGAM and rms

An assumption of the ordinal logistic regression is the proportional odds assumption. Using R and the 2 packages mentioned I have 2 ways to check that but I have questions in each one.

1) Using the rms package

Given the next commands

library(rms)
sf <- function(y)
c('Y>=1'=qlogis(mean(y >= 1)),'Y>=2'=qlogis(mean(y >= 2)),'Y>=3'=qlogis(mean(y >= 3)))
plot(s,which=1:3,pch=1:3,xlab='logit',main='',xlim=c(-2.5,2.5))

I have

lrm(formula = grade ~ Ki67 + Cyclin_E)

Frequencies of Missing Values Due to Each Variable
0        0        3

Model Likelihood     Discrimination    Rank Discrim.
Ratio Test            Indexes          Indexes

Obs            42    LR chi2     11.38    R2       0.268    C       0.728
1             11    d.f.            2    g        1.279    Dxy     0.456
2             15    Pr(> chi2) 0.0034    gr       3.592    gamma   0.458
3             16                         gp       0.192    tau-a   0.308
max |deriv| 1e-07                         Brier    0.166

Coef    S.E.   Wald Z Pr(>|Z|)
y>=2     -0.1895 0.8427 -0.22  0.8221
y>=3     -2.0690 0.9109 -2.27  0.0231
Ki67      0.0971 0.0330  2.94  0.0033
Cyclin_E -0.0076 0.0227 -0.33  0.7387

The s table gives: (unfortunately I don't know how to upload a graph made in R)

+--------+-------+--+----+---------+----------+
|        |       |N |Y>=1|Y>=2     |Y>=3      |
+--------+-------+--+----+---------+----------+
|Ki67    |[ 2, 9)|12|Inf |0.6931472|-1.0986123|
|        |[ 9,16)|12|Inf |0.3364722|-2.3978953|
|        |[16,24)|10|Inf |2.1972246| 0.0000000|
|        |[24,44]|11|Inf |2.3025851| 1.5040774|
+--------+-------+--+----+---------+----------+
|Cyclin_E|[ 3,16)|15|Inf |1.0116009|-0.1335314|
|        |[16,22)| 7|Inf |1.7917595|-0.9162907|
|        |[22,33)|10|Inf |1.3862944|-0.8472979|
|        |[33,80]|10|Inf |0.4054651|-0.4054651|
|        |Missing| 3|Inf |      Inf| 0.6931472|
+--------+-------+--+----+---------+----------+
|Overall |       |45|Inf |1.1284653|-0.4054651|
+--------+-------+--+----+---------+----------+

Where for the Ki67 I see that 3 out of the 4 differences logit(P[Y> = 2])-logit(P[Y> = 3]) are close to 2. Only the last one is quite lower (around 0.8). But here Ki67 is continuous and not categorical so I don't know if the results of the table are correct and there isn't any p-value to decide. By the way I run the above in SPSS and I didn't reject the assumption.

2) Using the VGAM package

Here using the next commands I have the model under the assumption of proportional odds

library(VGAM)
summary(fit1)

And the results

Coefficients:
Estimate Std. Error  z value
(Intercept):1  0.1894723   0.820442  0.23094
(Intercept):2  2.0690395   0.886732  2.33333
Ki67          -0.0970972   0.032423 -2.99467
Cyclin_E       0.0075887   0.021521  0.35261

Number of linear predictors:  2

Names of linear predictors: logit(P[Y< = 1]), logit(P[Y< = 2])

Dispersion Parameter for cumulative family:   1

Residual deviance: 79.86801 on 80 degrees of freedom

Log-likelihood: -39.93401 on 80 degrees of freedom

Number of iterations: 5

While using the next commands I have the model without the assumption of proportional odds

where unfortunately i receice the next message

Warning message: In vglm.fitter(x = x, y = y, w = w, offset = offset, Xm2 = Xm2, : convergence not obtained in 30 iterations

However if I type summary(fit2) I get results but again I don't know if they are correct. My intention was to use the next commands and get the answer but know I doubt if this is correct (by the way if I do it I get p-value=0.6.

pchisq(deviance(fit1)-deviance(fit2),
df=df.residual(fit1)-df.residual(fit2),lower.tail=FALSE)

So, regarding the methods mentioned above does anyone knows whether the results I get are valid or in the case of the VGAM package is there any way to increase the number of itterations?Is there any other way to check it?

-

I recommend partial residual plots using the rms package's lrm and residuals.lrm functions. You can also fit a series of binary models using different cutoffs for $Y$ and plot the log odds ratios vs. cutoff.

-

Are there robust standard errors available in the package? or will the sandwich package calculate them for you?

The dirty secret about model checking is there is almost never good statistical power to detect meaningful violations of your model. You can still check them, but more as post-hoc sensitivity analysis.

The best is to use methods that are agnostic as to whether your the assumptions under which your model is 'optimal' are correct. E.g. proportional hazards are almost always violated in cox models, and linear models likely never correspond to some true underlying linear structure. Yet, with robust standard errors, it doesn't matter, and we can still make statements about trends.

-
Thanks for your answer. Unfortunately I didn't find any robust standard errors but I found a solution.In vglm() there is an option to increase the number of iterations and so I achieved the convergence needed. However,could you or someone else tell me how I can upload a diagram made in R because I want to ask something regarding the plot.xmean.ordinaly() command but I don't know how to do it. – Nick Apr 8 '12 at 14:48
Although depending on your graphical interface to R there could be an easier way, and while you can always take a screenshot and prune it, the standard way is to create a device (e.g. png(filename), then plot() your graphics, then call dev.off() and you should have a nice image suited for uploading here). – miura Oct 1 '12 at 8:29
In the R rms package robust covariance matrix estimates are obtained by running the fit through robcov() after creating it with a fitting function such as orm. – Frank Harrell Oct 20 '15 at 0:44