Visualizing the Effects of Proportional-Odds Logistic Regression

Proportional-odds logistic regression is often used to model an ordered categorical response. By “ordered”, we mean categories that have a natural ordering, such as “Disagree”, “Neutral”, “Agree”, or “Everyday”, “Some days”, “Rarely”, “Never”. For a primer on proportional-odds logistic regression, see our post, Fitting and Interpreting a Proportional Odds Model. In this post we demonstrate how to visualize a proportional-odds model in R.

To begin, we load the effects package. The effects package provides functions for visualizing regression models. This post is essentially a tutorial for using the effects package with proportional-odds models. We also load the car, MASS and splines packages for particular functions, which we’ll explain as we encounter them. If you don’t have the effects or car packages, uncomment the lines below and run them in R. MASS and splines are recommended packages that come with R.

> # install.packages("effects")
> # install.packages("car")

> library(effects)
> library(car)
> library(MASS)
> library(splines) 

To demonstrate how to visualize a proportional-odds model we’ll use data from the World Values Surveys (1995-1997) for Australia, Norway, Sweden, and the United States. This dataset, WVS, comes with the effects package. Once we load the effects package, the data is ready to access.

> head(WVS)
      poverty religion degree country age gender
1  Too Little      yes     no     USA  44   male
2 About Right      yes     no     USA  40 female
3  Too Little      yes     no     USA  36 female
4    Too Much      yes    yes     USA  25 female
5  Too Little      yes    yes     USA  39   male
6 About Right      yes     no     USA  80 female

The response variable of interest is poverty, which is a 3-level ordered categorical variable. It contains the answer to the question “Do you think that what the government is doing for people in poverty in this country is about the right amount, too much, or too little?” The answers are an ordered categorical variable with levels “Too Little”, “About Right”, and “Too Much”. The other variables will serve as our predictors. These include country, gender, religion (belong to a religion?), education (hold a university degree?), and age in years. The data contains 5381 records.

Before we get started we need to note this example comes from an article on the effects package in the Journal of Statistical Software by John Fox and Jangman Hong, the authors of the effects package. You should definitely take the time to read through that article and cite it if you plan to use the effects package for your own research. What we seek to do in this blog post is elaborate on the example and provide some additional details.

Before we can visualize a proportional odds model we need to fit it. For this we use the polr function from the MASS package. The first model we fit models poverty as a function of country interacted with gender, religion, degree and age. The interaction allows the effects of the predictors to vary with each country.

> wvs.1 <- polr(poverty ~ country*(gender + religion + degree + age), data = WVS)
> summary(wvs.1)

Call:
polr(formula = poverty ~ country * (gender + religion + degree + 
    age), data = WVS)

Coefficients:
                               Value Std. Error  t value
countryNorway              0.5308176   0.286989  1.84961
countrySweden              0.5446552   0.546029  0.99748
countryUSA                -0.0347317   0.248059 -0.14001
gendermale                 0.0696120   0.090212  0.77165
religionyes                0.0094685   0.112476  0.08418
degreeyes                 -0.1242920   0.167603 -0.74158
age                        0.0155849   0.002597  6.00185
countryNorway:gendermale   0.1873611   0.144503  1.29659
countrySweden:gendermale   0.0563508   0.154414  0.36493
countryUSA:gendermale      0.2119735   0.139513  1.51938
countryNorway:religionyes -0.2186724   0.216256 -1.01118
countrySweden:religionyes -0.8789724   0.513263 -1.71252
countryUSA:religionyes     0.6002277   0.174433  3.44101
countryNorway:degreeyes    0.0558595   0.208202  0.26829
countrySweden:degreeyes    0.6281743   0.214295  2.93136
countryUSA:degreeyes       0.3030866   0.206394  1.46848
countryNorway:age         -0.0157142   0.004367 -3.59846
countrySweden:age         -0.0092122   0.004657 -1.97826
countryUSA:age             0.0005419   0.003975  0.13635

Intercepts:
                       Value   Std. Error t value
Too Little|About Right  0.7161  0.1535     4.6644
About Right|Too Much    2.5355  0.1578    16.0666

Residual Deviance: 10347.07 
AIC: 10389.07 

