Ordinal Models

Author
Affiliation

Matti Vuorre

Tilburg University

Introduction

What are ordinal data

Figure 1: A mood item
  • Ordinal data are common in psychology
  • Most common are Likert items
  • But also e.g. item above; school grades; number of forks you own; discrete temporal data

Methods for analysis

  • Metric models
    • Models that assume outcomes have a continuous distribution, e.g. t-test
    • Overestimate information in data; common & “simple”
  • Nonparametric statistics
    • e.g. analyses of signed ranks (R: ?wilcox.test, etc.)
    • Underestimate information in data; don’t scale well
  • Ordinal models
    • A zoo of models that treat outcomes appropriately as ordered categories
    • Let’s learn more!
    • “Ordinal Regression Models in Psychology: A Tutorial” (Bürkner and Vuorre 2019)

“Analyzing ordinal data with metric models: What could possibly go wrong?”

  • Liddell and Kruschke surveyed 68 Psychology articles that analysed ordinal data, and found that every article used metric models (2018)
  • Metric models on ordinal data can lead to false alarms, failures to detect true effects, distorted effect size estimates, and inversions of effects
  • Three main shortcomings of metric models:
    • Response categories may not be (e.g. psychologically) equidistant
    • Responses can be non-normally distributed
    • Can treat differences in variances of underlying variable inappropriately
  • I don’t mean to be an alarmist, or ignore practical considerations. We don’t know the empirical rate of differences. But…

“Analyzing ordinal data with metric models: What could possibly go wrong?”

Code 
movies <- read_rds(here("data/movie-data.rds"))
movies510 <- filter(movies, movie %in% c(5, 10))
movies510 <- pivot_longer(
  movies510,
  cols = -c(movie, title),
  names_to = "Rating",
  names_transform = ~str_remove(.x, "n") |> as.integer(),
  values_to = "Count"
) |>
  uncount(Count)
movies510 |>
  ggplot(aes(Rating)) +
  geom_histogram(binwidth = .5, center = 0) +
  scale_y_continuous(expand = expansion(c(0, .1))) +
  facet_wrap("movie", nrow = 1, labeller = label_both)
Figure 2: IMDB ratings of two movies.
Code 
y <- t.test(scale(Rating) ~ movie, data = movies510)
y <- tidy(y)
# Reverse because t-test subtracts latter level from former
y <- mutate(y, across(where(is.numeric), ~ -round(., 2)))

y2 <- clm(
  ordered(Rating) ~ movie,
  ~movie,
  link = "probit",
  data = movies510
)
y2 <- tidy(y2, conf.int = TRUE)[5, ]
y2 <- mutate(y2, across(where(is.numeric), ~ round(., 2)))
  • Welch’s t-test: Movie 10’s mean rating was significantly greater (Standardized difference = 0.14 [0.23, 0.04])
  • Cumulative probit model: Movie 10’s mean rating was significantly smaller (Difference = -0.2 [-0.31, -0.09])
  • I cherry-picked this example but it exists

Cumulative model

Ordinal models

  • There are many different ordinal models
  • We focus on the cumulative model (CM)
    • Generally the most useful / widely applicable model
  • IRT? SDT?
    • Also known as graded response model, SDT model for ratings, etc…

Example data

Code 
dat <- read_rds(here("data/forstmann.rds"))
  • 1,225 festivalgoers were asked about their mood, substance use, degree of experiencing a “transformative experience”
  • The mood rating item, \(Y\), had \(K + 1 = 6\) categories \(1, 2, ..., 6\)
Code 
head(dat) |>
  kable()
mood te gender age survey h24
5 3 1 30 Event3 0
5 2 2 50 Event5 0
6 7 1 32 Event4 0
6 7 2 33 Event2 0
5 1 1 22 Event3 0
5 7 2 22 Event4 0
Table 1: First six rows of data from Forstmann et al. (2020).

