Hi all
Ironically I’ve been working on something like this just today, just to build some intuition around ordinal regressions again as a refresher. Like @ceebeelee I am trying to build up the model from bare-bones without relying on OrderedLogistic
. I’m following the model McElreath describes in Statistical Rethinking 2nd ed, page 385, using the ordinal data from Kruschke’s DBDA2 ('OrdinalProbitData-1grp-1'
). To that data I’ve just added a column of 1’s for the intercept.
I think I have this right, prior predictive has the correct dimensions as far as I can see, but when sampling I encounter a Rewrite failure due to: local_IncSubtensor_serialize
AssertionError. I am a bit lost as to how to solve this one! Here is my code - I know this is a bit inefficient but just helps to cement knowledge for now:
c = {'K': range(d['y'].nunique() - 1),
'Kt': range(d['y'].nunique()),
'N': range(len(d))
}
with pm.Model(coords=c) as ordinal:
# Priors for cutpoints
cut = pm.Normal('cut',
mu=np.linspace(0, 5, num=len(c['K'])),
sigma=1,
transform=pm.distributions.transforms.univariate_ordered,
dims='K')
# Terms of the linear model
β = pm.Normal('β', mu=0, sigma=1)
# linear predictor
η = pm.Deterministic('η', β * d['intercept'].values, dims='N') # just the column of 1's I added
η1 = pm.Deterministic('η1', cut - η[:, None], dims=('N', 'K'))
Η = pm.Deterministic('Η', pm.math.invlogit(η1), dims=('N', 'K'))
# Calculate the category probabilities from this
p = pt.concatenate(
[
1 - Η[:, [0]],
Η[:, :-1] - Η[:, 1:],
Η[:, [-1]]
],
axis=1
)
# Likelihood
pm.Categorical('y', p=p, observed=d['Y'].values - 1, dims='N')
# inferencebutton
idata = pm.sample()
I’ve tried compiling this with pymc.sampling_jax.sample_numpyro_nuts
, and I get an incompatible shape for broadcasting error - one tensor has dims (1, 6) and requested shape is (100, 6). I think the problem seems to stem from the concatenation line where I compute the category probabilities directly, but I’ve used this in some Item Response Theory models which have basically the same logic as ordinal regression (graded response model), and it works just fine! Any tips really welcome here, I’m not sure if this is just me missing something.
Edit: I am on PyMC 5.6.