The summary output is imposing. In addition to 19 coefficients we have 2 intercepts. Larger coefficients with large t-values are indicative of important predictors, but with so many interactions it’s hard to to see what’s happening or what the model “says”. To evaluate whether the interactions are significant, we use the Anova function from the car package. By default the Anova function returns Type II tests, which tests each term after all others, save interactions. (The base R anova function performs Type I tests which tests each term sequentially.)

> Anova(wvs.1)
Analysis of Deviance Table (Type II tests)

Response: poverty
                 LR Chisq Df Pr(>Chisq)    
country           250.881  3  < 2.2e-16 ***
gender             10.749  1  0.0010435 ** 
religion            4.132  1  0.0420698 *  
degree              4.284  1  0.0384725 *  
age                49.950  1  1.577e-12 ***
country:gender      3.049  3  0.3841657    
country:religion   21.143  3  9.833e-05 ***
country:degree     12.861  3  0.0049476 ** 
country:age        17.529  3  0.0005501 ***

The Anova result shows all interactions except country:gender are significant. But what do the interactions mean? How do, say, country and age interact?

This is where the effects package enters. The effects package allows us to easily create effect displays. What are effect displays? The documentation for the effects package explains it this way:

“To create an effect display, predictors in a term are allowed to range over their combinations of values, while other predictors in the model are held to typical values.”

In other words, we take our model and use it to calculate predicted values for various combinations of certain “focal” predictors while holding other predictors at fixed values. Then we plot our predicted values versus the “focal” predictors to see how the response changes. Let’s demonstrate.

The two primary functions are Effect and plot. Effect generates the predictions and plot creates the display. Let’s say we’re interested in the age and country interaction. We want to visualize how age affects views on poverty for each country. Since our model includes an interaction, which was significant, we expect to see different trajectories for each country.

The following code generates the predicted values. The first argument, focal.predictors, is where we list the predictors we’re interested in. Notice it requires a vector, which is why we use the c() function. The second argument is the fitted model.

> Effect(focal.predictors = c("age","country"), wvs.1)

age*country effect (probability) for Too Little
    country
age  Australia    Norway    Sweden       USA
  20 0.5958889 0.5632140 0.6496015 0.4330782
  30 0.5578683 0.5635318 0.6349615 0.3939898
  40 0.5191570 0.5638496 0.6200674 0.3562127
  50 0.4802144 0.5641674 0.6049438 0.3201442
  60 0.4415107 0.5644851 0.5896166 0.2861049
  70 0.4035049 0.5648027 0.5741134 0.2543308
  80 0.3666230 0.5651203 0.5584631 0.2249734
  90 0.3312402 0.5654379 0.5426958 0.1981045

age*country effect (probability) for About Right
    country
age  Australia    Norway    Sweden       USA
  20 0.3050532 0.3250955 0.2699789 0.3918455
  30 0.3282699 0.3249058 0.2797785 0.4064109
  40 0.3502854 0.3247160 0.2895692 0.4171736
  50 0.3704966 0.3245261 0.2993161 0.4237414
  60 0.3883075 0.3243362 0.3089823 0.4258700
  70 0.4031610 0.3241462 0.3185295 0.4234793
  80 0.4145714 0.3239561 0.3279182 0.4166592
  90 0.4221528 0.3237659 0.3371079 0.4056632

age*country effect (probability) for Too Much
    country
age   Australia    Norway     Sweden       USA
  20 0.09905788 0.1116905 0.08041952 0.1750762
  30 0.11386182 0.1115624 0.08526005 0.1995993
  40 0.13055755 0.1114343 0.09036331 0.2266137
  50 0.14928897 0.1113065 0.09574006 0.2561143
  60 0.17018180 0.1111787 0.10140107 0.2880252
  70 0.19333407 0.1110511 0.10735706 0.3221899
  80 0.21880555 0.1109236 0.11361867 0.3583674
  90 0.24660694 0.1107962 0.12019630 0.3962323

Notice the output lists three sections of probabilities corresponding to each level of the response. The first section lists predicted probabilities for answering “Too Little” for each country for ages ranging from 20 to 90 in increments of 10. The predicted probability a 20-year-old from the USA answers “Too Little” is about 0.43. The second section is for “About right”. The predicted probability a 20-year-old from the USA answers “About Right” is about 0.39. Finally, the third section is for “Too Much”. The predicted probability a 20-year-old from the USA answers “Too Much” is about 0.18.