Cumulative model

  • CM assumes that the observed categorical variable \(Y\) is based on the categorization of an unobserved (“latent”) variable \(\tilde{Y}\) with \(K\) thresholds \(\tau = (\tau_1, \dots, \tau_k)\).
  • In this example, \(\tilde{Y}\) has a natural interpretation as current mood
  • We assume that \(\tilde{Y}\) has a normal distribution, but other choices are possible, such as (default) logistic
  • Describe the ordered distribution of responses using thresholds
    • \(Y = k \Leftrightarrow \tau_{k-1} < \tilde{Y} \leq \tau_k\)
  • These thresholds give the probability of each response category
    • \(Pr(Y = k) = \Phi(\tau_k) - \Phi(\tau_{k-1})\)
  • \(\tilde{Y}\) is amenable to regression (without intercept)
    • \(\tilde{Y} \sim N(\eta, \sigma = 1); \ \eta = b_1x_1 +...\)
    • \(Pr(Y = k \vert \eta) = \Phi(\tau_k - \eta) - \Phi(\tau_{k-1} - \eta)\)

Cumulative model

Code 
tab <- count(dat, mood, name = "Count") |>
  mutate(
    p = Count / sum(Count),
    cp = cumsum(p),
    z = qnorm(cp)
  )

p0 <- tab |>
  ggplot(aes(mood)) +
  geom_col(aes(y = Count)) +
  labs(x = "Mood") +
  scale_y_continuous(expand = expansion(c(0, .1)))
p0
Figure 3: Counts of responses in six mood categories.

Cumulative model

Code 
x <- tidy(ordinal::clm(ordered(mood) ~ 1, link = "probit", data = dat))
thresholds <- pull(x, estimate)
x <- tibble(
  x = seq(-4, 4, by = .01),
  y = dnorm(x)
)

p1 <- x |>
  ggplot(aes(x, y)) +
  geom_line(size = 1) +
  scale_y_continuous(expand = expansion(c(0, .3))) +
  scale_x_continuous(
    expression(tilde(Y))
  ) +
  theme(
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )
p1
Figure 4: Distribution of latent normal mood.

Cumulative model

Code 
p1 +
  scale_x_continuous(
    expression(tilde(Y)),
    breaks = thresholds,
    labels = c(expression(tau[1]), ~ tau[2], ~ tau[3], ~ tau[4], ~ tau[5])
  ) +
  geom_vline(xintercept = thresholds, size = .25)
Figure 5: Distribution of latent normal mood with estimated thresholds.

Cumulative model

Code 
tab <- count(dat, mood) |>
  mutate(
    p = n / sum(n),
    cp = cumsum(p),
    z = qnorm(cp)
  )
Code 
kable(tab, digits = 2)
mood n p cp z
1 5 0.00 0.00 -2.63
2 10 0.01 0.01 -2.23
3 23 0.02 0.03 -1.84
4 151 0.13 0.16 -0.98
5 661 0.57 0.73 0.62
6 310 0.27 1.00 Inf
Table 2: Calculating thresholds by hand.

Cumulative model

Code 
tab <- count(dat, mood) |>
  mutate(
    p = n / sum(n),
    cp = cumsum(p),
    z = qnorm(cp)
  )
Code 
p2 <- tab |>
  ggplot(aes(mood, cp)) +
  geom_line() +
  xlab("Mood") +
  geom_point(shape = 21, fill = "white")
p3 <- tab[-6, ] |>
  ggplot(aes(mood, z)) +
  geom_line() +
  xlab("Threshold") +
  geom_point(shape = 21, fill = "white")
(p0 | p2 | p3) +
  scale_x_continuous(breaks = 1:6)
Figure 6: Illustration of how thresholds are calculated from data.

In practice

Code 
library(brms) # Bayesian, slower, more flexible
library(ordinal) # Frequentist, fast, less flexible
  • So far we have described a weird link function + intercepts
  • Write your regressions in R (brms) modelling syntax
  • Effects on \(\tilde{Y}\) are directly interpretable

My first cumulative model

Code 
dat <- mutate(
  dat, 
  h24 = factor(h24, labels = c("No", "Yes"))
)
fit1 <- brm(
  mood ~ h24,
  family = cumulative("probit"),
  data = dat,
  file = here("models/ordinal-1")
)

  • family = cumulative(): CM

  • "probit": \(\tilde{Y} \sim {N}(\eta, \sigma = 1)\)

  • mood ~ h24: \(\eta = b_1\text{h24}\)

  • \(b_1\) is the degree to which mood is greater in people who used hallucinogens in the past 24 hours, compared to people who didn’t use

  • Scale of the latent variable (standard deviations)

