# Why are emmip( "response") y axis numbers not probabilities for ordinal regression?

I used emmeans functions (with help from this site) to obtain pairwise comparisons for different levels of variables in a model with interactions. Interpreting an interaction term in the context of study

The model int was estimated withlrm() from the rms package. The plot below was plotted with emmip(int, ses~fam, type = "response", CIs = T, levels = 0.95) and I also tested type = "linear.predictor". These two options returned exactly the same plot.

Using the emmip function with response mode to plot the effects, I expected to get the probabilities for the levels of ordinal response, but the y axis has values between-3 and -1.

I also printed the plotted means. The means and their were called with emmeans(model, specs = pairwise ~ ses*fam ) and I also extracted \$contrasts adding %>%confint(adjust = "tukey").

I cannot work out what these numbers are and why they are negative. What do they reflect and how do they relate to the pairwise comparisons from emmeans? I expected them to be probabilities for each outcome adding to 1.

In response to the comments: So the negative numbers come from the "Latent Variable" that underlies the ordinal regression model. A vignette for emmeans describes different ways ("modes") to express the results. With respect to rms models, it says:

rms models have an additional mode. With mode = "middle" (this is the default), the middle intercept is used, comparable to the default for rms::Predict(). This is quite similar in concept to mode = "latent", where all intercepts are averaged together.

As the intercepts of my model are:

 cat1|cat2 2.198193       cat2|cat3   3.158535


Could I ask how these means relate to the contrasts that can be obtained for the combinations of all the levels of ses and family type (from emmeans contrasts)?

For example in contrasts there is a significant difference between low ses family and and low ses family c, but the means of predicted response are not different as their confidence intervals overlap.

And is it possible and does it make sense to "remove" the average intercepts from the plot?

• I removed all the codes and will post in StackExchange. I asked an additional question about the difference between the meaning of differences from marginal means and contracts, with a kind request for an explanation.
– Milo
Commented May 13, 2023 at 16:45
• It would help if you could edit this question to show the specific code you used for the model (e.g., polr() versus lrm() versus orm()) and the call you made to emmeans() to get the table of contrasts that you show.
– EdM
Commented May 13, 2023 at 17:02
• I added the functions. I also read the discussion, Based on the linked discussion my understanding is that I would not be able to make any assumptions as the sample sizes for groups are very different, especially for the overlapping CI for two means. This is the first time I am trying interpreting more than the main model output for interaction so I know my questions are basic and appreciate the answers a lot
– Milo
Commented May 13, 2023 at 17:54
• I've updated the answer in a way that I should address your questions about the underlying statistical issues. It might also explain the problem that you had with trying to display probability estimates.
– EdM
Commented May 14, 2023 at 13:42

There are several ways to report the estimated outcomes from an ordinal regression model. With the emmeans package you can choose among them with a mode argument submitted to emmeans() or ref_grid().

Latent variable

According to a vignette page, the default is mode = "latent", and

With mode = "latent", the reference-grid predictions are made on the scale of the latent variable implied by the model. The scale and location of this latent variable are arbitrary...

The results you show are based on that underlying default "latent variable" model, an assumed continuous response (arbitrarily centered and scaled) that is converted to ordered categories as its value passes associated thresholds. The idea is illustrated in Figure 6.4 of these course notes. That "latent variable" model is fundamental and can be used to obtain other types of outcome estimates.

