# Bootstrapped coefficients for ordinal logistic regression with R

I am trying to get the bootstrapped confidence intervals of the coefficients for an ordinal logistic regression. Here below, my R code on fake data (reproducible example here below). This one does not work properly.

I suppose I need to enter a list of data with one line for each of the 20 subjects (this is the most simple way to proceed). Then the bootstrap with randomly select 20 rows (using sampling with replacement) to generate a new data set with 20 rows. That data set is converted into a new table of counts and a coefficient value is computed from that new “bootstrapped” table. This is repeated for each bootstrap sample. I can't get it! Thanks for your help.

####################
library(rms)
x=c(1,2,3,2,3,1,2,3,3,3,2,2,1,2,1,2,3,2,1,2)
y=c("math","eco","eco","lit","lit","eco","eco","math","math","lit","lit","math","eco","eco","math","lit","lit","math","eco","math")
Dataset<-data.frame(x,y)
h <- orm(x ~ y)
h

# calculate coefficients using bootstrap
library(boot)
logit.bootstrap <- function(data, indices) {
d<-data[indices,]
fit<-orm(x ~ y, data=data[indices,])
return(coefficients(fit))
}

# bootstrapping with 1000 replications
logit.boot <- boot(data=Dataset, statistic=logit.bootstrap,R=1000)

# view results
logit.boot
plot(logit.boot)

# get 95% confidence interval
boot.ci(logit.boot, type="all")
############################


I'm not really sure, but here are a few thoughts:

### Function not generating all statistics

My suspicion is that some of the bootstrap samples are not generating all the coefficients. This might make sense given that you are running ordinal regression. Thus, there might be something about needing to have enough data for each level of a predictor or an outcome variable.

See for example, this question on stackoverflow.

For example, when I ran the above code with only a few bootstrap samples (e.g., R = 100), quite often the code worked. Thus, presumably its only the occasional bootstrap sample that is causing problems. Presumably, if this is the issue, there would be some other bootstrap scheme that would resolve this issue.

### Try extracting each coefficient one at a time

Once again, I'm not sure if this is the answer, but I tried extracting each coefficient one at a time, and it seemed to get a bootstrap confidence interval for each coefficient. I'm not that familiar with ordinal regression, so treat with caution:

# calculate coefficients using bootstrap
library(boot)
logit.bootstrap <- function(data, indices, coefname = "y>=2") {
d<-data[indices,]
fit<-orm(x ~ y, data=d)
return(coefficients(fit)[coefname])
}

# get names of coefficients
z <- coef(orm(x ~ y, Dataset))
names(z)

# bootstrapping with 1000 replications
# loop over each coefficient and save each element to a list
lb <- lapply(names(z), function(X)
boot(data=Dataset, statistic= logit.bootstrap, R=1000, coefname= X))
names(lb) <- names(z)

# view results
sapply(lb, plot)

# get 95% confidence interval
lapply(lb, function(X) boot.ci(X, type="all"))