The binary model will give probabilities that depend on the frequency distribution of individual outcome categories within the coarsened category combinations. And it depends on the single threshold being ‘magical’, i.e., one every expert would
The partial prop. Odds model would be better than either, but may require knowing ahead of time which predictors are likely to not act in prop. Odds.
This study tested the effect of Fluoxetine on functional recovery after stroke, which was measured at 6 months using the modified Rankin Scale (mRS; 0-5). The authors used ordinal (PO) logistic regression for the analysis. Their model was the following:
\ mRS_{6m} \text ~ \ Group \ + \ Time \ + \ Predicted \ Outcome \ + \ Motor \ Deficit \ + \ Aphasia
where,
Time = time since stroke (2-8 days vs 9-15 days)
Predicted Outcome = probability of mRS of 0-2 at 6 months (≤0.15 vs >0.15), obtained from a prediction model
Motor Deficit = presence of motor deficit at randomization (yes vs no)
Aphasia = presence of aphasia at randomization (yes vs no)
I see the dichotomization of Time and Predicted Outcome as a problem. The very use of a prediction/ probability as a covariate, in place of the baseline mRS for instance, is problematic in my view. It seems that they didn’t have the baseline mRS available to use in the model. Also, they didn’t assess whether the PO assumption was met/ appropriate.
Any thoughts/ opinions?
To do a study of mRS without collecting mRS at baseline is a key design flaw. And the PO assumption is the least of their worries. The flaws you pointed out, especially dichotomization of a continuous risk measure, are huge. A spline in logit risk should have been used. Dichtomizing time is pretty bad also.
I saw that you mentioned on a different post that when we have an ordinal covariate (mRS) it is better to model it as nominal instead of as ordinal. Does that apply here too?
I didn’t mean to say that. I was saying that if there’s only say 4 levels it’s not worth the trouble of specially modeling it as ordinal.
I’m fitting this model but the bootstrap is not working:
m <- lrm(y ~ age + sex + lead73, x = T, y = T, data = lead)
b <- bootcov(m, B = 500, coef.reps = T, group = lead$y)
Warning: at least one resample excluded highest Y values, invalidating bootstrap covariance matrix estimate
The model has 39 intercepts.
@f2harrell, any suggestion?
Also, would the code below be correct for getting the contrast on the mean scale?
M <- Mean(m)
c <- contrast(b,
list(lead73 = quantile(lead$lead73, 0.75)),
list(lead73 = quantile(lead$lead73, 0.25)))
print(c, fun = M)
Thank you
fun= goes with contrast not print. And contrast will use the bootcov output to get approximate CLs on differences in means. The problem with some bootstrap samples not including some Y levels is a very tricky one. Sometimes it’s OK to give bootcov a group=Y type of argument to balance the bootstrap but performance of this has not been studied. At other times a bit of rounding of Y will suffice. This is another reason to go Bayesian.
Dear sir,
I had reworked the model in r so as to get better clarity. But i still am not able to understand if there are issues to the model. I had run the brant test so as to ensure that the proportional odds assumption holds, but i am getting the result as
In brant(modeltest) :
24111 combinations in table(dv,ivs) do not occur. Because of that, the test results might be invalid.
Does this mean that i shall have to go for the partial proportional odds model or generalized ordered logit model?
The Brant test is not good as it is anti-conservative. Instead see https://fharrell.org/post/impactpo
Dear sir,
I am not able to open this link as it might not be accessible here in India. Is there any alternative option I can take into consideration, as a pdf of the page or something similar.
Thanking you in advance
Sunena
I think it’s just an incorrect link. Try https://www.fharrell.com/post/impactpo/
sir,
Thank you for the link, it was really helpful. Is there any specific method in R that helps in checking whether the observations obtained from the sample population are unbiased while building the model, apart from the sampling technique chosen alone?
It’s all in the study design and how measurements are ascertained.
Dear Professor Harrell,
I’m running an ordinal regression model (orm) with the following formula, where HAM-D is the Hamilton Rating Scale for Depression:
orm(HAM-D ~ rcs(age) + sex + other covariates)
I’m having difficulty explaining the model output to some collaborators without framing it in terms of cut-offs of the HAM-D score. Could I ask you a couple of questions?
Can you direct me to a resource that explains how to calculate the estimated difference in exceedance probability and SEs for various cut-offs? For example, calculating Pr(\text{HAM-D} > 7 | \text{age} = 30) - Pr(\text{HAM-D} > 7 | \text{age} = 20).
Given the proportional odds (PO) assumption, would it be accurate to report the beta coefficients as “an increase in 1 unit of covariate A is associated with a exp(Beta)-fold increase in the odds of a HAM-D score > 7” instead of “odds of having a higher HAM-D score or HAM-D > j?”
Any suggestions or corrections would be greatly appreciated.
Thank you so much
Great questions. This has an example showing how to get the whole predicted Y | X distribution. And you can use the ExProb function in rms to derive a function to compute any desired exceedance probability like \Pr(Y \geq y | X). You might compute that for a variety of y.
What is perhaps more helpful to collaborators is to get predicted mean Y | X which the Mean function (full name Mean.orm) makes easy to do. In all our work in orthopedics where various functional status/disability scales are used we routinely presented partial effect plots when the estimated mean Y on the y-axes.
On your second question, an accurate way to state an odds ratio, when proportional odds is being assumed as you are doing, is for example:
Assume the age range you want to use is the inter-quartile range and that corresponds to ages 20 and 40: The age 40 : age 20 odds ratio for HAM-D \geq any specific number is … with a footnote that the “specific number” is any HAM-D value other than the lowest possible value. Then you can add:
e.g. OR for HAM-D \geq 7 which equals OR for HAMD-D \geq 10
But start with mean HAM-D.
I would like some advice a project I am currently working on that uses proportional odds models. The investigators for this project very much want to have an estimate of the typical change score between timepoints in a study with three time points comparing two randomization arms. They especially are interested in the change score because other studies have developed a minimal important difference (MID) for the instrument. While I could use a typical linear mixed effects model that assumes normality to estimate mean changes and differences in those changes between the arms, I thought it is better to use a proportional odds model with a random intercept to compare changes over time between the randomization arms. I am able to use this model to provide estimates of median scores but I don’t see how I could use a proportional odds to estimate the median change for both arms. My current solution has been to advocate for simply reporting the estimated medians and then using odds ratios scaled to the MID. Do you have any suggestions for how you might handle this situation?
As a follow-up question, I am curious about how to think about incorporating non-linear relationships between baseline and follow-up values when there are multiple time points. How should one handle the non-linear relationship between subsequent values in a study with multiple time points? I would think that a random intercept only accounts for a linear relationship on the link scale between subsequent values. One could use other random coefficients, but fitting such models could become challenging.
Thank you for any input!
Rob
Welcome to datamethods Rob.
As discussed here framing MCID in terms of change in a patient’s status is problematic, and as discussed here change from baseline is problematic in general. When using ordinal models, we need the response variable to be ordinal, and a difference between two ordinal variables (baseline and follow-up) is not necessarily ordinal. For example a change from 1 to 2 on an ordinal scale that runs from 1-10 may not be comparable to a change from 8 to 9.
One other general problem is which model is likely to fit the correlation structure of your data. As argued here, it is not necessarily natural to think of random effects to handle within-person correlation. Random intercepts induces a compound symmetric correlation structure, e.g., it assumes that the within-person correlation in responses measured at 2 months and 3 months is the same as the correlation between 2 and 24 months. This is very unlikely to fit the data. Ordinal Markov longitudinal models are excellent choices for handling the more likely-to-fit serial correlations we see in practice. See this and this.
To get at your question of how to use an ordinal model to get estimates on the original scale, see one of the replies posted earlier on this page where I talk about estimating the mean or differences in means. To estimate quantiles instead (which make a strong assumption that Y is continuous, i.e., has few ties in the data) use the Quantile function instead of the Mean function. The contrast function (rms::contrast.rms) can come into play here. If doing Bayesian ordinal modeling with rmsb::blrm it is easier to deal with uncertainty intervals of complex nonlinear transformations of the original ordinal model parameters (on the log-odds scale).
For simple examples of using ordinal models to estimate means and quantiles see this.
On your question about nonlinear baseline covariate effects, we try to think of these as general effects that have the same shape over time, i.e., hope for no covariate by time interaction (though that can be extended, if the sample size allows). So having for example a restricted cubic spline function of baseline with 4 knots should pose no problem.
It’s good to recall the reasons for going to this trouble.
Dear Dr Harrell,
Apologies if this is obvious or has been covered elsewhere, but I couldn’t see anything.
Following on from the discussions at the RMS25 course last week, I was thinking about how to apply an ordinal model to count data where we would usually use a Negative Binomial or Poisson model - for example, for the number of exacerbations in 1 year in respiratory diseases. In particular, it seems like this could accommodate deaths during follow-up. In a Negative Binomial or Poisson model, where a patient does not have complete follow-up, including an offset term would be one option to deal with this (under a MAR assumption).
It looks like an offset term can be included in an ordinal model. However, I’m not convinced this would have the same function as in a parametric count data model? Would the inclusion of a death category change this? The only other option I can think of is multiple imputation of the incomplete follow-up period data - do you think this the only/best option?
Any thoughts would be appreciated. Many thanks!
Consider this as count data makes an underlying assumption of a constant rate and makes it very hard to handle deaths and missing data. I don’t know about the offset question. A state transition model solves all these problems and makes it easy to model changing transition probabilities over time. I don’t get how you would handle deaths with the setup you described.
Apologies Dr Harrell, I confused myself by thinking about 2 problems at once, and explaining neither of them. The state transition model sounds like exactly what I need for the clinical trial data if there are deaths and incomplete follow-up - thank you!
The reason I was looking at an ordinal model was for a prediction model - I’ve been asked about predicting probability of at least 1 exacerbation. I’m not expecting many deaths at all, but I was just trying to think whether I could/should include any deaths by ranking them worse than the maximum observed exacerbations per year. This is from a registry with annual follow-up visits, so no partial annual follow-up.
The offset issue was an afterthought that I was curious about, which I should have asked separately.
Sorry again for my muddled question - I’ll try to state things clearer in future.