As you discovered, the "latent variable" for a scenario here is just the linear predictor $$X' \beta$$, where $$X$$ is the vector of predictor values and $$\beta$$ is the corresponding vector of coefficient estimates. When you ask emmeans for estimates of marginal means, it takes the grid of values of variables that you specify (in this case, the combinations of the predictors ses and fam) and returns the model's estimates for those values (by default, at the average values of variables that you don't specify).

As the vignette notes, however, "The scale and location of this latent variable are arbitrary." In the parameterization used by the lrm() and orm() functions in the rms package, the linear predictor is associated with the class probabilities:

$$\Pr(Y \ge j|X)=\frac{1}{1+\exp\left(-(\alpha_j + X' \beta)\right)}$$

where you have classes of outcome $$Y$$ labeled with increasing values of $$j$$ and the $$\alpha_j$$ are the corresponding intercepts estimated by the model. Clearly, you can shift the linear predictor $$X' \beta$$ by any offset and get the same outcome probability estimates if you make corresponding changes to the values of the $$\alpha_j$$ intercepts.

From that perspective, you can display the latent-variable estimates with any offset shift that you want. The mode="middle" choice will tend to center the estimates around 0, if you want to display the results in terms of the latent variable. It doesn't matter that the model estimates you displayed have negative values. If you ask for probability estimates, the corresponding values of $$\alpha_j$$ (your two intercept estimates) will put everything together properly.

Probability estimates

If you want outcome-class probabilities, you need to specify a different mode. There is a mode="prob" argument that can report class probability estimates for supported model types.

There's a "gotcha" with that mode, however, explained in the section of the vignette on multinomial responses, of which ordinal models are a special case:

Please note that, because the probabilities sum to 1 (and the latent values sum to 0) over the multivariate-response levels, all sensible results from emmeans() must involve that response as one of the factors. For example, if resp is a response with k levels, emmeans(model, ~ resp | trt) will yield the estimated multinomial distribution for each trt; but emmeans(model, ~ trt) will just yield the average probability of 1/k for each trt.

If you want to display class probabilities, you need to include the outcome among the variables that you specify to emmeans(). For example, to get a plot of outcome-class probabilities as a function of both ses and fam, you could specify:

emmip(int, outcome ~ fam|ses, mode = "prob", CIs = TRUE)


That will display estimates of outcome probabilities as a function of fam, in a separate facet for each ses level. The outcome probabilities for each combination of fam and ses will sum to 0, as expected.

A warning about implementation: in the older version of emmeans I'm using on this computer (emmeans_1.6.2-1), the mode="prob" argument didn't seem to work on an lrm() model, while it worked OK with an orm() model (which provided the same coefficient estimates). I don't know if that's still the case with newer versions.

Contrasts and statistical significance

A useful estimate of the "significance" of your predictors with respect to outcome includes all the terms involving each predictor. If you use a model generated by the rms package, the anova() function applied to the model provides a very useful display of overall significance and the significance of interaction and nonlinear terms in the model, via Wald tests.

At the other extreme of complexity, the pairwise comparisons among all predictor combinations take into account the corresponding coefficient estimates $$\beta$$, their variances, and their covariances. It applies the formula for the variance of a weighted sum of variables to determine if any of those coefficient combinations differs significantly from 0, with multiple-comparison corrections as specified.

Of course, the correction for multiple comparisons becomes more restrictive as the number of comparisons increases. If there are particular comparisons of major interest, it can be better to restrict your analysis to those rather than evaluating all pairwise comparisons.

For those pairwise differences, any constant offset shift in the linear predictors don't matter, as the offset will cancel out when you take the differences. That's another reason not to worry about the negative values displayed for the latent variables.

The point estimates of both of the contrasts that you show agree with the differences of the corresponding emmean values in the table you show. The standard error of low fam_a - low fam_b, however, is more than twice as large as for low fam_a - low fam_c, presumably because of a small number of cases in low fam_b. Although the point estimate of low fam_a - low fam_b is larger than for low fam_a - low fam_c, the larger standard error means that the first contrast isn't significant at $$p < 0.05$$ while the second is.

A final warning: overlapping confidence intervals can be consistent with a significant difference between two scenarios. The emmeans FAQ emphasizes that point. This answer goes into detail with respect to t-tests in simple scenarios, in which non-overlap of 95% CI is approximately equivalent to $$p < 0.005$$ for the difference between means.