Cumulative model

Code 
summary(fit1)
#>  Family: cumulative 
#>   Links: mu = probit; disc = identity 
#> Formula: mood ~ h24 
#>    Data: dat (Number of observations: 1160) 
#>   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup draws = 4000
#> 
#> Regression Coefficients:
#>              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept[1]    -2.63      0.15    -2.95    -2.34 1.00     2324     2419
#> Intercept[2]    -2.20      0.10    -2.40    -2.01 1.00     3627     3111
#> Intercept[3]    -1.80      0.07    -1.94    -1.66 1.00     4202     2939
#> Intercept[4]    -0.93      0.05    -1.02    -0.84 1.00     3977     3305
#> Intercept[5]     0.69      0.04     0.61     0.77 1.00     3898     3065
#> h24Yes           0.39      0.09     0.21     0.56 1.00     4301     2504
#> 
#> Further Distributional Parameters:
#>      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> disc     1.00      0.00     1.00     1.00   NA       NA       NA
#> 
#> Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).

Cumulative model

Code 
plot(fit1, pars = "b_", N = 6)

Cumulative model

Code 
thresholds <- fixef(fit1)[1:5, 1]
beta1 <- fixef(fit1)[6, 1]
x <- tibble(
  x = seq(-4, 4, by = .01),
  No = dnorm(x),
  Yes = dnorm(x, beta1)
)
x |>
  pivot_longer(c(No, Yes), names_to = "h24") |>
  ggplot(aes(x, value, col = h24)) +
  geom_vline(xintercept = thresholds, size = .25) +
  geom_line(size = 1) +
  scale_color_brewer(
    "Hallucinogens past 24 hours",
    palette = "Set1"
  ) +
  scale_y_continuous(expand = expansion(c(0, .3))) +
  scale_x_continuous(
    expression(tilde(Y)),
    breaks = thresholds,
    labels = c(expression(tau[1]), ~ tau[2], ~ tau[3], ~ tau[4], ~ tau[5])
  ) +
  theme(
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )

Cumulative model

Code 
conditional_effects(fit1, categorical = TRUE)

Cumulative model

Code 
pp_check(fit1, "bars_grouped", group = "h24")

Cumulative model

  • It is considered SOP to not assume equal variances when we do e.g. t tests
  • Metric models can deal terribly with different variances in \(\tilde{Y}\)
  • We can predict the variance (without intercept–must be fixed for baseline)
  • \(Pr(Y = k \vert \eta, disc) = \Phi(disc \times (\tau_{k+1} - \eta)) - \Phi(disc \times (\tau_{k} - \eta))\)
  • disc?
    • IRT: Discrimination parameter (slope of response function)
    • Predicted on the log scale \(disc = exp(\eta_{disc})\)
    • \(\sigma\) = \(1 / exp(disc)\)
  • \(\tilde{Y} \sim N(\eta, 1/exp(\eta_{disc})); \ \eta = b_1x_1 +...; \eta_{disc} = g_1x_2 + ...\)

Cumulative model

Code 
fit2 <- brm(
  bf(mood ~ h24) +
    lf(disc ~ 0 + h24, cmc = FALSE),
  family = cumulative("probit"),
  data = dat,
  file = here("models/ordinal-2")
)

Cumulative model

Code 
summary(fit2)
#>  Family: cumulative 
#>   Links: mu = probit; disc = log 
#> Formula: mood ~ h24 
#>          disc ~ 0 + h24
#>    Data: dat (Number of observations: 1160) 
#>   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup draws = 4000
#> 
#> Regression Coefficients:
#>              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept[1]    -2.62      0.16    -2.95    -2.33 1.00     2586     2250
#> Intercept[2]    -2.18      0.10    -2.40    -1.99 1.00     3440     2578
#> Intercept[3]    -1.78      0.07    -1.93    -1.64 1.00     4079     2950
#> Intercept[4]    -0.92      0.05    -1.01    -0.82 1.00     4161     3300
#> Intercept[5]     0.68      0.04     0.60     0.77 1.00     3863     3198
#> h24Yes           0.37      0.09     0.20     0.55 1.00     4322     2518
#> disc_h24Yes      0.08      0.08    -0.08     0.24 1.00     3838     3084
#> 
#> Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).

