Gregg Powell
2025-Apr-21 16:25 UTC
[R] 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>
Christofer Bogaso
2025-Apr-21 18:25 UTC
[R] Estimating regression with constraints in model coefficients
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