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
signature.asc
Description: OpenPGP digital signature
______________________________________________ R-help@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.
