Regression Models with Ordered Categorical Outcomes: clarification of Example

1

Something in the example at Regression Models with Ordered Categorical Outcomes — PyMC example gallery has thrown me off.
This is showing how to use a Dirichlet distribution for priors over cutpoints in an ordered response model.
What’s confusing to me is that the Dirichlet-version cutpoints are constrained to lie in the [0,1] interval, yet there is (just as in the other, sorted normals approach) no scale nor offset (constant) allowed for eta. Surely one normally needs either flexible cutpoints or flexible scale/offset??

2

Maybe you didn’t scroll far enough to the right in the display? It takes a Dirichlet, applies cumulative sum, then scales by multiplying by max - min and then adding min.

The other way to do this directly is by taking a base value c[0] and then defining c[1] = c[0] + delta[1], where you constrain delta > 0, for example by drawing it from a lognormal distribution. Lognormal also helps the boundaries collapsing to each other by making sure delta > 0 (the lognormal density goes to zero as the variate goes to zero).

If you’re OK with making them all positive, you can just take c_\text{diffs} \sim \textrm{lognormal}(...) and then setting c = \textrm{cumulativeSum}(c_\text{diffs}). Or you can just model the diffs as lognormal and give c[0] some kind of unconstrained prior.

3

Yet another option is to use the ordered transform:

        cp_scale = pm.HalfNormal(f"cutpoint_sd", sigma=2)
        cp_len = len(obs.columns)-1
        
        cutpoints = pm.ZeroSumNormal(
            "cutpoints",
            sigma=cp_scale,
            transform=pm.distributions.transforms.ordered,
            initval=np.linspace(-1,1,cp_len),
            shape=cp_len
        )
4

I remember the ordered not behaving nicely with the zero-sum-constraint. In general it’s better to encode non-default constraints in the generative graph, so that the forward model matches with the model logp

5

Thanks, Bob (et al). I think I understand what you’re suggesting. My query was just because min and max seemed to be fixed, not parameters to be estimated. At the moment I had adapted the gallery example to get this:

def unconstrainedUniform(K, itsname='cutpoints'):                                                    
    # Offset and scale for cutpoints                                                                                          
    c0 = pm.Normal(itsname+"_c0", 0, 3)                                                                                     
    s = pm.LogNormal(itsname+"_scale", )                                                                     
    return pm.Deterministic(                                                                                                  
        itsname,                                                                                                              
        c0 + s * pt.concatenate(                                                                                              
            [                                                                                                                 
                np.ones(1)*0,                                                                                                                    
                pt.extra_ops.cumsum(pm.Dirichlet(itsname+"_interior", a=np.ones(K - 2))) ,                                
            ]                                                                                                                 
        ),                                                                                                                    
    )
6

Cutpoint models have a multiplicative non-identifiability that needs to be identified. You have covariates x \in \mathbb{R}^{N \times K} and a regression coefficient vector \beta \in \mathbb{R}^K and you compare the linear predictor x \cdot \beta to the cutpoints e.g., x \cdot \beta < c_n. you can multiply \beta by a constant and the outpoint vector c by a constant and get the same result. So you need to identify the scale in some way. One way to do that is to fix max and min. You could just leave it as a simplex without a max and min and that would identify the scale.