Why do cutpoints and effects in ordinal regression become inflated when nested observations are used?

Hi - this is a bit of a more stats than brms/stan specific question, but I haven’t received a response on typical stats forums, and I also wonder if this is something other people using brms have encountered this and understand what is going on:

I’ve been trying to conduct some power analyses/calibration for a hierarchical, ordinal regression. In a simple case a group is measured at Time 1, treated, and measured again afterwards at Time 2. I simulate data using mvnorm function in R, it simply lets you set a correlation between paired observations, and so you can specify e.g., that Time2 should be correlated with Time1, but .5 higher on a standard normal scale. To get the data into and ordinal form I place cutpoints along the normal distribution to get 5 responses (<-1.5, <-.5, <.5, <1.5, >= 1.5).

If I run a cumulative ordinal regression with a probit link simply as Time 1 vs. Time 2, it perfectly recovers these cutpoints and the effect size of .5. However, if I specify in the model that observations are nested within subjects - i.e., each time 1 and 2 observation is paired - there is an inflation in the estimates of both the cutpoints and the effect size (e.g., rather than .5, the effect size might be .7 or so), and this inflation increases with increasing correlations between time 1 and time 2.

Why is this occurring? Notably, it doesn’t occur if I run the same regression models not as an ordinal regression but simply on the normally distributed data, so I think it is something specific to the ordinal model, and not something generally for nested models.

Why does nesting the observations within subjects inflate the estimates? I’ve attached screenshots of the output with vs. without the within ppt nesting below.


@JimBob hello, might this be related to the implementation of the hierarchical intercept in the cumulative model? Because the model strictly has no ‘Intercept’ term (as defined in the linear case), the (1|ppn) is included simultaneously on all the thresholds (they are all equally shifted). I’m thinking this would increase the total variation of the latent variable, but of course the thresholds can then adjust to accommodate it.

Can your model be expressed in terms of the predicted ordinal category probabilities (you can get this from ‘fitted()’)? I imagine those predictions should be equivalent, even if they differ on the latent scale.

Thanks for your response @AWoodward - when I get a chance in the next day or two, I will try this out and see. Is it correct though, how I coded the formula, to indicate that it is repeated measures for subjects, or is this done differently when you have a cumulative model?

Yes this specification is analogous to a random intercept model in the GLM context.

I think there is also scope in brms to have threshold-specific variability terms, but I haven’t personally explored that. Seems like the sort of thing that would need very informative data.

1 Like