We’ve been working on an R package for Markov ordinal transition models over the past year and it’s finally in a state that’s pretty usable.
I called it mostr (MOST was already taken) and you can find it here on github.
It’s not quite ready yet for CRAN…
My main goal was to lower the barrier for people to use MOST models, by using very similar workflow as the marginaleffects package. I focused mainly on marginal effects but mostr also supports predictions, contrasts, and variance estimation for particular patient profiles. Most of the functionality is shown in the vignettes.
Analytic variance estimation for second-order models, partial-PO models, and estimands like time benefit is still missing (it gets complicated very quickly), but optimized resampling approaches support any type of estimand.
NB: Especially for making our code faster (porting to C++), I relied quite a lot on Claude.
An about-to-be-released major update of the orm function adds two types of random effects that are pertinent for MOST models. The first is basic random intercepts which I think of as principled GEE to make treat effect standard errors more accurate. The second is an adjustment for an increasing-over-time variance of Y problem caused by conditioning on Y previous rather than on residuals from Y previous. I am thinking that the first type of random effect, at least, should be standard in MOST to make standard errors more robust.
For estimation I condition on random effect = 0 so that no code changes are needed for prediction.
Random effects are implemented so that scaling to large numbers of distinct Y values is efficient. For example, the ordinal package cannot efficiently handle 100 distinct Y values when there are random effects, but orm can efficiently handle thousands.
Making the Hmisc package function for computing state occupancy probabilities easy to use and making bootstrap confidence intervals for such derived quantities easy to get would be nice extensions of your package if you don’t already have those. Keep up the great work.
We have implemented our own bootstrapping approaches (see here and here), either NPM or fractional weighted bootstrap. We wrote our own implementation because we wanted to process all patients at once and as efficiently as possible (IIRC your Hmisc function only allows processing of one patient profile at a time).
I saw that you were working on this! We have adressed the intra-patient correlation using robust (sandwich) variance estimation, which controls frequentist type I error very well in all kinds of scenarios in our simulations.
I’m not yet sure how to handle random effects when calculating the variance for the ATE analytically, because I’m quite sure that just setting them to 0 will lead to undercoverage for derived estimands like difference in time in state 1 between groups, when actually randomly sampling from a superopulation. (Ofc we don’t randomly sample in reality, but one might still want to show coverage for that case)
I remember that when using Predict() using a blrm() model, RE are also set to 0. We built our own prediction function so the user can include or exclude the random effects, but we don’t support completely integrating them out.
The sandwich estimator will give fairly valid standard errors but will not correctly estimate the (remaining) Markov structure (dependence on Yprev) when the random effect variance is large. The sandwich estimator will also not be able to make the tailored dual random effects variance adjustment. The link I provided above explains this second form.
The sandwich estimator may also have trouble when cluster size variation is huge.
Tbh, I would expect Claude Opus 5.5 to pretty much one shot a Python version of this package.
I haven’t thought much about performance metrics yet. What did you have in mind for calibration. Currently more interested in MSE when some assumptions are violated / we borrow accross the scale when we shouldn’t.
The beauty of first order Markov models is that the previous Y value comes an ordinary covariate. So you can model transitions with any software or even ML. The only trick is to convert transition probabilities into state occupancy probabilities, to handle certain random effects patterns, and to allow the exposure variable to have different effects over time (e.g., one state occurs early, another state occurs late; non-proportional odds in time).