R help - Apr 2025 - Estimating regression with constraints in model coefficients

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/20250421/b7f8cea4/attachment.sig>

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 at
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