Cumulative model

Code 
as_draws_df(fit2, variable = "h24", regex = TRUE) |>
  mutate(sigma_h24 = 1 / exp(b_disc_h24Yes)) |>
  pivot_longer(contains("h24")) |>
  ggplot(aes(value, name)) +
  stat_histinterval()

Cumulative model

Code 
thresholds <- fixef(fit2)[1:5, 1]
beta1 <- fixef(fit2)[6, 1]
disc1 <- 1 / exp(fixef(fit2)[7, 1])
x <- tibble(
  x = seq(-4, 4, by = .01),
  No = dnorm(x),
  Yes = dnorm(x, beta1, disc1)
)
x |>
  pivot_longer(c(No, Yes), names_to = "h24") |>
  ggplot(aes(x, value, col = h24)) +
  geom_vline(xintercept = thresholds, size = .25) +
  geom_line(size = 1) +
  scale_y_continuous(expand = expansion(c(0, .3))) +
  scale_x_continuous(
    expression(tilde(Y)),
    breaks = thresholds,
    labels = c(expression(tau[1]), ~ tau[2], ~ tau[3], ~ tau[4], ~ tau[5])
  ) +
  theme(
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )

More ordinal models

Category specific effects

  • Cannot use CM*
  • Adjacent category model: predict decisions between categories
Figure 7: Adjacent category model.

Category specific effects

  • There are two cells with no observations
Code 
table(dat$h24, dat$mood)
#>      
#>         1   2   3   4   5   6
#>   No    5  10  22 136 561 242
#>   Yes   0   0   1  15 100  68
  • Those category-specific effects won’t be identified
  • If only there was a way to inject information to the model…
  • Bayes to the rescue!
Code 
weakly_informative_prior <- prior(normal(0, 1.5), class = "b")
fit3 <- brm(
  bf(mood ~ cs(h24)),
  family = acat("probit"),
  prior = weakly_informative_prior,
  data = dat,
  control = list(adapt_delta = .95),
  file = here("models/ordinal-3")
)

Category specific effects

Code 
summary(fit3)
#>  Family: acat 
#>   Links: mu = probit; disc = identity 
#> Formula: mood ~ cs(h24) 
#>    Data: dat (Number of observations: 1160) 
#>   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup draws = 4000
#> 
#> Regression Coefficients:
#>              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept[1]    -0.44      0.34    -1.11     0.20 1.00     3732     2856
#> Intercept[2]    -0.52      0.23    -0.97    -0.08 1.00     3159     2726
#> Intercept[3]    -1.09      0.12    -1.33    -0.84 1.00     3142     2779
#> Intercept[4]    -0.86      0.06    -0.97    -0.75 1.00     3525     2854
#> Intercept[5]     0.52      0.05     0.43     0.61 1.00     3997     3026
#> h24Yes[1]        0.38      1.28    -2.00     3.02 1.00     4609     2614
#> h24Yes[2]        1.05      1.13    -0.98     3.45 1.00     4261     3092
#> h24Yes[3]        0.61      0.51    -0.29     1.70 1.00     4065     2661
#> h24Yes[4]        0.28      0.16    -0.02     0.60 1.00     3709     2692
#> h24Yes[5]        0.28      0.11     0.06     0.49 1.00     4091     3115
#> 
#> Further Distributional Parameters:
#>      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> disc     1.00      0.00     1.00     1.00   NA       NA       NA
#> 
#> Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).

Category specific effects

Code 
conditional_effects(fit3, categorical = TRUE)

Category specific effects

Code 
x <- conditional_effects(fit1, categorical = TRUE)[[1]]
p1 <- x |>
  ggplot(aes(cats__, estimate__, col = h24)) +
  geom_pointrange(
    aes(ymin = lower__, ymax = upper__),
    position = position_dodge(.25)
  ) +
  labs(
    subtitle = "Cumulative model",
    x = "Response category (mood)",
    y = "Probability"
  )
x <- conditional_effects(fit3, categorical = TRUE)[[1]]
p2 <- x |>
  ggplot(aes(cats__, estimate__, col = h24)) +
  geom_pointrange(
    aes(ymin = lower__, ymax = upper__),
    position = position_dodge(.25)
  ) +
  labs(
    subtitle = "Adjacent category model (CS)",
    x = "Response category (mood)",
    y = "Probability"
  )