This information is much easier to digest as an effect display. Simply insert the original Effect function call into the plot function to create the effect display.

> plot(Effect(focal.predictors = c("age","country"), wvs.1), rug = FALSE)

Now we see how the model “works”. For example, in the upper right plot, we see that in the USA, the probability of answering “Too Much” increases rather dramatically with age, while the probabilities for answering the same in Norway and Sweden stay low and constant. Likewise we see that the probability of USA respondents answering “Too Little” decreases with age while the probabilities for Norway and Sweden stay rather high and constant. The effect display shows us where the interactions are happening and to what degree. (Note: setting rug = FALSE turns off the rug plot. When set to TRUE, the default, the marginal distribution of the predictor is displayed on the x axis. In this case we don’t find it very helpful since we have so much data.)

Recall that we need to use all predictors to generate these plots. Country and age are the focal predictors, so they are varied. This means gender, religion and degree are held fixed. We can find out the values they’re fixed at by saving the result of the Effect function and viewing the model matrix.

> e.out <- Effect(focal.predictors = c("age","country"), wvs.1)
> e.out$model.matrix[1,c("gendermale","religionyes","degreeyes")]
 gendermale religionyes   degreeyes 
  0.4935886   0.8539305   0.2124140 

This says gender was set to 0.4935886, religion to 0.8539305, and degree to 0.2124140. In the original data these are indicator variables that take values of 0 or 1 corresponding to No and Yes. Does it make sense to plug in decimals instead of 0s or 1s? It does if you think of modeling a population that is about 49% men, 85% religious, and 21% with a college degree. In fact this describes our sample. By default, the effects package will take the mean of numeric variables that are held fixed. We can verify that’s what it’s doing:

> mean(WVS$gender=="male")
[1] 0.4935886
> mean(WVS$religion=="yes")
[1] 0.8539305
> mean(WVS$degree=="yes")
[1] 0.212414

If we want to change these values, we can use the given.values argument. It needs to be a named vector that uses the terms as listed in the output summary. For example, to create an effect plot for religious men without a college degree:

> e.out <- Effect(focal.predictors = c("country","age"), mod = wvs.1, 
                given.values = c(gendermale = 1, religionyes = 1, degreeyes = 0))
> plot(e.out, rug = FALSE)

We see that the overall “story” of the display does not change; the same changes are happening in the same plots. But the overall probabilities have increased by evidence of the y axis now topping out at 0.7.

We can also change the range and number of focal predictors using the xlevels argument. For example, we can set age to range from 20 to 80 in steps of 10. Notice it needs to be a named list.

> e.out <- Effect(focal.predictors = c("country","age"), mod = wvs.1,
                 xlevels = list(age = seq(20,80,10)))
> plot(e.out, rug = FALSE)

We can investigate other interactions using the same syntax. Below is the effect display for the religion and country interaction. In this case, age is no longer a focal predictor and is held fixed at its mean (45.04).

> plot(Effect(focal.predictors = c("religion","country"), wvs.1), rug = FALSE)

We notice, for example, that religious people in the USA have a higher probability of answering “Too Much” compared to their counterparts in the other countries surveyed. We also notice there is much more uncertainty about estimates in Sweden for people who answered “No” to the religion question. This is due to the small numbers of respondents in those categories, as we can see with the xtabs function.

> xtabs(~ religion + country + poverty, data = WVS, country == "Sweden", 
        drop.unused.levels = TRUE)
, , poverty = Too Little

        country
religion Sweden
     no       7
     yes    597

, , poverty = About Right

        country
religion Sweden
     no       5
     yes    369

, , poverty = Too Much

        country
religion Sweden
     no       3
     yes     22

Following the example in Fox’s article, let’s fit another model that relaxes the linearity assumption for age. We can do this by generating what’s called a basis matrix for natural cubic splines. Instead of fitting a regular polynomial such as age + age^2, we fit piecewise cubic polynomials over the range of age separated by a certain number of intervals, or knots. The ns function in the splines package makes this easy to do. Below we use it in the model formula and specify 4 knots. Harrell (2001) suggests 3-5 knots is usually a good choice (p. 23), so 4 seems wise in this case.

