# Simulating data for an ordered logit model

I am currently working on an project where my dependent variable is available in the form of a Likert scale from 1-5. But before I will run my model on real data I want to test it on artificial data, but I am struggling how to simulate data for this kind of data.

I want to build a Bayesian ordered logit model (later I want to introduce random effects, but that doesn't matter at this point) in the spirit of Polson, Scott and Windle (2013) and their package "BayesLogit" in R.

I hope anyone can help!

Imagine you have one ordered categorical dependent variable, which is a Likert scale from 1 to 5 (I'll code it as y), and two independent variables that are normally and uniformly distributed (x1 and x2, respectively).

What a basic ordered logistic regression does is predict $k - 1$ log odds, where $k$ is how many levels the dependent variable has (in this case, 4):

First, the log of the sum of probabilities of being in 2, 3, 4, and 5 categories over the probability of being in the 1 category. Then, categories 3, 4, and 5 over the sum of 1 and 2, etc. These look like:

$\text{ln}({\hat{\pi}_2 + \hat{\pi}_3 + \hat{\pi}_4 + \hat{\pi}_5 \over \hat{\pi}_1})$, $\text{ln}({\hat{\pi}_3 + \hat{\pi}_4 + \hat{\pi}_5 \over \hat{\pi}_1 + \hat{\pi}_2})$, $\text{ln}({\hat{\pi}_4 + \hat{\pi}_5 \over \hat{\pi}_1 + \hat{\pi}_2 + \hat{\pi}_3})$, $\text{ln}({\hat{\pi}_5 \over \hat{\pi}_1 + \hat{\pi}_2 + \hat{\pi}_3 + \hat{\pi}_4})$

For ease, let's call these log-odds 1, 2, 3, and 4, respectively.

So how do we get these four log odds predictions? Each of them has their own equation:

$\text{log-odds}_1 = b_{01} + b_1X_1 + b_2X_2$

$\text{log-odds}_2 = b_{02} + b_1X_1 + b_2X_2$

$\text{log-odds}_3 = b_{03} + b_1X_1 + b_2X_2$

$\text{log-odds}_4 = b_{04} + b_1X_1 + b_2X_2$

You will note that the only thing different about these equations is that they have a different intercept. This is the proportional odds or parallel regressions assumption. Now, we can simulate data based on this model. Since they are ordered, $b_{01} > b_{02} > b_{03} > b_{04}$ (i.e., because the odds of a, b, c, and d happening, for example, cannot be less than the odds of a, b, and c happening).

First, let's make x1 and x2:

set.seed(1839)
n <- 150
x1 <- rnorm(n)
x2 <- runif(n)


Now, let's set the population parameters that predict the log-odds:

b01 <- 1
b02 <- 0.05
b03 <- -0.05
b04 <- -1
b1 <- 0.8
b2 <- 0.5


Now, we can predict each of the log-odds, given the equations above:

logodds1 <- b01 + b1 * x1 + b2 * x2
logodds2 <- b02 + b1 * x1 + b2 * x2
logodds3 <- b03 + b1 * x1 + b2 * x2
logodds4 <- b04 + b1 * x1 + b2 * x2


Now we can use $\text{logit}^{-1}$ to get the probability of the numerator for each person:

inv_logit <- function(logit) exp(logit) / (1 + exp(logit))
prob_2to5 <- inv_logit(logodds1)
prob_3to5 <- inv_logit(logodds2)
prob_4to5 <- inv_logit(logodds3)
prob_5 <- inv_logit(logodds4)


Now, we can do subtraction to get the probability of each category itself:

prob_1 <- 1 - prob_2to5
prob_2 <- prob_2to5 - prob_3to5
prob_3 <- prob_3to5 - prob_4to5
prob_4 <- prob_4to5 - prob_5


Just to make sure, are all probabilities above zero?

> table(c(prob_1, prob_2, prob_3, prob_4, prob_5) > 0)

TRUE
750


Cool. Now, we can use these probabilities to sample from numbers 1 through 5 (the categories of interest) to get our dependent variable, y:

y <- c()
for (i in 1:n) {
y[i] <- sample(
x = c(1:5),
size = 1,
prob = c(prob_1[i], prob_2[i], prob_3[i], prob_4[i], prob_5[i])
)
}


I could have done an apply statement there, but I'm being lazy.

Now, we can put these into a data.frame and run an ordinal regression:

> dat <- data.frame(x1, x2, y = factor(y))
>
> library(ordinal)
> summary(clm(y ~ x1 + x2, data = dat))
formula: y ~ x1 + x2
data:    dat

logit flexible  150  -204.23 420.47 8(1)  2.74e-07 2.2e+02

Coefficients:
Estimate Std. Error z value Pr(>|z|)
x1   0.8933     0.1786   5.003 5.65e-07 ***
x2   0.5810     0.5595   1.038    0.299
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Threshold coefficients:
Estimate Std. Error z value
1|2 -1.25461    0.35426  -3.541
2|3  0.02907    0.33321   0.087
3|4  0.25541    0.33261   0.768
4|5  0.91630    0.34070   2.689


Note that we recovered the parameters OK, but not exactly. Why? The sample function has random error. You could write that into a function, do it a number of times, and see how many times you recover your parameters, etc.

The code is not perfect, but it should get you started. It's also more verbose for instructional purposes.

Since you're doing it in a Bayesian framework, I would really appreciate an answer to my CrossValidated question on choosing semi-informative priors for ordinal regression.