[R] Estimating regression with constraints in model coefficients

[Gregg Powell]

Hello again Christofer,
Thanks for your thoughtful note — I’m glad the outline was helpful. Apologies for the long delay getting back to you. Had to do a small bit of research…

Recommended Text on Ordinal Log-Likelihoods:
You're right that most online sources jump straight to code or canned functions. For a solid theoretical treatment of ordinal models and their likelihoods, consider:
"Categorical Data Analysis" by Alan Agresti
– Especially Chapters 7 and 8 on ordinal logistic regression.
– Covers proportional odds models, cumulative logits, adjacent-category logits, and the derivation of likelihood functions.
– Provides not only equations but also intuition behind the model structure.
It’s a standard reference in the field and explains how to build log-likelihoods from first principles — including how the cumulative probabilities arise and why the difference of CDFs represents a category-specific probability.
Also helpful:
"An Introduction to Categorical Data Analysis" by Agresti (2nd ed) – A bit more accessible, and still covers the basics of ordinal models.
________________________________________

If You Want to Proceed Practically:
In parallel with theory, if you'd like a working R example coded from scratch — with:
• a custom likelihood for an ordinal (cumulative logit) model,
• fixed thresholds / no intercept,
• coefficient bounds,
• and a sum constraint on β

I’d be happy to attempt that using nloptr() or constrOptim(). You’d be able to walk through it and tweak it as necessary (no guarantee that I will get it right 😊)

Just let me know:
1. Whether you want the cutpoints fixed (and to what values), or if you want them estimated separately (with identifiability managed some other way);
2. What your bounds on the coefficients are (lower/upper vectors),
3. What value the sum of coefficients should equal (e.g., sum(β) = 1, or something else);
4. And whether you're working with the IDRE example data, or something else.

I can use the UCLA ologit.dta dataset as a basis if that's easiest to demo on, or if you have another dataset you’d prefer – again, let me know.

All the best!

gregg






On Monday, April 21st, 2025 at 11:25 AM, Christofer Bogaso <bogaso.christofer using gmail.com> wrote:

> 

> 

> Hi Gregg,
> 

> I am sincerely thankful for this workout.
> 

> Could you please suggest any text book on how to create log-likelihood
> for an ordinal model like this? Most of my online search point me
> directly to some R function etc, but a theoretical discussion on this
> subject would be really helpful to construct the same.
> 

> Thanks and regards,
> 

> On Mon, Apr 21, 2025 at 9:55 PM Gregg Powell g.a.powell using protonmail.com wrote:
> 

> > Christofer,
> > 

> > Given the constraints you mentioned—bounded parameters, no intercept, and a sum constraint—you're outside the capabilities of most off-the-shelf ordinal logistic regression functions in R like vglm or polr.
> > 

> > The most flexible recommendation at this point is to implement custom likelihood optimization using constrOptim() or nloptr::nloptr() with a manually coded log-likelihood function for the cumulative logit model.
> > 

> > Since you need:
> > 

> > Coefficient bounds (e.g., lb ≤ β ≤ ub),
> > 

> > No intercept,
> > 

> > And a constraint like sum(β) = C,
> > 

> > …you'll need to code your own objective function. Here's something of a high-level outline of the approach:
> > 

> > A. Model Setup
> > Let your design matrix X be n x p, and let Y take ordered values 1, 2, ..., K.
> > 

> > B. Parameters
> > Assume the thresholds (θ_k) are fixed (or removed entirely), and you’re estimating only the coefficient vector β. Your constraints are:
> > 

> > Box constraints: lb ≤ β ≤ ub
> > 

> > Equality constraint: sum(β) = C
> > 

> > C. Likelihood
> > The probability of category k is given by:
> > 

> > P(Y = k) = logistic(θ_k - Xβ) - logistic(θ_{k-1} - Xβ)
> > 

> > Without intercepts, this becomes more like:
> > 

> > P(Y ≤ k) = 1 / (1 + exp(-(c_k - Xβ)))
> > 

> > …where c_k are fixed thresholds.
> > 

> > Implementation using nloptr
> > This example shows a rough sketch using nloptr, which allows both equality and bound constraints:
> > 

> > > library(nloptr)
> > > 

> > > # Custom negative log-likelihood function
> > > negLL <- function(beta, X, y, K, cutpoints) {
> > > eta <- X %*% beta
> > > loglik <- 0
> > > for (k in 1:K) {
> > > pk <- plogis(cutpoints[k] - eta) - plogis(cutpoints[k - 1] - eta)
> > > loglik <- loglik + sum(log(pk[y == k]))
> > > }
> > > return(-loglik)
> > > }
> > > 

> > > # Constraint: sum(beta) = C
> > > sum_constraint <- function(beta, C) {
> > > return(sum(beta) - C)
> > > }
> > > 

> > > # Define objective and constraint wrapper
> > > objective <- function(beta) negLL(beta, X, y, K, cutpoints)
> > > eq_constraint <- function(beta) sum_constraint(beta, C = 2) # example >C
> > > 

> > > # Run optimization
> > > res <- nloptr(
> > > x0 = rep(0, ncol(X)),
> > > eval_f = objective,
> > > lb = lower_bounds,
> > > ub = upper_bounds,
> > > eval_g_eq = eq_constraint,
> > > opts = list(algorithm = "NLOPT_LD_SLSQP", xtol_rel = 1e-8)
> > > )
> > 

> > The next step would be writing the actual log-likelihood for your specific problem or verifying constraint implementation.
> > 

> > Manually coding the log-likelihood for an ordinal model is nontrivial... so a bit of a challenge :)
> > 

> > All the best,
> > gregg powell
> > Sierra Vista, AZ
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 603 bytes
Desc: OpenPGP digital signature
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20250423/d78ca286/attachment.sig>