Consequences of analysing ordinal data with metric methods

If the distances between the values of a rating scale cannot be assumed to be constant, then even the interpretation of a simple mean breaks down. Figure 1 shows how a 7-point proficiency scale may be mentally represented if participants avoid extreme responses and perceive crossing the midpoint as particularly impactful. In this case, the average proficiency expressed by the hypothetical responses {‘2’, ‘3’, ‘4’} is not identical to that expressed by {‘3’, ‘3’, ‘3’}, but actually reflects greater average proficiency. From this we see that when equidistance is violated, the very same metric mean (here, 3) may not express the same underlying quantity; rather, what the mean signifies depends on the particular distribution of responses. This lack of consistency between an underlying quantity and the metric mean naturally extends to effects: that is, to differences between means and regression slopes. Thus, at the very least, estimating ordinal effects with metric methods brings about a serious problem of interpretability.

Figure 1. The hypothetical mental representation of a 7-point scale with unequally-spaced intervals.

Importantly, as demonstrated by Liddell and Kruschke (Reference Liddell and Kruschke2018) with both simulated and real data, metric methods may produce serious inferential errors when applied to rating scales: namely, the detection of effects that do not exist; the failure to detect effects that do exist; and distortions of effect-size estimates. It is even possible for differences to be inverted: that is, for effects in one direction in the underlying construct (e.g., group A has larger proficiency than group B) to turn into effects in the opposite direction when mapped to metric means (e.g., group B > group A) (this can happen when the underlying distributions have different variances; see Liddell & Kruschke, Reference Liddell and Kruschke2018).

A (very) wrong model for ordinal data

The negative consequences described above arise because, in addition to unequal distances between responses, other properties of ordinal data are also fundamentally incompatible with the statistical models commonly used to analyse them. To see how this is the case, consider the following dataset.Footnote ¹ Schlenter (Reference Schlenter2019) conducted a study in which bilingual speakers listened to 56 German sentences with canonical and non-canonical orders of thematic arguments. In order to determine the acceptability of the two variants, these sentences were first rated by a group of 42 native speakers, using a 7-point scale (‘1 = not acceptable’…‘7 = very acceptable’).

A typical analysis would involve computing means in the two conditions (canonical: 6.35, non-canonical: 5.76), and then fitting a linear model: for example, a regression with condition as predictor. The results indicate that non-canonical sentences are less acceptable (b=−0.59, t=−9.64). As described above, it is hard to interpret such differences: it is not clear what it means for a construction to have 0.59 ‘units of acceptability’ less, and moreover, if equidistance cannot be assumed, the interpretation of this difference depends on the distribution of responses.

In addition, it can be shown that the assumptions of linear regression are severely violated. This can be assessed with a predictive check, in which instances of ‘predicted data’ are generated from the fitted statistical model and then compared to the observed data that the models were fitted on; this check is shown in Figure 2, panel a. A comparison of the lighter lines (predicted) against the darker line (observed) reveals that the linear model: (a) grossly underestimates the proportions of ‘6's and ‘7's; (b) predicts impossible non-integer responses (e.g., 6.27); and (c) predicts a sizeable proportion of responses outside the 7-point scale.

Figure 2. Predictive checks for a mixed-effects linear regression model (panel a), and for an ordinal model with flexible threshold locations (panel b). These assessments of model quality involve comparing the probability/counts of different responses (given in a 7-point acceptability scale) in data predicted by the fitted models (lighter lines in panel a, and dots in panel b) against the observed data that the models were fitted on (darker line in panel a, and bars in panel b). The two panels depict predicted and observed data differently because the linear model treats responses as continuous and the ordinal model appropriately treats them as discrete.

The mismatches happen because linear models are inherently continuous and assume that errors are normally distributed. In other words, the best inference we can draw from this model is that the observed data in each condition comes from a normally-distributed population. However, because these assumptions yield impossible data, we know that this inference is fundamentally flawed. To be sure, the assumptions of statistical models are never fully satisfied (“all models are wrong”; Box, Reference Box1976). However, linear models with normally-distributed errors are particularly wrong when applied to ordinal data.

A pervasive problem in bilingualism

