Suppose we have a 1–7 satisfaction or agreement scale and want to know whether one group shows more disagreement than another. A natural approach is to calculate a standard deviation for each unit and regress it on group indicators or other covariates.
The problem is that the variance of a bounded outcome depends mechanically on its mean. The Bhatia–Davis (2000) inequality says that, for a variable bounded in \([a,b]\) with mean \(\mu\),
\[\text{Var}(X) \leq (\mu - a)(b - \mu) \tag{37.1}\]
For a 1–7 item, the maximum possible standard deviation is therefore
\[\text{SD}_{\max}(\mu) = \sqrt{(\mu - 1)(7 - \mu)}\]
mu_grid <- seq(1, 7, length.out = 200)
sd_max <- sqrt(pmax((mu_grid - 1) * (7 - mu_grid), 0))
ggplot(data.frame(mu = mu_grid, sd_max = sd_max), aes(mu, sd_max)) +
geom_line(linewidth = 1, color = "#2d5fa3") +
labs(x = "Mean Rating (μ)", y = expression(SD[max](mu)),
title = "Maximum Possible SD as a Function of the Mean") +
theme_minimal(base_size = 12)
The maximum is 3 at the middle of the scale, where \(\mu=4\), and falls to zero at either endpoint. Satisfaction and performance ratings are often close to the upper endpoint. A group with a higher mean can then have a smaller observed standard deviation even when its underlying dispersion is the same.
Dividing the standard deviation by the mean does not remove the dependence. That adjustment assumes that variance grows proportionally with the mean, while variance on a bounded scale is compressed at both endpoints. In this example, dividing by \(\mu\) only moves the peak from \(\mu=4\) to \(\mu=1.75\).
A Gaussian location-scale model has a related problem. A model such as sigma ~ group treats the response as unbounded, so it can interpret compression near an endpoint as a real difference in residual variation.
One way to handle the bounds is to treat the Likert response as a discretized latent variable \(Y^*\):
\[Y_i^* = \eta_i + \sigma_i \epsilon_i, \qquad \epsilon_i \sim N(0, 1)\]
where \(\eta_i\) is the location predictor and \(\sigma_i\) is the latent residual standard deviation. The observed rating \(Y_i \in \{1, \dots, K\}\) is determined by ordered cutpoints \(\theta_1 < \dots < \theta_{K-1}\):
\[Y_i = k \iff \theta_{k-1} < Y_i^* \leq \theta_k\]
Under a probit link, the cumulative probability is:
\[P(Y_i \leq k) = \Phi\left(\frac{\theta_k - \eta_i}{\sigma_i}\right) = \Phi\bigl(\text{disc}_i (\theta_k - \eta_i)\bigr)\]
where discrimination is defined as the inverse of latent dispersion: \(\text{disc}_i \equiv 1/\sigma_i\).
In brms, we model \(\log(\text{disc}_i)\) directly. Higher discrimination means smaller latent dispersion, or more agreement. Lower discrimination means greater spread. The cutpoints are estimated on an unbounded latent scale, so an observed response near 7 does not by itself imply smaller latent dispersion. This interpretation assumes that the groups use the response categories in the same way.
Let’s simulate 100 subjects in each group, with six ratings for each subject. Both groups have latent standard deviation \(\sigma=1\). Group A has a latent mean near the ceiling (\(\mu_A=2.4\)), while Group B is closer to the middle (\(\mu_B=0.6\)).
set.seed(42)
n_per_group <- 100
n_raters <- 6
sigma_true <- 1.0 # Identical latent dispersion across groups
mu_A <- 2.4 # Group A: high mean (near ceiling)
mu_B <- 0.6 # Group B: mid-scale mean
subjects <- data.frame(
subject = 1:(2 * n_per_group),
group = rep(c("A", "B"), each = n_per_group),
mu = rep(c(mu_A, mu_B), each = n_per_group)
)
df <- subjects[rep(seq_len(nrow(subjects)), each = n_raters), ]
df$latent <- rnorm(nrow(df), mean = df$mu, sd = sigma_true)
# Discretize using 6 fixed cutpoints into categories 1 through 7
cuts <- c(-2, -1.2, -0.5, 0.2, 1.0, 1.8)
df$rating <- as.integer(cut(df$latent, breaks = c(-Inf, cuts, Inf), labels = 1:7))
df$rating_ord <- factor(df$rating, ordered = TRUE, levels = 1:7)subj_summary <- df |>
group_by(subject, group) |>
summarise(mean_rating = mean(rating), sd_rating = sd(rating), .groups = "drop")
subj_summary |>
group_by(group) |>
summarise(mean_rating = round(mean(mean_rating), 2),
mean_sd = round(mean(sd_rating), 2)) |>
knitr::kable(caption = "Observed summary statistics by group.")| group | mean_rating | mean_sd |
|---|---|---|
| A | 6.64 | 0.58 |
| B | 4.89 | 1.22 |
ggplot(subj_summary, aes(mean_rating, sd_rating, color = group)) +
geom_jitter(width = 0.05, height = 0, alpha = 0.5) +
geom_line(data = data.frame(mean_rating = mu_grid, sd_rating = sd_max),
aes(mean_rating, sd_rating), color = "grey40", inherit.aes = FALSE) +
labs(x = "Subject Mean Rating", y = "Subject SD",
title = "Boundary Effect: Group A's SD is Mechanically Compressed") +
theme_minimal(base_size = 12)
The observed standard deviation averages 0.58 in Group A and 1.22 in Group B. If we looked only at these numbers, we would conclude that raters agree much more about Group A, even though the two latent standard deviations are identical by construction.
I fit two models in brms. The Gaussian model makes \(\log(\sigma)\) a function of group. The ordinal cumulative-probit model instead makes \(\log(\text{disc})\) a function of group, with Group A used as the reference (\(\text{disc}_A=1\)).
# 1. Gaussian location-scale model
fit_gaussian <- brm(
bf(rating ~ group + (1 | subject), sigma ~ group),
data = df, family = gaussian(),
chains = 2, iter = 2000, warmup = 1000, cores = 2, seed = 1,
refresh = 0, silent = 2
)
# 2. Ordinal cumulative-probit location-scale model
fit_ordinal <- brm(
bf(rating_ord ~ group + (1 | subject), disc ~ 0 + group),
data = df, family = cumulative("probit"), init = 0,
prior = c(
prior(constant(0), class = "b", coef = "groupA", dpar = "disc"),
prior(normal(0, 1), class = "b", coef = "groupB", dpar = "disc")
),
chains = 2, iter = 2000, warmup = 1000, cores = 2, seed = 1,
refresh = 0, silent = 2
)| Estimate | Est.Error | Q2.5 | Q97.5 | |
|---|---|---|---|---|
| Intercept | 6.64 | 0.03 | 6.58 | 6.69 |
| sigma_Intercept | -0.41 | 0.03 | -0.46 | -0.35 |
| groupB | -1.75 | 0.06 | -1.86 | -1.63 |
| sigma_groupB | 0.68 | 0.04 | 0.60 | 0.76 |
| Estimate | Est.Error | Q2.5 | Q97.5 | |
|---|---|---|---|---|
| Intercept[1] | -4.39 | 0.32 | -5.06 | -3.81 |
| Intercept[2] | -3.65 | 0.23 | -4.12 | -3.22 |
| Intercept[3] | -2.98 | 0.17 | -3.34 | -2.66 |
| Intercept[4] | -2.27 | 0.12 | -2.52 | -2.04 |
| Intercept[5] | -1.44 | 0.07 | -1.59 | -1.30 |
| Intercept[6] | -0.63 | 0.06 | -0.74 | -0.52 |
| groupB | -1.89 | 0.10 | -2.11 | -1.69 |
| disc_groupA | 0.00 | 0.00 | 0.00 | 0.00 |
| disc_groupB | -0.01 | 0.08 | -0.17 | 0.15 |
The Gaussian coefficient on sigma_groupB is \(0.68\). It estimates that Group B’s residual standard deviation is \(\exp(0.68)\approx1.97\) times Group A’s (\(95\%\text{ CI}: [1.82,2.14]\)). This apparent difference comes from the scale ceiling.
The ordinal coefficient on disc_groupB is \(-0.01\) (\(95\%\text{ CI}: [-0.17,0.15]\)), which implies a latent standard-deviation ratio of \(\exp(0.01)\approx1.01\). The ordinal model recovers the equal dispersion used to generate the data.
The data were generated from the ordinal model’s own functional form: a latent normal variable, a constant \(\sigma\), and one cuts vector shared by both groups. It is therefore not surprising that the ordinal model fits well. What the example shows is narrower: boundary compression by itself does not create a difference in the latent dispersion parameter, while it does create one in the Gaussian model.
The example does not show that the model can distinguish dispersion from threshold differences, because I did not simulate any threshold difference. Common thresholds across groups, or a model that allows them to vary, are needed for identification. If groups use the categories differently, threshold shifts and disc can trade off: the cutpoints can absorb a real dispersion difference, or a threshold difference can appear as a change in disc. On real data, we need information outside the basic model, such as anchoring vignettes or a partial-invariance specification, to check this assumption.
A simpler descriptive measure is the ratio of the observed variance to the maximum possible variance at the same mean:
\[R = \frac{\text{Var}(Y)}{\text{Var}_{\max}(\mu)} \in [0, 1]\]
where \(R = 0\) indicates perfect agreement and \(R = 1\) indicates maximum polarization.
With only six ratings per subject, \(R\) is noisy and can be undefined when every rater chooses 7, giving \(0/0\). We can instead calculate it after pooling the ratings within each group:
| group | mean_rating | var_obs | var_max | R_pooled |
|---|---|---|---|---|
| A | 6.638333 | 0.4475306 | 2.039197 | 0.219 |
| B | 4.893333 | 1.7286222 | 8.201956 | 0.211 |
Both groups give \(R\approx0.21\)–\(0.22\). Dividing by the maximum attainable variance makes the descriptive measure more comparable across groups with different means. It is still only a normalization, not a model-based correction, and it has no latent dispersion parameter to interpret.
Raw standard deviations are hard to compare when group means differ on a bounded scale. If the goal is a model of underlying response dispersion, I would use the ordinal location-scale model and examine whether common thresholds are plausible. If the goal is only a descriptive comparison, pooled \(R=\text{Var}/\text{Var}_{\max}\) is easier to explain, as long as it is not interpreted as a latent parameter.