(p1 | p2) + plot_layout(guides = "collect")

Model comparison

Code 
fit4 <- brm(
  bf(mood ~ h24),
  family = acat("probit"),
  prior = weakly_informative_prior,
  data = dat,
  file = here("models/ordinal-4")
)

# Add LOO criteria to all models
fit1 <- add_criterion(fit1, "loo")
fit2 <- add_criterion(fit2, "loo")
fit3 <- add_criterion(fit3, "loo")
fit4 <- add_criterion(fit4, "loo")
Code 
loo_compare(
  fit1, fit2, fit3, fit4
)
#>      elpd_diff se_diff
#> fit4  0.0       0.0   
#> fit1 -0.7       0.6   
#> fit2 -1.2       0.5   
#> fit3 -2.9       2.1

Objections and counterarguments

Its too difficult

  • …
  • There are, of course, practical considerations
  • The weird link function and intercepts were difficult
  • Effects on latent variable are interpretable just like your betas in lm()

The results are the same anyway

  • How do you know?
  • Did you fit an ordinal model to confirm?
  • The prevalence of problems in metric models applied to ordinal data is an empirical questions, and results probably vary greatly between types of data & measures
  • Fit ordinal models whenever you can
  • Afford more nuanced interpretation of what’s going on in your data

Multilevel model

So far…

  • We did not consider variability in mood beyond hallucinogen use
  • Modeled H as fixed effect
  • Age? Survey (different events)? Interactions??
  • Multilevel model is a kind of interaction model

Multilevel model

Code 
fit5 <- brm(
  bf(mood ~ h24 + (h24 | survey)),
  family = cumulative("probit"),
  data = dat,
  file = here("models/ordinal-5")
)

Events

Code 
summary(fit5)
#>  Family: cumulative 
#>   Links: mu = probit; disc = identity 
#> Formula: mood ~ h24 + (h24 | survey) 
#>    Data: dat (Number of observations: 1160) 
#>   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup draws = 4000
#> 
#> Multilevel Hyperparameters:
#> ~survey (Number of levels: 6) 
#>                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept)           0.30      0.20     0.03     0.81 1.00      942      738
#> sd(h241)                0.33      0.28     0.01     1.06 1.00      966     1466
#> cor(Intercept,h241)    -0.08      0.57    -0.96     0.92 1.00     2555     2615
#> 
#> Regression Coefficients:
#>              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept[1]    -2.87      0.23    -3.30    -2.41 1.00     2176     1990
#> Intercept[2]    -2.41      0.19    -2.76    -2.02 1.00     2076     2064
#> Intercept[3]    -1.98      0.17    -2.30    -1.59 1.00     1913     2104
#> Intercept[4]    -1.08      0.16    -1.37    -0.71 1.00     1758     1960
#> Intercept[5]     0.57      0.16     0.29     0.94 1.00     1729     1936
#> h241             0.33      0.22    -0.09     0.84 1.00     1347     1144
#> 
#> Further Distributional Parameters:
#>      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> disc     1.00      0.00     1.00     1.00   NA       NA       NA
#> 
#> Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).

Another way to vary intercepts

This is not a great idea but in theory works.

