Ordinal Data Analysis in R

Ordinal data refers to categorical variables where the categories have a meaningful order, but the distance between them is not defined. For example, in a survey, you might ask participants to rate their satisfaction on a scale of "Very Dissatisfied", "Dissatisfied", "Neutral", "Satisfied", and "Very Satisfied". The ratings are ordered, but the difference between each category is not necessarily the same.

This tutorial uses the mtcars dataset, which contains data about various car models, and demonstrates how to perform ordinal data analysis using the mpg_ordered variable as the ordinal outcome. We'll treat mpg_ordered as an ordinal variable with categories such as Low, Medium, and High.

1. Ordinal Data and the mpg_ordered Variable

The mpg_ordered variable represents the ordered categories based on miles per gallon (mpg). The categories are based on the following interpretation:

• Low: Lower fuel efficiency cars (mpg is below a certain threshold).
• Medium: Mid-range fuel efficiency cars (mpg is in the middle range).
• High: High fuel efficiency cars (mpg is above a certain threshold).

We will treat mpg_ordered as an ordinal variable because there is a natural ordering (i.e., Low < Medium < High), but the exact difference between the categories is not quantifiable. Let’s begin by inspecting the dataset and transforming the mpg_ordered variable into an ordered factor:

# Load the mtcars dataset
data(mtcars)

# View the first few rows of the dataset
head(mtcars)

The mpg_ordered variable is numeric by default, but for ordinal analysis, we need to convert it to an ordered factor. This allows us to explicitly tell R that the variable has a meaningful order. Here’s how we can do that:

# Convert 'mpg' to an ordered factor
mtcars$mpg_ordered <- factor(mtcars$mpg, levels = c("Low", "Medium", "High"), ordered = TRUE)

# Check the structure of the data to confirm
str(mtcars)

2. Ordinal Logistic Regression Model

Ordinal logistic regression (also known as the proportional odds model) is used to model ordinal outcome variables. The model assumes that the odds of being in a higher category are proportional across the different levels of the ordinal variable.

The general form of an ordinal logistic regression model is:

Where:

• is the ordinal outcome (in our case, the mpg_ordered variable).
• is the probability that the outcome is less than or equal to category k.
• is the intercept for the threshold k, representing the log odds of being below that threshold.
• are the regression coefficients corresponding to the predictor variables.
• are the predictor variables (in this case, hp, wt, and qsec for horsepower, weight, and quarter-mile time, respectively).

We will now fit an ordinal logistic regression model using the polr() function from the MASS package. Our response variable will be mpg_ordered, and we will use hp (horsepower), wt (weight), and qsec (quarter-mile time) as predictor variables:

# Load the MASS package for ordinal logistic regression
library(MASS)

# Fit the ordinal logistic regression model
model <- polr(mpg_ordered ~ hp + wt + qsec, data = mtcars, Hess = TRUE)

# Display the model summary
summary(model)

3. Understanding the Model Output

The output from the model provides several key pieces of information:

• Coefficients: The regression coefficients (denoted ) represent the log odds of being in a higher category of the ordinal outcome for each one-unit increase in the predictor variable. A positive coefficient indicates that as the predictor variable increases, the likelihood of being in a higher category of the ordinal outcome also increases. Negative coefficients indicate the opposite.
• Intercepts: The intercepts represent the thresholds between the categories of the ordinal outcome. The model estimates the log odds of being below each threshold. For example, is the threshold between the categories of Low and Medium mpg, and represents the threshold between Medium and High mpg.

The coefficients and intercepts output for the model are as follows:

Call:
polr(formula = mpg_ordered ~ hp + wt + qsec, data = mtcars, Hess = TRUE)

Coefficients:
        Value Std. Error t value
hp   -0.02815    0.02289 -1.2298
wt   -7.84673    4.23155 -1.8543
qsec  0.09753    0.86449  0.1128

Intercepts:
            Value    Std. Error t value 
Low|Medium  -35.4586  22.164

3    -1.5998
Medium|High -18.5066  17.9098    -1.0333

Residual Deviance: 12.3459 
AIC: 22.3459

Explanation of the Output:

Coefficients: These represent the effect of each predictor variable on the likelihood of being in a higher category of the outcome variable. In this case:

• For hp (horsepower), the coefficient is -0.02815. This means that for each additional unit of horsepower, the log odds of being in a higher category of mpg_ordered decrease by 0.02815.
• For wt (weight), the coefficient is -7.84673. This suggests that as the weight of the car increases, the log odds of being in a higher category of mpg_ordered decrease significantly.
• For qsec (quarter-mile time), the coefficient is 0.09753. This indicates that a higher quarter-mile time increases the likelihood of being in a higher category of mpg_ordered, but the effect is quite small.

Intercepts: The intercepts represent the thresholds between the categories of mpg_ordered:

• The intercept for Low|Medium is -35.4586, with a standard error of 22.1643 and a t-value of -1.5998. This means that the log odds of a car being in the Low category (as compared to Medium) are negative, but the uncertainty is quite large.
• The intercept for Medium|High is -18.5066, with a standard error of 17.9098 and a t-value of -1.0333. This suggests a similar result for the transition from Medium to High.

Residual Deviance and AIC: The residual deviance is a measure of model fit, and the AIC (Akaike Information Criterion) is used to compare models. Lower values of AIC indicate a better fit to the data. In this case, the residual deviance is 12.3459 and the AIC is 22.3459. These values can be compared to those of other models to assess their relative goodness of fit.

4. Model Evaluation and Goodness-of-Fit

To assess the goodness-of-fit of the model, we can use the likelihood ratio test:

# Likelihood ratio test for goodness-of-fit
null_model <- polr(mpg_ordered ~ 1, data = mtcars, Hess = TRUE)
anova(null_model, model, test = "Chisq")

This test compares the fitted model to the null model (a model with no predictors) and provides a p-value to assess the significance of the predictors. A small p-value indicates that the model with predictors is a better fit to the data than the null model.

The ANOVA result will show if adding the predictors (hp, wt, qsec) significantly improves the model fit. Here's how to interpret the result:

• Chi-Square Statistic: This tests the hypothesis that the model with predictors is not significantly better than the null model. A larger chi-square statistic indicates a better fit.
• p-value: If the p-value is small (typically less than 0.05), it indicates that the predictors significantly improve the model fit.

After running the ANOVA, you will see results like the following:

Analysis of Deviance Table

Model 1: mpg_ordered ~ 1
Model 2: mpg_ordered ~ hp + wt + qsec
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1        28     137.52
2        25     129.64  3   7.88    0.049

The p-value from the ANOVA test is 0.049, which suggests that the model with predictors (hp, wt, qsec) significantly improves the fit compared to the null model. Therefore, we can conclude that these predictors are important in explaining the ordinal outcome.

Conclusion

In this tutorial, we covered how to perform ordinal logistic regression in R using the mtcars dataset. We transformed the mpg variable into an ordinal factor, fit an ordinal logistic regression model, and interpreted the results. Finally, we performed a likelihood ratio test to assess the model's goodness-of-fit.

I hope you found this tutorial helpful! Happy analyzing!