23.2.1 Implementation in JAGS

To estimate the parameters, we need to express Equations 23.1-23.4 in JAGS (or Stan), and establish prior distributions. The first part of the JAGS model specification says that each ordinal datum comes from a categorical distribution that has probabilities computed as in Equations 23.1-23.3. In the JAGS code below, Ntotal is the total number of data points, nYlevels is the number of outcome levels, and the matrix pr[i,k] holds the predicted probability that datum y[i] has value k. In other words, pr[i,k] is p(y_i = k|μ, σ, {θ_j}). $R$ecall that the cumulative normal function in JAGS (and $R$) is called pnorm. Now, see if you can find Equations 23.1-23.3 in the first part of the model specification:

model {
 for (i in 1:N total) {
 y[i] ~ dcat(pr[i,1:nYlevels]) (
 pr[i,1] < - pnorm(thresh[1], mu, 1/sigma ^ 2)
 for (k in 2:(nYlevels-1)) {
 pr[i,k] < - max(0, pnorm(thresh[k], mu, 1/sigma ^ 2)
 - pnorm(thresh[k−1], mu, 1/sigma ^ 2))
 }
 pr[i,nYlevels] < ‒ 1 - pnorm(thresh[nYlevels-1], mu, 1/sigma ^ 2)
 }

In the JAGS model specification above, the cumulative normal probabilities that are written as

could instead be coded as

This alternative form (above) better matches the form in Equation 23.1 and may be used instead if it suits your aesthetic preferences. I avoid the form with 0 and 1 because it can be dangerous, from a programming perspective, to use unlabeled constants.

You may be wondering why the expression for pr[i,k], above, uses the maximum of zero and the difference between cumulative normal probabilities. The reason is that the threshold values are randomly generated by the MCMC chain, and it is remotely possible that the value of thresh[k] would be less than the value of thresh[k−1], which would produce a negative difference of cumulative normal probabilitieswhich makes no sense. By setting the difference to zerothe likelihood of y = k is zero, and the candidate threshold values are not accepted.

The model specification continues with the prior distribution on the parameters. As usual, in the absence of specific prior knowledge, we use a prior that is broad on the scale of data. The mean and standard deviation should be somewhere in the vicinity of the data, which can only range from 1 to nYlevels on the ordinal outcome scale. Thus, the prior is declared as:

mu ~ dnorm((1+nYlevels)/2, 1/(nYlevels)~2)
sigma ~ dunif(nYlevels/1000, nYlevels*10)

The prior on the thresholds centers θ_k at k + 0.5 but allows a considerable range:

for (k in 2:(nYlevels-2)){ # 1 and nYlevels-1 are fixed, not stochastic
 thresh[k] ~ dnorm(k+0.5, 1/2~2)
}
} # end of model

Notice above that only the non-end thresholds are given priors, because only they are estimated. The end thresholds are fixed as specified in Equation 23.4. But how is Equation 23.4 implemented in JAGS? This requires specifying the fixed components of the threshold vector in the data list that is delivered to JAGS separately. The other components of the threshold vector are estimated (i.e., stochastically set by the MCMC process). Estimated components are given the value NA, which JAGS interprets as meaning that JAGS should fill them in. Fixed components are given whatever constant values are desired. Thus, in $R$, we state:

thresh = rep(NA, nYlevels-1)   # start by putting NA in all components
thresh[1]=1+0.5   # overwrite 1st component with 1+0.5
thresh[nYlevels-1] = nYlevels ‒ 0.5   # overwrite last component with nYlevels-0.5
dataList = list(  # package thresh into the dataList
 thresh = thresh,
 # … other data and constants
)

The full specification is in the program Jags-Yord-Xnom1grp-Mnormal.R, which is called from the high-level script Jags-Yord-Xnom1grp-Mnormal-Example.R.

23.2.2 Examples: Bayesian estimation recovers true parameter values

Figures 23.2 and 23.3 show results of two examples. In both cases, the data were generated from the model with particular known parameter values, and the goal of the examples is to demonstrate that Bayesian estimation accurately recovers the parameter values even when μ is extreme and when the thresholds are not evenly spaced.

Figure 23.2 shows an example in which the ordinal data happen to be piled up on one end of the scale. The data are shown as the histogram in the top-right subpanel of the figure. (A real case of this sort of distribution will be provided later.) The generating parameter values are μ = 1.0 and σ = 2.5, with the thresholds equally spaced at θ₁ = 1.5, θ₂ = 2.5, θ₃ = 3.5, θ₄ = 4.5, θ₅ = 5.5, and θ₆ = 6.5. As you can see in the marginal posterior distributions of μ and σ, the Bayesian estimation very accurately recovers the generating parameters.

The top-right subpanel of Figure 23.2 superimposes posterior predictive probabilities of the outcomes. At each outcome value, a dot plots the median posterior predictive probability and a vertical segment indicates the 95% HDI of posterior predictive probabilities. A key aspect of this thresholded cumulative-normal model is that it accurately describes the ordinal data distribution.

The posterior distribution of an effect size is shown in the second row, right column of Figure 23.2. Effect size for a single group must be defined with respect to a comparison value C; for purposes of illustration, the comparison value here is chosen to be C = 2.0. Effect size is defined as the difference between the mean and the comparison value, relative to the standard deviation: (μ ‒ C)/σ. (This type of effect size was introduced in the case of a single group of metric values, as displayed in Figure 16.3, p. 457.) The posterior distribution leaves no doubt that the effect size is different from zero. A similar conclusion is reached by examining the posterior distribution of μ relative to the comparison value C = 2: We see that credible values of μ are clearly different from 2. Later, we will see that a conventional t test, treating the ordinal data as if they were metric, comes to a very different conclusion.

The bottom panel of Figure 23.2 shows the posterior distribution on the thresholds. The plot does not show the marginal distributions of individual thresholds because the threshold values are strongly correlated with each other, and therefore, the marginal distributions can be misleading. $R$emember that, at each step in the MCMC chain, the parameter values are jointly credible. Consider what happens at a particular step in the chain. If θ₂ is randomly chosen to be a bit higher than usual, the upward shift of the threshold not only increases the likelihood of outcome 2 but also decreases the likelihood of outcome 3. To compensate and maintain the likelihood of outcome 3, threshold θ₃ also is chosen to be a bit higher than usual. This domino effect continues through all the thresholds. Consequently, if one threshold is higher than usual, all thresholds tend to be higher than usual, and if one threshold is lower than usual, then all thresholds tend to be lower than usual. In other words, the thresholds are strongly correlated across steps in the chain. To display this linkage of the thresholds, I have chosen to plot the threshold values against the mean of the thresholds. At a particular step s in the MCMC chain, the jointly credible threshold values are θ₁ (s),…, θ_K-1(s), and the mean threshold at that step in the chain is $\overline{\theta }\left(s\right)\text{?}=\text{?}{\displaystyle {\sum }_k^{K-1}{\theta }_k\left(s\right)/\left(K-1\right)}$. The threshold values are then plotted as the points $\langle {\theta }_1\left(s\right),\text{?}\overline{\theta }\left(s\right)\rangle ,\dots ,\langle {\theta }_{K-1}\left(s\right),\text{?}\overline{\theta }\left(s\right)\rangle$. Notice that these points are all at the same height on the vertical axis. At different steps in the chain, the mean threshold will be higher or lower. You can see in the bottom panel of Figure 23.2 that this plotting method produces distinct clusters of points for each threshold value. You can see that the clusters for θ₁ and θ_K−1 have fixed values at 1.5 and K ‒ 0.5, as they should. Each cluster has a vertical dashed line that indicates its mean. The mean values for the estimated thresholds are very near the values that generated the data. You can also see that the higher thresholds have greater uncertainty than the lower thresholds, which makes sense because there are fewer high data points than low data points.

A second example is shown in Figure 23.3. This example illustrates a case in which the thresholds are not equally spaced, such that the penultimate ordinal values are emphasized and the overall data distribution is bimodal. The caption of the figure reports the true values of the parameters. You can see from the posterior distribution that the true parameter values are accurately estimated, including the thresholds. Importantly, the posterior predictive distribution in the top-right subpanel accurately describes the bimodal distribution of the outcomes.

23.3.1 Implementation in JAGS

The model specification in JAGS is a simple extension of the specification for one group that was explained in the previous sections. The extension merely uses two means and two standard deviations (one for each group). It uses the same thresholds for both groups. As in the one-group case, the number of outcome levels is denoted nYlevels and the total number data points across both groups is denoted Ntotal. As was the case for two groups of metric data explored in Section 16.3 (p. 468), the group index of respondent i is denoted x[i] (and can have values 1 or 2). The JAGS model specification is virtually identical to the case of one group, except what was mu is now mu[x[i]] and what was sigma is now sigma[x[i]]:

model {
 for(i in 1 : Ntotal) {
 y[i] ~ dcat(pr[i,1:nYlevels])
 pr[i,1] <- pnorm(thresh[1], mu[x[i]], 1/sigma[x[i]]~2)
 for (k in 2:(nYlevels-1)) {
 pr[i,k] <- max(0, pnorm(thresh[k], mu[x[i]], 1/sigma[x[i]]~2)
 - pnorm(thresh[k−1], mu[x[i]], 1/sigma[x[i]] ~ 2))
 }
 pr[i,nYlevels] <- 1 - pnorm(thresh[nYlevels-1], mu[x[i]], 1/sigma[x[i]]~2)
 }
 for(j in 1:2){ # 2 groups
 mu[j] ~ dnorm((1+nYlevels)/2, 1/(nYlevels)~2)
 sigma[j] ~ dunif(nYlevels/1000, nYlevels*10)
 }
 for(kin 2:(nYlevels-2)) { # 1 and nYlevels-1 are fixed, not stochastic
 thresh[k] ~ dnorm(k+0.5, 1/2~2)
 }
}

The complete model specification and other functions are defined in Jags-Yord-Xnom2grp-MnormalHet.R, and the high-level script that calls the functions is named Jags-Yord-Xnom2grp-MnormalHet-Example.R.

23.3.2 Examples: Not funny

Figure 23.4 shows a case with artificial data to demonstrate accurate recovery of known parameters. In this case, the ordinal data emphasize the low end of the scale. The true parameters have different means but equal variances, and the thresholds are equally spaced. The figure caption provides exact values. The posterior distribution shows that the parameter values are accurately recovered, and the posterior predictive distribution accurately describes the distribution of data. However, the uncertainty is large enough that 95% HDI on the effect size includes zero.

When the data of Figure 23.4 are treated as if they were metric and submitted to NHST, the conclusions are different. In particular, a t test concludes that the means are significantly different and an F test concludes that the variances are also significantly different (see the figure caption for details). These conclusions are not to be trusted because their p values assume that the data are normally distributed, which is clearly not the case here.

Figure 23.5 shows a case with real ordinal data that hug one end of the response scale. The data are funniness ratings of two jokes.² For each of 25 jokes, people were asked “How funny did you find the joke?” on a rating scale from 1 to 7 with 1 marked “Not at all” and 7 marked “Very.” Figure 23.5 shows results from two of the jokes that were rated as relatively less funny than others.³ I won't report the jokes here because, after all, people did not think these particular jokes were very funny! As can be seen in the posterior distribution of Figure 23.5, Bayesian estimation shows clear differences in the underlying means and variances of the two jokes, while simultaneously describing the data distributions very well. However, if the data are treated as metric and submitted to NHST, no significant differences emerge, as detailed in the caption. Again, the conclusions from NHST are not to be trusted because the p values assume that the data are normally distributed, which is clearly not the case here.

23.4.1 Implementation in JAGS

The model of Figure 23.6 is easy to implement in JAGS (or Stan). In the program, we generalize to multiple metric predictors (instead of only one metric predictor). There are only two changes from the model specification of the previous applications. First, μ is defined as a linear function of the predictors, just as in the previous programs for multiple regression (linear, logistic, softmax, or conditional logistic). Second, the predictors are standardized to improve MCMC efficiency, again as in the programs for multiple regression. Here is the JAGS model specification:

# Standardize the data:
data {
 for(j in 1 : Nx) {
 xm[j] <- mean(x[,j])
 xsd [j] < - sd(x[,j])
 for(i in 1 : Ntotal) {
 zx[i,j] <- (x[i,j] - xm[j]) / xsd [j]
 }
 }
}
# Specify the model for standardized data:
model {
 for (i in 1 : Ntotal) {
 y[i] ~ dcat(pr[i,1:nYlevels])
 pr[i,1] <- pnorm(thresh[1], mu[i], 1/sigma~2)
 for (k in 2:(nYlevels-1)) {
 pr[i,k] <- max(0, pnorm(thresh[k], mu[i], 1/sigma~2)
 - pnorm(thresh[k−1], mu[i], 1/sigma~2))
 }
 pr[i,nYlevels] <- 1 - pnorm(thresh[nYlevels-1], mu[i], 1/sigma~2)
 mu[i] <- z beta0 + sum(zbeta[1 : Nx] * zx[i, 1 : Nx])
 }
 # Priors vague on anchored ordinal scale:
 zbeta0 ~ dnorm((1+nYlevels)/2, 1/(nYlevels)~2)
 for(j in 1 : Nx) {
 zbeta[j] ~ dnorm(0, 1/(nYlevels)~2)
 }
 zsigma ~ dunif(nYlevels/1000, nYlevels*10)
 for(kin 2:(nYlevels-2)) { # 1 and nYlevels-1 are fixed
 thresh[k] ~ dnorm(k+0.5, 1/2^2)
 }
 # Transform to original scale:
 beta[1: Nx] <- (zbeta[1 : Nx] / xsd[1 : Nx])
 beta0 < - zbeta0 - sum(zbeta[1 : Nx] * xm[1 : Nx] / x sd[1 : Nx])
 sigma <- zsigma
}

You will have noticed in the code, above, that the predicted values were not standardized. Why? Because the predicted values are ordinal and do not have meaningful metric intervals between them. The ordinal values are just ordered categories, with probabilities described by the categorical distribution dcat in JAGS, which is a distribution over integer indices. Because the predicted values are not standardized, the transformation from standardized parameters to original-scale parameters is exactly like what was used for logistic regression, back in Equation 21.1 (p. 625). Full details are in the files Jags-Yord-XmetMulti-Mnormal.R and Jags-Yord-XmetMulti-Mnormal-Example.R.

23.4.2 Example: Happiness and money

We start with examples that have a single metric predictor. First, a case with fictitious data generated by known parameter values, which demonstrates that Bayesian estimation accurately recovers the parameter values. Figure 23.7 shows the data in its upper panel. The ordinal predicted values range from 1 to 7 in this case. You can see that when the predictor value (x) is small, there are a preponderance of small y values, and when the predictor value is large, there are mostly large y values. The parameter values that were used to generated the data are specified in the figure caption.

Figure 23.7 displays aspects of the posterior predictive distribution superimposed on the data. A smattering of credible regression lines is displayed. It must be remembered, however, that the regression lines refer to the underlying metric predicted variable, not to the ordinal predicted variable. Thus, the regression lines are merely suggestive and should be used to get a visual impression of the uncertainty in the slope and intercept. More directly pertinent are the posterior predicted probabilities of each outcome at selected x values marked by vertical lines. At the selected x values, the posterior predicted probabilities of the outcomes are plotted as horizontal bars “rising” leftward away from the vertical baseline. To see the probability distribution, it is easiest to tilt your head to the left, and imagine the bars as if marking the probabilities of intervals under the normal curve. In this example, when x is small, the bars on small values of y are tall relative to the bars on large values of y. And when x is large, the bars on small values of y are short relative to the bars on large values of y.

The main point of the example in Figure 23.7, other than familiarizing you with the graphical conventions, is that Bayesian estimation accurately recovers the true parameter values that generated the data. This is in contrast with estimates derived by treating the ordinal data as if they were metric. The least-squares estimates of the slope, intercept, and residual standard deviation are reported in the caption of the figure, and they are much too small in this case.

Our second example involving a single metric predictor uses real data from the Chinese household income project of 2002, which surveyed personal income and other aspects of people in urban and rural areas of the People's $R$ epublic of China (Shi, 2009). In particular, one survey item was a person's total household assets, measured in Chinese yuan (which had an exchange rate in 2002 of about 8.27 yuan per 1.00 US dollar). Another survey item asked, “Generally speaking, do you feel happy?” with response options of 1 = not happy at all, 2 = not very happy, 3 = so-so, 4 = happy, and 5 = very happy (and “don't know,” which was given just 1% of the time and is excluded from this analysis).⁴ There were 6,835 data points, as plotted in Figure 23.8.

One thing to notice from the data in Figure 23.8 is that most happiness ratings are 4, with the next most prominent being 3. One way that this preponderance of “4” ratings is handled by the parameter estimates is by setting the thresholds for categories 2 and 3 to relatively low values, as shown in the bottom panel. This makes the interval for 4 relatively wide.

The posterior distribution in Figure 23.8 indicates that the mean happiness rating clearly increases as assets increase. But it's not a huge increase: If assets go up by 82,770 yuan (the equivalent of 10,000 US dollars in 2002), then mean happiness increases 0.34 points (on the underlying metric scale). Moreover, at any given level of assets, there is large variability in happiness ratings: The marginal posterior distribution of σ has a mode of about 0.85 (on a response scale that ranges from 1 to 5). Note that the value of σ is estimated with high certainty, in that its 95% HDI extends only from about 0.83 to 0.88. This is worth reiterating to be sure the meaning of σ is clear: σ indicates variance in the data, and it is big, but the posterior estimate of the variance is also precise because the sample size is large.

If we were to treat the ordinal ratings as if they were metric values, the least-squares estimates of the regression parameters are a bit less than the Bayesian estimates, as reported in the caption of Figure 23.8, despite this being a case in which the data fall mostly in the middle of the response scale. For an example of treating these happiness data as metric, see Jiang, Lu, and Sato (2009), who defended the practice by citing comparisons reported by Ferrer-i-Carbonell and Frijters (2004).

The general trends indicated by this analysis are typical findings in many studies. But the analysis given here is meant only for pedagogical purposes, not for drawing strong conclusions about the relation of happiness to money. In particular, this example used linear trend for simplicity, and more sophisticated analyses (perhaps with other data) could examine nonlinear trends. Of course, there are many influences on happiness other than money, and the link between money and happiness can be mediated by other factors (e.g., Blanchflower & Oswald, 2004; Johnson & Krueger, 2006, and references cited therein).

23.4.3 Example: Movies—They don't make 'em like they used to

We now consider a case with two predictors. A movie critic rated many movies on a 1-7 scale⁵ and the data about the movies included their length (i.e., duration in minutes) and year of release. These data were assembled by Moore (2006) and come from reviews by Maltin (1996). It's not necessary that there should be any relation between length or year of release and the reviewer's ratings, linear or otherwise. Nevertheless, we can apply thresholded cumulative-normal linear regression to explore possible trends.

Figure 23.9 shows the data and results of the analysis. The upper panel of Figure 23.9 shows a scatter plot in which each point is a movie, with each datum plotted as a numeral 1-7 that indicates the movie's rating, positioned on axes that indicate the length and year of the movie. You can see that there is considerable variation in the ratings that is not accounted for by the predictors, because movies of many different ratings are intermingled among each other.

Despite the “noise” in the data, the analysis reveals credibly nonzero relationships between the predictors and the reviewer's ratings. The marginal posterior distribution on the regression coefficients indicates that ratings tend to rise as length increases, but ratings tend to decline as year of release increases. In other words, the reviewer tends to like old, long movies more than recent, short movies.

To get a visual impression of the linear trend in the data, Figure 23.9 plots a few credible threshold lines. At a step in the MCMC chain, consider that jointly credible thresholds and regression coefficients. We can solve for the 〈 x₁, x₂〉 loci that fall exactly at threshold θ_k by starting with the definition of the predicted value μ and setting it equal to θ_k:

$\begin{array}{ll}\mu \hfill & ={\theta }_k\hfill \\ \hfill & ={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2\hfill \end{array}$

hence,

$x_2\text{?}=\text{?}\left(\frac{{\theta }_k-{\beta }_0}{{\beta }_2}\right)\text{?}+\text{?}\left(\frac{-{\beta }_1}{{\beta }_2}\right)\text{?}{x}_1$

The lines determined by Equation 23.5 are plotted in Figure 23.9 at a few different steps in the MCMC chain, with thresholds from the same step plotted with lines of the same type (solid, dashed, or dotted). Threshold lines from the same step in the chain must be parallel because the regression coefficients are constant at that step, but are different at another step. The threshold lines in Figure 23.9 are level contours on the underlying metric planar surface, and the lines reveal that the ratings increase toward the top left, that is, as x₁ decreases and x₂ increases.

Curiously, however, the extreme thresholds fall beyond where any data occur. For instance, the threshold line for transitioning from 6 to 7 appears at the extreme upper left of the plotting region and there are no data points at all in the region corresponding to 7 above the 6-to-7 threshold. The analogous remarks apply to the lower-right corner: There are no data points in the region corresponding to 1, below the 1-to-2 threshold. How can this be? The next section explains.

23.4.4 Why are some thresholds outside the data?

To motivate an understanding of the thresholds, consider the artificial data presented in Figure 23.10. These data imitate the movie ratings, but with two key differences. First and foremost, the artificial data have much smaller noise, with σ = 0.20 as opposed to σ ≈ 1.25 in the real data. Second, the artificial data have points that span the entire range of both predictors, unlike the real data which had points mostly in the central region.

Bayesian estimation was applied to the artificial movie data in the same way as the real movie ratings. Notice that the posterior predictive thresholds plotted in Figure 23.10 very clearly cleave parallel regions of distinct ordinal values. The example demonstrates that the threshold lines do make intuitive sense when there is low noise and a broad range of data.

This leaves us needing an intuitive explanation for why the thresholds do not seem so clear for the real data. Figure 23.11 is designed to drive your intuition. The figure shows a single predictor on the horizontal axis of the panels, with the probability of outcome on the vertical axis. Each outcome is marked by a different line style (solid, dashed, dotted, and so on). All that changes across panels is the amount of noise. In the top panel, the noise is small, with σ = 0.1. Notice in this case that each outcome clearly dominates its corresponding interval. For example, outcome 2 occurs with nearly 100% probability between θ₁ = 1.5 and θ₂ = 2.5, and outcome 2 rarely occurs outside that interval. The next panel down has a bit more noise, with σ = 0.5. You can see that each outcome has maximum probability within its corresponding interval, but there is considerable smearing of outcomes into adjacent intervals. This smearing is worse in the next panel where σ = 1.0. When the noise is large, as in the bottom panel where σ = 2.0, there is so much “smearing” of outcomes across thresholds that the most probable outcome within intervals is not necessarily the outcome centered on that interval. For example, between θ₁ = 1.5 and θ₂ = 2.5 the most probable outcome is 1, not 2. To understand why this happens, imagine the normal distribution (from Figure 23.6) centered on θ₁ = 1.5. A full 50% of the normal falls below θ₁ = 1.5, so the probability of outcome 1 is 50%, but the probability of outcome 2 is much less because only a small portion of the normal distribution falls within the interval from θ₁ = 1.5to θ₂ = 2.5.

Suppose that the noise is large and that the range of the predictor values in the data is narrow relative to the thresholds. For example, suppose that the predictor values in the data go only from 2 to 4 in the bottom panel of Figure 23.11. Because of the smearing of outcome probabilities, the data would probably have many outcomes of 1 and 5 despite the fact that the predictor values never are less than θ₁ or greater than θ₄. It is this situation that occurs for the real movie data: The most credible parameter values are a large value for σ with the extreme thresholds set a bit outside the predictor values that actually occur in the data.

The preceding discussion has referred to σ as “noise” merely for linguistic ease. Calling the outcomes “noisy” does not mean the underlying generator ofthe outcomes is inherently wildly random. The “noise” is merely variation in the outcome that cannot be accounted for by the particular model we have chosen with the particular predictors we have chosen. A different model and/or different predictors might account for the outcomes well with little residual noise. In this sense, the noise is in the model, not in the data.