How widespread is the analysis of ordinal data with metric methods in bilingualism research? Although there are examples of the application of appropriate methods (e.g., Kissling, Reference Kissling2018; Tare et al., Reference Tare, Golonka, Lancaster, Bonilla, Doughty, Belnap, Jackson, Sanz and Morales-Front2018), this is far from common. We quantified this tendency by searching for all articles published in Bilingualism: Language and Cognition in 2019–2020, containing the words ‘rating’ or related words. Forty-six articles analysed ordinal data, mostly proficiency ratings. Of these, only 3 appropriately summarised it with counts instead of means, but skipped inferential statistics or reported inappropriate ones; only 2 used a method that can be applied to ordinal data (Spearman's correlation), but still computed means. No article used both descriptive and inferential statistics that respect the properties of ordinal data. The remainder of this article describes and exemplifies a solution to this pervasive problem.

Modelling unequal distances

The modelling of unequal distances can be better understood by comparing models with and without equidistance. We will make use of a simulated dataset of proficiency ratings, consisting of the (randomly generated) responses of 1,000 participants, given on a 7-point scale. The dataset aims to represent a heterogeneous group of L2 speakers with a mean level at the midpoint of the scale. Importantly, the data was generated assuming a representation similar to that in Figure 1: that is, with unequal distances between values. The counts of each response are shown in Figure 4 below (as bars), and suggest an exaggerated tendency for choosing the midpoint of the scale.

Figure 4. Predictive checks for ordinal models with flexible thresholds (full dots) and equidistant thresholds (empty triangles). These assessments of model quality involve comparing the probability/counts of different responses (given in a 7-point acceptability scale) in data predicted by the fitted models (dots and triangles) against the observed data that the models were fitted on (shown as bars).

We use the brm() function to fit two (Bayesian) thresholded-cumulative models to this dataset. We use a probit link function, which means that we assume a normally-distributed latent variable. In the first model (m.equidistant), the distance between each threshold is assumed to be the same:

In the second model (m.flexible), we allow consecutive thresholds to be estimated at any distance from one another:

A summary of the flexible model is shown in Table 1 with estimates of each threshold's position and associated 95% credible intervals (this is the range within which a parameter falls with 95% probability). The estimates are expressed in standard deviation (SD) units; for example, the first threshold (labelled Intercept[1]) is 1.92 SDs below the mean of the latent distribution. The thresholds are more easily interpreted when visualised, as in Figure 3. In order to capture the observed proportions of each response (‘1’ to ‘7’), the thresholds are estimated to be closer together (e.g., second and third thresholds) or further apart.

Table 1 Summary of a model with flexible thresholds, with estimates and 95% credible intervals for the position of each threshold. For comparison, estimates from the equidistant model are also displayed, as well as the ‘true’ population values which were used to generate the data.

The arrows at the bottom of Figure 3 depict the ‘true’ population values, which in this case (and unlike with real data) are already known. The models can be assessed by whether they can recover these underlying parameters, and it can be seen that the flexible model does a very good job. The equidistant model fares worse (see Table 1), especially for the third and fourth thresholds, which are estimated as much closer together than the population values.

One way of assessing a model's goodness-of-fit is the previously discussed predictive check, in which data samples are predicted from the fitted model and compared to the observed data. For Bayesian models, these checks can be easily conducted with the pp_check() function; we use bar plots because the models are discrete rather than continuous:

The results are displayed in Figure 4 (in a single plot), and show that the equidistant model severely underestimates the number of ‘4’ responses and overestimates the number of ‘3’ and ‘5’ responses. By contrast, the predictions of the flexible model are very close to the observed data. Another way of assessing models is to formally compare them: for example, with cross-validation (Bürkner, Reference Bürkner2017; Vasishth et al., Reference Vasishth, Nicenboim, Beckman, Li and Kong2018); this also shows an advantage for the flexible model (see Appendix S1 in the Supplementary Material).

From the various comparisons, we would conclude that flexible thresholds are necessary to appropriately model this dataset, and by extension, that the psychological distances between response values are not constant across the scale. In particular, the estimated threshold locations show that ‘5’ responses actually express greater underlying proficiency than a metric 5 (and ‘3’ responses express lower proficiency than a 3), so that the difference between a ‘3’-participant and a ‘5’-participant should be interpreted as particularly large. Such inferences about how the underlying construct relates to the observed responses are missed when metric models are used, and are an important advantage of ordinal models.

Effects of categorical predictors

To illustrate how the effects of predictors are estimated in ordinal models, we go back to Schlenter (Reference Schlenter2019)'s dataset, in which canonical and non-canonical sentences were rated on a 7-point acceptability scale. In thresholded-cumulative models, effects are estimated in the same way as in linear regression (i.e., by ‘shifting’ the mean of a normal distribution), but they take place at the latent variable, rather than on metric responses.

We fit a (mixed-effects) thresholded-cumulative model with condition (canonical, non-canonical) as predictor. Because each participant and item are associated with multiple responses, participant and item are included as random effects:

The effect of condition on acceptability is estimated as -0.68 [-0.83, -0.54] (see model summary in Table S1 in the Supplementary Material). The effect is expressed in SD units, which in this case can be interpreted as a standardized effect size, similar to Cohen's d (Reference Cohen1988; Glass, McGaw, & Smith, Reference Glass, McGaw and Smith1981). Both the effect's magnitude and its narrow credible interval indicate a substantial difference between conditions, with lower acceptability for non-canonical sentences.

An important aid to interpretation is the examination of the model's conditional effects: that is, the predicted proportions of responses in each condition:

The predictions indicate that canonical sentences have very high acceptability, with a large proportion of ‘7 = very acceptable’ responses (64% [52, 75]); also see Figure S1 in the Supplementary Material. The proportion of ‘6's is lower but still substantial (28% [20, 34]), and other responses are less frequent (‘5’ responses: 6% [3, 10]). In turn, the lower acceptability of non-canonical sentences is expressed by a much smaller proportion of ‘7's (37% [26, 49]), and by more lower-acceptability responses (‘6’ responses: 38% [34, 42]; ‘5's: 14% [10, 19]).

Such predictions are finer than those obtained from the linear model above (cf. ‘0.59 less acceptability’), because they are expressed as probabilities of discrete responses and not in the inappropriate metric scale. Note also that these proportions are not calculated directly from the data. Rather, they constitute better, more generalisable inferences, because they come from a model that: (a) estimates the effect of condition across all responses; (b) takes into account the whole structure of the data in terms of its participants and items; and (c) can potentially include other sources of information, like adjustments for covariates.

Finally, we assess the model's goodness-of-fit with a predictive check, displayed in Figure 2, panel b:

Whereas the predictions of the linear model seriously mismatched the data (see above), the ordinal model fares remarkably well, and readily accommodates both the discreteness of responses and their non-normality.

Effects of continuous predictors

The usefulness of examining the model's conditional effects is particularly clear in the case of continuous predictors. To illustrate, we will make use of a dataset collected by Puebla (Reference Puebla2016) in which 55 second language German speakers self-rated their speaking proficiency on a 7-point scale (responses covered only the four highest categories, ‘functional’, ‘good’, ‘very good’, and ‘nativelike’). We will estimate the effect of age of acquisition (AoA) on self-ratings of proficiency (see, e.g., Hakuta et al., Reference Hakuta, Bialystok and Wiley2003).

Responses were coded as categories (e.g., ‘good’), so we first establish their order, and then fit a thresholded-cumulative model with AoA as a continuous predictor:

The effect of AoA on the latent proficiency variable is estimated as -0.15 [-0.22, -0.09] per year (see model summary in Table S2 in the Supplementary Material). Thus, the effect across the whole AoA range (4–24 years) is 3 SDs on the latent scale, indicating a large difference between early and late bilinguals, with lower proficiency at later AoAs.

We plot the conditional effects of this model (Figure 5):Footnote ³

Figure 5. Effect of AoA on self-rated speaking proficiency. Separate predictions are plotted for each response category (i.e., different probabilities for ‘functional’, ‘good’, ‘very good’, ‘nativelike’). Given that responses are mutually exclusive, their predicted proportions add up to 100% at each AoA. Note that even though the model predicts an AoA effect for each response category, all predictions arise from a single (linear) effect of AoA on the latent proficiency variable. Shaded bands indicate 95% credible intervals. Data from Puebla (Reference Puebla2016).

At the earlier AoAs there is a predominance of ‘nativelike’ responses. As AoA increases, ‘nativelike’ responses decline sharply and there is an increase in the proportions of the lower-proficiency categories (i.e., ‘nativelike’ responses are progressively replaced by ‘very good’ and ‘good’ responses). By an AoA of 13, all three top proficiency choices are likely responses, but there is a predominance of ‘very good's. As AoA increases further (after puberty), the most common response becomes ‘good’ and some ‘functional’ responses start emerging.

These inferences are much more informative than what could be afforded by a linear model, because predictions (and their uncertainty) are appropriately expressed in terms of the different response categories.