[R] Estimating regression with constraints in model coefficients

[Ben Bolker]

   Section 2 of the vignette for the ordinal package:

https://cran.r-project.org/web/packages/ordinal/vignettes/clm_article.pdf

   gives a reasonably complete, if short, definition/discussion of the 
log-likelihood framework for ordinal models. It's probably also 
discussed in Venables and Ripley (Modern Applied Statistics in S), since 
the polr function is in the MASS package ...

On 2025-04-21 2:25 p.m., Christofer Bogaso 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
> 
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide https://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

-- 
Dr. Benjamin Bolker
Professor, Mathematics & Statistics and Biology, McMaster University
Director, School of Computational Science and Engineering
 > E-mail is sent at my convenience; I don't expect replies outside of 
working hours.