## Code

```
# Function to curtail to [0.02, 0.98] before taking logit
<- function(p) qlogis(pmax(pmin(p, 1. - 0.02), 0.02))
lg # Function to curtail log ORs to [-6, 6]
<- function(x) pmin(pmax(x, -6), 6)
cu
# Function to quantify degree of non-proportional odds first by computing the
# Gini's mean difference of the difference between
# two logit of ECDFs. Quantifies variability of differences over y
# When ECDF is 0 or 1 replace by 0.02, 0.98 so can take logit
# Note that ecdf produces a function and when it is called with an
# x-value that outside the range of the data the value computed is 0 or 1
# Computes a second index of non-PO by getting a weighted standard deviation of
# all possible log ORs, where weights are inverse of variance of log OR
# Note that when a P(Y>=y) is 0 or 1 the weight is zero because variance is infinite.
<- function(y1, y2, pl=FALSE, xlim=range(ys),
npod ylim=range(lg(f1(r)), lg(f2(r))), axes=TRUE) {
<- ecdf(y1)
f1 <- ecdf(y2)
f2 <- c(y1, y2)
y <- range(y)
r <- sort(unique(y))
ys # There cannot be non-PO if only 2 levels of y, and if no overlap
# there is no way to assess non-PO
if(length(ys) <= 2 || max(y1) <= min(y2) || max(y2) <= min(y1))
<- npo2 <- 0.
npo1 else {
<- lg(f1(ys)) - lg(f2(ys))
dif <- GiniMd(dif)
npo1 <- w <- numeric(length(ys) - 1)
lor <- length(y1)
n1 <- length(y2)
n2 for(j in 2:length(ys)) {
<- ys[j]
y <- mean(y1 >= y)
p1 <- mean(y2 >= y)
p2 if(min(p1, p2) == 0 || max(p1, p2) == 1) lor[j-1] <- w[j-1] <- 0
else {
-1] <- log(p2 / (1. - p2)) - log(p1 / (1. - p1))
lor[j-1] <- 1. / ((1. / (n1 * p1 * (1. - p1))) +
w [j1. / (n2 * p2 * (1. - p2))))
(
}
}<- sqrt(wtd.var(lor, w, normwt=TRUE))
npo2
}if(pl) {
plot (ys, lg(f1(ys)), type='s', xlab='', ylab='',
xlim=xlim, ylim=ylim, axes=axes)
lines(ys, lg(f2(ys)), type='s', col='red')
text(xlim[1], ylim[2],
paste0('npo1=', round(npo1, 2), '\nnpo2=', round(npo2, 2)),
adj=c(0,1))
}c(npo1=npo1, npo2=npo2)
}
getRs('reptools.r', put='source') # defines kabl, htmlList
getRs('hashCheck.r', put='source') # defines runifChanged
<- 50; n0 <- n1 <- 25
n par(mfrow=c(3,4), mar=c(1.5,1.5,.5,.5), mgp=c(.5, .4, 0))
set.seed(368)
<- matrix(NA, nrow=11, ncol=2)
z for(i in 1 : 11) {
<- runif(n)
p0 # Note that sample uses prob as weights and they need not sum to 1
<- sample(1 : n, n0, prob=p0, replace=TRUE)
y0 <- runif(n)
p1 <- sample(1 : n, n1, prob=p1, replace=TRUE)
y1 <- npod(y0, y1, pl=TRUE, xlim=c(0, 50), ylim=c(-4, 4))
z[i, ]
}plot(z[, 1], z[, 2], xlab='npo1', ylab='npo2')
```