> wvs.2 <- polr(poverty ~ country*(gender + religion + degree + ns(age, 4)),data = WVS)
> summary(wvs.2)
> Anova(wvs.2)

Due to space considerations we don’t print the output of the summary and Anova functions. The model summary shows information for 31 coefficients and is very difficult to interpret. The Anova result is similar in substance to the first model, showing all interactions except country:gender significant. It’s worth noting the AIC of the second model is considerably lower than the first (10373.85 vs 10389.07), suggesting a better fit.

Let’s look at the country and age interaction while allowing age to range from 20 – 80:

> plot(Effect(focal.predictors = c("country","age"), mod = wvs.2, 
            xlevels = list(age = 20:80)), rug = FALSE)

We see that using a natural spline allows a nonlinear effect of age. For example we see the probability of answering “Too Little” in the USA decreases sharply from 20 to 30, increases from about age 30 to 45, and then decreases and levels out through age 80.

The effects package also allows us to create “stacked” effect displays for proportional-odds models. We do this by setting style="stacked" in the plot function.

> plot(Effect(focal.predictors = c("country","age"), mod = wvs.2, 
            xlevels = list(age = 20:80)), 
     rug = FALSE,
     style="stacked")

This plot is useful for allowing us to compare probabilities across the response categories. For example, in Norway and Sweden, people are most likely to answer “Too Little” regardless of age. The blue shaded regions dominate their graphs.

We can also create a “latent” version of the effect display. In this plot, the y axis is on the logit scale, which we interpret to be a latent, or hidden, scale from which the ordered categories are derived. We create it by setting latent = TRUE in the Effect function.

> plot(Effect(focal.predictors = c("country","age"), mod = wvs.2, 
            xlevels = list(age = 20:80),
            latent = TRUE),
     rug = FALSE,
     ylim = c(0,3.5))

This plot is useful when we’re more interested in classification than probability. The horizontal lines in the plots correspond to the intercepts in the summary output. We can think of these lines as threshholds that define where we crossover from one category to the next on the latent scale. The TL-AR line indicates the boundary between the “Too Little” and “About Right” categories. The AR-TM line indicates the boundary between the “About Right” and “Too Much” categories. Like the “stacked” effect display, we see that someone from Norway and Sweden would be expected to answer “Too Little” regardless of age, though the confidence ribbon indicates this expectation is far from certain, especially for older and younger respondents. On the other hand, most USA respondents are expected to answer “About Right”.

Let’s see how the “latent” plot changes when we set the non-focal predictors to college-educated, non-religious female.

> plot(Effect(focal.predictors = c("country","age"), mod = wvs.2, 
            xlevels = list(age = 20:80),
            given.values = c(gendermale = 0, religionyes = 0, degreeyes = 1),
            latent = TRUE),
     rug = FALSE,
     ylim = c(0,3.5))

Notice now that predicted classification for Sweden is “About Right” over the age range but with increased uncertainty. We also see increased chances of answering “Too Little” for certain age ranges in the USA.

We can add gender as a focal predictor to compare plots for males versus females:

> plot(Effect(focal.predictors = c("country","age","gender"), mod = wvs.2, 
            xlevels = list(age = 20:80),
            latent = TRUE),
     rug = FALSE,
     ylim = c(0,3.5))

Since we didn’t fit a 3-way interaction between country, gender and age, the trajectories do not change between genders. They simply shift horizontally between the two levels of gender.

References

Fox, J. and J. Hong (2009). Effect displays in R for multinomial and proportional-odds logit models: Extensions to the effects package. Journal of Statistical Software 32:1, 1–24, <http://www.jstatsoft.org/v32/i01/>.

Harrell, Frank. Regression Modeling Strategies. Springer, 2001.

R Core Team (2017). R: A language and environment for statistical computing. R Foundation for Statistical
Computing, Vienna, Austria. URL https://www.R-project.org/.

For questions or clarifications regarding this article, contact the UVa Library StatLab: statlab@virginia.edu

Clay Ford
Statistical Research Consultant
University of Virginia Library
May 10, 2017