Code 
fit6 <- brm(
  bf(mood | resp_thres(gr = survey) ~ h24 + (0 + h24 | survey)),
  family = cumulative("probit"),
  data = dat,
  file = here("models/ordinal-6")
)
Code 
summary(fit6)
#>  Family: cumulative 
#>   Links: mu = probit; disc = identity 
#> Formula: mood | resp_thres(gr = survey) ~ h24 + (0 + h24 | survey) 
#>    Data: dat (Number of observations: 1160) 
#>   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup draws = 4000
#> 
#> Multilevel Hyperparameters:
#> ~survey (Number of levels: 6) 
#>                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(h24No)             0.89      0.73     0.03     2.69 1.01      709     1618
#> sd(h24Yes)            1.11      0.89     0.04     3.21 1.00      686     1809
#> cor(h24No,h24Yes)     0.55      0.52    -0.80     1.00 1.00     1218     1918
#> 
#> Regression Coefficients:
#>                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept[Event1,1]    -6.00      3.61   -15.40    -1.66 1.01     1367     1040
#> Intercept[Event1,2]    -3.31      1.73    -7.23    -0.63 1.00     1536     1660
#> Intercept[Event1,3]    -1.20      0.92    -2.61     0.91 1.00     1224     3040
#> Intercept[Event1,4]    -0.39      0.91    -1.75     1.74 1.00     1223     3260
#> Intercept[Event1,5]     1.14      0.91    -0.24     3.29 1.00     1214     3155
#> Intercept[Event2,1]    -5.02      2.92   -13.13    -1.41 1.00     1151      959
#> Intercept[Event2,2]    -1.87      0.82    -3.14     0.03 1.00     1188     3161
#> Intercept[Event2,3]    -1.48      0.81    -2.65     0.37 1.00     1165     3394
#> Intercept[Event2,4]    -0.50      0.80    -1.61     1.36 1.00     1151     3219
#> Intercept[Event2,5]     1.22      0.80     0.11     3.09 1.00     1136     3183
#> Intercept[Event3,1]    -6.01      3.01   -13.79    -2.01 1.00     1501      954
#> Intercept[Event3,2]    -3.56      1.63    -7.30    -1.02 1.00     1538     1243
#> Intercept[Event3,3]    -2.09      0.90    -3.48     0.02 1.00     1087     2748
#> Intercept[Event3,4]    -0.48      0.85    -1.60     1.57 1.00     1017     3002
#> Intercept[Event3,5]     1.16      0.85     0.05     3.22 1.00     1005     2911
#> Intercept[Event4,1]    -5.68      3.22   -14.18    -1.69 1.00     2057     1522
#> Intercept[Event4,2]    -3.21      1.52    -6.96    -0.67 1.00     2630     2764
#> Intercept[Event4,3]    -1.71      0.86    -3.13     0.19 1.00     1436     3028
#> Intercept[Event4,4]    -0.35      0.82    -1.53     1.55 1.00     1240     2942
#> Intercept[Event4,5]     1.49      0.82     0.32     3.40 1.00     1282     3273
#> Intercept[Event5,1]    -1.65      0.75    -2.92     0.04 1.00     1611     1805
#> Intercept[Event5,2]    -1.14      0.73    -2.38     0.49 1.00     1613     2456
#> Intercept[Event5,3]    -0.76      0.73    -1.99     0.89 1.00     1629     2350
#> Intercept[Event5,4]    -0.18      0.72    -1.39     1.45 1.00     1637     2360
#> Intercept[Event5,5]     1.46      0.73     0.28     3.09 1.00     1568     2333
#> Intercept[Event6,1]    -2.07      0.81    -3.57    -0.34 1.00     1275     2478
#> Intercept[Event6,2]    -1.78      0.78    -3.13    -0.06 1.00     1171     2301
#> Intercept[Event6,3]    -1.36      0.75    -2.57     0.33 1.00     1048     2427
#> Intercept[Event6,4]    -0.55      0.73    -1.66     1.12 1.00     1076     2288
#> Intercept[Event6,5]     0.79      0.73    -0.28     2.46 1.00     1111     2113
#> h241                    0.09      0.75    -1.69     1.64 1.00      950     1060
#> 
#> Further Distributional Parameters:
#>      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> disc     1.00      0.00     1.00     1.00   NA       NA       NA
#> 
#> Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).

References

Bürkner, Paul-Christian, and Matti Vuorre. 2019. “Ordinal Regression Models in Psychology: A Tutorial.” Advances in Methods and Practices in Psychological Science 2 (1): 77–101. https://doi.org/10.1177/2515245918823199.
Forstmann, Matthias, Daniel A. Yudkin, Annayah M. B. Prosser, S. Megan Heller, and Molly J. Crockett. 2020. “Transformative Experience and Social Connectedness Mediate the Mood-Enhancing Effects of Psychedelic Use in Naturalistic Settings.” Proceedings of the National Academy of Sciences 117 (5): 2338–46. https://doi.org/10.1073/pnas.1918477117.
Liddell, Torrin M., and John K. Kruschke. 2018. “Analyzing Ordinal Data with Metric Models: What Could Possibly Go Wrong?” Journal of Experimental Social Psychology 79 (November): 328–48. https://doi.org/10.1016/j.jesp.2018.08.009.