I have a dataset of patients with associated risk factors and outcomes. I am interested in the effect COPD (hxcopd), which is a binary variable, has on an ordinal outcome with 6 possible outcomes (outcome) increasing in severity (0 through 5). I will be writing out my train of thought. I am not very confident in this interpretation so I am looking for any criticisms of my thought process or ultimate intepretation.
In R (polr) the ordinal logistic regression model is parameterized as:
$$logit (P(Y \le j)) = \beta_{j0} – \eta_{1}x_1 – \cdots – \eta_{p} x_p$$
Due to the parallel lines assumption, even though I have six categories, the coefficient of COPD (hxcopd) stays the same across the five categories. The two equations for hxcopd = 1 and hxcopd = 0 are
$$ \begin{eqnarray} logit (P(Y \le j | x_1=1) & = & \beta_{j0} – \eta_{1} \\ logit (P(Y \le j | x_1=0) & = & \beta_{j0} \end{eqnarray} $$
Then $logit (P(Y \le j)|x_1=1) -logit (P(Y \le j)|x_1=0) = – \eta_{1}.$
To perform this R, I create an unadjusted fit using the following R function:
unadjfit <- MASS::polr(outcome ~ hxcopd, data = df)
To summarize the results:
summary(unadjfit)
Call:
MASS::polr(formula = outcome ~ hxcopd, data = df)
Coefficients:
Value Std. Error t value
hxcopdTRUE 0.331 0.1297 2.552
Intercepts:
Value Std. Error t value
0|1 3.5431 0.0985 35.9645
1|2 3.6880 0.1015 36.3219
2|3 4.1911 0.1150 36.4485
3|4 4.8785 0.1431 34.0918
4|5 5.9261 0.2178 27.2052
Residual Deviance: 2964.116
AIC: 2976.116
The coefficient in the summary is $\eta_1$.
1) The above statement is true, correct? I can not find a clear answer in the polr documentation.
The output shows that for patients with COPD, the log odds of being category 0 (versus category > 0) is actually $-\hat{\eta}_1=-0.331$ or $0.331$ points lower than patients without COPD.
2) Could I also say "that for patients with COPD, the log odds of being in a category $\leq J$ vs. $>J$ is $0.331$ points lower than patients without COPD"?
Since the coefficient $ – \eta_{1}$ represents a one unit change in the log odds when moving between the two exposures (i.e., hxcopd = 1 to hxcopd = 0), we can write:
$$logit (P(Y \le j|x_1=1) -logit (P(Y \le j|x_1=0) = – \eta_{1}.$$
Exponentiate both sides of this equation and use $log(b)-log(a) = log(b/a)$:
$$\frac{P(Y \le j |x_1=1)}{P(Y>j|x_1=1)} / \frac{P(Y \le j |x_1=0)}{P(Y>j|x_1=0)} = exp( -\eta_{1}).$$
Which by the proportional odds assumption can be simplified:
$$\frac{P(Y \le j |x_1=1)}{P(Y>j|x_1=1)} = p_1 / (1-p_1) $$
$$\frac{P(Y \le j |x_1=0)}{P(Y>j|x_1=0)} = p_0 / (1-p_0)$$
The odds ratio is defined as:
$$\frac{p_1 / (1-p_1) }{p_0 / (1-p_0)} = exp( -\eta_{1}).$$
But since R does not give us $-\eta$, but rather $\eta$, exp(coef(unadjfit)) gives a different OR:
exp(coef(unadjfit))
hxcopdTRUE
1.392349
Since $exp(-\eta_{1}) = \frac{1}{exp(\eta_{1})}$,
$$exp(\eta_{1}) = \frac{p_0 / (1-p_0) }{p_1 / (1-p_1)}.$$
From the output, $\hat{\eta}_1=0.331$,the odds ratio $exp(\hat{\eta}_1)=1.39$ is actually $\frac{p_0 / (1-p_0) }{p_1 / (1-p_1)}$.
3) So a correct interpretation of the polr output is that people without COPD have higher odds of being in a category $\leq J$ vs. $>J$ when compared to patients with COPD?
If I wanted to interpret differently, I could flip the odds ratio around:
$$ \begin{eqnarray} exp(-\eta_{1}) & = & \frac{p_1 / (1-p_1)}{p_0/(1-p_0)} \\ & = & \frac{p_1 (1-p_0)}{p_0(1-p_1)} \\ & = & \frac{(1-p_0)/p_0}{(1-p_1)/p_1} \\ & = & \frac{P (Y >j | x=0)/P(Y \le j|x=0)}{P(Y > j | x=1)/P(Y \le j | x=1)}. \end{eqnarray} $$
Since $exp(-\eta_{1}) = \frac{1}{exp(\eta_{1})}$,
$$\frac{P (Y >j | x=1)/P(Y \le j|x=1)}{P(Y > j | x=0)/P(Y \le j | x=0)} = exp(\eta).$$
Instead of interpreting the odds of being in category $\leq J$, we can interpret the odds of being in category $>J$.
4) Is it appropriate to interpret $exp(\hat{\eta}) = exp(0.331) = 1.39$ as "people with COPD have 1.39 times the odds of being in category $>J$ compared to people without COPD"?
5) If I wanted to report the results of the ordinal regression in an academic manuscript, what is usually expected? I assume the OR and 95% CI. Is it customary to also report the intercepts for each level of the outcome?
Credit to UCLA Statistical Consulting for their great walkthrough of ordinal logistic regression.
Best Answer
Nicely laid out question, Dylan. I'll take a stab at answering it but will keep my answer practical (i.e., without using mathematical equations).
Will you change the sign of the hxcopd coefficient for reporting purposes?
The first thing you need to determine when looking at the Coefficients output produced by polr is whether you are going to change the sign of the reported coefficient for the purposes of your interpretation or not. In your case, are you going to interpret the coefficient of hxcopdTRUE directly (i.e., 0.331) without changing its sign or are you going to interpret the changed-sign coefficient of -0.331?
What groupings of values for your response variable are you really interested in comparing?
If you are NOT going to change the sign of the reported coefficient by multiplying that coefficient by -1 (i.e., if you are going to interpret 0.331), the ensuing interpretation will allow you to compare these groupings of values for your response variable in terms of log odds:
If you ARE going to change the sign of the reported coefficient by multiplying that coefficient by -1 (i.e., if you are going to interpret -0.331), then your interpretation will involve the following groupings of values for the response variable:
In the latter case, you are comparing more versus less severity; in the former, you are comparing less versus more severity. Thus, you have to be careful which case you choose so that your interpretation appropriately conveys the underlying comparisons.
No change in sign for coefficient of hxcopd
Say you choose to interpret the coefficient of hxcopdTRUE of 0.331 without changing its sign. That coefficient tells you the following:
The odds of having a severity rating of 5 rather than 0, 1, 2, 3 or 4 are estimated to be 1.39 times higher (or 39% higher) for those with COPD than for those without COPD;
The odds of having a severity rating of 4 or 5 rather than 0, 1, 2 or 3 are estimated to be 1.39 times higher (or 39% higher) for those with COPD than for those without COPD;
The odds of having a severity rating of 3, 4 or 5 rather than 0, 1 or 2 are estimated to be 1.39 times higher (or 39% higher) for those with COPD than for those without COPD;
The odds of having a severity rating of 2, 3, 4 or 5 rather than 0 or 1 are estimated to be 1.39 times higher (or 39% higher) for those with COPD than for those without COPD;
The odds of having a severity rating of 1, 2, 3, 4 or 5 rather than 0 are estimated to be 1.39 times higher (or 39% higher) for those with COPD than for those without COPD.
Other language you may see people use in this context would be "the odds are 1.39-fold higher" or "the odds are higher by a multiplicative factor of 1.39".
The above interpretations are repetitive so you would most likely want to consolidate them in a single statement along these lines (or whatever makes sense in your specific setting):
The odds of having a higher rather than a lower severity rating (e.g., 1, 2, 3, 4 or 5 rather than 0; ) are estimated to be 1.39 times higher (or 39% higher) for those with COPD than for those without COPD.
Change in sign for coefficient of hxcopd
Now, if you DO change the sign of your coefficient for hxcopd, your interpretation will also change since you have to interpret -0.331 or exp(-0.331) instead of 0.331 or exp(0.331).
On the log odds scale, you would have this type of interpretation:
The log odds of having a severity rating of 0 rather than 1, 2, 3 or 4 are estimated to be 0.331 points lower for those with COPD than for those without COPD;
The log odds of having a severity rating of 0 or 1 rather than 2, 3, 4 or 5 are estimated to be 0.331 points lower for those with COPD than for those without COPD;
The log odds of having a severity rating of 0, 1 or 2 rather than 3, 4 or 5 are estimated to be 0.331 times lower for those with COPD than for those without COPD;
The log odds of having a severity rating of 0, 1, 2 or 3 rather than 4 or 5 are estimated to be 0.331 points lower for those with COPD than for those without COPD;
The log odds of having a severity rating of 0, 1, 2, 3 or 4 rather than 5 are estimated to be 0.331 times lower for those with COPD than for those without COPD.
On the odds scale, you would have to say things like the ones below, since exp(-0.331) = 0.72 and (0.72-1)x100% = -28%:
The odds of having a severity rating of 0 rather than 1, 2, 3 or 4 are 0.72 times lower (or 28% lower) for those with COPD than for those without COPD;
The odds of having a severity rating of 0 or 1 rather than 2, 3, 4 or 5 are 0.72 times lower (or 28% lower) for those with COPD than for those without COPD;
The odds of having a severity rating of 0, 1 or 2 rather than 3, 4 or 5 are 0.72 times lower (or 28% lower) for those with COPD than for those without COPD;
The odds of having a severity rating of 0, 1, 2 or 3 rather than 4 or 5 are 0.72 times lower (or 28% lower) for those with COPD than for those without COPD;
The odds of having a severity rating of 0, 1, 2, 3 or 4 rather than 5 are 0.72 times lower (or 28% lower) for those with COPD than for those without COPD.
The consolidated statement for this last case could look like this:
The odds of having a lower rather than a higher severity rating (e.g., 0, 1, 2, 3 or 4 rather than 5) are estimated to be 0.72 times lower (or 28% lower) for those with COPD than for those without COPD.
In a manuscript, you would most likely have to report a consolidated statement and add the 95% confidence intervals to the reported points (on the log scale) or to the reported odds ratios (on the odds scale). You would also have to explain that you checked whether the proportional odds assumptions holds for your data. Finally, you would need to be clear about what groupings of values for your response variable you are reporting, as explained above.
I assume you have already read this post: https://stats.idre.ucla.edu/r/faq/ologit-coefficients/. It's worth going through it to convince yourself that you are indeed reporting the appropriate quantities in your case. In particular, after fitting your model, look at the following:
The phat object will report the probability that your response variable takes on a particular value among 0, 1, 2, 3, 4 or 5, separately for those without COPD and those with COPD.
Then, if you want to compute the odds of having a rating of 5 rather than 0,1,2,3 or 4, say, among those with COPD, you would just divide the reported probability for a rating of 5 in the "with COPD row" (i.e., the second row of phat) by the sum of the reported probabilities for ratings of 0, 1, 2, 3 or 4 in the same row. The same odds among those without COPD would be derived by dividing the reported probability for a rating of 5 in the "without COPD row" (i.e., the first row of phat) by the sum of the reported probabilities for ratings of 0, 1, 2, 3 or 4 in the same row. The ratio of the two odds will give you the odds ratio of having a rating of 5 rather than 0,1,2,3 or 4 for those with COPD relative to those without COPD. If this coincides with what comes out of R via the interpretation process described above, you are on the right path!
Addendum
Brant's Wald test is used by some to verify the reasonableness of the proportional odds assumption for each predictor variable in your model and for all of them together (as explained, for instance, in this article by Richard Williams on Understanding and interpreting generalized ordered logit models:
https://www3.nd.edu/~rwilliam/gologit2/UnderStandingGologit2016.pdf.
R has a brant package for this: https://medium.com/evangelinelee/brant-test-for-proportional-odds-in-r-b0b373a93aa2.
There is also the possibility of using a likelihood ratio test for testing the proportionality of odds assumption, as mentioned for instance in this article: Assessing proportionality assumption in the adjacent category logistic regression model by Dolgun et al.: https://www.intlpress.com/site/pub/files/_fulltext/journals/sii/2014/0007/0002/SII-2014-0007-0002-a012.pdf. The likelihood ratio test is an omnibus test of proportionality of odds (hence it considers all predictor variables together). See here, for example: https://stat.ethz.ch/pipermail/r-help/2014-November/423706.html.
You can also check this assumption visually in addition to using formal statistical tests.
One thing you may find helpful in addition to checking assumptions is visualizing the results of your modelling using the effects package in R, as explained here in the post Visualizing the Effects of Proportional-Odds Logistic Regression: https://data.library.virginia.edu/visualizing-the-effects-of-proportional-odds-logistic-regression/.