# RMS Semiparametric Ordinal Longitudinal Model

**URL:** <https://discourse.datamethods.org/t/rms-semiparametric-ordinal-longitudinal-model/4819>\
**Category:** models\
**Tags:** ordinal, longitudinal\
**Created:** [September 17, 2021, 1:12am UTC](https://discourse.datamethods.org/t/rms-semiparametric-ordinal-longitudinal-model/4819 "2021-09-17T01:12:41Z")\
**Posts on this page:** 1\
**Showing post:** 110

<div class="post-metadata">

**Author:** ![Johannes\_Schwenke](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/johannes_schwenke/32/5052_2.png) [@Johannes\_Schwenke](https://discourse.datamethods.org/u/Johannes_Schwenke)\
**Post date:** [February 16, 2026, 11:32am UTC](https://discourse.datamethods.org/t/rms-semiparametric-ordinal-longitudinal-model/4819/110 "2026-02-16T11:32:17Z")

</div>

I’ve been toiling away on an R Package for easier modeling of longitudinal ordinal outcomes with VGAM. If one wants to target the average treatment effect (ATE), interval estimation can be very time-intensive with nonparametric bootstrap refitting: for a trial with about 1,000 patients, it can take around 30 minutes without parallelization, even with fairly optimized/vectorized code.

I therefore explored alternatives, including drawing coefficient vectors from a (cluster-robust) multivariate normal (MVN). That is much faster (often by at least an order of magnitude).

However, after reading [Ye et al. 2023](https://pmc.ncbi.nlm.nih.gov/articles/PMC10665030/) and discussion on [Bluesky](https://bsky.app/profile/schwenkej.bsky.social/post/3melnzmk5n22a), my  
understanding is that MVN simulation mainly propagates coefficient uncertainty while conditioning on the observed covariate distribution X (i.e., closer to a sample average treatment effect (?)). In preliminary simulations under a superpopulation setup, I see slight CI undercoverage with this approach, which could worry some stakeholders.

If I understand correctly, nonparametric bootstrap should capture variability in X, as we resample X, but it just takes too long when simulating many scenarios, even on 80+ cores.

So I’ve been experimenting with a one-step score-based approach related to [Score Based Approach to Wild Bootstrap Inference](https://eml.berkeley.edu/~pkline/papers/ScoreFinal_web.pdf), and I’d appreciate feedback on whether this is theoretically reasonable.

My approach is as follows: We approximate uncertainty by simulation.

- Draw cluster (patient) weights w\_g \sim \mathrm{Exp}(1) (exponential multiplier weights, analogous to Bayesian bootstrap weighting).

- Compute centered multipliers u\_g = w\_g - 1, and form a perturbed score at \hat\beta: \sum\_g u\_g S\_g(\hat\beta), where S\_g is the patient-aggregated score contribution.

- Use one Newton step with the model’s variance-covariance information to get approximate perturbed coefficients instead of fully refitting.

- Repeat this many times to generate \beta^{\*}, then run the state-occupancy calculation pipeline.

- For marginalization, reuse the same draw’s cluster weights (mapped to baseline patients and normalized), so variability in X is propagated along with coefficient uncertainty.

I know that the Bayesian Bootstrap doesn’t actually have the goal of frequentist coverage… But from some initial simulations things do look better than when simulating from an MVN. However, I usually assume that there is some theoretical arguments against most things one could be doing. So, is there a clear theoretical reason this approach should be avoided for frequentist inference on superpopulation-targeted marginal effects? Thanks!

---

_[View the full topic](https://discourse.datamethods.org/t/rms-semiparametric-ordinal-longitudinal-model/4819)._